Changeset 236 for trunk/src/pivlab.m
- Timestamp:
- Apr 12, 2011, 12:12:19 AM (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/pivlab.m
r233 r236 19 19 % mask: =[] for no mask 20 20 % roi: 4 element vector defining a region of interest: x position, y position, width, height, (in image indices), for the whole image, roi=[]; 21 function [xtable ytable utable vtable ctable typevectorresult_conv errormsg] = pivlab (image1,image2,ibx2,iby2,isx2,isy2,shiftx,shifty, GridIndices, subpixfinder,mask)21 function [xtable ytable utable vtable ctable F result_conv errormsg] = pivlab (image1,image2,ibx2,iby2,isx2,isy2,shiftx,shifty, GridIndices, subpixfinder,mask) 22 22 %this funtion performs the DCC PIV analysis. Recent window-deformation 23 23 %methods perform better and will maybe be implemented in the future. 24 errormsg='';25 warning off %MATLAB:log:logOfZero26 [npy_ima npx_ima]=size(image1);27 if ~isequal(size(image1),size(image2))28 errormsg='image pair with unequal size';29 return30 end31 xroi=0;32 yroi=0;33 image1_roi=double(image1);34 image2_roi=double(image2);35 if numel(mask)>036 cellmask=mask;37 mask=zeros(size(image1));38 for i=1:size(cellmask,1);39 masklayerx=cellmask{i,1};40 masklayery=cellmask{i,2};41 mask = mask + poly2mask(masklayerx-xroi,masklayery-yroi,npy_ima,npx_ima); %kleineres eingangsbild und maske geshiftet42 end43 else44 mask=zeros(size(image1));45 end46 mask(mask>1)=1;47 48 % ibx=2*ibx2-1;%ibx and iby odd, reduced by 1 if even49 % iby=2*iby2-1;50 % miniy=1+iby251 % minix=1+ibx252 % maxiy=step*(floor(size(image1_roi,1)/step))-(iby-1)+iby2 %statt size deltax von ROI nehmen53 % maxix=step*(floor(size(image1_roi,2)/step))-(ibx-1)+ibx254 % numelementsy=floor((maxiy-miniy)/step+1);55 % numelementsx=floor((maxix-minix)/step+1);56 %57 % LAy=miniy;58 % LAx=minix;59 % LUy=size(image1_roi,1)-maxiy;60 % LUx=size(image1_roi,2)-maxix;61 % shift4centery=round((LUy-LAy)/2);62 % shift4centerx=round((LUx-LAx)/2);63 % if shift4centery<0 %shift4center will be negative if in the unshifted case the left border is bigger than the right border. the vectormatrix is hence not centered on the image. the matrix cannot be shifted more towards the left border because then image2_crop would have a negative index. The only way to center the matrix would be to remove a column of vectors on the right side. but then we weould have less data....64 % shift4centery=0;65 % end66 % if shift4centerx<0 %shift4center will be negative if in the unshifted case the left border is bigger than the right border. the vectormatrix is hence not centered on the image. the matrix cannot be shifted more towards the left border because then image2_crop would have a negative index. The only way to center the matrix would be to remove a column of vectors on the right side. but then we weould have less data....67 % shift4centerx=0;68 % end69 % miniy=miniy+shift4centery;70 % minix=minix+shift4centerx;71 % maxix=maxix+shift4centerx;72 % maxiy=maxiy+shift4centery;73 74 image1_roi=padarray(image1_roi,[iby2 ibx2], min(min(image1_roi)));%add a border around the image with minimum image value75 image2_roi=padarray(image2_roi,[iby2 ibx2], min(min(image1_roi)));76 mask=padarray(mask,[iby2 ibx2],0);77 %SubPixOffset=0.5;%odd values chosen for ibx and iby78 79 24 nbvec=size(GridIndices,1); 80 25 xtable=zeros(nbvec,1); … … 82 27 utable=xtable; 83 28 vtable=xtable; 84 u2table=xtable; 85 v2table=xtable; 86 s2n=xtable; 87 typevector=ones(size(xtable)); 29 ctable=xtable; 30 F=xtable; 31 result_conv=[]; 32 errormsg=''; 33 %warning off %MATLAB:log:logOfZero 34 [npy_ima npx_ima]=size(image1); 35 if ~isequal(size(image2),[npy_ima npx_ima]) 36 errormsg='image pair with unequal size'; 37 return 38 end 88 39 89 nrx=0; 90 nrxreal=0; 91 nry=0; 92 increments=0; 40 %% mask 41 testmask=0; 42 if exist('mask','var') && ~isempty(mask) 43 testmask=1; 44 if ~isequal(size(mask),[npy_ima npx_ima]) 45 errormsg='mask must be an image with the same size as the images'; 46 return 47 end 48 % Convention for mask 49 % mask >200 : velocity calculated 50 % 200 >=mask>150;velocity not calculated, interpolation allowed (bad spots) 51 % 150>=mask >100: velocity not calculated, nor interpolated 52 % 100>=mask> 20: velocity not calculated, impermeable (no flux through mask boundaries) 53 % 20>=mask: velocity=0 54 test_noflux=(mask<=100) ; 55 test_undefined=(mask<=200 & mask>100 ); 56 image1(test_undefined)=min(min(image1))*ones(size(image1));% put image to zero in the undefined area 57 image2(test_undefined)=min(min(image1))*ones(size(image1));% put image to zero in the undefined area 58 end 59 image1=double(image1); 60 image2=double(image2); 93 61 94 %% MAINLOOP 62 %% calculate correlations: MAINLOOP 63 corrmax=0; 64 sum_square=1;% default 95 65 for ivec=1:nbvec 96 66 iref=GridIndices(ivec,1); 97 67 jref=GridIndices(ivec,2); 98 %jref=npy_ima-PointCoord(ivec,2)+1; 99 image1_crop=image1_roi(jref-iby2:jref+iby2,iref-ibx2:iref+ibx2); 100 image2_crop=image2_roi(jref+shifty-isy2:jref+shifty+isy2,iref+shiftx-isx2:iref+shiftx+isx2); 101 image1_crop=image1_crop-mean(mean(image1_crop)); 102 image2_crop=image2_crop-mean(mean(image2_crop)); 103 if mask(jref,iref)==0 104 %reference: Oliver Pust, PIV: Direct Cross-Correlation 105 % image2_crop: sub image with the size of the search area in image 2 106 % image1_crop: sub image of the correlation box in image 1 107 % %image2_crop is bigger than image1_crop. Zeropading is therefore not 108 result_conv= conv2(image2_crop,rot90(image1_crop,2),'valid'); 109 %necessary. 'Valid' makes sure that no zero padded content is returned. 110 corrmax= max(max(result_conv)); 111 result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255 112 % result_conv=flipdim(result_conv,2);%reverse x direction 113 %Find the 255 peak 114 [y,x] = find(result_conv==255); 115 if isnan(y)==0 & isnan(x)==0 116 try 117 if subpixfinder==1 118 [vector] = SUBPIXGAUSS (result_conv,x,y); 119 elseif subpixfinder==2 120 [vector] = SUBPIX2DGAUSS (result_conv,x,y); 121 end 122 catch ME 123 errormsg=ME.message 124 vector=[0 0]; %if something goes wrong with cross correlation..... 68 testmask_ij=0; 69 test0=0; 70 if testmask 71 if mask(jref,iref)<=20 72 vector=[0 0]; 73 test0=1; 74 else 75 mask_crop1=mask(jref-iby2:jref+iby2,iref-ibx2:iref+ibx2); 76 mask_crop2=mask(jref+shifty-isy2:jref+shifty+isy2,iref+shiftx-isx2:iref+shiftx+isx2); 77 if ~isempty(find(mask_crop1<=200 & mask_crop1>100,1)) || ~isempty(find(mask_crop2<=200 & mask_crop2>100,1)); 78 testmask_ij=1; 79 end 80 end 81 end 82 if ~test0 83 image1_crop=image1(jref-iby2:jref+iby2,iref-ibx2:iref+ibx2); 84 image2_crop=image2(jref+shifty-isy2:jref+shifty+isy2,iref+shiftx-isx2:iref+shiftx+isx2); 85 image1_crop=image1_crop-mean(mean(image1_crop)); 86 image2_crop=image2_crop-mean(mean(image2_crop)); 87 %reference: Oliver Pust, PIV: Direct Cross-Correlation 88 result_conv= conv2(image2_crop,flipdim(flipdim(image1_crop,2),1),'valid'); 89 corrmax= max(max(result_conv)); 90 result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255 91 %Find the correlation max, at 255 92 [y,x] = find(result_conv==255); 93 if ~isnan(y) & ~isnan(x) 94 try 95 if subpixfinder==1 96 [vector] = SUBPIXGAUSS (result_conv,x,y); 97 elseif subpixfinder==2 98 [vector] = SUBPIX2DGAUSS (result_conv,x,y); 125 99 end 126 else 100 sum_square=sum(sum(image1_crop.*image1_crop)); 101 ctable(ivec)=corrmax/sum_square;% correlation value 102 if vector(1)>shiftx+isx2-ibx2+subpixfinder || vector(2)>shifty+isy2-iby2+subpixfinder 103 F(ivec)=-2;%vector reaches the border of the search zone 104 end 105 catch ME 127 106 vector=[0 0]; %if something goes wrong with cross correlation..... 107 F(ivec)=3; 128 108 end 129 else %if mask was not 0 then130 vector=[0 0]; 131 typevector(ivec)=0;109 else 110 vector=[0 0]; %if something goes wrong with cross correlation..... 111 F(ivec)=3; 132 112 end 133 134 %Create the vector matrix x, y, u, v 135 xtable(ivec)=GridIndices(ivec,1); 136 ytable(ivec)=GridIndices(ivec,2); 137 utable(ivec)=vector(1); 138 vtable(ivec)=vector(2); 139 sum_square=sum(sum(image1_crop.*image1_crop)); 140 ctable(ivec)=corrmax/sum_square; 113 if testmask_ij 114 F(ivec)=3; 115 end 116 end 117 118 %Create the vector matrix x, y, u, v 119 xtable(ivec)=iref+vector(1)/2;% convec flow (velocity taken at the point middle from imgae1 and 2) 120 ytable(ivec)=jref+vector(2)/2; 121 utable(ivec)=vector(1); 122 vtable(ivec)=vector(2); 141 123 end 142 result_conv=result_conv /255;124 result_conv=result_conv*corrmax/(255*sum_square);% keep the last correlation matrix for output 143 125 144 126
Note: See TracChangeset
for help on using the changeset viewer.