Ignore:
Timestamp:
Jul 19, 2024, 5:57:46 AM (2 months ago)
Author:
sommeria
Message:

filter_tps modified to avoid possible bugs in case of few vectors, subpixel interpolation slightly modified in civ

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/series/stereo_civ.m

    r1152 r1161  
    3838%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
    3939
    40 function [Data,errormsg,result_conv]= stereo_civ(Param)
     40function [Data,errormsg,result_conv, XmlData]= stereo_civ(Param)
    4141Data=[];
    4242errormsg='';
    … …  
    8585WaitbarHandle=findobj(hseries,'Tag','Waitbar');%handle of waitbar in GUI series
    8686
     87%%%%%%%%%%%%%%%% Modif pour test %%%%%%%%%%%%%%%%%%
     88% CheckInputFile=isfield(Param,'InputTable');%= 1 in test use for TestCiv (no nc file involved)
     89CheckOutputFile = isfield(Param,'OutputSubDir') ; %= 1 in test use for TestPatch (no nc file produced)
     90%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     91
    8792%% input files and indexing
    8893MaxIndex_i=Param.IndexRange.MaxIndex_i;
    … …  
    95100time=[];
    96101for iview=1:size(Param.InputTable,1)
     102    %XmlFileName=find_imadoc(Param.InputTable{iview,1},Param.InputTable{iview,2},Param.InputTable{iview,3},Param.InputTable{iview,5});
    97103    XmlFileName=find_imadoc(Param.InputTable{iview,1},Param.InputTable{iview,2});
    98104    if isempty(XmlFileName)
    … …  
    235241
    236242%% Output directory
    237 OutputDir=[Param.OutputSubDir Param.OutputDirExt];
     243OutputDir='/Test';
     244if CheckOutputFile
     245    OutputDir=[Param.OutputSubDir Param.OutputDirExt];
     246end
     247
    238248
    239249ListGlobalAttribute={'Conventions','Program','CivStage'};
    … …  
    252262end
    253263if isempty(time)% time = index i  by default
    254     MaxIndex_i=max(i2_series_Civ1);
    255     MaxIndex_j=max(j2_series_Civ1);
     264    MaxIndex_i=max(i2_series_Civ1, [], 'all');
     265    MaxIndex_j=max(j2_series_Civ1, [], 'all');
    256266    time=(1:MaxIndex_i)'*ones(1,MaxIndex_j);
    257267    time=[zeros(1,MaxIndex_j);time];% insert a first line of zeros
    … …  
    312322
    313323       
    314         [A,Rangx,Rangy]=phys_ima(A,XmlData,1);
     324        [A,Rangx,Rangy]=phys_ima(A,XmlData,1);%transform image A in phys coordinates
    315325        [Npy,Npx]=size(A{1});
    316326        PhysImageA=fullfile_uvmat(RootPath_A,Civ1Dir,RootFile_A,'.png','_1a',i1_series_Civ1(ifield),[],1);
    … …  
    372382        Data.Civ1_F=reshape(F,[],1);
    373383
     384        %%%%%%%%%%%%%%%% ajout fonction test %%%%%%%
     385        if ~isfield (Param.ActionInput,'Civ2') && ~isfield (Param.ActionInput,'Civ3') && ~isfield(Param.ActionInput,'Patch1')
     386            [Data.Xmid, Data.Ymid, Data.Uphys, Data.Vphys, Data.Zphys, ...
     387                Data.Yphys, Data.Xphys, Data.Error] = getPhysValues(Rangx, ...
     388                Rangy, Npx, Npy, Data.Civ1_X, Data.Civ1_Y, Data.Civ1_U, ...
     389                Data.Civ1_V, XmlData);
     390   
     391            if ~isempty(errormsg)
     392                disp_uvmat('ERROR',errormsg,checkrun)
     393                return
     394            end
     395        end
     396        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     397
    374398    end
    375399   
    … …  
    414438            return
    415439        end
    416        
    417         Data.ListGlobalAttribute=[Data.ListGlobalAttribute {'Patch1_Rho','Patch1_Threshold','Patch1_SubDomain'}];
     440
     441        %%%%%%%%%%%%%%%%%%%% modif fonction test %%%%%%%%%%%%%%%%%%%%%%%
     442        try
     443            Data.ListGlobalAttribute=[Data.ListGlobalAttribute {'Patch1_Rho','Patch1_Threshold','Patch1_SubDomain'}];
     444        catch
     445            Data.ListGlobalAttribute={'Patch1_Rho','Patch1_Threshold','Patch1_SubDomain'};
     446        end
     447        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     448
     449           
    418450        Data.Patch1_FieldSmooth=Param.ActionInput.Patch1.FieldSmooth;
    419451        Data.Patch1_MaxDiff=Param.ActionInput.Patch1.MaxDiff;
    420452        Data.Patch1_SubDomainSize=Param.ActionInput.Patch1.SubDomainSize;
    421         nbvar=length(Data.ListVarName);
    422         Data.ListVarName=[Data.ListVarName {'Civ1_U_smooth','Civ1_V_smooth','Civ1_SubRange','Civ1_NbCentres','Civ1_Coord_tps','Civ1_U_tps','Civ1_V_tps'}];
    423         Data.VarDimName=[Data.VarDimName {'nb_vec_1','nb_vec_1',{'nb_coord','nb_bounds','nb_subdomain_1'},'nb_subdomain_1',...
    424             {'nb_tps_1','nb_coord','nb_subdomain_1'},{'nb_tps_1','nb_subdomain_1'},{'nb_tps_1','nb_subdomain_1'}}];
     453
     454        %%%%%%%%%%%%%%%%%%% modif fonction test %%%%%%%%%%%%%%%%%%%%%%%%
     455        try
     456            nbvar=length(Data.ListVarName);
     457            Data.ListVarName=[Data.ListVarName {'Civ1_U_smooth','Civ1_V_smooth','Civ1_SubRange','Civ1_NbCentres','Civ1_Coord_tps','Civ1_U_tps','Civ1_V_tps'}];
     458            Data.VarDimName=[Data.VarDimName {'nb_vec_1','nb_vec_1',{'nb_coord','nb_bounds','nb_subdomain_1'},'nb_subdomain_1',...
     459                {'nb_tps_1','nb_coord','nb_subdomain_1'},{'nb_tps_1','nb_subdomain_1'},{'nb_tps_1','nb_subdomain_1'}}];
     460        catch
     461            nbvar=0;
     462            Data.ListVarName={'Civ1_U_smooth','Civ1_V_smooth','Civ1_SubRange','Civ1_NbCentres','Civ1_Coord_tps','Civ1_U_tps','Civ1_V_tps'};
     463            Data.VarDimName={'nb_vec_1','nb_vec_1',{'nb_coord','nb_bounds','nb_subdomain_1'},'nb_subdomain_1',...
     464                {'nb_tps_1','nb_coord','nb_subdomain_1'},{'nb_tps_1','nb_subdomain_1'},{'nb_tps_1','nb_subdomain_1'}};
     465        end
     466       
     467        if ~isfield(Param.ActionInput,'Civ1')
     468            Data.Civ1_X = Param.Civ1_X ;
     469            Data.Civ1_Y = Param.Civ1_Y ;
     470            Data.Civ1_U = Param.Civ1_U ;
     471            Data.Civ1_V = Param.Civ1_V ;
     472            Data.Civ1_FF = Param.Civ1_FF ;
     473        end
     474        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     475       
    425476        Data.VarAttribute{nbvar+1}.Role='vector_x';
    426477        Data.VarAttribute{nbvar+2}.Role='vector_y';
    … …  
    436487        end
    437488        [Data.Civ1_SubRange,Data.Civ1_NbCentres,Data.Civ1_Coord_tps,Data.Civ1_U_tps,Data.Civ1_V_tps,tild,Ures, Vres,tild,FFres]=...
    438             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);
     489            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); 
     490           
    439491        Data.Civ1_U_smooth(ind_good)=Ures;
    440492        Data.Civ1_V_smooth(ind_good)=Vres;
    441493        Data.Civ1_FF(ind_good)=FFres;
    442         Data.CivStage=3;
     494        Data.CivStage=3;
     495
     496        if ~isfield (Param.ActionInput,'Civ2') && ~isfield (Param.ActionInput,'Civ3')
     497
     498            %%%%%%%%% modif fonction test %%%%%%%%%%%%%%%%%
     499            try
     500                [Data.Xmid, Data.Ymid, Data.Uphys, Data.Vphys, Data.Zphys, ...
     501                    Data.Yphys, Data.Xphys, Data.Error] = getPhysValues(Rangx, ...
     502                    Rangy, Npx, Npy, Data.Civ1_X, Data.Civ1_Y, Data.Civ1_U_smooth, ...
     503                    Data.Civ1_V_smooth, XmlData);
     504            catch % si on a pas Rangx, Rangy, Npx, Npy on ouvre les images pour les recalculer
     505                try
     506                    ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(ifield),[],j1_series_Civ1(ifield));
     507                    [A{1},VideoObject_A] = read_image(ImageName_A,FileType_A,VideoObject_A,FrameIndex_A_Civ1(ifield));
     508                    ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(ifield),[],j2_series_Civ1(ifield));
     509                    [A{2},VideoObject_B] = read_image(ImageName_B,FileType_B,VideoObject_B,FrameIndex_B_Civ1(ifield));
     510                catch ME
     511                    if ~isempty(ME.message)
     512                        disp_uvmat('ERROR', ['error reading input image: ' ME.message],checkrun)
     513                        return
     514                    end
     515                end
     516                [A,Rangx,Rangy]=phys_ima(A,XmlData,1);%transform image A in phys coordinates
     517                [Npy,Npx]=size(A{1});
     518                [Data.Xmid, Data.Ymid, Data.Uphys, Data.Vphys, Data.Zphys, ...
     519                    Data.Yphys, Data.Xphys, Data.Error] = getPhysValues(Rangx, ...
     520                    Rangy, Npx, Npy, Data.Civ1_X, Data.Civ1_Y, Data.Civ1_U_smooth, ...
     521                    Data.Civ1_V_smooth, XmlData);
     522            end             
     523            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     524   
     525            if ~isempty(errormsg)
     526                disp_uvmat('ERROR',errormsg,checkrun)
     527                return
     528            end
     529        end
    443530               
    444531    end
     532
     533   
     534   
    445535   
    446536    %% Civ2
    … …  
    553643        Data.Civ2_F=reshape(F,[],1);
    554644        Data.CivStage=Data.CivStage+1;
     645
     646        %%%%%%%%%%%%%%%% ajout fonction test %%%%%%%
     647        if ~isfield (Param.ActionInput,'Civ3') && ~isfield(Param.ActionInput,'Patch2')
     648            [Data.Xmid, Data.Ymid, Data.Uphys, Data.Vphys, Data.Zphys, ...
     649                Data.Yphys, Data.Xphys, Data.Error] = getPhysValues(Rangx, ...
     650                Rangy, Npx, Npy, Data.Civ2_X, Data.Civ2_Y, Data.Civ2_U, ...
     651                Data.Civ2_V, XmlData);
     652   
     653            if ~isempty(errormsg)
     654                disp_uvmat('ERROR',errormsg,checkrun)
     655                return
     656            end
     657        end
     658        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    555659    end
    556660   
    … …  
    612716        Data.Civ2_FF(ind_good)=FFres;
    613717        Data.CivStage=Data.CivStage+1;
     718
     719
     720        if ~isfield (Param.ActionInput,'Civ3')
     721            [Data.Xmid, Data.Ymid, Data.Uphys, Data.Vphys, Data.Zphys, ...
     722                Data.Yphys, Data.Xphys, Data.Error] = getPhysValues(Rangx, ...
     723                Rangy, Npx, Npy, Data.Civ2_X, Data.Civ2_Y, Data.Civ2_U_smooth, ...
     724                Data.Civ2_V_smooth, XmlData);
     725   
     726            if ~isempty(errormsg)
     727                disp_uvmat('ERROR',errormsg,checkrun)
     728                return
     729            end
     730        end
     731
    614732    end
    615733   
    … …  
    769887        Data.CivStage=Data.CivStage+1;
    770888           
    771          % get z from u and v (displacements)     
    772         Data.Xmid=Rangx(1)+(Rangx(2)-Rangx(1))*(Data.Civ3_X-0.5)/(Npx-1);%temporary coordinate (velocity taken at the point middle from imgae 1 and 2)
    773         Data.Ymid=Rangy(2)+(Rangy(1)-Rangy(2))*(Data.Civ3_Y-0.5)/(Npy-1);%temporary coordinate (velocity taken at the point middle from imgae 1 and 2)
    774         Data.Uphys=Data.Civ3_U_smooth*(Rangx(2)-Rangx(1))/(Npx-1);
    775         Data.Vphys=Data.Civ3_V_smooth*(Rangy(1)-Rangy(2))/(Npy-1);
    776         [Data.Zphys,Data.Xphys,Data.Yphys,Data.Error]=shift2z(Data.Xmid,Data.Ymid,Data.Uphys,Data.Vphys,XmlData); %Data.Xphys and Data.Xphys are real coordinate (geometric correction more accurate than xtemp/ytemp)
     889       
     890        [Data.Xmid, Data.Ymid, Data.Uphys, Data.Vphys, Data.Zphys, ...
     891            Data.Yphys, Data.Xphys, Data.Error] = getPhysValues(Rangx, ...
     892            Rangy, Npx, Npy, Data.Civ3_X, Data.Civ3_Y, Data.Civ3_U_smooth, ...
     893            Data.Civ3_V_smooth, XmlData);
     894
    777895        if ~isempty(errormsg)
    778896            disp_uvmat('ERROR',errormsg,checkrun)
    … …  
    784902   
    785903    %% write result in a netcdf file if requested
     904
    786905    if LSM ~= 1 % store all data
    787         if exist('ncfile','var')
    788             errormsg=struct2nc(ncfile,Data);
     906        if exist('ncfile2','var')
     907            errormsg=struct2nc(ncfile2,Data);
    789908            if isempty(errormsg)
    790                 disp([ncfile ' written'])
     909                disp([ncfile2 ' written'])
    791910            else
    792911                disp(errormsg)
    … …  
    794913        end
    795914    else
     915       
     916       
    796917       % store only phys data
    797         Data_light.ListVarName={'Xphys','Yphys','Zphys','Civ3_C','DX','DY','Error'};
     918        Data_light.ListVarName={'Xphys','Yphys','Zphys','Civ_C','DX','DY','Error'};
    798919        Data_light.VarDimName={'nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3'};
    799920        Data_light.VarAttribute{1}.Role='coord_x';
    … …  
    802923         Data_light.VarAttribute{5}.Role='vector_x';
    803924         Data_light.VarAttribute{6}.Role='vector_y';
    804         ind_good=find(Data.Civ3_FF==0);
    805         Data_light.Zphys=Data.Zphys(ind_good);
    806         Data_light.Yphys=Data.Yphys(ind_good);
    807         Data_light.Xphys=Data.Xphys(ind_good);
    808         Data_light.Civ3_C=Data.Civ3_C(ind_good);
    809         Data_light.DX=Data.Uphys(ind_good);
    810         Data_light.DY=Data.Vphys(ind_good);
    811         Data_light.Error=Data.Error(ind_good);
     925       
     926         % vrifie la dernire passe effectue
     927        if isfield (Param.ActionInput,'Civ3')
     928            try
     929                ind_good=find(Data.Civ3_FF==0);
     930                Data_light.Civ_C=Data.Civ3_C(ind_good);
     931            catch
     932                Data_light.Civ_C=Data.Civ3_C;
     933            end
     934        elseif isfield (Param.ActionInput,'Civ2')
     935            try
     936                ind_good=find(Data.Civ2_FF==0);
     937                Data_light.Civ_C=Data.Civ2_C(ind_good);
     938            catch
     939                Data_light.Civ_C=Data.Civ2_C;
     940            end
     941        elseif isfield (Param.ActionInput,'Civ1')
     942            try
     943                ind_good=find(Data.Civ1_FF==0);
     944                Data_light.Civ_C=Data.Civ1_C(ind_good);
     945            catch
     946                Data_light.Civ_C=Data.Civ1_C;
     947            end
     948        end
     949       
     950        try
     951            Data_light.Zphys=Data.Zphys(ind_good);
     952            Data_light.Yphys=Data.Yphys(ind_good);
     953            Data_light.Xphys=Data.Xphys(ind_good);
     954            Data_light.DX=Data.Uphys(ind_good);
     955            Data_light.DY=Data.Vphys(ind_good);
     956            Data_light.Error=Data.Error(ind_good);
     957        catch
     958            Data_light.Zphys=Data.Zphys;
     959            Data_light.Yphys=Data.Yphys;
     960            Data_light.Xphys=Data.Xphys;
     961            Data_light.DX=Data.Uphys;
     962            Data_light.DY=Data.Vphys;
     963            Data_light.Error=Data.Error;
     964        end
     965
     966
    812967       if exist('ncfile2','var')
    813968            errormsg=struct2nc(ncfile2,Data_light);
    … …  
    827982end
    828983
    829 
     984end
    830985
    831986% 'civ': function piv.m adapted from PIVlab http://pivlab.blogspot.com/
    … …  
    10681223end
    10691224result_conv=result_conv*corrmax/(255*sum_square);% keep the last correlation matrix for output
     1225end
    10701226
    10711227%------------------------------------------------------------------------
    … …  
    10981254end
    10991255vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
     1256end
    11001257
    11011258%------------------------------------------------------------------------
    … …  
    11411298end
    11421299vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
     1300end
    11431301
    11441302%'RUN_FIX': function for fixing velocity fields:
    … …  
    11891347end
    11901348
     1349end
    11911350
    11921351%------------------------------------------------------------------------
    … …  
    12501409        j2_series=ind2*ones(1,size(i_series,2));
    12511410        check_bounds=zeros(size(i1_series));% no limitations due to min-max indices
     1411    end
    12521412end
    12531413
    … …  
    12581418% ymid+ v/2: set of apparent phys y coordinates in the ref plane, image B
    12591419% XmlData: content of the xml files containing geometric calibration parameters
     1420
     1421
    12601422function [z,Xphy,Yphy,Error]=shift2z(xmid, ymid, u, v,XmlData)
    12611423z=0;
    … …  
    13201482Xphy=mean(xnew,1);
    13211483Yphy=mean(ynew,1);
    1322 
    1323 
    1324 
    1325 
    1326 
     1484end
     1485
     1486function [Xmid, Ymid, Uphys, Vphys, Zphys, Yphys, Xphys, Error] ...
     1487    = getPhysValues(Rangx, Rangy, Npx, Npy, Data_Civ_X, Data_Civ_Y, ...
     1488    Data_Civ_U_smooth, Data_Civ_V_smooth, XmlData)
     1489   
     1490    % get z from u and v (displacements)     
     1491    Xmid=Rangx(1)+(Rangx(2)-Rangx(1))*(Data_Civ_X-0.5)/(Npx-1);%temporary coordinate (velocity taken at the point middle from imgae 1 and 2)
     1492    Ymid=Rangy(2)+(Rangy(1)-Rangy(2))*(Data_Civ_Y-0.5)/(Npy-1);%temporary coordinate (velocity taken at the point middle from imgae 1 and 2)
     1493    Uphys=Data_Civ_U_smooth*(Rangx(2)-Rangx(1))/(Npx-1);
     1494    Vphys=Data_Civ_V_smooth*(Rangy(1)-Rangy(2))/(Npy-1);
     1495    [Zphys,Xphys,Yphys,Error]=shift2z(Xmid,Ymid,Uphys,Vphys,XmlData); %Data.Xphys and Data.Xphys are real coordinate (geometric correction more accurate than xtemp/ytempy
     1496
     1497end
Note: See TracChangeset for help on using the changeset viewer.