Changeset 1160 for trunk/src

Timestamp:

Jul 16, 2024, 9:57:54 AM (2 months ago)

Author:

sommeria

Message:

script_delete_PCO added

Location:

trunk/src

Files:

• script_delete_PCO.m (added)
• script_delete_rdvision.m (modified) (2 diffs)
• series/civ_3D.m (modified) (17 diffs)
• series/civ_series.m (modified) (10 diffs)

trunk/src/script_delete_rdvision.m

@@ -85 +85 @@
                 if exist(ExtractFolder,'dir') && exist(PngFolder,'dir')
                     filename_seq=fullfile(ExtractFolder,'im.seq');
+                    NumberOfFrames=0;
                     try
-                    s=ini2struct(filename_seq);
-                    FileInfo=s.sequenceSettings;
-                    if isfield(s.sequenceSettings,'numberoffiles')
-                        NumberOfFrames=str2double(s.sequenceSettings.numberoffiles);
-                    else
-                        status='bad seq file';
-                    end
+                        s=ini2struct(filename_seq);
+                        FileInfo=s.sequenceSettings;
+
+                        if isfield(s.sequenceSettings,'numberoffiles')
+                            NumberOfFrames=str2double(s.sequenceSettings.numberoffiles);
+                        else
+                            status='bad seq file';
+                        end
                     catch ME
                         disp(['error in ' filename_seq])
                     end
-                    DirPng=dir(PngFolder);
+                    DirPng=dir(PngFolder);         
+
                     if numel(DirPng)==NumberOfFrames+2
-                        Checkdelete=1;
+                        PngCells=struct2cell(DirPng);
+                        mm=cell2mat(PngCells(4,:));% check the sizes of extracted images
+                        sizemax=max(mm);
+                        sizemin=min(mm(3:end));
+                        disp(['max size(Mbytes)=' num2str(sizemax/1000000)]);
+                        if min(mm(3:end))<0.9*sizemax
+                            status=['WARNING' 'min size(Mbytes)=' num2str(sizemin/1000000)];
+                        else
+                            Checkdelete=1;% approve deletion of the source bin files
+                        end
                     else
                         status=['extraction not finished,' num2str(numel(DirPng)-2) ' images extracted'];
@@ -141 +153 @@
     end
 end
+'END'

trunk/src/series/civ_3D.m

@@ -16 +16 @@
 %           if absent, the fct looks for input data in Param.ActionInput     (test mode)
 %     Param.OutputSubDir: sets the folder name of output file(s,
-%           if absent no file is produced, result in the output structure Data (test mode)
 %     Param.ActionInput: substructure with the parameters provided by the GUI civ_input
 %                      .Civ1: parameters for civ1cc
@@ -185 +184 @@
 Data.VarAttribute{1}.Role='coord_z';
 Data.VarAttribute{2}.Role='coord_x';
-ata.VarAttribute{3}.Role='coord_y';
+Data.VarAttribute{3}.Role='coord_y';
 Data.VarAttribute{4}.Role='vector_x';
 Data.VarAttribute{5}.Role='vector_y';
@@ -310 +309 @@
         % calculate velocity data at the first position z, corresponding to j_index=par_civ1.Dz
         [Data.Civ1_X(1,:,:),Data.Civ1_Y(1,:,:), Data.Civ1_U(1,:,:), Data.Civ1_V(1,:,:),Data.Civ1_W(1,:,:),...
-            Data.Civ1_C(1,:,:), Data.Civ1_FF(1,:,:), result_conv, errormsg] = civ3D (par_civ1);
+            Data.Civ1_C(1,:,:), Data.Civ1_FF(1,:,:), ~, errormsg] = civ3D (par_civ1);
         if ~isempty(errormsg)
             disp_uvmat('ERROR',errormsg,checkrun)
             return
         end
