# source:trunk/src/pivlab.m@227

Last change on this file since 227 was 227, checked in by sommeria, 13 years ago

add function sub_field_series to apply the sub_field operation to a series of fileds (for instance subtracting a background to an image series)

File size: 8.9 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,ibx,iby 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% ibx,iby: size of the correlation box along x and y (in px)
17% step: mesh of the measurement points (in px)
18% subpixfinder=1 or 2 controls the curve fitting of the image correlation
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,ibx2,iby2,isx2,isy2,shiftx,shifty, PointCoord, subpixfinder,mask)
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
25% 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    xroi=0;
34    yroi=0;
35    image1_roi=double(image1);
36    image2_roi=double(image2);
37% end
45    end
46else
48end
50% ibx=2*ibx2-1%ibx and iby odd, reduced by 1 if even
51% iby=2*iby2-1
52% miniy=1+iby2
53% minix=1+ibx2
54% maxiy=step*(floor(size(image1_roi,1)/step))-(iby-1)+iby2 %statt size deltax von ROI nehmen
55% maxix=step*(floor(size(image1_roi,2)/step))-(ibx-1)+ibx2
56% numelementsy=floor((maxiy-miniy)/step+1);
57% numelementsx=floor((maxix-minix)/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((LUy-LAy)/2);
64% shift4centerx=round((LUx-LAx)/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
76image1_roi=padarray(image1_roi,[iby2 ibx2], min(min(image1_roi)));%add a border around the image with minimum image value
79SubPixOffset=0.5;%odd values chosen for ibx and iby
80
81nbvec=size(PointCoord,1);
82xtable=zeros(nbvec,1);
83ytable=xtable;
84utable=xtable;
85vtable=xtable;
86u2table=xtable;
87v2table=xtable;
88s2n=xtable;
89typevector=ones(size(xtable));
90
91nrx=0;
92nrxreal=0;
93nry=0;
94increments=0;
95
96%% MAINLOOP
97for ivec=1:nbvec
98    iref=PointCoord(ivec,1);
99    jref=PointCoord(ivec,2);
100    image1_crop=image1_roi(jref-iby2:jref+iby2,iref-ibx2:iref+ibx2);
101    image2_crop=image2_roi(jref+shifty-isy2:jref+shifty+isy2,iref+shiftx-isx2:iref+shiftx+isx2);
102%     n=round((j-miniy)/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+iby-1, i:i+ibx-1);
112%         image2_crop=image2_roi(ceil(j-iby/2):ceil(j+1.5*iby-1), ceil(i-ibx/2):ceil(i+1.5*ibx-1));
113%         corrmax=0;
115           %reference: Oliver Pust, PIV: Direct Cross-Correlation
116           % image2_crop: sub image with the size of the search area in image 2
117           % image1_crop: sub image of the correlation box in image 1
118           % %image2_crop is bigger than image1_crop. Zeropading is therefore not
119            result_conv= conv2(image2_crop,rot90(image1_crop,2),'valid');
120            %necessary. 'Valid' makes sure that no zero padded content is returned.
121            corrmax= max(max(result_conv));
122            result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255
123            %Find the 255 peak
124            [y,x] = find(result_conv==255);
125            if isnan(y)==0 & isnan(x)==0
126                try
127                    if subpixfinder==1
128                        [vector] = SUBPIXGAUSS (result_conv,ibx,iby,x,y,SubPixOffset);
129                    elseif subpixfinder==2
130                        [vector] = SUBPIX2DGAUSS (result_conv,ibx,iby,x,y,SubPixOffset);
131                    end
132                catch
133                    vector=[NaN NaN]; %if something goes wrong with cross correlation.....
134                end
135            else
136                vector=[NaN NaN]; %if something goes wrong with cross correlation.....
137            end
138        else %if mask was not 0 then
139            vector=[NaN NaN];
140            typevector(ivec)=0;
141        end
142
143        %Create the vector matrix x, y, u, v
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));
149end
150% xtable=xtable-ibx2;
151% ytable=ytable-iby2;
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
161function [vector] = SUBPIXGAUSS (result_conv,ibx,iby,x,y,SubPixOffset)
162if size(x,1)>1 %if there are more than 1 peaks just take the first
163    x=x(1:1);
164end
165if size(y,1)>1 %if there are more than 1 peaks just take the first
166    y=y(1:1);
167end
168if (x <= (size(result_conv,1)-1)) && (y <= (size(result_conv,1)-1)) && (x >= 1) && (y >= 1)
169    %the following 8 lines are copyright (c) 1998, Uri Shavit, Roi Gurka, Alex Liberzon, Technion  Israel Institute of Technology
170    %http://urapiv.wordpress.com
171    f0 = log(result_conv(y,x));
172    f1 = log(result_conv(y-1,x));
173    f2 = log(result_conv(y+1,x));
174    peaky = y+ (f1-f2)/(2*f1-4*f0+2*f2);
175    f0 = log(result_conv(y,x));
176    f1 = log(result_conv(y,x-1));
177    f2 = log(result_conv(y,x+1));
178    peakx = x+ (f1-f2)/(2*f1-4*f0+2*f2);
179    %
180    SubpixelX=peakx-(ibx/2)-SubPixOffset;
181    SubpixelY=peaky-(iby/2)-SubPixOffset;
182    vector=[SubpixelX, SubpixelY];
183else
184    vector=[NaN NaN];
185end
186
187function [vector] = SUBPIX2DGAUSS (result_conv,ibx,iby,x,y,SubPixOffset)
188if size(x,1)>1 %if there are more than 1 peaks just take the first
189    x=x(1:1);
190end
191if size(y,1)>1 %if there are more than 1 peaks just take the first
192    y=y(1:1);
193end
194if (x <= (size(result_conv,1)-1)) && (y <= (size(result_conv,1)-1)) && (x >= 1) && (y >= 1)
195    for i=-1:1
196        for j=-1:1
197            %following 15 lines based on
198            %H. Nobach Æ M. Honkanen (2005)
199            %Two-dimensional Gaussian regression for sub-pixel displacement
200            %estimation in particle image velocimetry or particle position
201            %estimation in particle tracking velocimetry
202            %Experiments in Fluids (2005) 38: 511515
203            c10(j+2,i+2)=i*log(result_conv(y+j, x+i));
204            c01(j+2,i+2)=j*log(result_conv(y+j, x+i));
205            c11(j+2,i+2)=i*j*log(result_conv(y+j, x+i));
206            c20(j+2,i+2)=(3*i^2-2)*log(result_conv(y+j, x+i));
207            c02(j+2,i+2)=(3*j^2-2)*log(result_conv(y+j, x+i));
208            %c00(j+2,i+2)=(5-3*i^2-3*j^2)*log(result_conv_norm(maxY+j, maxX+i));
209        end
210    end
211    c10=(1/6)*sum(sum(c10));
212    c01=(1/6)*sum(sum(c01));
213    c11=(1/4)*sum(sum(c11));
214    c20=(1/6)*sum(sum(c20));
215    c02=(1/6)*sum(sum(c02));
216    %c00=(1/9)*sum(sum(c00));
217
218    deltax=(c11*c01-2*c10*c02)/(4*c20*c02-c11^2);
219    deltay=(c11*c10-2*c01*c20)/(4*c20*c02-c11^2);
220    peakx=x+deltax;
221    peaky=y+deltay;
222
223    SubpixelX=peakx-(ibx/2)-SubPixOffset;
224    SubpixelY=peaky-(iby/2)-SubPixOffset;
225    vector=[SubpixelX, SubpixelY];
226else
227    vector=[NaN NaN];
228end
Note: See TracBrowser for help on using the repository browser.