Changeset 236 for trunk/src/pivlab.m


Ignore:
Timestamp:
Apr 12, 2011, 12:12:19 AM (13 years ago)
Author:
sommeria
Message:

correct Matlab PIV, remove call to image tool box. Improve menu of uvmat VelType? (replacement of buttons)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/pivlab.m

    r233 r236  
    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 result_conv errormsg] = pivlab (image1,image2,ibx2,iby2,isx2,isy2,shiftx,shifty, GridIndices, subpixfinder,mask)
     21function [xtable ytable utable vtable ctable F result_conv errormsg] = pivlab (image1,image2,ibx2,iby2,isx2,isy2,shiftx,shifty, GridIndices, 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.
    24 errormsg='';
    25 warning off %MATLAB:log:logOfZero
    26 [npy_ima npx_ima]=size(image1);
    27 if ~isequal(size(image1),size(image2))
    28     errormsg='image pair with unequal size';
    29     return
    30 end
    31     xroi=0;
    32     yroi=0;
    33     image1_roi=double(image1);
    34     image2_roi=double(image2);
    35 if numel(mask)>0
    36     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 geshiftet
    42     end
    43 else
    44     mask=zeros(size(image1));
    45 end
    46 mask(mask>1)=1;
    47 
    48 % ibx=2*ibx2-1;%ibx and iby odd, reduced by 1 if even
    49 % iby=2*iby2-1;
    50 % miniy=1+iby2
    51 % minix=1+ibx2
    52 % maxiy=step*(floor(size(image1_roi,1)/step))-(iby-1)+iby2 %statt size deltax von ROI nehmen
    53 % maxix=step*(floor(size(image1_roi,2)/step))-(ibx-1)+ibx2
    54 % 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 % end
    66 % 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 % end
    69 % 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 value
    75 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 iby
    78 
    7924nbvec=size(GridIndices,1);
    8025xtable=zeros(nbvec,1);
    … …  
    8227utable=xtable;
    8328vtable=xtable;
    84 u2table=xtable;
    85 v2table=xtable;
    86 s2n=xtable;
    87 typevector=ones(size(xtable));
     29ctable=xtable;
     30F=xtable;
     31result_conv=[];
     32errormsg='';
     33%warning off %MATLAB:log:logOfZero
     34[npy_ima npx_ima]=size(image1);
     35if ~isequal(size(image2),[npy_ima npx_ima])
     36    errormsg='image pair with unequal size';
     37    return
     38end
    8839
    89 nrx=0;
    90 nrxreal=0;
    91 nry=0;
    92 increments=0;
     40%% mask
     41testmask=0;
     42if 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
     58end
     59image1=double(image1);
     60image2=double(image2);
    9361
    94 %% MAINLOOP
     62%% calculate correlations: MAINLOOP
     63corrmax=0;
     64sum_square=1;% default
    9565for ivec=1:nbvec
    9666    iref=GridIndices(ivec,1);
    9767    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);
    12599                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
    127106                vector=[0 0]; %if something goes wrong with cross correlation.....
     107                F(ivec)=3;
    128108            end
    129         else %if mask was not 0 then
    130             vector=[0 0];
    131             typevector(ivec)=0;
     109        else
     110            vector=[0 0]; %if something goes wrong with cross correlation.....
     111            F(ivec)=3;
    132112        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);
    141123end
    142 result_conv=result_conv/255;
     124result_conv=result_conv*corrmax/(255*sum_square);% keep the last correlation matrix for output
    143125
    144126
Note: See TracChangeset for help on using the changeset viewer.