-        for z_index=2:floor((NbSlice-1)/par_civ1.Dz)% loop on slices
+        %for z_index=2:floor((NbSlice-1)/par_civ1.Dz)% loop on slices
+        for z_index=2:floor((NbSlice-SearchRange_z)/par_civ1.Dz)% loop on slices
             par_civ1.ImageA=circshift(par_civ1.ImageA,-par_civ1.Dz);%shift the indices in the image block upward by par_civ1.Dz
             par_civ1.ImageB=circshift(par_civ1.ImageA,-par_civ1.Dz);
@@ -336 +336 @@
             end
         end
-        Data.Civ1_C=double(Data.Civ1_C);
-        Data.Civ1_C(Data.Civ1_C==Inf)=0;
+        Data.Civ1_V=-Data.Civ1_V;%reverse v
+        Data.Civ1_Y=npy-Data.Civ1_Y+1;%reverse y
         [npz,npy,npx]=size(Data.Civ1_X);
         Data.Coord_z=SearchRange_z+1:par_civ1.Dz:NbSlice-1;
-        Data.Coord_y=1:npy;
-        Data.Coord_x=1:npx;
+      %  Data.Coord_y=flip(1:npy);
+       % Data.Coord_x=1:npx;
         % case of mask TO ADAPT
         if par_civ1.CheckMask&&~isempty(par_civ1.Mask)
@@ -445 +445 @@
         % perform Patch calculation using the UVMAT fct 'filter_tps'
 
-        [Data.Civ1_SubRange,Data.Civ1_NbCentres,Data.Civ1_Coord_tps,Data.Civ1_U_tps,Data.Civ1_V_tps,tild,Ures, Vres,tild,FFres]=...
+        [Data.Civ1_SubRange,Data.Civ1_NbCentres,Data.Civ1_Coord_tps,Data.Civ1_U_tps,Data.Civ1_V_tps,~,Ures, Vres,~,FFres]=...
             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);
         Data.Civ1_U_smooth(ind_good)=Ures;% take the interpolated (smoothed) velocity values for good vectors, keep civ1 data for the other
@@ -659 +659 @@
         if exist('ncfile','var')
             CivFile=ncfile;
-            [Data,tild,tild,errormsg]=nc2struct(CivFile);%read civ1 and detect_false1 data in the existing netcdf file
+            [Data,~,~,errormsg]=nc2struct(CivFile);%read civ1 and detect_false1 data in the existing netcdf file
             if ~isempty(errormsg)
                 disp_uvmat('ERROR',errormsg,checkrun)
                 return
             end
-            %         elseif isfield(Param,'Civ2_X')% use Civ2 data as input in Param (test mode)
-            %             Data.ListGlobalAttribute={};
-            %             Data.ListVarName={};
-            %             Data.VarDimName={};
-            %             Data.Civ2_X=Param.Civ2_X;
-            %             Data.Civ2_Y=Param.Civ2_Y;
-            %             Data.Civ2_U=Param.Civ2_U;
-            %             Data.Civ2_V=Param.Civ2_V;
-            %             Data.Civ2_FF=Param.Civ2_FF;
         end
     end
@@ -787 +778 @@
 
 function [xtable,ytable,utable,vtable,wtable,ctable,FF,result_conv,errormsg] = civ3D (par_civ)
+%% check image sizes
+[npz,npy_ima,npx_ima]=size(par_civ.ImageA);% npz = 
+if ~isequal(size(par_civ.ImageB),[npz npy_ima npx_ima])
+    errormsg='image pair with unequal size';
+    return
+end
 
 %% prepare measurement grid
-nbinterv_x=floor((par_civ.ImageWidth-1)/par_civ.Dx);
+nbinterv_x=floor((npx_ima-1)/par_civ.Dx);
 gridlength_x=nbinterv_x*par_civ.Dx;
-minix=ceil((par_civ.ImageWidth-gridlength_x)/2);
-nbinterv_y=floor((par_civ.ImageHeight-1)/par_civ.Dy);
+minix=ceil((npx_ima-gridlength_x)/2);
+nbinterv_y=floor((npy_ima-1)/par_civ.Dy);
 gridlength_y=nbinterv_y*par_civ.Dy;
