Ignore:
Timestamp:
Jul 11, 2024, 4:13:03 PM (3 months ago)
Author:
sommeria
Message:

filter_time added

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/series/civ_3D.m

    r1156 r1157  
    130130end
    131131
    132 [tild,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
     132[~,i1_series,~,j1_series,~]=get_file_series(Param);
    133133iview_A=0;% series index (iview) for the first image series
    134134iview_B=0;% series index (iview) for the second image series (only non zero for option 'shift' comparing two image series )
    … …  
    156156PairCiv1=Param.ActionInput.PairIndices.ListPairCiv1;
    157157
    158 [i1_series_Civ1,i2_series_Civ1,j1_series_Civ1,j2_series_Civ1,check_bounds,NomTypeNc]=...
     158[i1_series_Civ1,i2_series_Civ1,j1_series_Civ1,j2_series_Civ1,~,NomTypeNc]=...
    159159                        find_pair_indices(PairCiv1,i1_series{1},j1_series{1},MinIndex_i,MaxIndex_i,MinIndex_j,MaxIndex_j);
    160160                   
    … …  
    168168%% prepare output Data
    169169ListGlobalAttribute={'Conventions','Program','CivStage'};
    170 Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes
    171 Data.Program='civ_series';
     170Data.Conventions='uvmat/civdata_3D';% states the conventions used for the description of field variables and attributes
     171Data.Program='civ_3D';
     172if isfield(Param,'UvmatRevision')
     173    Data.Program=[Data.Program ', uvmat r' Param.UvmatRevision];
     174end
    172175Data.CivStage=0;%default
    173176list_param=(fieldnames(Param.ActionInput.Civ1))';
    174     list_param(strcmp('TestCiv1',list_param))=[];% remove the parameter TestCiv1 from the list
    175     Civ1_param=regexprep(list_param,'^.+','Civ1_$0');% insert 'Civ1_' before  each string in list_param
    176     Civ1_param=[{'Civ1_ImageA','Civ1_ImageB','Civ1_Time','Civ1_Dt'} Civ1_param]; %insert the names of the two input images
     177list_param(strcmp('TestCiv1',list_param))=[];% remove the parameter TestCiv1 from the list
     178Civ1_param=regexprep(list_param,'^.+','Civ1_$0');% insert 'Civ1_' before  each string in list_param
     179Civ1_param=[{'Civ1_ImageA','Civ1_ImageB','Civ1_Time','Civ1_Dt'} Civ1_param]; %insert the names of the two input images
    177180Data.ListGlobalAttribute=[ListGlobalAttribute Civ1_param];
    178             % set the list of variables
    179         Data.ListVarName={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_FF'};%  cell array containing the names of the fields to record
    180         Data.VarDimName={'nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1'};
    181         Data.VarAttribute{1}.Role='coord_x';
    182         Data.VarAttribute{2}.Role='coord_y';
    183         Data.VarAttribute{3}.Role='vector_x';
    184         Data.VarAttribute{4}.Role='vector_y';
    185         Data.VarAttribute{5}.Role='vector_z';
    186         Data.VarAttribute{6}.Role='ancillary';
    187         Data.VarAttribute{7}.Role='errorflag';
     181% set the list of variables
     182Data.ListVarName={'Coord_z','Coord_y','Coord_x','Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_FF'};%  cell array containing the names of the fields to record
     183Data.VarDimName={'npz','npy','npx',{'npz','npy','npx'},{'npz','npy','npx'},{'npz','npy','npx'},{'npz','npy','npx'},...
     184    {'npz','npy','npx'},{'npz','npy','npx'},{'npz','npy','npx'}};
     185Data.VarAttribute{1}.Role='coord_x';
     186Data.VarAttribute{2}.Role='coord_y';
     187Data.VarAttribute{3}.Role='vector_x';
     188Data.VarAttribute{4}.Role='vector_y';
     189Data.VarAttribute{5}.Role='vector_z';
     190Data.VarAttribute{6}.Role='ancillary';
     191Data.VarAttribute{7}.Role='errorflag';
     192Data.Coord_z=j1_series_Civ1;
     193% path for output nc file
     194OutputPath=fullfile(Param.OutputPath,num2str(Param.Experiment),num2str(Param.Device),[Param.OutputSubDir Param.OutputDirExt]);
    188195
    189196%% get timing from the ImaDoc file or input video
    … …  
    227234    par_civ1.ImageA=zeros(2*SearchRange_z+1,npy,npx);
    228235par_civ1.ImageB=zeros(2*SearchRange_z+1,npy,npx);
     236
     237
    229238
    230239%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    … …  
    268277   
    269278    Data.CivStage=1;
    270     % output nc file
    271     ncfile_out=fullfile_uvmat(OutputPath,Param.InputTable{1,2},Param.InputTable{1,3},Param.InputTable{1,5},...
    272         NomTypeNc,i1_series_Civ1(ifield),i2_series_Civ1(ifield),j1_series_Civ1(ifield),j2_series_Civ1(ifield));
     279 
     280    ncfile_out=fullfile_uvmat(OutputPath,'',Param.InputTable{1,3},'.nc',...
     281                '_1-1',i1_series_Civ1(ifield),i2_series_Civ1(ifield)); %output file name
    273282
    274283    if ~CheckOverwrite && exist(ncfile_out,'file')
    … …  
    284293        disp('civ1 started')
    285294
    286 
    287         % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1.ImageNameA and B)
    288 
    289         par_civ1.ImageA=zeros(2*SearchRange_z+1,npy,npx);
     295        % read input images by vertical blocks with nbre of images equal to 2*SearchRange+1,
     296        par_civ1.ImageA=zeros(2*SearchRange_z+1,npy,npx);%image block initiation
    290297        par_civ1.ImageB=zeros(2*SearchRange_z+1,npy,npx);
    291         %first vertical block centered at image index islice=par_civ1.Dz
    292         islice=par_civ1.Dz;
    293          for iz=1:par_civ1.Dz+SearchRange_z
    294               ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1,ifield),[],j1_series_Civ1(iz,1));%
     298       
     299        z_index=1;%first vertical block centered at image index z_index=par_civ1.Dz
     300        for iz=1:par_civ1.Dz+SearchRange_z
     301            j_image_index=j1_series_Civ1(iz,1)% j index of the current image
     302            ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1,ifield),[],j_image_index);%
    295303            A= read_image(ImageName_A,FileType_A);
    296             ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(1,ifield),[],j1_series_Civ1(iz,1));
     304            ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(1,ifield),[],j_image_index);
    297305            B= read_image(ImageName_B,FileType_B);
    298306            par_civ1.ImageA(iz+par_civ1.Dz-1,:,:) = A;
    299307            par_civ1.ImageB(iz+par_civ1.Dz-1,:,:) = B;
    300          end
    301          % caluclate velocity data (y and v in indices, reverse to y component)
    302             [Data.Civ1_X(islice,:,:),Data.Civ1_Y(islice,:,:), utable, vtable,wtable, ctable, FF, result_conv, errormsg] = civ3D (par_civ1);
     308        end
     309        % calculate velocity data at the first position z, corresponding to j_index=par_civ1.Dz
     310        [Data.Civ1_X(1,:,:),Data.Civ1_Y(1,:,:), Data.Civ1_U(1,:,:), Data.Civ1_V(1,:,:),Data.Civ1_W(1,:,:),...
     311            Data.Civ1_C(1,:,:), Data.Civ1_FF(1,:,:), result_conv, errormsg] = civ3D (par_civ1);
     312        if ~isempty(errormsg)
     313            disp_uvmat('ERROR',errormsg,checkrun)
     314            return
     315        end
     316        for z_index=2:floor((NbSlice-1)/par_civ1.Dz)% loop on slices
     317            par_civ1.ImageA=circshift(par_civ1.ImageA,-par_civ1.Dz);%shift the indices in the image block upward by par_civ1.Dz
     318            par_civ1.ImageB=circshift(par_civ1.ImageA,-par_civ1.Dz);
     319            for iz=1:par_civ1.Dz %read the new images at the end of the image block
     320                j_image_index=j1_series_Civ1(z_index*par_civ1.Dz+SearchRange_z-par_civ1.Dz+iz,1)
     321                ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1,ifield),[],j_image_index);%
     322                A= read_image(ImageName_A,FileType_A);
     323                ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(1,ifield),[],j_image_index);
     324                B= read_image(ImageName_B,FileType_B);
     325                par_civ1.ImageA(2*SearchRange_z+1-par_civ1.Dz+iz,:,:) = A;
     326                par_civ1.ImageB(2*SearchRange_z+1-par_civ1.Dz+iz,:,:) = B;
     327            end
     328            % caluclate velocity data (y and v in indices, reverse to y component)
     329            [Data.Civ1_X(z_index,:,:),Data.Civ1_Y(z_index,:,:), Data.Civ1_U(z_index,:,:), Data.Civ1_V(z_index,:,:),Data.Civ1_W(z_index,:,:),...
     330                Data.Civ1_C(z_index,:,:), Data.Civ1_FF(z_index,:,:), result_conv, errormsg] = civ3D (par_civ1);
     331
    303332            if ~isempty(errormsg)
    304333                disp_uvmat('ERROR',errormsg,checkrun)
    305334                return
    306335            end
    307         for islice=2*par_civ1.Dz:NbSlice% loop on slices for the first image in the pair
    308             par_civ1.ImageA=circshift(par_civ1.ImageA,-par_civ1.Dz);%shift the indces in the block upward by par_civ1.Dz
    309             par_civ1.ImageB=circshift(par_civ1.ImageA,-par_civ1.Dz);
    310               for iz=1:par_civ1.Dz
    311             ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1,ifield),[],j1_series_Civ1(islice+SearchRange_z-par_civ1.Dz+iz,1));%
    312             A= read_image(ImageName_A,FileType_A);
    313             ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(1,ifield),[],j1_series_Civ1(islice+SearchRange_z-par_civ1.Dz+iz,1));
    314             B= read_image(ImageName_B,FileType_B);
    315             par_civ1.ImageA(2*SearchRange_z+1-par_civ1.Dz+iz,:,:) = A;
    316             par_civ1.ImageB(2*SearchRange_z+1-par_civ1.Dz+iz,:,:) = B;
    317             end
    318             % caluclate velocity data (y and v in indices, reverse to y component)
    319             [Data.Civ1_X(islice,:,:),Data.Civ1_Y(islice,:,:), utable, vtable,wtable, ctable, FF, result_conv, errormsg] = civ3D (par_civ1);
    320             if ~isempty(errormsg)
    321                 disp_uvmat('ERROR',errormsg,checkrun)
    322                 return
    323             end
    324 
    325         end
     336        end
     337        Data.Civ1_C=double(Data.Civ1_C);
     338        Data.Civ1_C(Data.Civ1_C==Inf)=0;
     339        [npz,npy,npx]=size(Data.Civ1_X);
     340        Data.Coord_z=1:npz;
     341        Data.Coord_y=1:npy;
     342        Data.Coord_x=1:npx;
    326343        % case of mask TO ADAPT
    327344        if par_civ1.CheckMask&&~isempty(par_civ1.Mask)
    … …  
    668685        end
    669686        Data.ListGlobalAttribute=[Data.ListGlobalAttribute Fix2_param];
    670         Data.Civ2_FF=double(detect_false(Param.ActionInput.Fix2,Data.Civ2_C,Data.Civ2_U,Data.Civ2_V,Data.Civ2_FF));
     687        Data.Civ2_FF=detect_false(Param.ActionInput.Fix2,Data.Civ2_C,Data.Civ2_U,Data.Civ2_V,Data.Civ2_FF);
    671688        Data.CivStage=Data.CivStage+1;
    672689    end
    … …  
    728745    time_total=toc(tstart);
    729746    disp(['ellapsed time ' num2str(time_total/60,2) ' minutes'])
    730     disp(['time image reading ' num2str(time_input,2) ' s'])
    731747    disp(['time civ1 ' num2str(time_civ1,2) ' s'])
    732748    disp(['time patch1 ' num2str(time_patch1,2) ' s'])
    733749    disp(['time civ2 ' num2str(time_civ2,2) ' s'])
    734750    disp(['time patch2 ' num2str(time_patch2,2) ' s'])
    735     disp(['time other ' num2str((time_total-time_input-time_civ1-time_patch1-time_civ2-time_patch2),2) ' s'])
    736     % end
     751
    737752end
    738753
    … …  
    860875    for ivec_x=1:nbvec_x
    861876        for ivec_y=1:nbvec_y
    862             ivec_y
    863             iref=round(par_civ.Grid(ivec_y,ivec_x,1)+0.5)% xindex on the image A for the middle of the correlation box
    864             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
     877            iref=round(par_civ.Grid(ivec_y,ivec_x,1)+0.5);% xindex on the image A for the middle of the correlation box
     878            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
    865879            subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage
    866880            subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage
    … …  
    918932
    919933               
    920                     %reference: Oliver Pust, PIV: Direct Cross-Correlation
     934                   
    921935                    for kz=1:par_civ.SearchBoxSize(3)
    922936                        subima2=squeeze(image2_crop(kz,:,:));
    … …  
    928942               
    929943                    end
    930                     [corrmax,z]=max(max_xy);
     944                    [corrmax,zmax]=max(max_xy);
    931945               
    932                     x=xk(z);
    933                     y=yk(z);
     946                    x=xk(zmax);
     947                    y=yk(zmax);
    934948                    result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255
    935                     subimage2_crop=squeeze(image2_crop(z,y:y+2*iby2/mesh,x:x+2*ibx2/mesh));%subimage of image 2 corresponding to the optimum displacement of first image
    936                     sum_square=sum(sum(squeeze(image1_crop(z,:,:).*image1_crop(z,:,:))));
     949                    subimage2_crop=squeeze(image2_crop(zmax,y:y+2*iby2/mesh,x:x+2*ibx2/mesh));%subimage of image 2 corresponding to the optimum displacement of first image
     950                    sum_square=sum(sum(squeeze(image1_crop(zmax,:,:).*image1_crop(zmax,:,:))));
    937951                    sum_square=sum_square*sum(sum(subimage2_crop.*subimage2_crop));% product of variances of image 1 and 2
    938952                    sum_square=sqrt(sum_square);% srt of the variance product to normalise correlation
    … …  
    940954       
    941955                            if par_civ.CorrSmooth==1
    942                                 [vector,FF(ivec_y,ivec_x)] = SUBPIXGAUSS (result_conv(z,:,:),x,y);%TODO: improve by max optimisation along z
     956                                [vector,FF(ivec_y,ivec_x)] = SUBPIXGAUSS (squeeze(result_conv(zmax,:,:)),x,y);%TODO: improve by max optimisation along z
    943957                            elseif par_civ.CorrSmooth==2
    944                                 [vector,FF(ivec_y,ivec_x)] = SUBPIX2DGAUSS (result_conv(z,:,:),x,y);
     958                                [vector,FF(ivec_y,ivec_x)] = SUBPIX2DGAUSS (squeeze(result_conv(zmax,:,:)),x,y);
    945959                            else
    946                                 [vector,FF(ivec_y,ivec_x)] = quadr_fit(result_conv(z,:,:),x,y);
     960                                [vector,FF(ivec_y,ivec_x)] = quadr_fit(squeeze(result_conv(zmax,:,:)),x,y);
    947961                            end
    948962                            utable(ivec_y,ivec_x)=vector(1)*mesh+shiftx(ivec_y,ivec_x);
    … …  
    952966                            iref=round(xtable(ivec_y,ivec_x)+0.5);% nearest image index for the middle of the vector
    953967                            jref=round(ytable(ivec_y,ivec_x)+0.5);
    954                             wtable(ivec_y,ivec_x)=z-kref;
     968                            wtable(ivec_y,ivec_x)=zmax-kref;
    955969                            % eliminate vectors located in the mask
    956970                            if  checkmask && (iref<1 || jref<1 ||iref>npx_ima || jref>npy_ima ||( par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100))
    … …  
    962976
    963977                    else
    964                         FF(ivec_y,ivec_x)=1;
     978                        FF(ivec_y,ivec_x)=true;
    965979                    end
    966980                end
    … …  
    975989% OUPUT:
    976990% vector = optimum displacement vector with subpixel correction
    977 % F =flag: =0 OK
    978 %           =-2 , warning: max too close to the edge of the search box (1 pixel margin)
     991% FF =flag: ='true' max too close to the edge of the search box (1 pixel margin)
    979992% INPUT:
    980993% x,y: position of the maximum correlation at integer values
    981994
    982 function [vector,F] = SUBPIXGAUSS (result_conv,x,y)
     995function [vector,FF] = SUBPIXGAUSS (result_conv,x,y)
    983996%------------------------------------------------------------------------
    984 % vector=[0 0]; %default
    985 F=0;
     997vector=[0 0]; %default
    986998[npy,npx]=size(result_conv);
    987999result_conv(result_conv<1)=1; %set to 1 correlation values smaller than 1  (=0 by discretisation, to avoid divergence in the log)
    9881000%the following 8 lines are copyright (c) 1998, Uri Shavit, Roi Gurka, Alex Liberzon, Technion ??? Israel Institute of Technology
    9891001%http://urapiv.wordpress.com
    990 peaky = y;
    991 if y < npy && y > 1
    992     f0 = log(result_conv(y,x));
    993     f1 = log(result_conv(y-1,x));
    994     f2 = log(result_conv(y+1,x));
    995     peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
    996 else
    997     F=1; % warning flag for vector truncated by the limited search box
    998 end
    999 peakx=x;
    1000 if x < npx-1 && x > 1
    1001     f0 = log(result_conv(y,x));
    1002     f1 = log(result_conv(y,x-1));
    1003     f2 = log(result_conv(y,x+1));
    1004     peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
    1005 else
    1006     F=1; % warning flag for vector truncated by the limited search box
    1007 end
    1008 vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
     1002FF=false;
     1003    peaky = y;
     1004    if y < npy && y > 1
     1005        f0 = log(result_conv(y,x));
     1006        f1 = log(result_conv(y-1,x));
     1007        f2 = log(result_conv(y+1,x));
     1008        peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
     1009    else
     1010        FF=true; % error flag for vector truncated by the limited search box in y
     1011    end
     1012    peakx=x;
     1013    if x < npx-1 && x > 1
     1014        f0 = log(result_conv(y,x));
     1015        f1 = log(result_conv(y,x-1));
     1016        f2 = log(result_conv(y,x+1));
     1017        peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
     1018    else
     1019        FF=true; % warning flag for vector truncated by the limited search box in x
     1020    end
     1021    vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
    10091022
    10101023%------------------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.