Ignore:
Timestamp:
Jan 26, 2015, 12:37:56 AM (6 years ago)
Author:
sommeria
Message:

civtest_implemented_civ2

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/series/civ_series.m

    r855 r856  
    5858    Data.AllowInputSort='off';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default)
    5959    Data.WholeIndexRange='off';% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default)
    60     if isfield(Data,'ActionInput') && strcmp(Data.ActionInput.PairIndices.ListPairMode,'pair j1-j2')
     60    if isfield(Data,'ActionInput') && isfield(Data.ActionInput,'PairIndices')&& strcmp(Data.ActionInput.PairIndices.ListPairMode,'pair j1-j2')
    6161        Data.Desable_j_index='on';% hide the j index in series (set by the pair choice in civ_input)
    6262    end
     
    174174                j1_series_Civ2=j1_series_Civ1;
    175175                j2_series_Civ2=j2_series_Civ1;
    176                 NomTypeNc=NomType;
     176                NomTypeNc=Param.InputTable{1,4};
    177177            case 'PIV volume'
    178178                % TODO, TODO
     
    329329
    330330%%%%% MAIN LOOP %%%%%%
     331maskname='';% initiate the mask name
    331332for ifield=1:NbField
    332333    if ~isempty(RUNHandle)% update the waitbar in interactive mode with GUI series  (checkrun=1)
     
    360361    if isfield (Param.ActionInput,'Civ1')
    361362        par_civ1=Param.ActionInput.Civ1;
    362         if CheckInputFile % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1)
     363        if CheckInputFile % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1.ImageNameA and B)
    363364            try
    364365                ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(ifield),[],j1_series_Civ1(ifield));
     
    422423        Data.VarAttribute{4}.Role='vector_y';
    423424        Data.VarAttribute{5}.Role='warnflag';
    424        
     425        if par_civ1.CheckMask&&~isempty(par_civ1.Mask)
     426            if strcmp(maskname,par_civ1.Mask)% mask exist, not already read in civ1
     427                par_civ1.Mask=mask; %use mask already opened
     428            else
     429                par_civ1.Mask=imread(par_civ1.Mask);%update the mask, an store it for future use
     430                mask=par_civ1.Mask;
     431                maskname=par_civ1.Mask;
     432            end
     433        end
    425434        if strcmp(Param.ActionInput.ListCompareMode, 'PIV volume')
    426435            Data.ListVarName=[Data.ListVarName 'Civ1_Z'];
     
    429438            for ivol=1:NbSlice
    430439                % caluclate velocity data (y and v in indices, reverse to y component)
    431                 [xtable ytable utable vtable ctable F result_conv errormsg] = civ (par_civ1);
     440                [xtable, ytable, utable, vtable, ctable, F, result_conv, errormsg] = civ (par_civ1);
    432441                if ~isempty(errormsg)
    433442                     disp_uvmat('ERROR',errormsg,checkrun)
     
    444453        else %usual PIV
    445454            % caluclate velocity data (y and v in indices, reverse to y component)
    446             [xtable ytable utable vtable ctable F result_conv errormsg] = civ (par_civ1);
     455            [xtable, ytable, utable, vtable, ctable, F, result_conv, errormsg] = civ (par_civ1);
    447456            if ~isempty(errormsg)
    448457                disp_uvmat('ERROR',errormsg,checkrun)
     
    473482                [Data,tild,tild,errormsg]=nc2struct(CivFile);%read civ1 and fix1 data in the existing netcdf file
    474483            end
    475         else
     484        elseif isfield(Param,'Civ1_X')
    476485            Data.ListGlobalAttribute={};
    477486            Data.ListVarName={};
     
    572581    if isfield (Param.ActionInput,'Civ2')
    573582        par_civ2=Param.ActionInput.Civ2;
    574         par_civ2.ImageA=[];
    575         par_civ2.ImageB=[];
    576         i1=i1_series_Civ2(ifield);
    577         i2=i1;
    578         if ~isempty(i2_series_Civ2)
    579             i2=i2_series_Civ2(ifield);
    580         end
    581         j1=1;
    582         if ~isempty(j1_series_Civ2)
    583             j1=j1_series_Civ2(ifield);
    584         end
    585         j2=j1;
    586         if ~isempty(j2_series_Civ2)
    587             j2=j2_series_Civ2(ifield);
    588         end
    589         ImageName_A_Civ2=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1,[],j1);
    590 
    591         if strcmp(ImageName_A_Civ2,ImageName_A) && isequal(FrameIndex_A_Civ1(ifield),FrameIndex_A_Civ2(ifield))
    592             par_civ2.ImageA=par_civ1.ImageA;
    593         else
    594             [par_civ2.ImageA,VideoObject_A] = read_image(ImageName_A_Civ2,FileType_A,VideoObject_A,FrameIndex_A_Civ2(ifield));
    595         end
    596         ImageName_B_Civ2=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2,[],j2);
    597         if strcmp(ImageName_B_Civ2,ImageName_B) && isequal(FrameIndex_B_Civ1(ifield),FrameIndex_B_Civ2)
    598             par_civ2.ImageB=par_civ1.ImageB;
    599         else
    600             [par_civ2.ImageB,VideoObject_B] = read_image(ImageName_B_Civ2,FileType_B,VideoObject_B,FrameIndex_B_Civ2(ifield));
    601         end     
    602        
    603         ncfile=fullfile_uvmat(RootPath_A,OutputDir,RootFile_A,'.nc',NomTypeNc,i1,i2,j1,j2);
    604         par_civ2.ImageWidth=FileInfo_A.Width;
    605         par_civ2.ImageHeight=FileInfo_A.Height;
    606        
    607         if isfield(par_civ2,'Grid')% grid points set as input file
    608             if ischar(par_civ2.Grid)%read the grid file if the input is a file name
    609                 par_civ2.Grid=dlmread(par_civ2.Grid);
    610                 par_civ2.Grid(1,:)=[];%the first line must be removed (heading in the grid file)
    611             end
    612         else% automatic grid
    613             minix=floor(par_civ2.Dx/2)-0.5;
    614             maxix=minix+par_civ2.Dx*floor((par_civ2.ImageWidth-1)/par_civ2.Dx);
    615             miniy=floor(par_civ2.Dy/2)-0.5;
    616             maxiy=minix+par_civ2.Dy*floor((par_civ2.ImageHeight-1)/par_civ2.Dy);
    617             [GridX,GridY]=meshgrid(minix:par_civ2.Dx:maxix,miniy:par_civ2.Dy:maxiy);
    618             par_civ2.Grid(:,1)=reshape(GridX,[],1);
    619             par_civ2.Grid(:,2)=reshape(GridY,[],1);
    620         end
    621         Shiftx=zeros(size(par_civ2.Grid,1),1);% shift expected from civ1 data
    622         Shifty=zeros(size(par_civ2.Grid,1),1);
    623         nbval=zeros(size(par_civ2.Grid,1),1);
     583        if CheckInputFile % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1.ImageNameA and B)
     584            par_civ2.ImageA=[];
     585            par_civ2.ImageB=[];
     586            i1=i1_series_Civ2(ifield);
     587            i2=i1;
     588            if ~isempty(i2_series_Civ2)
     589                i2=i2_series_Civ2(ifield);
     590            end
     591            j1=1;
     592            if ~isempty(j1_series_Civ2)
     593                j1=j1_series_Civ2(ifield);
     594            end
     595            j2=j1;
     596            if ~isempty(j2_series_Civ2)
     597                j2=j2_series_Civ2(ifield);
     598            end
     599            ImageName_A_Civ2=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1,[],j1);
     600           
     601            if strcmp(ImageName_A_Civ2,ImageName_A) && isequal(FrameIndex_A_Civ1(ifield),FrameIndex_A_Civ2(ifield))
     602                par_civ2.ImageA=par_civ1.ImageA;
     603            else
     604                [par_civ2.ImageA,VideoObject_A] = read_image(ImageName_A_Civ2,FileType_A,VideoObject_A,FrameIndex_A_Civ2(ifield));
     605            end
     606            ImageName_B_Civ2=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2,[],j2);
     607            if strcmp(ImageName_B_Civ2,ImageName_B) && isequal(FrameIndex_B_Civ1(ifield),FrameIndex_B_Civ2)
     608                par_civ2.ImageB=par_civ1.ImageB;
     609            else
     610                [par_civ2.ImageB,VideoObject_B] = read_image(ImageName_B_Civ2,FileType_B,VideoObject_B,FrameIndex_B_Civ2(ifield));
     611            end
     612            if strcmp(Param.ActionInput.ListCompareMode,'PIV')
     613                ncfile=fullfile_uvmat(RootPath_A,OutputDir,RootFile_A,'.nc',NomTypeNc,i1,i2,j1,j2);
     614            else % displacement
     615                ncfile=fullfile_uvmat(RootPath_A,OutputDir,RootFile_A,'.nc',NomTypeNc,i2,[],j2);
     616            end
     617            par_civ2.ImageWidth=FileInfo_A.Width;
     618            par_civ2.ImageHeight=FileInfo_A.Height;
     619            if isfield(par_civ2,'Grid')% grid points set as input file
     620                if ischar(par_civ2.Grid)%read the grid file if the input is a file name
     621                    par_civ2.Grid=dlmread(par_civ2.Grid);
     622                    par_civ2.Grid(1,:)=[];%the first line must be removed (heading in the grid file)
     623                end
     624            else% automatic grid
     625                minix=floor(par_civ2.Dx/2)-0.5;
     626                maxix=minix+par_civ2.Dx*floor((par_civ2.ImageWidth-1)/par_civ2.Dx);
     627                miniy=floor(par_civ2.Dy/2)-0.5;
     628                maxiy=minix+par_civ2.Dy*floor((par_civ2.ImageHeight-1)/par_civ2.Dy);
     629                [GridX,GridY]=meshgrid(minix:par_civ2.Dx:maxix,miniy:par_civ2.Dy:maxiy);
     630                par_civ2.Grid(:,1)=reshape(GridX,[],1);
     631                par_civ2.Grid(:,2)=reshape(GridY,[],1);
     632            end
     633        end
    624634        if par_civ2.CheckDeformation
    625635            DUDX=zeros(size(par_civ2.Grid,1),1);
     
    629639        end
    630640       
    631          % get the guess from patch1 or patch2 (case 'CheckCiv3')
    632         if isfield (par_civ2,'CheckCiv3') && par_civ2.CheckCiv3 %get the guess from  patch2
    633            SubRange= Data.Civ2_SubRange;
    634            NbCentres=Data.Civ2_NbCentres;
    635            Coord_tps=Data.Civ2_Coord_tps;
    636            U_tps=Data.Civ2_U_tps;
    637            V_tps=Data.Civ2_V_tps;
    638            CivStage=Data.CivStage;
    639            Civ1_Dt=Data.Civ2_Dt;
    640            Data=[];%reinitialise the result structure Data
    641            Data.ListGlobalAttribute={'Conventions','Program','CivStage'};
    642            Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes
    643            Data.Program='civ_series';
    644            Data.CivStage=CivStage;
    645            Data.ListVarName={};
    646            Data.VarDimName={};
    647         else % get the guess from patch1
    648            SubRange= Data.Civ1_SubRange;
    649            NbCentres=Data.Civ1_NbCentres;
    650            Coord_tps=Data.Civ1_Coord_tps;
    651            U_tps=Data.Civ1_U_tps;
    652            V_tps=Data.Civ1_V_tps;
    653            Civ1_Dt=Data.Civ1_Dt;
    654         end
    655         NbSubDomain=size(SubRange,3);       
     641        % get the guess from patch1 or patch2 (case 'CheckCiv3')
     642        if CheckInputFile % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1.ImageNameA and B)
     643            if isfield (par_civ2,'CheckCiv3') && par_civ2.CheckCiv3 %get the guess from  patch2
     644                SubRange= Data.Civ2_SubRange;
     645                NbCentres=Data.Civ2_NbCentres;
     646                Coord_tps=Data.Civ2_Coord_tps;
     647                U_tps=Data.Civ2_U_tps;
     648                V_tps=Data.Civ2_V_tps;
     649                CivStage=Data.CivStage;
     650                Civ1_Dt=Data.Civ2_Dt;
     651                Data=[];%reinitialise the result structure Data
     652                Data.ListGlobalAttribute={'Conventions','Program','CivStage'};
     653                Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes
     654                Data.Program='civ_series';
     655                Data.CivStage=CivStage;
     656                Data.ListVarName={};
     657                Data.VarDimName={};
     658            else % get the guess from patch1
     659                SubRange= Data.Civ1_SubRange;
     660                NbCentres=Data.Civ1_NbCentres;
     661                Coord_tps=Data.Civ1_Coord_tps;
     662                U_tps=Data.Civ1_U_tps;
     663                V_tps=Data.Civ1_V_tps;
     664                Civ1_Dt=Data.Civ1_Dt;
     665            end
     666        else
     667            SubRange= par_civ2.Civ1_SubRange;
     668            NbCentres=par_civ2.Civ1_NbCentres;
     669            Coord_tps=par_civ2.Civ1_Coord_tps;
     670            U_tps=par_civ2.Civ1_U_tps;
     671            V_tps=par_civ2.Civ1_V_tps;
     672            Civ1_Dt=par_civ2.Civ1_Dt;
     673                         Data.ListVarName={};
     674                Data.VarDimName={};
     675        end
     676        Shiftx=zeros(size(par_civ2.Grid,1),1);% shift expected from civ1 data
     677        Shifty=zeros(size(par_civ2.Grid,1),1);
     678        nbval=zeros(size(par_civ2.Grid,1),1);% nbre of interpolated values at each grid point (from the different patch subdomains)
     679        NbSubDomain=size(SubRange,3);
    656680        for isub=1:NbSubDomain% for each sub-domain of Patch1
    657681            nbvec_sub=NbCentres(isub);% nbre of Civ vectors in the subdomain
     
    660684            epoints = par_civ2.Grid(ind_sel,:);% coordinates of interpolation sites
    661685            ctrs=Coord_tps(1:nbvec_sub,:,isub) ;%(=initial points) ctrs
    662             nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
    663             EM = tps_eval(epoints,ctrs);               
     686            nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for each interpolation point (in case of subdomain overlap)
     687            EM = tps_eval(epoints,ctrs);
    664688            Shiftx(ind_sel)=Shiftx(ind_sel)+EM*U_tps(1:nbvec_sub+3,isub);
    665689            Shifty(ind_sel)=Shifty(ind_sel)+EM*V_tps(1:nbvec_sub+3,isub);
     
    672696            end
    673697        end
    674         mask='';
    675         if par_civ2.CheckMask&&~isempty(par_civ2.Mask)&& ~strcmp(maskname,par_civ2.Mask)% mask exist, not already read in civ1
    676             mask=imread(par_civ2.Mask);
     698        if par_civ2.CheckMask&&~isempty(par_civ2.Mask)
     699            if strcmp(maskname,par_civ2.Mask)% mask exist, not already read in civ1
     700                par_civ2.Mask=mask; %use mask already opened
     701            else
     702                par_civ2.Mask=imread(par_civ2.Mask);%update the mask, and store it for future use
     703                mask=par_civ2.Mask;
     704                maskname=par_civ2.Mask;
     705            end
    677706        end
    678707        ibx2=ceil(par_civ2.CorrBoxSize(1)/2);
     
    680709        par_civ2.SearchBoxSize(1)=2*ibx2+9;% search ara +-4 pixels around the guess
    681710        par_civ2.SearchBoxSize(2)=2*iby2+9;
    682         Civ2_Dt=time(i2+1,j2+1)-time(i1+1,j1+1);
     711        if strcmp(Param.ActionInput.ListCompareMode,'displacement')
     712            Civ1_Dt=1;
     713            Civ2_Dt=1;
     714        elseif exist('time','var')
     715            Civ2_Dt=time(i2+1,j2+1)-time(i1+1,j1+1);
     716        else
     717            Civ2_Dt=1;%TO CHECK, TEST
     718        end
    683719        par_civ2.SearchBoxShift=(Civ2_Dt/Civ1_Dt)*[Shiftx(nbval>=1)./nbval(nbval>=1) Shifty(nbval>=1)./nbval(nbval>=1)];
    684720        par_civ2.Grid=[par_civ2.Grid(nbval>=1,1)-par_civ2.SearchBoxShift(:,1)/2 par_civ2.Grid(nbval>=1,2)-par_civ2.SearchBoxShift(:,2)/2];% grid taken at the extrapolated origin of the displacement vectors
     
    690726        end
    691727        % calculate velocity data (y and v in indices, reverse to y component)
    692         [xtable ytable utable vtable ctable F] = civ (par_civ2);
     728        [xtable, ytable, utable, vtable, ctable, F,result_conv,errormsg] = civ (par_civ2);
     729       
    693730        list_param=(fieldnames(Param.ActionInput.Civ2))';
    694731        Civ2_param=regexprep(list_param,'^.+','Civ2_$0');% insert 'Civ2_' before  each string in list_param
    695732        Civ2_param=[{'Civ2_ImageA','Civ2_ImageB','Civ2_Time','Civ2_Dt'} Civ2_param]; %insert the names of the two input images
    696         %indicate the values of all the global attributes in the output data
     733        %indicate the values of all the global attributes in the output data
     734        if exist('ImageName_A','var')
    697735        Data.Civ2_ImageA=ImageName_A;
    698736        Data.Civ2_ImageB=ImageName_B;
    699737        Data.Civ2_Time=(time(i2+1,j2+1)+time(i1+1,j1+1))/2;
    700738        Data.Civ2_Dt=Civ2_Dt;
    701 %         Data.Civ2_Time=1;
    702 %         Data.Civ2_Dt=1;
     739        end
     740        %         Data.Civ2_Time=1;
     741        %         Data.Civ2_Dt=1;
    703742        for ilist=1:length(list_param)
    704743            Data.(Civ2_param{4+ilist})=Param.ActionInput.Civ2.(list_param{ilist});
     
    725764    %% Fix2
    726765    if isfield (Param.ActionInput,'Fix2')
    727                 list_param=fieldnames(Param.ActionInput.Fix2)';
     766        list_param=fieldnames(Param.ActionInput.Fix2)';
    728767        Fix2_param=regexprep(list_param,'^.+','Fix2_$0');% insert 'Fix1_' before  each string in ListFixParam
    729768        %indicate the values of all the global attributes in the output data
     
    824863%INPUT:
    825864% par_civ: structure of input parameters, with fields:
    826 %  .CorrBoxSize
    827 %  .SearchBoxSize
    828 %  .SearchBoxShift
    829 %  .ImageHeight
    830 %  .ImageWidth
    831 %  .Dx, Dy
    832 %  .Grid
    833 %  .Mask
    834 %  .MinIma
     865%  .ImageA: first image for correlation (matrix)
     866%  .ImageB: second image for correlation(matrix)
     867%  .CorrBoxSize: 1,2 vector giving the size of the correlation box in x and y
     868%  .SearchBoxSize:  1,2 vector giving the size of the search box in x and y
     869%  .SearchBoxShift: 1,2 vector or 2 column matrix (for civ2) giving the shift of the search box in x and y
     870%  .CorrSmooth: =1 or 2 determines the choice of the sub-pixel determination of the correlation max
     871%  .ImageWidth: nb of pixels of the image in x
     872%  .Dx, Dy: mesh for the PIV calculation
     873%  .Grid: grid giving the PIV calculation points (alternative to .Dx .Dy)
     874%  .Mask: name of a mask file or mask image matrix itself
     875%  .MinIma: thresholds for image luminosity
    835876%  .MaxIma
    836 %  .image1:first image (matrix)
    837 % image2: second image (matrix)
    838 % ibx2,iby2: half size of the correlation box along x and y, in px (size=(2*iby2+1,2*ibx2+1)
    839 % isx2,isy2: half size of the search box along x and y, in px (size=(2*isy2+1,2*isx2+1)
    840 % shiftx, shifty: shift of the search box (in pixel index, yshift reversed)
    841 % step: mesh of the measurement points (in px)
    842 % subpixfinder=1 or 2 controls the curve fitting of the image correlation
    843 % mask: =[] for no mask
    844 % roi: 4 element vector defining a region of interest: x position, y position, width, height, (in image indices), for the whole image, roi=[];
     877%  .CheckDeformation=1 for subpixel interpolation and image deformation (linear transform)
     878%  .DUDX: matrix of deformation obtained from patch at each grid point
     879%  .DUDY
     880%  .DVDX:
     881%  .DVDY
     882
    845883function [xtable,ytable,utable,vtable,ctable,F,result_conv,errormsg] = civ (par_civ)
    846884
    847885%% prepare measurement grid
    848886if isfield(par_civ,'Grid')% grid points set as input
    849     if ischar(par_civ.Grid)%read the drid file if the input is a file name
     887    if ischar(par_civ.Grid)%read the grid file if the input is a file name
    850888        par_civ.Grid=dlmread(par_civ.Grid);
    851889        par_civ.Grid(1,:)=[];%the first line must be removed (heading in the grid file)
    852890    end
     891    % else par_civ.Grid is already an array
    853892else% automatic grid
    854893    minix=floor(par_civ.Dx/2)-0.5;
     
    889928        return    % get the grid only, no civ calculation
    890929    elseif ischar(par_civ.Mask)
    891         par_civ.Mask=imread(par_civ.Mask);
     930        par_civ.Mask=imread(par_civ.Mask);% read the mask if not allready done
    892931    end
    893932end
     
    948987sum_square=1;% default
    949988mesh=1;% default
    950 CheckDecimal=isfield(par_civ,'CheckDecimal')&& par_civ.CheckDecimal==1;
    951 if CheckDecimal
    952     mesh=0.2;%mesh in pixels for subpixel image interpolation
    953     CheckDeformation=isfield(par_civ,'CheckDeformation')&& par_civ.CheckDeformation==1;
    954 end
    955 % vector=[0 0];%default
    956 for ivec=1:nbvec
    957     iref=round(par_civ.Grid(ivec,1)+0.5);% xindex on the image A for the middle of the correlation box
    958     jref=round(par_civ.ImageHeight-par_civ.Grid(ivec,2)+0.5);% yindex on the image B for the middle of the correlation box
    959     %if ~(checkmask && par_civ.Mask(jref,iref)<=20) %velocity not set to zero by the black mask
    960     %         if jref-iby2<1 || jref+iby2>par_civ.ImageHeight|| iref-ibx2<1 || iref+ibx2>par_civ.ImageWidth||...
    961     %               jref+shifty(ivec)-isy2<1||jref+shifty(ivec)+isy2>par_civ.ImageHeight|| iref+shiftx(ivec)-isx2<1 || iref+shiftx(ivec)+isx2>par_civ.ImageWidth  % we are outside the image
    962     %             F(ivec)=3;
    963     %         else
    964     F(ivec)=0;
    965     subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage
    966     subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage
    967     subrange2_x=iref+shiftx(ivec)-isx2:iref+shiftx(ivec)+isx2;%x indices defining the second subimage
    968     subrange2_y=jref+shifty(ivec)-isy2:jref+shifty(ivec)+isy2;%y indices defining the second subimage
    969     image1_crop=MinA*ones(numel(subrange1_y),numel(subrange1_x));% default value=min of image A
    970     image2_crop=MinA*ones(numel(subrange2_y),numel(subrange2_x));% default value=min of image A
    971     check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1
    972     check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight;
    973     check2_x=subrange2_x>=1 & subrange2_x<=par_civ.ImageWidth;% check which points in the subimage 2 are contained in the initial image 2
    974     check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight;
    975    
    976     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
    977     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
    978     image1_mean=mean(mean(image1_crop));
    979     image2_mean=mean(mean(image2_crop));
    980     %threshold on image minimum
    981     if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma)
    982         F(ivec)=3;
    983     end
    984     %threshold on image maximum
    985     if check_MaxIma && (image1_mean > par_civ.MaxIma || image2_mean > par_civ.MaxIma)
    986         F(ivec)=3;
    987     end
    988     %         end
    989     if F(ivec)~=3
    990         image1_crop=image1_crop-image1_mean;%substract the mean
    991         image2_crop=image2_crop-image2_mean;
    992         if CheckDecimal
    993             xi=(1:mesh:size(image1_crop,2));
    994             yi=(1:mesh:size(image1_crop,1))';
     989CheckDeformation=isfield(par_civ,'CheckDeformation')&& par_civ.CheckDeformation==1;
     990if CheckDeformation
     991    mesh=0.25;%mesh in pixels for subpixel image interpolation (x 4 in each direction)
     992end
     993if par_civ.CorrSmooth~=0 % par_civ.CorrSmooth=0 implies no civ computation (just input image and grid points viven)
     994    for ivec=1:nbvec
     995        iref=round(par_civ.Grid(ivec,1)+0.5);% xindex on the image A for the middle of the correlation box
     996        jref=round(par_civ.ImageHeight-par_civ.Grid(ivec,2)+0.5);% yindex on the image B for the middle of the correlation box
     997        F(ivec)=0;
     998        subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage
     999        subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage
     1000        subrange2_x=iref+shiftx(ivec)-isx2:iref+shiftx(ivec)+isx2;%x indices defining the second subimage
     1001        subrange2_y=jref+shifty(ivec)-isy2:jref+shifty(ivec)+isy2;%y indices defining the second subimage
     1002        image1_crop=MinA*ones(numel(subrange1_y),numel(subrange1_x));% default value=min of image A
     1003        image2_crop=MinA*ones(numel(subrange2_y),numel(subrange2_x));% default value=min of image A
     1004        check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1
     1005        check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight;
     1006        check2_x=subrange2_x>=1 & subrange2_x<=par_civ.ImageWidth;% check which points in the subimage 2 are contained in the initial image 2
     1007        check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight;     
     1008        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
     1009        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
     1010        image1_mean=mean(mean(image1_crop));
     1011        image2_mean=mean(mean(image2_crop));
     1012        %threshold on image minimum
     1013        if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma)
     1014            F(ivec)=3;
     1015        end
     1016        %threshold on image maximum
     1017        if check_MaxIma && (image1_mean > par_civ.MaxIma || image2_mean > par_civ.MaxIma)
     1018            F(ivec)=3;
     1019        end
     1020       
     1021        if F(ivec)~=3
     1022            image1_crop=image1_crop-image1_mean;%substract the mean
     1023            image2_crop=image2_crop-image2_mean;
    9951024            if CheckDeformation
     1025                xi=(1:mesh:size(image1_crop,2));
     1026                yi=(1:mesh:size(image1_crop,1))';
    9961027                [XI,YI]=meshgrid(xi-ceil(size(image1_crop,2)/2),yi-ceil(size(image1_crop,1)/2));
    9971028                XIant=XI-par_civ.DUDX(ivec)*XI-par_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2);
    9981029                YIant=YI-par_civ.DVDX(ivec)*XI-par_civ.DVDY(ivec)*YI+ceil(size(image1_crop,1)/2);
    9991030                image1_crop=interp2(image1_crop,XIant,YIant);
    1000             else
    1001                 image1_crop=interp2(image1_crop,xi,yi);
    1002             end
    1003             xi=(1:mesh:size(image2_crop,2));
    1004             yi=(1:mesh:size(image2_crop,1))';
    1005             image2_crop=interp2(image2_crop,xi,yi);
    1006         end
    1007         sum_square=sum(sum(image1_crop.*image1_crop));
    1008         %reference: Oliver Pust, PIV: Direct Cross-Correlation
    1009         result_conv= conv2(image2_crop,flipdim(flipdim(image1_crop,2),1),'valid');
    1010         corrmax= max(max(result_conv));
    1011         result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255
    1012         %Find the correlation max, at 255
    1013         [y,x] = find(result_conv==255,1);
    1014         if ~isempty(y) && ~isempty(x)
    1015             try
    1016                 if par_civ.CorrSmooth==1
    1017                     [vector,F(ivec)] = SUBPIXGAUSS (result_conv,x,y);
    1018                 elseif par_civ.CorrSmooth==2
    1019                     [vector,F(ivec)] = SUBPIX2DGAUSS (result_conv,x,y);
    1020                 end
    1021 %                 if ~isfield(par_civ,'CheckDeformation')
    1022                 utable(ivec)=vector(1)*mesh+shiftx(ivec);
    1023                 vtable(ivec)=vector(2)*mesh+shifty(ivec);
    1024 %                 else
    1025 %                                 utable(ivec)=shiftx(ivec);% TEST TEST !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1026 %                 vtable(ivec)=shifty(ivec);
    1027 %                 end
    1028                 xtable(ivec)=iref+utable(ivec)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2)
    1029                 ytable(ivec)=jref+vtable(ivec)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge)
    1030                 iref=round(xtable(ivec));% image index for the middle of the vector
    1031                 jref=round(ytable(ivec));
    1032                 if checkmask && par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100
    1033                     utable(ivec)=0;
    1034                     vtable(ivec)=0;
     1031                xi=(1:mesh:size(image2_crop,2));
     1032                yi=(1:mesh:size(image2_crop,1))';
     1033                image2_crop=interp2(image2_crop,xi,yi);
     1034            end
     1035            sum_square=sum(sum(image1_crop.*image1_crop));
     1036            %reference: Oliver Pust, PIV: Direct Cross-Correlation
     1037            result_conv= conv2(image2_crop,flipdim(flipdim(image1_crop,2),1),'valid');
     1038            corrmax= max(max(result_conv));
     1039            result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255
     1040            %Find the correlation max, at 255
     1041            [y,x] = find(result_conv==255,1);
     1042            if ~isempty(y) && ~isempty(x)
     1043                try
     1044                    if par_civ.CorrSmooth==1
     1045                        [vector,F(ivec)] = SUBPIXGAUSS (result_conv,x,y);
     1046                    elseif par_civ.CorrSmooth==2
     1047                        [vector,F(ivec)] = SUBPIX2DGAUSS (result_conv,x,y);
     1048                    end
     1049                    utable(ivec)=vector(1)*mesh+shiftx(ivec);
     1050                    vtable(ivec)=vector(2)*mesh+shifty(ivec);
     1051                    xtable(ivec)=iref+utable(ivec)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2)
     1052                    ytable(ivec)=jref+vtable(ivec)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge)
     1053                    iref=round(xtable(ivec));% image index for the middle of the vector
     1054                    jref=round(ytable(ivec));
     1055                    if checkmask && par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100
     1056                        utable(ivec)=0;
     1057                        vtable(ivec)=0;
     1058                        F(ivec)=3;
     1059                    end
     1060                    ctable(ivec)=corrmax/sum_square;% correlation value
     1061                catch ME
    10351062                    F(ivec)=3;
    10361063                end
    1037                 ctable(ivec)=corrmax/sum_square;% correlation value
    1038             catch ME
     1064            else
    10391065                F(ivec)=3;
    10401066            end
    1041         else
    1042             F(ivec)=3;
    10431067        end
    10441068    end
Note: See TracChangeset for help on using the changeset viewer.