-miniy=ceil((par_civ.ImageHeight-gridlength_y)/2);
-[GridX,GridY]=meshgrid(minix:par_civ.Dx:par_civ.ImageWidth-1,miniy:par_civ.Dy:par_civ.ImageHeight-1);
-% minix:par_civ.Dx:par_civ.ImageWidth-1
-% miniy:par_civ.Dy:par_civ.ImageHeight-1
-% minix=floor(par_civ.Dx/2)-0.5;
-% maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx);
-% miniy=floor(par_civ.Dy/2)-0.5;% first automatic grid point at half the mesh Dy
-% maxiy=miniy+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy);
-% [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy);
+miniy=ceil((npy_ima-gridlength_y)/2);
+[GridX,GridY]=meshgrid(minix:par_civ.Dx:npx_ima-1,miniy:par_civ.Dy:npy_ima-1);
 par_civ.Grid(:,:,1)=GridX;
 par_civ.Grid(:,:,2)=GridY;% increases with array index,
 [nbvec_y,nbvec_x,~]=size(par_civ.Grid);
-% 
-% 
-% minix=floor(par_civ.Dx/2)-0.5;
-%     maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx);
-%     miniy=floor(par_civ.Dy/2)-0.5;% first automatic grid point at half the mesh Dy
-%     maxiy=miniy+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy);
-%     [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy);
-%     par_civ.Grid(:,1)=reshape(GridX,[],1);
-%     par_civ.Grid(:,2)=reshape(GridY,[],1);% increases with array index
+
 
 %% prepare correlation and search boxes
@@ -826 +808 @@
 shifty=-round(par_civ.SearchBoxShift(:,2));% sign minus because image j index increases when y decreases
 if numel(shiftx)==1% case of a unique shift for the whole field( civ1)
-    shiftx=shiftx*ones(nbvec_y,nbvec_x,1);
-    shifty=shifty*ones(nbvec_y,nbvec_x,1);
+    shiftx=shiftx*ones(nbvec_y,nbvec_x);
+    shifty=shifty*ones(nbvec_y,nbvec_x);
 end
 
 %% Array initialisation and default output  if par_civ.CorrSmooth=0 (just the grid calculated, no civ computation)
 xtable=round(par_civ.Grid(:,:,1)+0.5)-0.5;
-ytable=round(par_civ.ImageHeight-par_civ.Grid(:,:,2)+0.5)-0.5;% y index corresponding to the position in image coordiantes
-utable=shiftx;%zeros(nbvec,1);
-vtable=shifty;%zeros(nbvec,1);
+%ytable=round(npy_ima-par_civ.Grid(:,:,2)+0.5)-0.5;% y index corresponding to the position in image coordiantes
+ytable=round(par_civ.Grid(:,:,2)+0.5)-0.5;
+utable=shiftx;
+vtable=shifty;
 wtable=zeros(size(utable));
-ctable=zeros(nbvec_y,nbvec_x,1);
-FF=zeros(nbvec_y,nbvec_x,1);
+ctable=zeros(nbvec_y,nbvec_x);
+FF=false(nbvec_y,nbvec_x);
 result_conv=[];
 errormsg='';
@@ -844 +827 @@
 check_MinIma=isfield(par_civ,'MinIma');% test for image luminosity threshold
 check_MaxIma=isfield(par_civ,'MaxIma') && ~isempty(par_civ.MaxIma);
-
-[npz,npy_ima,npx_ima]=size(par_civ.ImageA);
-if ~isequal(size(par_civ.ImageB),[npz npy_ima npx_ima])
-    errormsg='image pair with unequal size';
-    return
-end
 
 %% Apply mask
@@ -860 +837 @@
 checkmask=0;
 MinA=min(min(min(par_civ.ImageA)));
-%MinB=min(min(par_civ.ImageB));
-%check_undefined=false(size(par_civ.ImageA));
 if isfield(par_civ,'Mask') && ~isempty(par_civ.Mask)
     checkmask=1;
