source: trunk/src/pivlab.m @ 225

Last change on this file since 225 was 225, checked in by sommeria, 10 years ago

civ : bug in binary check repaired. Introduction of the option PIVlab for PIV (still under test)

File size: 9.5 KB
Line 
1% 'pivlab': function piv.m adapted from PIVlab http://pivlab.blogspot.com/
2%--------------------------------------------------------------------------
3% function [xtable ytable utable vtable typevector] = pivlab (image1,image2,interrogationarea, step, subpixfinder, mask, roi)
4%
5% OUTPUT:
6% xtable: set of x coordinates
7% ytable: set of y coordiantes
8% utable: set of u displacements (along x)
9% vtable: set of v displacements (along y)
10% ctable: max image correlation for each vector
11% typevector: set of flags, =1 for good, =0 for NaN vectors
12%
13%INPUT:
14% image1:first image (matrix)
15% image2: second image (matrix)
16% interrogationarea: size of the correlation box (in px)
17% step: mesh of the measurement points (in px)
18% subpixfinder=1 or 2 controls the curve fitting of the image correlation
19% mask: =[] for no mask
20% roi: 4 element vector defining a region of interest: x position, y position, width, height, (in image indices), for the whole image, roi=[];
21function [xtable ytable utable vtable ctable typevector] = pivlab (image1,image2,interrogationarea, step, subpixfinder, mask, roi)
22%this funtion performs the DCC PIV analysis. Recent window-deformation
23%methods perform better and will maybe be implemented in the future.
24warning off %MATLAB:log:logOfZero
25if 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));
32else
33    xroi=0;
34    yroi=0;
35    image1_roi=double(image1);
36    image2_roi=double(image2);
37end
38if numel(mask)>0
39    cellmask=mask;
40    mask=zeros(size(image1_roi));
41    for i=1:size(cellmask,1);
42        masklayerx=cellmask{i,1};
43        masklayery=cellmask{i,2};
44        mask = mask + poly2mask(masklayerx-xroi,masklayery-yroi,size(image1_roi,1),size(image1_roi,2)); %kleineres eingangsbild und maske geshiftet
45    end
46else
47    mask=zeros(size(image1_roi));
48end
49mask(mask>1)=1;
50
51miniy=1+(ceil(interrogationarea/2));
52minix=1+(ceil(interrogationarea/2));
53maxiy=step*(floor(size(image1_roi,1)/step))-(interrogationarea-1)+(ceil(interrogationarea/2)); %statt size deltax von ROI nehmen
54maxix=step*(floor(size(image1_roi,2)/step))-(interrogationarea-1)+(ceil(interrogationarea/2));
55
56numelementsy=floor((maxiy-miniy)/step+1);
57numelementsx=floor((maxix-minix)/step+1);
58
59LAy=miniy;
60LAx=minix;
61LUy=size(image1_roi,1)-maxiy;
62LUx=size(image1_roi,2)-maxix;
63shift4centery=round((LUy-LAy)/2);
64shift4centerx=round((LUx-LAx)/2);
65if 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;
67end
68if 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;
70end
71miniy=miniy+shift4centery;
72minix=minix+shift4centerx;
73maxix=maxix+shift4centerx;
74maxiy=maxiy+shift4centery;
75
76image1_roi=padarray(image1_roi,[ceil(interrogationarea/2) ceil(interrogationarea/2)], min(min(image1_roi)));
77image2_roi=padarray(image2_roi,[ceil(interrogationarea/2) ceil(interrogationarea/2)], min(min(image1_roi)));
78mask=padarray(mask,[ceil(interrogationarea/2) ceil(interrogationarea/2)],0);
79if (rem(interrogationarea,2) == 0) %for the subpixel displacement measurement
80    SubPixOffset=1;
81else
82    SubPixOffset=0.5;
83end
84xtable=zeros(numelementsy,numelementsx);
85ytable=xtable;
86utable=xtable;
87vtable=xtable;
88u2table=xtable;
89v2table=xtable;
90s2n=xtable;
91typevector=ones(numelementsy,numelementsx);;
92
93nrx=0;
94nrxreal=0;
95nry=0;
96increments=0;
97
98%% MAINLOOP
99%handles=guihandles(getappdata(0,'hgui'));
100for 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((j-miniy)/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+interrogationarea-1, i:i+interrogationarea-1);
118        image2_crop=image2_roi(ceil(j-interrogationarea/2):ceil(j+1.5*interrogationarea-1), ceil(i-interrogationarea/2):ceil(i+1.5*interrogationarea-1));
119        corrmax=0;
120        if mask(round(j+interrogationarea/2),round(i+interrogationarea/2))==0
121           %reference: Oliver Pust, PIV: Direct Cross-Correlation
122           % image2_crop: sub image with the size of the search area in image 2
123           % image1_crop: sub image of the correlation box in image 1
124           % %image2_crop is bigger than image1_crop. Zeropading is therefore not
125            result_conv= conv2(image2_crop,rot90(image1_crop,2),'valid');         
126            %necessary. 'Valid' makes sure that no zero padded content is returned.
127            corrmax= max(max(result_conv));
128            result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255
129            %Find the 255 peak
130            [y,x] = find(result_conv==255);
131            if isnan(y)==0 & isnan(x)==0
132                try
133                    if subpixfinder==1
134                        [vector] = SUBPIXGAUSS (result_conv,interrogationarea,x,y,SubPixOffset);
135                    elseif subpixfinder==2
136                        [vector] = SUBPIX2DGAUSS (result_conv,interrogationarea,x,y,SubPixOffset);
137                    end
138                catch
139                    vector=[NaN NaN]; %if something goes wrong with cross correlation.....
140                end
141            else
142                vector=[NaN NaN]; %if something goes wrong with cross correlation.....
143            end
144        else %if mask was not 0 then
145            vector=[NaN NaN];
146            typevector(nry,nrxreal)=0;
147        end
148
149        %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
157end
158xtable=xtable-ceil(interrogationarea/2);
159ytable=ytable-ceil(interrogationarea/2);
160
161xtable=xtable+xroi;
162ytable=ytable+yroi;
163
164utable(utable>interrogationarea/1.5)=NaN;
165vtable(utable>interrogationarea/1.5)=NaN;
166vtable(vtable>interrogationarea/1.5)=NaN;
167utable(vtable>interrogationarea/1.5)=NaN;
168
169function [vector] = SUBPIXGAUSS (result_conv,interrogationarea,x,y,SubPixOffset)
170if size(x,1)>1 %if there are more than 1 peaks just take the first
171    x=x(1:1);
172end
173if size(y,1)>1 %if there are more than 1 peaks just take the first
174    y=y(1:1);
175end
176if (x <= (size(result_conv,1)-1)) && (y <= (size(result_conv,1)-1)) && (x >= 1) && (y >= 1)
177    %the following 8 lines are copyright (c) 1998, Uri Shavit, Roi Gurka, Alex Liberzon, Technion – Israel Institute of Technology
178    %http://urapiv.wordpress.com
179    f0 = log(result_conv(y,x));
180    f1 = log(result_conv(y-1,x));
181    f2 = log(result_conv(y+1,x));
182    peaky = y+ (f1-f2)/(2*f1-4*f0+2*f2);
183    f0 = log(result_conv(y,x));
184    f1 = log(result_conv(y,x-1));
185    f2 = log(result_conv(y,x+1));
186    peakx = x+ (f1-f2)/(2*f1-4*f0+2*f2);
187    %
188    SubpixelX=peakx-(interrogationarea/2)-SubPixOffset;
189    SubpixelY=peaky-(interrogationarea/2)-SubPixOffset;
190    vector=[SubpixelX, SubpixelY];
191else
192    vector=[NaN NaN];
193end
194
195function [vector] = SUBPIX2DGAUSS (result_conv,interrogationarea,x,y,SubPixOffset)
196if size(x,1)>1 %if there are more than 1 peaks just take the first
197    x=x(1:1);
198end
199if size(y,1)>1 %if there are more than 1 peaks just take the first
200    y=y(1:1);
201end
202if (x <= (size(result_conv,1)-1)) && (y <= (size(result_conv,1)-1)) && (x >= 1) && (y >= 1)
203    for i=-1:1
204        for j=-1:1
205            %following 15 lines based on
206            %H. Nobach Æ M. Honkanen (2005)
207            %Two-dimensional Gaussian regression for sub-pixel displacement
208            %estimation in particle image velocimetry or particle position
209            %estimation in particle tracking velocimetry
210            %Experiments in Fluids (2005) 38: 511–515
211            c10(j+2,i+2)=i*log(result_conv(y+j, x+i));
212            c01(j+2,i+2)=j*log(result_conv(y+j, x+i));
213            c11(j+2,i+2)=i*j*log(result_conv(y+j, x+i));
214            c20(j+2,i+2)=(3*i^2-2)*log(result_conv(y+j, x+i));
215            c02(j+2,i+2)=(3*j^2-2)*log(result_conv(y+j, x+i));
216            %c00(j+2,i+2)=(5-3*i^2-3*j^2)*log(result_conv_norm(maxY+j, maxX+i));
217        end
218    end
219    c10=(1/6)*sum(sum(c10));
220    c01=(1/6)*sum(sum(c01));
221    c11=(1/4)*sum(sum(c11));
222    c20=(1/6)*sum(sum(c20));
223    c02=(1/6)*sum(sum(c02));
224    %c00=(1/9)*sum(sum(c00));
225
226    deltax=(c11*c01-2*c10*c02)/(4*c20*c02-c11^2);
227    deltay=(c11*c10-2*c01*c20)/(4*c20*c02-c11^2);
228    peakx=x+deltax;
229    peaky=y+deltay;
230
231    SubpixelX=peakx-(interrogationarea/2)-SubPixOffset;
232    SubpixelY=peaky-(interrogationarea/2)-SubPixOffset;
233    vector=[SubpixelX, SubpixelY];
234else
235    vector=[NaN NaN];
236end
Note: See TracBrowser for help on using the repository browser.