Changeset 576 for trunk/src/civ_matlab.m
- Timestamp:
- Mar 4, 2013, 8:13:53 AM (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/civ_matlab.m
r575 r576 51 51 check_civ1=0;%default 52 52 check_patch1=0;%default 53 % ImageFileA=Param.Civ1.ImageA; 54 % ImageFileB=Param.Civ1.ImageB; 55 % if isfield(Param,'Civ2') 56 % ImageFileA_civ2=Param.Civ2.ImageA; 57 % ImageFileB_civ2=Param.Civ2.ImageB; 58 % end 53 59 54 % case of input Param set by an xml file (batch mode) 60 55 if ischar(Param) 61 56 Param=xml2struct(Param); %if Param is the name of an xml file, read this file as a Matlab structure 62 % if isfield(Param,'Civ1')63 % if strcmp(Param.Civ1.FileTypeA,'video')64 % Param.Civ1.ImageA=VideoReader(regexprep(Param.Civ1.ImageA,'''','\'));% remove spurious ' appearing in Windows;65 % elseif strcmp(Param.Civ1.FileTypeA,'mmreader')66 % Param.Civ1.ImageA=mmreader(regexprep(Param.Civ1.ImageA,'''','\'));% remove spurious ' appearing in Windows67 % end68 % if strcmp(Param.Civ1.FileTypeB,'video')69 % Param.Civ1.ImageB=VideoReader(regexprep(Param.Civ1.ImageB,'''','\'));% remove spurious ' appearing in Windows70 % elseif strcmp(Param.Civ1.FileTypeB,'mmreader')71 % Param.Civ1.ImageB=mmreader(regexprep(Param.Civ1.ImageB,'''','\'));% remove spurious ' appearing in Windows72 % end73 % end74 % if isfield(Param,'Civ2')75 % if strcmp(Param.Civ2.FileTypeA,'video')76 % Param.Civ2.ImageA=VideoReader(regexprep(Param.Civ2.ImageA,'''','\'));% remove spurious ' appearing in Windows;77 % elseif strcmp(Param.Civ2.FileTypeA,'mmreader')78 % Param.Civ2.ImageA=mmreader(regexprep(Param.Civ2.ImageA,'''','\'));% remove spurious ' appearing in Windows79 % if strcmp(Param.Civ2.FileTypeB,'video')80 % Param.Civ2.ImageB=VideoReader(regexprep(Param.Civ2.ImageB,'''','\'));% remove spurious ' appearing in Windows81 % elseif strcmp(Param.Civ2.FileTypeB,'mmreader')82 % Param.Civ2.ImageB=mmreader(regexprep(Param.Civ2.ImageB,'''','\'));% remove spurious ' appearing in Windows83 % end;84 % end85 % end86 57 end 87 58 … … 116 87 % caluclate velocity data (y and v in indices, reverse to y component) 117 88 [xtable ytable utable vtable ctable F result_conv errormsg] = civ (par_civ1); 118 119 % to try the reverse_pair method, uncomment below120 % [xtable1 ytable1 utable1 vtable1 ctable1 F1 result_conv1 errormsg1] = civ (Param.Civ1);121 % Param.Civ1.reverse_pair=1;122 % [xtable2 ytable2 utable2 vtable2 ctable2 F2 result_conv2 errormsg2] = civ (Param.Civ1);123 % xtable=[xtable1; xtable2];124 % ytable=[ytable1; ytable2];125 % utable=[utable1; -utable2];126 % vtable=[vtable1; -vtable2];127 % ctable=[ctable1; ctable2];128 % F=[F1; F2];129 % result_conv=[result_conv1; result_conv2];130 % errormsg=[errormsg1; errormsg2];131 89 if ~isempty(errormsg) 132 90 return … … 271 229 end 272 230 end 273 ibx2=ceil(par_civ2.CorrBoxSize(1)/2); 274 iby2=ceil(par_civ2.CorrBoxSize(2)/2); 275 isx2=ibx2+4;% search ara +-4 pixels around the guess 276 isy2=iby2+4; 277 % shift from par_civ2.filename_nc1 278 % shiftx=velocity interpolated at position 279 miniy=max(1+isy2,1+iby2); 280 minix=max(1+isx2,1+ibx2); 281 maxiy=min(size(par_civ2.ImageA,1)-isy2,size(par_civ2.ImageA,1)-iby2); 282 maxix=min(size(par_civ2.ImageA,2)-isx2,size(par_civ2.ImageA,2)-ibx2); 283 [GridX,GridY]=meshgrid(minix:par_civ2.Dx:maxix,miniy:par_civ2.Dy:maxiy); 284 GridX=reshape(GridX,[],1); 285 GridY=reshape(GridY,[],1); 286 Shiftx=zeros(size(GridX));% shift expected from civ1 data 287 Shifty=zeros(size(GridX)); 288 nbval=zeros(size(GridX)); 231 % ibx2=ceil(par_civ2.CorrBoxSize(1)/2); 232 % iby2=ceil(par_civ2.CorrBoxSize(2)/2); 233 % isx2=ibx2+4;% search ara +-4 pixels around the guess 234 % isy2=iby2+4; 235 % % shift from par_civ2.filename_nc1 236 % % shiftx=velocity interpolated at position 237 % miniy=max(1+isy2,1+iby2); 238 % minix=max(1+isx2,1+ibx2); 239 % maxiy=min(size(par_civ2.ImageA,1)-isy2,size(par_civ2.ImageA,1)-iby2); 240 % maxix=min(size(par_civ2.ImageA,2)-isx2,size(par_civ2.ImageA,2)-ibx2); 241 % [GridX,GridY]=meshgrid(minix:par_civ2.Dx:maxix,miniy:par_civ2.Dy:maxiy); 242 % GridX=reshape(GridX,[],1); 243 % GridY=reshape(GridY,[],1); 244 if isfield(par_civ2,'Grid')% grid points set as input file 245 if ischar(par_civ2.Grid)%read the grid file if the input is a file name 246 par_civ2.Grid=dlmread(par_civ2.Grid); 247 par_civ2.Grid(1,:)=[];%the first line must be removed (heading in the grid file) 248 end 249 else% automatic grid 250 minix=floor(par_civ2.Dx/2)-0.5; 251 maxix=minix+par_civ2.Dx*floor((par_civ2.ImageWidth-1)/par_civ2.Dx); 252 miniy=floor(par_civ2.Dy/2)-0.5; 253 maxiy=minix+par_civ2.Dy*floor((par_civ2.ImageHeight-1)/par_civ2.Dy); 254 [GridX,GridY]=meshgrid(minix:par_civ2.Dx:maxix,miniy:par_civ2.Dy:maxiy); 255 par_civ2.Grid(:,1)=reshape(GridX,[],1); 256 par_civ2.Grid(:,2)=reshape(GridY,[],1); 257 end 258 Shiftx=zeros(size(par_civ2.Grid,1));% shift expected from civ1 data 259 Shifty=zeros(size(par_civ2.Grid,1)); 260 nbval=zeros(size(par_civ2.Grid,1)); 289 261 if par_civ2.CheckDeformation 290 DUDX=zeros(size( GridX));291 DUDY=zeros(size( GridX));292 DVDX=zeros(size( GridX));293 DVDY=zeros(size( GridX));262 DUDX=zeros(size(par_civ2.Grid,1)); 263 DUDY=zeros(size(par_civ2.Grid,1)); 264 DVDX=zeros(size(par_civ2.Grid,1)); 265 DVDY=zeros(size(par_civ2.Grid,1)); 294 266 end 295 267 NbSubDomain=size(Data.Civ1_SubRange,3); … … 316 288 mask=imread(par_civ2.Mask); 317 289 end 318 par_civ2.SearchBoxSize(1)=2*isx2+1; 319 par_civ2.SearchBoxSize(2)=2*isy2+1; 290 ibx2=ceil(par_civ2.CorrBoxSize(1)/2); 291 iby2=ceil(par_civ2.CorrBoxSize(2)/2); 292 % isx2=ibx2+4;% search ara +-4 pixels around the guess 293 % isy2=iby2+4; 294 par_civ2.SearchBoxSize(1)=2*ibx2+9;% search ara +-4 pixels around the guess 295 par_civ2.SearchBoxSize(2)=2*iby2+9; 296 %par_civ2.SearchBoxSize(1)=2*isx2+1; 297 %par_civ2.SearchBoxSize(2)=2*isy2+1; 320 298 par_civ2.SearchBoxShift=[Shiftx(nbval>=1)./nbval(nbval>=1) Shifty(nbval>=1)./nbval(nbval>=1)]; 321 % par_civ2.SearchBoxShift(2)=Shifty(nbval>=1)./nbval(nbval>=1);322 299 par_civ2.Grid=[GridX(nbval>=1)-par_civ2.SearchBoxShift(:,1)/2 GridY(nbval>=1)-par_civ2.SearchBoxShift(:,2)/2];% grid taken at the extrapolated origin of the displacement vectors 323 300 if par_civ2.CheckDeformation … … 329 306 % caluclate velocity data (y and v in indices, reverse to y component) 330 307 [xtable ytable utable vtable ctable F] = civ (par_civ2); 331 % diff_squared=(utable-par_civ2.Shiftx).*(utable-par_civ2.Shiftx)+(vtable+par_civ2.Shifty).*(vtable+par_civ2.Shifty);332 % F(diff_squared>=4)=4; %flag vectors whose distance to the guess exceeds 2 pixels333 308 list_param=(fieldnames(Param.Civ2))'; 334 309 list_remove={'pxcmx','pxcmy','npx','npy','gridflag','maskflag','term_a','term_b','T0'}; … … 350 325 end 351 326 Data.ListGlobalAttribute=[Data.ListGlobalAttribute Civ2_param {'Civ2_Time','Civ2_Dt'}]; 352 % Data.Civ2_Time=par_civ2.Time;353 % Data.Civ2_Dt=par_civ2.Dt;354 327 nbvar=numel(Data.ListVarName); 355 328 Data.ListVarName=[Data.ListVarName {'Civ2_X','Civ2_Y','Civ2_U','Civ2_V','Civ2_F','Civ2_C'}];% cell array containing the names of the fields to record … … 475 448 %methods perform better and will maybe be implemented in the future. 476 449 477 %% prepare grid 450 %% prepare measurement grid 451 if isfield(par_civ,'Grid')% grid points set as input 452 if ischar(par_civ.Grid)%read the drid file if the input is a file name 453 par_civ.Grid=dlmread(par_civ.Grid); 454 par_civ.Grid(1,:)=[];%the first line must be removed (heading in the grid file) 455 end 456 else% automatic grid 457 minix=floor(par_civ.Dx/2)-0.5; 458 maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx); 459 miniy=floor(par_civ.Dy/2)-0.5; 460 maxiy=minix+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy); 461 [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy); 462 par_civ.Grid(:,1)=reshape(GridX,[],1); 463 par_civ.Grid(:,2)=reshape(GridY,[],1); 464 end 465 nbvec=size(par_civ.Grid,1); 466 467 %% prepare correlation and search boxes 478 468 ibx2=ceil(par_civ.CorrBoxSize(1)/2); 479 469 iby2=ceil(par_civ.CorrBoxSize(2)/2); … … 482 472 shiftx=round(par_civ.SearchBoxShift(:,1)); 483 473 shifty=-round(par_civ.SearchBoxShift(:,2));% sign minus because image j index increases when y decreases 484 if isfield(par_civ,'Grid') 485 if ischar(par_civ.Grid)%read the drid file if the input is a file name 486 par_civ.Grid=dlmread(par_civ.Grid); 487 par_civ.Grid(1,:)=[];%the first line must be removed (heading in the grid file) 488 end 489 else% automatic measurement grid 490 % ibx2=ceil(par_civ.Bx/2); 491 % iby2=ceil(par_civ.By/2); 492 % isx2=ceil(par_civ.Searchx/2); 493 % isy2=ceil(par_civ.Searchy/2); 494 miniy=max(1+isy2+shifty,1+iby2); 495 minix=max(1+isx2-shiftx,1+ibx2); 496 maxiy=min(par_civ.ImageHeight-isy2+shifty,par_civ.ImageHeight-iby2); 497 maxix=min(par_civ.ImageWidth-isx2-shiftx,par_civ.ImageWidth-ibx2); 498 [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy); 499 par_civ.Grid(:,1)=reshape(GridX,[],1); 500 par_civ.Grid(:,2)=reshape(GridY,[],1); 501 end 502 nbvec=size(par_civ.Grid,1); 503 if numel(shiftx)==1 474 if numel(shiftx)==1% case of a unique shift for the whole field( civ1) 504 475 shiftx=shiftx*ones(nbvec,1); 505 476 shifty=shifty*ones(nbvec,1); 506 477 end 478 507 479 %% Default output 508 480 xtable=par_civ.Grid(:,1); … … 561 533 % 20>=mask: velocity=0 562 534 checkmask=0; 535 MinA=min(min(par_civ.ImageA)); 536 MinB=min(min(par_civ.ImageB)); 563 537 if isfield(par_civ,'Mask') && ~isempty(par_civ.Mask) 564 538 checkmask=1; … … 569 543 % check_noflux=(par_civ.Mask<100) ;%TODO: to implement 570 544 check_undefined=(par_civ.Mask<200 & par_civ.Mask>=20 ); 571 par_civ.ImageA(check_undefined)= min(min(par_civ.ImageA));% put image A to zero (i.e. the min image value) in the undefined area572 par_civ.ImageB(check_undefined)= min(min(par_civ.ImageB));% put image B to zero (i.e. the min image value) in the undefined area545 par_civ.ImageA(check_undefined)=MinA;% put image A to zero (i.e. the min image value) in the undefined area 546 par_civ.ImageB(check_undefined)=MinB;% put image B to zero (i.e. the min image value) in the undefined area 573 547 end 574 548 … … 587 561 jref=round(par_civ.ImageHeight-par_civ.Grid(ivec,2)+0.5);% yindex on the image B for the middle of the correlation box 588 562 if ~(checkmask && par_civ.Mask(jref,iref)<=20) %velocity not set to zero by the black mask 589 if jref-iby2<1 || jref+iby2>par_civ.ImageHeight|| iref-ibx2<1 || iref+ibx2>par_civ.ImageWidth||... 590 jref+shifty(ivec)-isy2<1||jref+shifty(ivec)+isy2>par_civ.ImageHeight|| iref+shiftx(ivec)-isx2<1 || iref+shiftx(ivec)+isx2>par_civ.ImageWidth % we are outside the image 591 F(ivec)=3; 592 else 593 image1_crop=par_civ.ImageA(jref-iby2:jref+iby2,iref-ibx2:iref+ibx2);%extract a subimage (correlation box) from image A 594 image2_crop=par_civ.ImageB(jref+shifty(ivec)-isy2:jref+shifty(ivec)+isy2,iref+shiftx(ivec)-isx2:iref+shiftx(ivec)+isx2);%extract a larger subimage (search box) from image B 563 % if jref-iby2<1 || jref+iby2>par_civ.ImageHeight|| iref-ibx2<1 || iref+ibx2>par_civ.ImageWidth||... 564 % jref+shifty(ivec)-isy2<1||jref+shifty(ivec)+isy2>par_civ.ImageHeight|| iref+shiftx(ivec)-isx2<1 || iref+shiftx(ivec)+isx2>par_civ.ImageWidth % we are outside the image 565 % F(ivec)=3; 566 % else 567 F(ivec)=0; 568 subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage 569 subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage 570 subrange2_x=iref+shiftx(ivec)-isx2:iref+shiftx(ivec)+isx2;%x indices defining the second subimage 571 subrange2_y=jref+shifty(ivec)-isy2:jref+shifty(ivec)+isy2;%y indices defining the second subimage 572 image1_crop=MinA*ones(numel(subrange1_y),numel(subrange1_x));% default value=min of image A 573 image2_crop=MinA*ones(numel(subrange2_y),numel(subrange2_x));% default value=min of image A 574 check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1 575 check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight; 576 check2_x=subrange2_x>=1 & subrange2_x<=par_civ.ImageWidth;% check which points in the subimage 2 are contained in the initial image 2 577 check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight; 578 image1_crop(check1_y,check1_x)=par_civ.ImageA(subrange1_y(check1_y),subrange1_x(check1_x));%extract a subimage (correlation box) from image A 579 image2_crop(check2_y,check2_x)=par_civ.ImageB(subrange2_y(check2_y),subrange2_x(check2_x));%extract a larger subimage (search box) from image B 595 580 image1_mean=mean(mean(image1_crop)); 596 581 image2_mean=mean(mean(image2_crop)); … … 603 588 F(ivec)=3; 604 589 end 605 end590 % end 606 591 if F(ivec)~=3 607 592 image1_crop=image1_crop-image1_mean;%substract the mean
Note: See TracChangeset
for help on using the changeset viewer.