@@ -868 +843 @@
         return
     end
-    check_undefined=(par_civ.Mask<200 & par_civ.Mask>=20 );
+    check_undefined=(par_civ.Mask<200 & par_civ.Mask>=20 );% logical of the size of the image block =true for masked pixels
 end
 
@@ -875 +850 @@
 sum_square=1;% default
 mesh=1;% default
-CheckDeformation=isfield(par_civ,'CheckDeformation')&& par_civ.CheckDeformation==1;
-if CheckDeformation
-    mesh=0.25;%mesh in pixels for subpixel image interpolation (x 4 in each direction)
-    par_civ.CorrSmooth=2;% use SUBPIX2DGAUSS (take into account more points near the max)
-end
+% CheckDeformation=isfield(par_civ,'CheckDeformation')&& par_civ.CheckDeformation==1;
+% if CheckDeformation
+%     mesh=0.25;%mesh in pixels for subpixel image interpolation (x 4 in each direction)
+%     par_civ.CorrSmooth=2;% use SUBPIX2DGAUSS (take into account more points near the max)
+% end
 
 if par_civ.CorrSmooth~=0 % par_civ.CorrSmooth=0 implies no civ computation (just input image and grid points given)
@@ -885 +860 @@
         for ivec_y=1:nbvec_y
             iref=round(par_civ.Grid(ivec_y,ivec_x,1)+0.5);% xindex on the image A for the middle of the correlation box
-            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
+            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
             subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage
             subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage
@@ -892 +867 @@
             image1_crop=MinA*ones(npz,numel(subrange1_y),numel(subrange1_x));% default value=min of image A
             image2_crop=MinA*ones(npz,numel(subrange2_y),numel(subrange2_x));% default value=min of image A
-            check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1
-            check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight;
-            check2_x=subrange2_x>=1 & subrange2_x<=par_civ.ImageWidth;% check which points in the subimage 2 are contained in the initial image 2
-            check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight;
+            check1_x=subrange1_x>=1 & subrange1_x<=npx_ima;% check which points in the subimage 1 are contained in the initial image 1
+            check1_y=subrange1_y>=1 & subrange1_y<=npy_ima;
+            check2_x=subrange2_x>=1 & subrange2_x<=npx_ima;% check which points in the subimage 2 are contained in the initial image 2
+            check2_y=subrange2_y>=1 & subrange2_y<=npy_ima;
             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
             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
-            if checkmask
-                mask1_crop=ones(numel(subrange1_y),numel(subrange1_x));% default value=1 for mask
-                mask2_crop=ones(numel(subrange2_y),numel(subrange2_x));% default value=1 for mask
-                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
-                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
-                sizemask=sum(sum(mask1_crop))/(numel(subrange1_y)*numel(subrange1_x));%size of the masked part relative to the correlation sub-image
+            if checkmask% case of mask images
+                mask1_crop=true(npz,2*iby2+1,2*ibx2+1);% default value=1 for mask
+                mask2_crop=true(npz,2*isy2+1,2*isx2+1);% default value=1 for mask
+                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
+                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
+                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
                 if sizemask > 1/2% eliminate point if more than half of the correlation box is masked
                     FF(ivec_y,ivec_x)=1; %
@@ -909 +884 @@
                     vtable(ivec_y,ivec_x)=NaN;
                 else
-                    FF(ivec_y,ivec_x)=0;
                     image1_crop=image1_crop.*~mask1_crop;% put to zero the masked pixels (mask1_crop='true'=1)
                     image2_crop=image2_crop.*~mask2_crop;
-                    image1_mean=mean(mean(image1_crop))/(1-sizemask);
-                    image2_mean=mean(mean(image2_crop))/(1-sizemask);
+                    image1_mean=mean(mean(image1_crop,2),3)/(1-sizemask);
+                    image2_mean=mean(mean(image2_crop,2),3)/(1-sizemask);
                 end
             else
