Changeset 227 for trunk/src/pivlab.m
 Timestamp:
 Mar 31, 2011, 1:42:51 PM (13 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/pivlab.m
r225 r227 1 1 % 'pivlab': function piv.m adapted from PIVlab http://pivlab.blogspot.com/ 2 2 % 3 % function [xtable ytable utable vtable typevector] = pivlab (image1,image2,i nterrogationarea,step, subpixfinder, mask, roi)3 % function [xtable ytable utable vtable typevector] = pivlab (image1,image2,ibx,iby step, subpixfinder, mask, roi) 4 4 % 5 5 % OUTPUT: … … 14 14 % image1:first image (matrix) 15 15 % image2: second image (matrix) 16 % i nterrogationarea: size of the correlation box(in px)16 % ibx,iby: size of the correlation box along x and y (in px) 17 17 % step: mesh of the measurement points (in px) 18 18 % subpixfinder=1 or 2 controls the curve fitting of the image correlation 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 typevector] = pivlab (image1,image2,i nterrogationarea, step, subpixfinder, mask, roi)21 function [xtable ytable utable vtable ctable typevector] = pivlab (image1,image2,ibx2,iby2,isx2,isy2,shiftx,shifty, PointCoord, subpixfinder,mask) 22 22 %this funtion performs the DCC PIV analysis. Recent windowdeformation 23 23 %methods perform better and will maybe be implemented in the future. 24 24 warning off %MATLAB:log:logOfZero 25 if numel(roi)>026 xroi=roi(1);27 yroi=roi(2);28 widthroi=roi(3);29 heightroi=roi(4);30 image1_roi=double(image1(yroi:yroi+heightroi,xroi:xroi+widthroi));31 image2_roi=double(image2(yroi:yroi+heightroi,xroi:xroi+widthroi));32 else25 % if numel(roi)>0 26 % xroi=roi(1); 27 % yroi=roi(2); 28 % widthroi=roi(3); 29 % heightroi=roi(4); 30 % image1_roi=double(image1(yroi:yroi+heightroi,xroi:xroi+widthroi)); 31 % image2_roi=double(image2(yroi:yroi+heightroi,xroi:xroi+widthroi)); 32 % else 33 33 xroi=0; 34 34 yroi=0; 35 35 image1_roi=double(image1); 36 36 image2_roi=double(image2); 37 end37 % end 38 38 if numel(mask)>0 39 39 cellmask=mask; … … 48 48 end 49 49 mask(mask>1)=1; 50 51 miniy=1+(ceil(interrogationarea/2)); 52 minix=1+(ceil(interrogationarea/2)); 53 maxiy=step*(floor(size(image1_roi,1)/step))(interrogationarea1)+(ceil(interrogationarea/2)); %statt size deltax von ROI nehmen 54 maxix=step*(floor(size(image1_roi,2)/step))(interrogationarea1)+(ceil(interrogationarea/2)); 55 56 numelementsy=floor((maxiyminiy)/step+1); 57 numelementsx=floor((maxixminix)/step+1); 58 59 LAy=miniy; 60 LAx=minix; 61 LUy=size(image1_roi,1)maxiy; 62 LUx=size(image1_roi,2)maxix; 63 shift4centery=round((LUyLAy)/2); 64 shift4centerx=round((LUxLAx)/2); 65 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.... 66 shift4centery=0; 67 end 68 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.... 69 shift4centerx=0; 70 end 71 miniy=miniy+shift4centery; 72 minix=minix+shift4centerx; 73 maxix=maxix+shift4centerx; 74 maxiy=maxiy+shift4centery; 75 76 image1_roi=padarray(image1_roi,[ceil(interrogationarea/2) ceil(interrogationarea/2)], min(min(image1_roi))); 77 image2_roi=padarray(image2_roi,[ceil(interrogationarea/2) ceil(interrogationarea/2)], min(min(image1_roi))); 78 mask=padarray(mask,[ceil(interrogationarea/2) ceil(interrogationarea/2)],0); 79 if (rem(interrogationarea,2) == 0) %for the subpixel displacement measurement 80 SubPixOffset=1; 81 else 82 SubPixOffset=0.5; 83 end 84 xtable=zeros(numelementsy,numelementsx); 50 % ibx=2*ibx21%ibx and iby odd, reduced by 1 if even 51 % iby=2*iby21 52 % miniy=1+iby2 53 % minix=1+ibx2 54 % maxiy=step*(floor(size(image1_roi,1)/step))(iby1)+iby2 %statt size deltax von ROI nehmen 55 % maxix=step*(floor(size(image1_roi,2)/step))(ibx1)+ibx2 56 % numelementsy=floor((maxiyminiy)/step+1); 57 % numelementsx=floor((maxixminix)/step+1); 58 % 59 % LAy=miniy; 60 % LAx=minix; 61 % LUy=size(image1_roi,1)maxiy; 62 % LUx=size(image1_roi,2)maxix; 63 % shift4centery=round((LUyLAy)/2); 64 % shift4centerx=round((LUxLAx)/2); 65 % 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.... 66 % shift4centery=0; 67 % end 68 % 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.... 69 % shift4centerx=0; 70 % end 71 % miniy=miniy+shift4centery; 72 % minix=minix+shift4centerx; 73 % maxix=maxix+shift4centerx; 74 % maxiy=maxiy+shift4centery; 75 76 image1_roi=padarray(image1_roi,[iby2 ibx2], min(min(image1_roi)));%add a border around the image with minimum image value 77 image2_roi=padarray(image2_roi,[iby2 ibx2], min(min(image1_roi))); 78 mask=padarray(mask,[iby2 ibx2],0); 79 SubPixOffset=0.5;%odd values chosen for ibx and iby 80 81 nbvec=size(PointCoord,1); 82 xtable=zeros(nbvec,1); 85 83 ytable=xtable; 86 84 utable=xtable; … … 89 87 v2table=xtable; 90 88 s2n=xtable; 91 typevector=ones( numelementsy,numelementsx);;89 typevector=ones(size(xtable)); 92 90 93 91 nrx=0; … … 97 95 98 96 %% MAINLOOP 99 %handles=guihandles(getappdata(0,'hgui')); 100 for j = miniy:step:maxiy %vertical loop 101 nry=nry+1; 102 % if increments<6 %reduced display refreshing rate 103 % increments=increments+1; 104 % else 105 % increments=1; 106 % %set(handles.progress, 'string' , ['Frame progress: ' int2str(j/maxiy*100) '%']);drawnow; 107 % end 108 n=round((jminiy)/maxiy*100); 109 for i = minix:step:maxix % horizontal loop 110 nrx=nrx+1;%used to determine the pos of the vector in resulting matrix 111 if nrxreal < numelementsx 112 nrxreal=nrxreal+1; 113 else 114 nrxreal=1; 115 end 116 startpoint=[i j]; 117 image1_crop=image1_roi(j:j+interrogationarea1, i:i+interrogationarea1); 118 image2_crop=image2_roi(ceil(jinterrogationarea/2):ceil(j+1.5*interrogationarea1), ceil(iinterrogationarea/2):ceil(i+1.5*interrogationarea1)); 119 corrmax=0; 120 if mask(round(j+interrogationarea/2),round(i+interrogationarea/2))==0 97 for ivec=1:nbvec 98 iref=PointCoord(ivec,1); 99 jref=PointCoord(ivec,2); 100 image1_crop=image1_roi(jrefiby2:jref+iby2,irefibx2:iref+ibx2); 101 image2_crop=image2_roi(jref+shiftyisy2:jref+shifty+isy2,iref+shiftxisx2:iref+shiftx+isx2); 102 % n=round((jminiy)/maxiy*100); 103 % for i = minix:step:maxix % horizontal loop 104 % nrx=nrx+1;%used to determine the pos of the vector in resulting matrix 105 % if nrxreal < numelementsx 106 % nrxreal=nrxreal+1; 107 % else 108 % nrxreal=1; 109 % end 110 % startpoint=[i j]; 111 % image1_crop=image1_roi(j:j+iby1, i:i+ibx1); 112 % image2_crop=image2_roi(ceil(jiby/2):ceil(j+1.5*iby1), ceil(iibx/2):ceil(i+1.5*ibx1)); 113 % corrmax=0; 114 if mask(jref,iref)==0 121 115 %reference: Oliver Pust, PIV: Direct CrossCorrelation 122 116 % image2_crop: sub image with the size of the search area in image 2 … … 132 126 try 133 127 if subpixfinder==1 134 [vector] = SUBPIXGAUSS (result_conv,i nterrogationarea,x,y,SubPixOffset);128 [vector] = SUBPIXGAUSS (result_conv,ibx,iby,x,y,SubPixOffset); 135 129 elseif subpixfinder==2 136 [vector] = SUBPIX2DGAUSS (result_conv,i nterrogationarea,x,y,SubPixOffset);130 [vector] = SUBPIX2DGAUSS (result_conv,ibx,iby,x,y,SubPixOffset); 137 131 end 138 132 catch … … 144 138 else %if mask was not 0 then 145 139 vector=[NaN NaN]; 146 typevector( nry,nrxreal)=0;140 typevector(ivec)=0; 147 141 end 148 142 149 143 %Create the vector matrix x, y, u, v 150 xtable(nry,nrxreal)=startpoint(1)+interrogationarea/2; 151 ytable(nry,:)=startpoint(1,2)+interrogationarea/2; 152 utable(nry,nrxreal)=vector(1); 153 vtable(nry,nrxreal)=vector(2); 154 ctable(nry,nrxreal)=corrmax/sum(sum(image1_crop.*image1_crop)); 155 end 156 157 end 158 xtable=xtableceil(interrogationarea/2); 159 ytable=ytableceil(interrogationarea/2); 160 161 xtable=xtable+xroi; 162 ytable=ytable+yroi; 163 164 utable(utable>interrogationarea/1.5)=NaN; 165 vtable(utable>interrogationarea/1.5)=NaN; 166 vtable(vtable>interrogationarea/1.5)=NaN; 167 utable(vtable>interrogationarea/1.5)=NaN; 168 169 function [vector] = SUBPIXGAUSS (result_conv,interrogationarea,x,y,SubPixOffset) 144 xtable(ivec)=PointCoord(ivec,2); 145 ytable(ivec)=PointCoord(ivec,1); 146 utable(ivec)=vector(1); 147 vtable(ivec)=vector(2); 148 ctable(ivec)=corrmax/sum(sum(image1_crop.*image1_crop)); 149 end 150 % xtable=xtableibx2; 151 % ytable=ytableiby2; 152 % 153 % xtable=xtable+xroi; 154 % ytable=ytable+yroi; 155 % 156 % utable(utable>ibx/1.5)=NaN; 157 % vtable(utable>ibx/1.5)=NaN; 158 % vtable(vtable>iby/1.5)=NaN; 159 % utable(vtable>iby/1.5)=NaN; 160 161 function [vector] = SUBPIXGAUSS (result_conv,ibx,iby,x,y,SubPixOffset) 170 162 if size(x,1)>1 %if there are more than 1 peaks just take the first 171 163 x=x(1:1); … … 186 178 peakx = x+ (f1f2)/(2*f14*f0+2*f2); 187 179 % 188 SubpixelX=peakx(i nterrogationarea/2)SubPixOffset;189 SubpixelY=peaky(i nterrogationarea/2)SubPixOffset;180 SubpixelX=peakx(ibx/2)SubPixOffset; 181 SubpixelY=peaky(iby/2)SubPixOffset; 190 182 vector=[SubpixelX, SubpixelY]; 191 183 else … … 193 185 end 194 186 195 function [vector] = SUBPIX2DGAUSS (result_conv,i nterrogationarea,x,y,SubPixOffset)187 function [vector] = SUBPIX2DGAUSS (result_conv,ibx,iby,x,y,SubPixOffset) 196 188 if size(x,1)>1 %if there are more than 1 peaks just take the first 197 189 x=x(1:1); … … 229 221 peaky=y+deltay; 230 222 231 SubpixelX=peakx(i nterrogationarea/2)SubPixOffset;232 SubpixelY=peaky(i nterrogationarea/2)SubPixOffset;223 SubpixelX=peakx(ibx/2)SubPixOffset; 224 SubpixelY=peaky(iby/2)SubPixOffset; 233 225 vector=[SubpixelX, SubpixelY]; 234 226 else
Note: See TracChangeset
for help on using the changeset viewer.