Changeset 1160 for trunk


Ignore:
Timestamp:
Jul 16, 2024, 9:57:54 AM (5 months ago)
Author:
sommeria
Message:

script_delete_PCO added

Location:
trunk/src
Files:
1 added
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/script_delete_rdvision.m

    r1159 r1160  
    8585                if exist(ExtractFolder,'dir') && exist(PngFolder,'dir')
    8686                    filename_seq=fullfile(ExtractFolder,'im.seq');
     87                    NumberOfFrames=0;
    8788                    try
    88                     s=ini2struct(filename_seq);
    89                     FileInfo=s.sequenceSettings;
    90                     if isfield(s.sequenceSettings,'numberoffiles')
    91                         NumberOfFrames=str2double(s.sequenceSettings.numberoffiles);
    92                     else
    93                         status='bad seq file';
    94                     end
     89                        s=ini2struct(filename_seq);
     90                        FileInfo=s.sequenceSettings;
     91
     92                        if isfield(s.sequenceSettings,'numberoffiles')
     93                            NumberOfFrames=str2double(s.sequenceSettings.numberoffiles);
     94                        else
     95                            status='bad seq file';
     96                        end
    9597                    catch ME
    9698                        disp(['error in ' filename_seq])
    9799                    end
    98                     DirPng=dir(PngFolder);
     100                    DirPng=dir(PngFolder);         
     101
    99102                    if numel(DirPng)==NumberOfFrames+2
    100                         Checkdelete=1;
     103                        PngCells=struct2cell(DirPng);
     104                        mm=cell2mat(PngCells(4,:));% check the sizes of extracted images
     105                        sizemax=max(mm);
     106                        sizemin=min(mm(3:end));
     107                        disp(['max size(Mbytes)=' num2str(sizemax/1000000)]);
     108                        if min(mm(3:end))<0.9*sizemax
     109                            status=['WARNING' 'min size(Mbytes)=' num2str(sizemin/1000000)];
     110                        else
     111                            Checkdelete=1;% approve deletion of the source bin files
     112                        end
    101113                    else
    102114                        status=['extraction not finished,' num2str(numel(DirPng)-2) ' images extracted'];
    … …  
    141153    end
    142154end
     155'END'

  • trunk/src/series/civ_3D.m

    r1158 r1160  
    1616%           if absent, the fct looks for input data in Param.ActionInput     (test mode)
    1717%     Param.OutputSubDir: sets the folder name of output file(s,
    18 %           if absent no file is produced, result in the output structure Data (test mode)
    1918%     Param.ActionInput: substructure with the parameters provided by the GUI civ_input
    2019%                      .Civ1: parameters for civ1cc
    … …  
    185184Data.VarAttribute{1}.Role='coord_z';
    186185Data.VarAttribute{2}.Role='coord_x';
    187 ata.VarAttribute{3}.Role='coord_y';
     186Data.VarAttribute{3}.Role='coord_y';
    188187Data.VarAttribute{4}.Role='vector_x';
    189188Data.VarAttribute{5}.Role='vector_y';
    … …  
    310309        % calculate velocity data at the first position z, corresponding to j_index=par_civ1.Dz
    311310        [Data.Civ1_X(1,:,:),Data.Civ1_Y(1,:,:), Data.Civ1_U(1,:,:), Data.Civ1_V(1,:,:),Data.Civ1_W(1,:,:),...
    312             Data.Civ1_C(1,:,:), Data.Civ1_FF(1,:,:), result_conv, errormsg] = civ3D (par_civ1);
     311            Data.Civ1_C(1,:,:), Data.Civ1_FF(1,:,:), ~, errormsg] = civ3D (par_civ1);
    313312        if ~isempty(errormsg)
    314313            disp_uvmat('ERROR',errormsg,checkrun)
    315314            return
    316315        end
    317         for z_index=2:floor((NbSlice-1)/par_civ1.Dz)% loop on slices
     316        %for z_index=2:floor((NbSlice-1)/par_civ1.Dz)% loop on slices
     317        for z_index=2:floor((NbSlice-SearchRange_z)/par_civ1.Dz)% loop on slices
    318318            par_civ1.ImageA=circshift(par_civ1.ImageA,-par_civ1.Dz);%shift the indices in the image block upward by par_civ1.Dz
    319319            par_civ1.ImageB=circshift(par_civ1.ImageA,-par_civ1.Dz);
    … …  
    336336            end
    337337        end
    338         Data.Civ1_C=double(Data.Civ1_C);
    339         Data.Civ1_C(Data.Civ1_C==Inf)=0;
     338        Data.Civ1_V=-Data.Civ1_V;%reverse v
     339        Data.Civ1_Y=npy-Data.Civ1_Y+1;%reverse y
    340340        [npz,npy,npx]=size(Data.Civ1_X);
    341341        Data.Coord_z=SearchRange_z+1:par_civ1.Dz:NbSlice-1;
    342         Data.Coord_y=1:npy;
    343         Data.Coord_x=1:npx;
     342      %  Data.Coord_y=flip(1:npy);
     343       % Data.Coord_x=1:npx;
    344344        % case of mask TO ADAPT
    345345        if par_civ1.CheckMask&&~isempty(par_civ1.Mask)
    … …  
    445445        % perform Patch calculation using the UVMAT fct 'filter_tps'
    446446
    447         [Data.Civ1_SubRange,Data.Civ1_NbCentres,Data.Civ1_Coord_tps,Data.Civ1_U_tps,Data.Civ1_V_tps,tild,Ures, Vres,tild,FFres]=...
     447        [Data.Civ1_SubRange,Data.Civ1_NbCentres,Data.Civ1_Coord_tps,Data.Civ1_U_tps,Data.Civ1_V_tps,~,Ures, Vres,~,FFres]=...
    448448            filter_tps([Data.Civ1_X(ind_good) Data.Civ1_Y(ind_good)],Data.Civ1_U(ind_good),Data.Civ1_V(ind_good),[],Data.Patch1_SubDomainSize,Data.Patch1_FieldSmooth,Data.Patch1_MaxDiff);
    449449        Data.Civ1_U_smooth(ind_good)=Ures;% take the interpolated (smoothed) velocity values for good vectors, keep civ1 data for the other
    … …  
    659659        if exist('ncfile','var')
    660660            CivFile=ncfile;
    661             [Data,tild,tild,errormsg]=nc2struct(CivFile);%read civ1 and detect_false1 data in the existing netcdf file
     661            [Data,~,~,errormsg]=nc2struct(CivFile);%read civ1 and detect_false1 data in the existing netcdf file
    662662            if ~isempty(errormsg)
    663663                disp_uvmat('ERROR',errormsg,checkrun)
    664664                return
    665665            end
    666             %         elseif isfield(Param,'Civ2_X')% use Civ2 data as input in Param (test mode)
    667             %             Data.ListGlobalAttribute={};
    668             %             Data.ListVarName={};
    669             %             Data.VarDimName={};
    670             %             Data.Civ2_X=Param.Civ2_X;
    671             %             Data.Civ2_Y=Param.Civ2_Y;
    672             %             Data.Civ2_U=Param.Civ2_U;
    673             %             Data.Civ2_V=Param.Civ2_V;
    674             %             Data.Civ2_FF=Param.Civ2_FF;
    675666        end
    676667    end
    … …  
    787778
    788779function [xtable,ytable,utable,vtable,wtable,ctable,FF,result_conv,errormsg] = civ3D (par_civ)
     780%% check image sizes
     781[npz,npy_ima,npx_ima]=size(par_civ.ImageA);% npz =
     782if ~isequal(size(par_civ.ImageB),[npz npy_ima npx_ima])
     783    errormsg='image pair with unequal size';
     784    return
     785end
    789786
    790787%% prepare measurement grid
    791 nbinterv_x=floor((par_civ.ImageWidth-1)/par_civ.Dx);
     788nbinterv_x=floor((npx_ima-1)/par_civ.Dx);
    792789gridlength_x=nbinterv_x*par_civ.Dx;
    793 minix=ceil((par_civ.ImageWidth-gridlength_x)/2);
    794 nbinterv_y=floor((par_civ.ImageHeight-1)/par_civ.Dy);
     790minix=ceil((npx_ima-gridlength_x)/2);
     791nbinterv_y=floor((npy_ima-1)/par_civ.Dy);
    795792gridlength_y=nbinterv_y*par_civ.Dy;
    796 miniy=ceil((par_civ.ImageHeight-gridlength_y)/2);
    797 [GridX,GridY]=meshgrid(minix:par_civ.Dx:par_civ.ImageWidth-1,miniy:par_civ.Dy:par_civ.ImageHeight-1);
    798 % minix:par_civ.Dx:par_civ.ImageWidth-1
    799 % miniy:par_civ.Dy:par_civ.ImageHeight-1
    800 % minix=floor(par_civ.Dx/2)-0.5;
    801 % maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx);
    802 % miniy=floor(par_civ.Dy/2)-0.5;% first automatic grid point at half the mesh Dy
    803 % maxiy=miniy+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy);
    804 % [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy);
     793miniy=ceil((npy_ima-gridlength_y)/2);
     794[GridX,GridY]=meshgrid(minix:par_civ.Dx:npx_ima-1,miniy:par_civ.Dy:npy_ima-1);
    805795par_civ.Grid(:,:,1)=GridX;
    806796par_civ.Grid(:,:,2)=GridY;% increases with array index,
    807797[nbvec_y,nbvec_x,~]=size(par_civ.Grid);
    808 %
    809 %
    810 % minix=floor(par_civ.Dx/2)-0.5;
    811 %     maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx);
    812 %     miniy=floor(par_civ.Dy/2)-0.5;% first automatic grid point at half the mesh Dy
    813 %     maxiy=miniy+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy);
    814 %     [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy);
    815 %     par_civ.Grid(:,1)=reshape(GridX,[],1);
    816 %     par_civ.Grid(:,2)=reshape(GridY,[],1);% increases with array index
     798
    817799
    818800%% prepare correlation and search boxes
    … …  
    826808shifty=-round(par_civ.SearchBoxShift(:,2));% sign minus because image j index increases when y decreases
    827809if numel(shiftx)==1% case of a unique shift for the whole field( civ1)
    828     shiftx=shiftx*ones(nbvec_y,nbvec_x,1);
    829     shifty=shifty*ones(nbvec_y,nbvec_x,1);
     810    shiftx=shiftx*ones(nbvec_y,nbvec_x);
     811    shifty=shifty*ones(nbvec_y,nbvec_x);
    830812end
    831813
    832814%% Array initialisation and default output  if par_civ.CorrSmooth=0 (just the grid calculated, no civ computation)
    833815xtable=round(par_civ.Grid(:,:,1)+0.5)-0.5;
    834 ytable=round(par_civ.ImageHeight-par_civ.Grid(:,:,2)+0.5)-0.5;% y index corresponding to the position in image coordiantes
    835 utable=shiftx;%zeros(nbvec,1);
    836 vtable=shifty;%zeros(nbvec,1);
     816%ytable=round(npy_ima-par_civ.Grid(:,:,2)+0.5)-0.5;% y index corresponding to the position in image coordiantes
     817ytable=round(par_civ.Grid(:,:,2)+0.5)-0.5;
     818utable=shiftx;
     819vtable=shifty;
    837820wtable=zeros(size(utable));
    838 ctable=zeros(nbvec_y,nbvec_x,1);
    839 FF=zeros(nbvec_y,nbvec_x,1);
     821ctable=zeros(nbvec_y,nbvec_x);
     822FF=false(nbvec_y,nbvec_x);
    840823result_conv=[];
    841824errormsg='';
    … …  
    844827check_MinIma=isfield(par_civ,'MinIma');% test for image luminosity threshold
    845828check_MaxIma=isfield(par_civ,'MaxIma') && ~isempty(par_civ.MaxIma);
    846 
    847 [npz,npy_ima,npx_ima]=size(par_civ.ImageA);
    848 if ~isequal(size(par_civ.ImageB),[npz npy_ima npx_ima])
    849     errormsg='image pair with unequal size';
    850     return
    851 end
    852829
    853830%% Apply mask
    … …  
    860837checkmask=0;
    861838MinA=min(min(min(par_civ.ImageA)));
    862 %MinB=min(min(par_civ.ImageB));
    863 %check_undefined=false(size(par_civ.ImageA));
    864839if isfield(par_civ,'Mask') && ~isempty(par_civ.Mask)
    865840    checkmask=1;
    … …  
    868843        return
    869844    end
    870     check_undefined=(par_civ.Mask<200 & par_civ.Mask>=20 );
     845    check_undefined=(par_civ.Mask<200 & par_civ.Mask>=20 );% logical of the size of the image block =true for masked pixels
    871846end
    872847
    … …  
    875850sum_square=1;% default
    876851mesh=1;% default
    877 CheckDeformation=isfield(par_civ,'CheckDeformation')&& par_civ.CheckDeformation==1;
    878 if CheckDeformation
    879     mesh=0.25;%mesh in pixels for subpixel image interpolation (x 4 in each direction)
    880     par_civ.CorrSmooth=2;% use SUBPIX2DGAUSS (take into account more points near the max)
    881 end
     852% CheckDeformation=isfield(par_civ,'CheckDeformation')&& par_civ.CheckDeformation==1;
     853% if CheckDeformation
     854%     mesh=0.25;%mesh in pixels for subpixel image interpolation (x 4 in each direction)
     855%     par_civ.CorrSmooth=2;% use SUBPIX2DGAUSS (take into account more points near the max)
     856% end
    882857
    883858if par_civ.CorrSmooth~=0 % par_civ.CorrSmooth=0 implies no civ computation (just input image and grid points given)
    … …  
    885860        for ivec_y=1:nbvec_y
    886861            iref=round(par_civ.Grid(ivec_y,ivec_x,1)+0.5);% xindex on the image A for the middle of the correlation box
    887             jref=round(par_civ.ImageHeight-par_civ.Grid(ivec_y,ivec_x,2)+0.5);%  j index  for the middle of the correlation box in the image A
     862            jref=round(npy_ima-par_civ.Grid(ivec_y,ivec_x,2)+0.5);%  j index  for the middle of the correlation box in the image A
    888863            subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage
    889864            subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage
    … …  
    892867            image1_crop=MinA*ones(npz,numel(subrange1_y),numel(subrange1_x));% default value=min of image A
    893868            image2_crop=MinA*ones(npz,numel(subrange2_y),numel(subrange2_x));% default value=min of image A
    894             check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1
    895             check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight;
    896             check2_x=subrange2_x>=1 & subrange2_x<=par_civ.ImageWidth;% check which points in the subimage 2 are contained in the initial image 2
    897             check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight;
     869            check1_x=subrange1_x>=1 & subrange1_x<=npx_ima;% check which points in the subimage 1 are contained in the initial image 1
     870            check1_y=subrange1_y>=1 & subrange1_y<=npy_ima;
     871            check2_x=subrange2_x>=1 & subrange2_x<=npx_ima;% check which points in the subimage 2 are contained in the initial image 2
     872            check2_y=subrange2_y>=1 & subrange2_y<=npy_ima;
    898873            image1_crop(:,check1_y,check1_x)=par_civ.ImageA(:,subrange1_y(check1_y),subrange1_x(check1_x));%extract a subimage (correlation box) from image A
    899874            image2_crop(:,check2_y,check2_x)=par_civ.ImageB(:,subrange2_y(check2_y),subrange2_x(check2_x));%extract a larger subimage (search box) from image B
    900             if checkmask
    901                 mask1_crop=ones(numel(subrange1_y),numel(subrange1_x));% default value=1 for mask
    902                 mask2_crop=ones(numel(subrange2_y),numel(subrange2_x));% default value=1 for mask
    903                 mask1_crop(check1_y,check1_x)=check_undefined(subrange1_y(check1_y),subrange1_x(check1_x));%extract a mask subimage (correlation box) from image A
    904                 mask2_crop(check2_y,check2_x)=check_undefined(subrange2_y(check2_y),subrange2_x(check2_x));%extract a mask subimage (search box) from image B
    905                 sizemask=sum(sum(mask1_crop))/(numel(subrange1_y)*numel(subrange1_x));%size of the masked part relative to the correlation sub-image
     875            if checkmask% case of mask images
     876                mask1_crop=true(npz,2*iby2+1,2*ibx2+1);% default value=1 for mask
     877                mask2_crop=true(npz,2*isy2+1,2*isx2+1);% default value=1 for mask
     878                mask1_crop(npz,check1_y,check1_x)=check_undefined(npz,subrange1_y(check1_y),subrange1_x(check1_x));%extract a mask subimage (correlation box) from image A
     879                mask2_crop(npz,check2_y,check2_x)=check_undefined(npz,subrange2_y(check2_y),subrange2_x(check2_x));%extract a mask subimage (search box) from image B
     880                sizemask=sum(sum(mask1_crop,2),3)/(npz*(2*ibx2+1)*(2*iby2+1));%size of the masked part relative to the correlation sub-image
    906881                if sizemask > 1/2% eliminate point if more than half of the correlation box is masked
    907882                    FF(ivec_y,ivec_x)=1; %
    … …  
    909884                    vtable(ivec_y,ivec_x)=NaN;
    910885                else
    911                     FF(ivec_y,ivec_x)=0;
    912886                    image1_crop=image1_crop.*~mask1_crop;% put to zero the masked pixels (mask1_crop='true'=1)
    913887                    image2_crop=image2_crop.*~mask2_crop;
    914                     image1_mean=mean(mean(image1_crop))/(1-sizemask);
    915                     image2_mean=mean(mean(image2_crop))/(1-sizemask);
     888                    image1_mean=mean(mean(image1_crop,2),3)/(1-sizemask);
     889                    image2_mean=mean(mean(image2_crop,2),3)/(1-sizemask);
    916890                end
    917891            else
    918                 image1_mean=mean(mean(image1_crop));
    919                 image2_mean=mean(mean(image2_crop));
    920             end
    921             %threshold on image minimum
     892                image1_mean=mean(mean(image1_crop,2),3);
     893                image2_mean=mean(mean(image2_crop,2),3);
     894            end
     895     
    922896            if FF(ivec_y,ivec_x)==0
    923897                if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma)
    924                     FF(ivec_y,ivec_x)=1;
    925                     %threshold on image maximum
     898                    FF(ivec_y,ivec_x)=1;       %threshold on image minimum
     899               
    926900                elseif check_MaxIma && (image1_mean > par_civ.MaxIma || image2_mean > par_civ.MaxIma)
    927                     FF(ivec_y,ivec_x)=1;
     901                    FF(ivec_y,ivec_x)=1;    %threshold on image maximum
    928902                end
    929903                if FF(ivec_y,ivec_x)==1
    930                     utable(ivec_y,ivec_x)=NaN;
     904                    utable(ivec_y,ivec_x)=NaN; 
    931905                    vtable(ivec_y,ivec_x)=NaN;
    932906                else
    933                     %mask
    934                     if checkmask
    935                         image1_crop=(image1_crop-image1_mean).*~mask1_crop;%substract the mean, put to zero the masked parts
    936                         image2_crop=(image2_crop-image2_mean).*~mask2_crop;
    937                     else
    938                         image1_crop=(image1_crop-image1_mean);
    939                         image2_crop=(image2_crop-image2_mean);
    940                     end
     907                    % case of mask
     908                    % if checkmask
     909                    %     image1_crop=(image1_crop-image1_mean).*~mask1_crop;%substract the mean, put to zero the masked parts
     910                    %     image2_crop=(image2_crop-image2_mean).*~mask2_crop;% TO MODIFY !!!!!!!!!!!!!!!!!!!!!
     911                    % else
     912                    %     % image1_crop=(image1_crop-image1_mean);
     913                    %     % image2_crop=(image2_crop-image2_mean);
     914                    % end
    941915
    942916                    npxycorr=par_civ.SearchBoxSize(1:2)-par_civ.CorrBoxSize(1:2)+1;
    943917                    result_conv=zeros([par_civ.SearchBoxSize(3) npxycorr]);%initialise the conv product
    944918                    max_xy=zeros(par_civ.SearchBoxSize(3),1);%initialise the max correlation vs z
    945                     xk=ones(par_civ.SearchBoxSize(3),1);%initialise the x displacement corresponding to the max corr
    946                     yk=ones(par_civ.SearchBoxSize(3),1);%initialise the y displacement corresponding to the max corr
    947                  
    948                     for kz=kref:kref%par_civ.SearchBoxSize(3)
    949                         subima2=squeeze(image2_crop(kz,:,:));
    950                         subima1=squeeze(image1_crop(kref,:,:));
     919                    xk=ones(npz,1);%initialise the x displacement corresponding to the max corr
     920                    yk=ones(npz,1);%initialise the y displacement corresponding to the max corr
     921                    subima1=squeeze(image1_crop(kref,:,:))-image1_mean(kref);
     922                    for kz=1:npz
     923                        subima2=squeeze(image2_crop(kz,:,:))-image2_mean(kz);   
    951924                        correl_xy=conv2(subima2,flip(flip(subima1,2),1),'valid');
    952                         result_conv(kz,:,:)= 255*correl_xy;
     925                        result_conv(kz,:,:)=correl_xy;
    953926                        max_xy(kz)=max(max(correl_xy));
    954                         % max_xy(kz)=max(max(correl_xy));
    955                          [yk(kz),xk(kz)]=find(correl_xy==max_xy(kz),1);
     927                         [yk(kz),xk(kz)]=find(correl_xy>=0.999*max_xy(kz),1);%find the indices of the maximum, use 0.999 to avoid rounding errors
    956928                    end
    957929                    [corrmax,dz]=max(max_xy);
     930                    result_conv=result_conv*255/corrmax;% normalise with a max =255
    958931                    dx=xk(dz);
    959932                    dy=yk(dz);
    960                  
    961                     %
    962                     % result_conv=255*result_conv/corrmax; %normalize, peak=always 255
    963                     % for kz=1:par_civ.SearchBoxSize(3)
    964                     %     [dy,dx] = find(result_conv(kz,:,:)==255,1)
    965                     %     if ~isempty(dy)
    966                     %         dz=kz;
    967                     %         break
    968                     %     end
    969                     % end
    970                     subimage2_crop=squeeze(image2_crop(dz,dy:dy+2*iby2/mesh,dx:dx+2*ibx2/mesh));%subimage of image 2 corresponding to the optimum displacement of first image
    971                     sum_square=sum(sum(squeeze(image1_crop(dz,:,:).*image1_crop(dz,:,:))));
     933                    subimage2_crop=squeeze(image2_crop(dz,dy:dy+2*iby2/mesh,dx:dx+2*ibx2/mesh))-image2_mean(dz);%subimage of image 2 corresponding to the optimum displacement of first image
     934                    sum_square=sum(sum(subima1.*subima1));
    972935                    sum_square=sum_square*sum(sum(subimage2_crop.*subimage2_crop));% product of variances of image 1 and 2
    973936                    sum_square=sqrt(sum_square);% srt of the variance product to normalise correlation
     937                    % if ivec_x==3 && ivec_y==4
     938                    %     'test'
     939                    % end
    974940                    if ~isempty(dz)&& ~isempty(dy) && ~isempty(dx)
    975941                        if par_civ.CorrSmooth==1
    … …  
    987953                        iref=round(xtable(ivec_y,ivec_x)+0.5);% nearest image index for the middle of the vector
    988954                        jref=round(ytable(ivec_y,ivec_x)+0.5);
    989                  
     955                       % if utable(ivec_y,ivec_x)==0 && vtable(ivec_y,ivec_x)==0
     956                       %     'test'
     957                       % end
    990958                        % eliminate vectors located in the mask
    991959                        if  checkmask && (iref<1 || jref<1 ||iref>npx_ima || jref>npy_ima ||( par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100))
    … …  
    1017985
    1018986[npz,npy,npx]=size(result_conv);
    1019 result_conv(result_conv<1)=1; %set to 1 correlation values smaller than 1  (=0 by discretisation, to avoid divergence in the log)
     987result_conv(result_conv<1)=1; %set to 1 correlation values smaller than 1  (to avoid divergence in the log)
    1020988%the following 8 lines are copyright (c) 1998, Uri Shavit, Roi Gurka, Alex Liberzon, Technion ??? Israel Institute of Technology
    1021989%http://urapiv.wordpress.com
    1022990FF=false;
    1023991peakz=z;peaky=y;peakx=x;
    1024 if z < npz && z > 1 && y < npy && y > 1 && x < npx && x > 1
    1025     f0 = log(result_conv(z,y,x));
    1026     f1 = log(result_conv(z-1,y,x));
    1027     f2 = log(result_conv(z+1,y,x));
    1028     peakz = peakz+ (f1-f2)/(2*f1-4*f0+2*f2);
    1029 
    1030     f1 = log(result_conv(z,y-1,x));
    1031     f2 = log(result_conv(z,y+1,x));
    1032     peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
    1033 
    1034     f1 = log(result_conv(z,y,x-1));
    1035     f2 = log(result_conv(z,y,x+1));
    1036     peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
    1037 else
    1038     FF=true;
    1039 end
     992% if z < npz && z > 1 && y < npy && y > 1 && x < npx && x > 1
     993%     f0 = log(result_conv(z,y,x));
     994%     f1 = log(result_conv(z-1,y,x));
     995%     f2 = log(result_conv(z+1,y,x));
     996%     peakz = peakz+ (f1-f2)/(2*f1-4*f0+2*f2);
     997%
     998%     f1 = log(result_conv(z,y-1,x));
     999%     f2 = log(result_conv(z,y+1,x));
     1000%     peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
     1001%
     1002%     f1 = log(result_conv(z,y,x-1));
     1003%     f2 = log(result_conv(z,y,x+1));
     1004%     peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
     1005% else
     1006%     FF=true;
     1007% end
    10401008
    10411009vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1 peakz-floor(npz/2)-1];

  • trunk/src/series/civ_series.m

    r1158 r1160  
    952952%  .DVDY
    953953
    954 function [xtable,ytable,utable,vtable,ctable,F,result_conv,errormsg] = civ (par_civ)
     954function [xtable,ytable,utable,vtable,ctable,FF,result_conv,errormsg] = civ (par_civ)
    955955
    956956%% prepare measurement grid
    … …  
    972972    par_civ.Grid(:,2)=reshape(GridY,[],1);% increases with array index
    973973    %
    974     %
    975974    % minix=floor(par_civ.Dx/2)-0.5;
    976975    % maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx);
    … …  
    10011000vtable=shifty;%zeros(nbvec,1);
    10021001ctable=zeros(nbvec,1);
    1003 F=zeros(nbvec,1);
     1002FF=zeros(nbvec,1);
    10041003result_conv=[];
    10051004errormsg='';
    … …  
    10581057        iref=round(par_civ.Grid(ivec,1)+0.5);% xindex on the image A for the middle of the correlation box
    10591058        jref=round(par_civ.ImageHeight-par_civ.Grid(ivec,2)+0.5);%  j index  for the middle of the correlation box in the image A
    1060         F(ivec)=0;
     1059        FF(ivec)=0;
    10611060        subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage
    10621061        subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage
    … …  
    10781077            sizemask=sum(sum(mask1_crop))/(numel(subrange1_y)*numel(subrange1_x));%size of the masked part relative to the correlation sub-image
    10791078            if sizemask > 1/2% eliminate point if more than half of the correlation box is masked
    1080                 F(ivec)=1; %
     1079                FF(ivec)=1; %
    10811080                utable(ivec)=NaN;
    10821081                vtable(ivec)=NaN;
    … …  
    10921091        end
    10931092        %threshold on image minimum
    1094         if F(ivec)~=1
     1093        if FF(ivec)==0
    10951094            if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma)
    1096                 F(ivec)=1;
     1095                FF(ivec)=1;
    10971096                %threshold on image maximum
    10981097            elseif check_MaxIma && (image1_mean > par_civ.MaxIma || image2_mean > par_civ.MaxIma)
    1099                 F(ivec)=1;
    1100             end
    1101             if F(ivec)==1
     1098                FF(ivec)=1;
     1099            end
     1100            if FF(ivec)==1
    11021101                utable(ivec)=NaN;
    11031102                vtable(ivec)=NaN;
    … …  
    11271126                sum_square=sum(sum(image1_crop.*image1_crop));
    11281127                %reference: Oliver Pust, PIV: Direct Cross-Correlation
     1128                %%%%%% correlation calculation
    11291129                result_conv= conv2(image2_crop,flip(flip(image1_crop,2),1),'valid');
     1130                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    11301131                corrmax= max(max(result_conv));
    11311132                result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255
    … …  
    11381139                    try
    11391140                        if par_civ.CorrSmooth==1
    1140                             [vector,F(ivec)] = SUBPIXGAUSS (result_conv,x,y);
     1141                            [vector,FF(ivec)] = SUBPIXGAUSS (result_conv,x,y);
    11411142                        elseif par_civ.CorrSmooth==2
    1142                             [vector,F(ivec)] = SUBPIX2DGAUSS (result_conv,x,y);
     1143                            [vector,FF(ivec)] = SUBPIX2DGAUSS (result_conv,x,y);
    11431144                        else
    1144                             [vector,F(ivec)] = quadr_fit(result_conv,x,y);
     1145                            [vector,FF(ivec)] = quadr_fit(result_conv,x,y);
    11451146                        end
    11461147                        utable(ivec)=vector(1)*mesh+shiftx(ivec);
    … …  
    11541155                            utable(ivec)=0;
    11551156                            vtable(ivec)=0;
    1156                             F(ivec)=1;
     1157                            FF(ivec)=1;
    11571158                        end
    11581159                        ctable(ivec)=corrmax/sum_square;% correlation value
    11591160                    catch ME
    1160                         F(ivec)=1;
     1161                        FF(ivec)=1;
    11611162                        disp(ME.message)
    11621163                    end
    11631164                else
    1164                     F(ivec)=1;
     1165                    FF(ivec)=1;
    11651166                end
    11661167            end
    … …  
    11881189%http://urapiv.wordpress.com
    11891190peaky = y;peakx=x;
    1190 if y < npy && y > 1 && x < npx-1 && x > 1
    1191     f0 = log(result_conv(y,x));
    1192     f1 = log(result_conv(y-1,x));
    1193     f2 = log(result_conv(y+1,x));
    1194     peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
    1195     f1 = log(result_conv(y,x-1));
    1196     f2 = log(result_conv(y,x+1));
    1197     peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
    1198 else
    1199     F=1; % warning flag for vector truncated by the limited search box
    1200 end
     1191% if y < npy && y > 1 && x < npx-1 && x > 1
     1192%     f0 = log(result_conv(y,x));
     1193%     f1 = log(result_conv(y-1,x));
     1194%     f2 = log(result_conv(y+1,x));
     1195%     peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
     1196%     f1 = log(result_conv(y,x-1));
     1197%     f2 = log(result_conv(y,x+1));
     1198%     peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
     1199% else
     1200%     F=1; % warning flag for vector truncated by the limited search box
     1201% end
    12011202
    12021203vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
Note: See TracChangeset for help on using the changeset viewer.