-                image1_mean=mean(mean(image1_crop));
-                image2_mean=mean(mean(image2_crop));
-            end
-            %threshold on image minimum
+                image1_mean=mean(mean(image1_crop,2),3);
+                image2_mean=mean(mean(image2_crop,2),3);
+            end
+     
             if FF(ivec_y,ivec_x)==0
                 if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma)
-                    FF(ivec_y,ivec_x)=1;
-                    %threshold on image maximum
+                    FF(ivec_y,ivec_x)=1;       %threshold on image minimum
+                
                 elseif check_MaxIma && (image1_mean > par_civ.MaxIma || image2_mean > par_civ.MaxIma)
-                    FF(ivec_y,ivec_x)=1;
+                    FF(ivec_y,ivec_x)=1;    %threshold on image maximum
                 end
                 if FF(ivec_y,ivec_x)==1
-                    utable(ivec_y,ivec_x)=NaN;
+                    utable(ivec_y,ivec_x)=NaN;  
                     vtable(ivec_y,ivec_x)=NaN;
                 else
-                    %mask
-                    if checkmask
-                        image1_crop=(image1_crop-image1_mean).*~mask1_crop;%substract the mean, put to zero the masked parts
-                        image2_crop=(image2_crop-image2_mean).*~mask2_crop;
-                    else
-                        image1_crop=(image1_crop-image1_mean);
-                        image2_crop=(image2_crop-image2_mean);
-                    end
+                    % case of mask
+                    % if checkmask
+                    %     image1_crop=(image1_crop-image1_mean).*~mask1_crop;%substract the mean, put to zero the masked parts
+                    %     image2_crop=(image2_crop-image2_mean).*~mask2_crop;% TO MODIFY !!!!!!!!!!!!!!!!!!!!!
+                    % else
+                    %     % image1_crop=(image1_crop-image1_mean);
+                    %     % image2_crop=(image2_crop-image2_mean);
+                    % end
 
                     npxycorr=par_civ.SearchBoxSize(1:2)-par_civ.CorrBoxSize(1:2)+1;
                     result_conv=zeros([par_civ.SearchBoxSize(3) npxycorr]);%initialise the conv product
                     max_xy=zeros(par_civ.SearchBoxSize(3),1);%initialise the max correlation vs z
-                    xk=ones(par_civ.SearchBoxSize(3),1);%initialise the x displacement corresponding to the max corr 
-                    yk=ones(par_civ.SearchBoxSize(3),1);%initialise the y displacement corresponding to the max corr 
-                  
-                    for kz=kref:kref%par_civ.SearchBoxSize(3)
-                        subima2=squeeze(image2_crop(kz,:,:));
-                        subima1=squeeze(image1_crop(kref,:,:));
+                    xk=ones(npz,1);%initialise the x displacement corresponding to the max corr 
+                    yk=ones(npz,1);%initialise the y displacement corresponding to the max corr 
+                    subima1=squeeze(image1_crop(kref,:,:))-image1_mean(kref);
+                    for kz=1:npz
+                        subima2=squeeze(image2_crop(kz,:,:))-image2_mean(kz);   
                         correl_xy=conv2(subima2,flip(flip(subima1,2),1),'valid');
-                        result_conv(kz,:,:)= 255*correl_xy;
+                        result_conv(kz,:,:)=correl_xy;
                         max_xy(kz)=max(max(correl_xy));
