Changeset 227 for trunk/src/pivlab.m


Ignore:
Timestamp:
Mar 31, 2011, 1:42:51 PM (13 years ago)
Author:
sommeria
Message:

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:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/pivlab.m

    r225 r227  
    11% 'pivlab': function piv.m adapted from PIVlab http://pivlab.blogspot.com/
    22%--------------------------------------------------------------------------
    3 % function [xtable ytable utable vtable typevector] = pivlab (image1,image2,interrogationarea, step, subpixfinder, mask, roi)
     3% function [xtable ytable utable vtable typevector] = pivlab (image1,image2,ibx,iby step, subpixfinder, mask, roi)
    44%
    55% OUTPUT:
    … …  
    1414% image1:first image (matrix)
    1515% image2: second image (matrix)
    16 % interrogationarea: size of the correlation box (in px)
     16% ibx,iby: size of the correlation box along x and y (in px)
    1717% step: mesh of the measurement points (in px)
    1818% subpixfinder=1 or 2 controls the curve fitting of the image correlation
    1919% mask: =[] for no mask
    2020% 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,interrogationarea, step, subpixfinder, mask, roi)
     21function [xtable ytable utable vtable ctable typevector] = pivlab (image1,image2,ibx2,iby2,isx2,isy2,shiftx,shifty, PointCoord, subpixfinder,mask)
    2222%this funtion performs the DCC PIV analysis. Recent window-deformation
    2323%methods perform better and will maybe be implemented in the future.
    2424warning 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
     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
    3333    xroi=0;
    3434    yroi=0;
    3535    image1_roi=double(image1);
    3636    image2_roi=double(image2);
    37 end
     37% end
    3838if numel(mask)>0
    3939    cellmask=mask;
    … …  
    4848end
    4949mask(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))-(interrogationarea-1)+(ceil(interrogationarea/2)); %statt size deltax von ROI nehmen
    54 maxix=step*(floor(size(image1_roi,2)/step))-(interrogationarea-1)+(ceil(interrogationarea/2));
    55 
    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 
    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*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
     77image2_roi=padarray(image2_roi,[iby2 ibx2], min(min(image1_roi)));
     78mask=padarray(mask,[iby2 ibx2],0);
     79SubPixOffset=0.5;%odd values chosen for ibx and iby
     80
     81nbvec=size(PointCoord,1);
     82xtable=zeros(nbvec,1);
    8583ytable=xtable;
    8684utable=xtable;
    … …  
    8987v2table=xtable;
    9088s2n=xtable;
    91 typevector=ones(numelementsy,numelementsx);;
     89typevector=ones(size(xtable));
    9290
    9391nrx=0;
    … …  
    9795
    9896%% 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((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
     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;
     114        if mask(jref,iref)==0
    121115           %reference: Oliver Pust, PIV: Direct Cross-Correlation
    122116           % image2_crop: sub image with the size of the search area in image 2
    … …  
    132126                try
    133127                    if subpixfinder==1
    134                         [vector] = SUBPIXGAUSS (result_conv,interrogationarea,x,y,SubPixOffset);
     128                        [vector] = SUBPIXGAUSS (result_conv,ibx,iby,x,y,SubPixOffset);
    135129                    elseif subpixfinder==2
    136                         [vector] = SUBPIX2DGAUSS (result_conv,interrogationarea,x,y,SubPixOffset);
     130                        [vector] = SUBPIX2DGAUSS (result_conv,ibx,iby,x,y,SubPixOffset);
    137131                    end
    138132                catch
    … …  
    144138        else %if mask was not 0 then
    145139            vector=[NaN NaN];
    146             typevector(nry,nrxreal)=0;
     140            typevector(ivec)=0;
    147141        end
    148142
    149143        %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=xtable-ceil(interrogationarea/2);
    159 ytable=ytable-ceil(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));
     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)
    170162if size(x,1)>1 %if there are more than 1 peaks just take the first
    171163    x=x(1:1);
    … …  
    186178    peakx = x+ (f1-f2)/(2*f1-4*f0+2*f2);
    187179    %
    188     SubpixelX=peakx-(interrogationarea/2)-SubPixOffset;
    189     SubpixelY=peaky-(interrogationarea/2)-SubPixOffset;
     180    SubpixelX=peakx-(ibx/2)-SubPixOffset;
     181    SubpixelY=peaky-(iby/2)-SubPixOffset;
    190182    vector=[SubpixelX, SubpixelY];
    191183else
    … …  
    193185end
    194186
    195 function [vector] = SUBPIX2DGAUSS (result_conv,interrogationarea,x,y,SubPixOffset)
     187function [vector] = SUBPIX2DGAUSS (result_conv,ibx,iby,x,y,SubPixOffset)
    196188if size(x,1)>1 %if there are more than 1 peaks just take the first
    197189    x=x(1:1);
    … …  
    229221    peaky=y+deltay;
    230222
    231     SubpixelX=peakx-(interrogationarea/2)-SubPixOffset;
    232     SubpixelY=peaky-(interrogationarea/2)-SubPixOffset;
     223    SubpixelX=peakx-(ibx/2)-SubPixOffset;
     224    SubpixelY=peaky-(iby/2)-SubPixOffset;
    233225    vector=[SubpixelX, SubpixelY];
    234226else
Note: See TracChangeset for help on using the changeset viewer.