-                        % max_xy(kz)=max(max(correl_xy));
-                         [yk(kz),xk(kz)]=find(correl_xy==max_xy(kz),1);
+                         [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
                     end
                     [corrmax,dz]=max(max_xy);
+                    result_conv=result_conv*255/corrmax;% normalise with a max =255
                     dx=xk(dz);
                     dy=yk(dz);
-                 
-                    % 
-                    % result_conv=255*result_conv/corrmax; %normalize, peak=always 255
-                    % for kz=1:par_civ.SearchBoxSize(3)
-                    %     [dy,dx] = find(result_conv(kz,:,:)==255,1)
-                    %     if ~isempty(dy)
-                    %         dz=kz;
-                    %         break
-                    %     end
-                    % end
-                    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
-                    sum_square=sum(sum(squeeze(image1_crop(dz,:,:).*image1_crop(dz,:,:))));
+                    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
+                    sum_square=sum(sum(subima1.*subima1));
                     sum_square=sum_square*sum(sum(subimage2_crop.*subimage2_crop));% product of variances of image 1 and 2
                     sum_square=sqrt(sum_square);% srt of the variance product to normalise correlation
+                    % if ivec_x==3 && ivec_y==4
+                    %     'test'
+                    % end
                     if ~isempty(dz)&& ~isempty(dy) && ~isempty(dx)
                         if par_civ.CorrSmooth==1
@@ -987 +953 @@
                         iref=round(xtable(ivec_y,ivec_x)+0.5);% nearest image index for the middle of the vector
                         jref=round(ytable(ivec_y,ivec_x)+0.5);
-                  
+                       % if utable(ivec_y,ivec_x)==0 && vtable(ivec_y,ivec_x)==0
+                       %     'test'
+                       % end
                         % eliminate vectors located in the mask
                         if  checkmask && (iref<1 || jref<1 ||iref>npx_ima || jref>npy_ima ||( par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100))
@@ -1017 +985 @@
 
 [npz,npy,npx]=size(result_conv);
-result_conv(result_conv<1)=1; %set to 1 correlation values smaller than 1  (=0 by discretisation, to avoid divergence in the log)
+result_conv(result_conv<1)=1; %set to 1 correlation values smaller than 1  (to avoid divergence in the log)
 %the following 8 lines are copyright (c) 1998, Uri Shavit, Roi Gurka, Alex Liberzon, Technion ??? Israel Institute of Technology
 %http://urapiv.wordpress.com
 FF=false;
 peakz=z;peaky=y;peakx=x;
-if z < npz && z > 1 && y < npy && y > 1 && x < npx && x > 1
-    f0 = log(result_conv(z,y,x));
-    f1 = log(result_conv(z-1,y,x));
-    f2 = log(result_conv(z+1,y,x));
-    peakz = peakz+ (f1-f2)/(2*f1-4*f0+2*f2);
-
-    f1 = log(result_conv(z,y-1,x));
-    f2 = log(result_conv(z,y+1,x));
-    peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
-
-    f1 = log(result_conv(z,y,x-1));
-    f2 = log(result_conv(z,y,x+1));
-    peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
-else
-    FF=true;
-end
+% if z < npz && z > 1 && y < npy && y > 1 && x < npx && x > 1
+%     f0 = log(result_conv(z,y,x));
+%     f1 = log(result_conv(z-1,y,x));
+%     f2 = log(result_conv(z+1,y,x));
+%     peakz = peakz+ (f1-f2)/(2*f1-4*f0+2*f2);
+% 
+%     f1 = log(result_conv(z,y-1,x));
+%     f2 = log(result_conv(z,y+1,x));
+%     peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
+% 
+%     f1 = log(result_conv(z,y,x-1));
+%     f2 = log(result_conv(z,y,x+1));
+%     peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
+% else
+%     FF=true;
+% end
 
 vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1 peakz-floor(npz/2)-1];

trunk/src/series/civ_series.m

@@ -952 +952 @@
 %  .DVDY
 
-function [xtable,ytable,utable,vtable,ctable,F,result_conv,errormsg] = civ (par_civ)
+function [xtable,ytable,utable,vtable,ctable,FF,result_conv,errormsg] = civ (par_civ)
 
 %% prepare measurement grid
@@ -972 +972 @@
     par_civ.Grid(:,2)=reshape(GridY,[],1);% increases with array index
     % 
-    % 
     % minix=floor(par_civ.Dx/2)-0.5;
     % maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx);
@@ -1001 +1000 @@
 vtable=shifty;%zeros(nbvec,1);
 ctable=zeros(nbvec,1);
-F=zeros(nbvec,1);
+FF=zeros(nbvec,1);
 result_conv=[];
 errormsg='';
@@ -1058 +1057 @@
         iref=round(par_civ.Grid(ivec,1)+0.5);% xindex on the image A for the middle of the correlation box
         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
-        F(ivec)=0;
+        FF(ivec)=0;
         subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage
         subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage
@@ -1078 +1077 @@
             sizemask=sum(sum(mask1_crop))/(numel(subrange1_y)*numel(subrange1_x));%size of the masked part relative to the correlation sub-image
             if sizemask > 1/2% eliminate point if more than half of the correlation box is masked
-                F(ivec)=1; %
+                FF(ivec)=1; %
                 utable(ivec)=NaN;
                 vtable(ivec)=NaN;
@@ -1092 +1091 @@
         end
         %threshold on image minimum
-        if F(ivec)~=1
+        if FF(ivec)==0
             if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma)
-                F(ivec)=1;
+                FF(ivec)=1;
                 %threshold on image maximum
             elseif check_MaxIma && (image1_mean > par_civ.MaxIma || image2_mean > par_civ.MaxIma)
-                F(ivec)=1;
-            end
-            if F(ivec)==1
+                FF(ivec)=1;
+            end
+            if FF(ivec)==1
                 utable(ivec)=NaN;
                 vtable(ivec)=NaN;
@@ -1127 +1126 @@
                 sum_square=sum(sum(image1_crop.*image1_crop));
                 %reference: Oliver Pust, PIV: Direct Cross-Correlation
+                %%%%%% correlation calculation
                 result_conv= conv2(image2_crop,flip(flip(image1_crop,2),1),'valid');
+                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 corrmax= max(max(result_conv));
                 result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255
@@ -1138 +1139 @@
                     try
                         if par_civ.CorrSmooth==1
-                            [vector,F(ivec)] = SUBPIXGAUSS (result_conv,x,y);
+                            [vector,FF(ivec)] = SUBPIXGAUSS (result_conv,x,y);
                         elseif par_civ.CorrSmooth==2
-                            [vector,F(ivec)] = SUBPIX2DGAUSS (result_conv,x,y);
+                            [vector,FF(ivec)] = SUBPIX2DGAUSS (result_conv,x,y);
                         else
-                            [vector,F(ivec)] = quadr_fit(result_conv,x,y);
+                            [vector,FF(ivec)] = quadr_fit(result_conv,x,y);
                         end
                         utable(ivec)=vector(1)*mesh+shiftx(ivec);
@@ -1154 +1155 @@
                             utable(ivec)=0;
                             vtable(ivec)=0;
-                            F(ivec)=1;
+                            FF(ivec)=1;
                         end
                         ctable(ivec)=corrmax/sum_square;% correlation value
                     catch ME
-                        F(ivec)=1;
+                        FF(ivec)=1;
                         disp(ME.message)
                     end
                 else
-                    F(ivec)=1;
+                    FF(ivec)=1;
                 end
             end
@@ -1188 +1189 @@
 %http://urapiv.wordpress.com
 peaky = y;peakx=x;
-if y < npy && y > 1 && x < npx-1 && x > 1
-    f0 = log(result_conv(y,x));
-    f1 = log(result_conv(y-1,x));
-    f2 = log(result_conv(y+1,x));
-    peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
-    f1 = log(result_conv(y,x-1));
-    f2 = log(result_conv(y,x+1));
-    peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
-else
-    F=1; % warning flag for vector truncated by the limited search box
-end
+% if y < npy && y > 1 && x < npx-1 && x > 1
+%     f0 = log(result_conv(y,x));
+%     f1 = log(result_conv(y-1,x));
+%     f2 = log(result_conv(y+1,x));
+%     peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
+%     f1 = log(result_conv(y,x-1));
+%     f2 = log(result_conv(y,x+1));
+%     peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
+% else
+%     F=1; % warning flag for vector truncated by the limited search box
+% end
 
 vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];