Ignore:
Timestamp:
Feb 6, 2015, 7:55:10 PM (6 years ago)
Author:
sommeria
Message:

bug_repaired

File:
1 edited

Legend:

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

    r862 r864  
    4141Data=[];
    4242errormsg='';
    43 
    4443%% set the input elements needed on the GUI series when the action is selected in the menu ActionName or InputTable refreshed
    4544if isstruct(Param) && isequal(Param.Action.RUN,0)% function activated from the GUI series but not RUN
     
    5049    path_series=fileparts(which('series'));
    5150    addpath(fullfile(path_series,'series'))
    52    % Data=civ_input(Param);% introduce the civ parameters using the GUI civ_input
    53    Data=stereo_input(Param)
     51    Data=stereo_input(Param);% introduce the civ parameters using the GUI civ_input
    5452    if isempty(Data)
    5553        Data=Param;% if  civ_input has been cancelled, keep previous parameters
     
    263261end
    264262
    265 
     263tic
    266264%%%%% MAIN LOOP %%%%%%
    267265for ifield=1:NbField
     
    271269        break
    272270    end
     271    % variable for light saving or not.
     272       LSM=Param.ActionInput.CheckLSM;
     273       
    273274    Civ1Dir=OutputDir;
    274275
    275         ncfile=fullfile_uvmat(RootPath_A,Civ1Dir,RootFile_A,'.nc',NomTypeNc,i2_series_Civ1(ifield),[],...
     276        ncfile=fullfile_uvmat(RootPath_A,Civ1Dir,[RootFile_A,'_All'],'.nc',NomTypeNc,i2_series_Civ1(ifield),[],...
    276277            j1_series_Civ1(ifield),j2_series_Civ1(ifield));
     278       
     279       
     280        ncfile2=fullfile_uvmat(RootPath_A,Civ1Dir,[RootFile_A,'_Light'],'.nc',NomTypeNc,i2_series_Civ1(ifield),[],...
     281            j1_series_Civ1(ifield),j2_series_Civ1(ifield));
     282       
    277283    %% Civ1
     284   
    278285    % if Civ1 computation is requested
    279286    if isfield (Param.ActionInput,'Civ1')
     
    294301        PhysImageA=fullfile_uvmat(RootPath_A,Civ1Dir,RootFile_A,'.png','_1a',i1_series_Civ1(ifield),[],1);
    295302        PhysImageB=fullfile_uvmat(RootPath_A,Civ1Dir,RootFile_A,'.png','_1a',i1_series_Civ1(ifield),[],2);
     303        if LSM ~= 1
    296304        imwrite(A{1},PhysImageA)
    297305        imwrite(A{2},PhysImageB)
     306        end
     307       
    298308        par_civ1.ImageA=A{1};
    299309        par_civ1.ImageB=A{2};
     
    338348       
    339349        % calculate velocity data (y and v in indices, reverse to y component)
    340         [xtable ytable utable vtable ctable F result_conv errormsg] = civ (par_civ1);
     350        [xtable, ytable, utable, vtable, ctable, F, result_conv, errormsg] = civ (par_civ1);
    341351        Data.Civ1_X=reshape(xtable,[],1);
    342352        Data.Civ1_Y=reshape(par_civ1.ImageHeight-ytable+1,[],1);
     
    346356        Data.Civ1_C=reshape(ctable,[],1);
    347357        Data.Civ1_F=reshape(F,[],1);
    348 %         Data.Xphys=Rangx(1)+(Rangx(2)-Rangx(1))*(Data.Civ1_X-0.5)/(Npx-1);
    349 %         Data.Yphys=Rangy(1)+(Rangy(2)-Rangy(1))*(Data.Civ1_Y-0.5)/(Npy-1);
    350 %         
    351 %         U=Data.Civ1_U*(Rangx(2)-Rangx(1))/(Npx-1);
    352 %         V=Data.Civ1_V*(Rangy(2)-Rangy(1))/(Npy-1);
    353 %         [Data.Zphys,Data.Civ1_E]=shift2z(Data.Xphys,Data.Yphys,U,V,XmlData);
    354 %         if ~isempty(errormsg)
    355 %             disp_uvmat('ERROR',errormsg,checkrun)
    356 %             return
    357 %         end
     358
    358359    end
    359360   
     
    425426        Data.Civ1_FF(ind_good)=FFres;
    426427        Data.CivStage=3;
     428             
     429       
     430       
     431       
     432%               
     433%          % get z from u and v (displacements)
     434%       
     435%         tempXmid=Rangx(1)+(Rangx(2)-Rangx(1))*(Data.Civ1_X-0.5)/(Npx-1);%temporary coordinate (velocity taken at the point middle from imgae 1 and 2)
     436%         tempYmid=Rangy(1)+(Rangy(2)-Rangy(1))*(Data.Civ1_Y-0.5)/(Npy-1);%temporary coordinate (velocity taken at the point middle from imgae 1 and 2)
     437%         tempUphys=Data.Civ1_U_smooth*(Rangx(2)-Rangx(1))/(Npx-1);
     438%         tempVphys=Data.Civ1_V_smooth*(Rangy(2)-Rangy(1))/(Npy-1);
     439%         [tempZphys,tempXphys,tempYphys,tempCiv3_E]=shift2z(tempXmid,tempYmid,tempUphys,tempVphys,XmlData); %Data.Xphys and Data.Xphys are real coordinate (geometric correction more accurate than xtemp/ytemp)
     440%         temp=find(Data.Civ1_FF~=0);
     441%         tempXmid(temp)=[];
     442%         tempYmid(temp)=[];
     443%         tempZphys(temp)=[];
     444%             
    427445    end
    428446   
     
    446464            j2=j2_series_Civ2(ifield);
    447465        end
    448         par_civ2.ImageWidth=size(par_civ2.ImageA,2);%FileInfo_A.Width;
    449         par_civ2.ImageHeight=size(par_civ2.ImageA,1);%FileInfo_A.Height;
     466        par_civ2.ImageWidth=size(par_civ2.ImageA,2);
     467        par_civ2.ImageHeight=size(par_civ2.ImageA,1);
     468       
    450469        if isfield(par_civ2,'Grid')% grid points set as input file
    451470            if ischar(par_civ2.Grid)%read the grid file if the input is a file name
     
    461480            par_civ2.Grid(:,1)=reshape(GridX,[],1);
    462481            par_civ2.Grid(:,2)=reshape(GridY,[],1);
     482           
     483           
    463484        end
    464485        Shiftx=zeros(size(par_civ2.Grid,1),1);% shift expected from civ1 data
     
    508529        end
    509530        % calculate velocity data (y and v in indices, reverse to y component)
    510         [xtable ytable utable vtable ctable F] = civ (par_civ2);
     531        [xtable, ytable, utable, vtable, ctable, F] = civ (par_civ2);
    511532        list_param=(fieldnames(Param.ActionInput.Civ2))';
    512533        Civ2_param=regexprep(list_param,'^.+','Civ2_$0');% insert 'Civ2_' before  each string in list_param
     
    523544       
    524545        nbvar=numel(Data.ListVarName);
    525         Data.ListVarName=[Data.ListVarName {'Civ2_X','Civ2_Y','Civ2_U','Civ2_V','Civ2_F','Civ2_C','Xphys','Yphys','Zphys','Civ2_E'}];%  cell array containing the names of the fields to record
    526         Data.VarDimName=[Data.VarDimName {'nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2'}];
     546        Data.ListVarName=[Data.ListVarName {'Civ2_X','Civ2_Y','Civ2_U','Civ2_V','Civ2_F','Civ2_C'}];%  cell array containing the names of the fields to record
     547        Data.VarDimName=[Data.VarDimName {'nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2'}];
    527548        Data.VarAttribute{nbvar+1}.Role='coord_x';
    528549        Data.VarAttribute{nbvar+2}.Role='coord_y';
     
    598619       
    599620           
     621 
     622       
     623    end
     624   
     625       
     626    %% Civ3
     627   
     628    if isfield (Param.ActionInput,'Civ3')
     629        par_civ3=Param.ActionInput.Civ3;
     630        par_civ3.ImageA=par_civ1.ImageA;
     631        par_civ3.ImageB=par_civ1.ImageB;
     632        %         if ~isfield(Param.Civ1,'ImageA')
     633%         i1=i1_series_Civ3(ifield);
     634%         i2=i1;
     635%         if ~isempty(i2_series_Civ3)
     636%             i2=i2_series_Civ3(ifield);
     637%         end
     638%         j1=1;
     639%         if ~isempty(j1_series_Civ3)
     640%             j1=j1_series_Civ3(ifield);
     641%         end
     642%         j2=j1;
     643%         if ~isempty(j2_series_Civ3)
     644%             j2=j2_series_Civ3(ifield);
     645%         end
     646        par_civ3.ImageWidth=size(par_civ3.ImageA,2);
     647        par_civ3.ImageHeight=size(par_civ3.ImageA,1);
     648       
     649        if isfield(par_civ3,'Grid')% grid points set as input file
     650            if ischar(par_civ3.Grid)%read the grid file if the input is a file name
     651                par_civ3.Grid=dlmread(par_civ3.Grid);
     652                par_civ3.Grid(1,:)=[];%the first line must be removed (heading in the grid file)
     653            end
     654        else% automatic grid
     655            minix=floor(par_civ3.Dx/2)-0.5;
     656            maxix=minix+par_civ3.Dx*floor((par_civ3.ImageWidth-1)/par_civ3.Dx);
     657            miniy=floor(par_civ3.Dy/2)-0.5;
     658            maxiy=minix+par_civ3.Dy*floor((par_civ3.ImageHeight-1)/par_civ3.Dy);
     659            [GridX,GridY]=meshgrid(minix:par_civ3.Dx:maxix,miniy:par_civ3.Dy:maxiy);
     660            par_civ3.Grid(:,1)=reshape(GridX,[],1);
     661            par_civ3.Grid(:,2)=reshape(GridY,[],1);
     662           
     663           
     664        end
     665        Shiftx=zeros(size(par_civ3.Grid,1),1);% shift expected from civ2 data
     666        Shifty=zeros(size(par_civ3.Grid,1),1);
     667        nbval=zeros(size(par_civ3.Grid,1),1);
     668        if par_civ3.CheckDeformation
     669            DUDX=zeros(size(par_civ3.Grid,1),1);
     670            DUDY=zeros(size(par_civ3.Grid,1),1);
     671            DVDX=zeros(size(par_civ3.Grid,1),1);
     672            DVDY=zeros(size(par_civ3.Grid,1),1);
     673        end
     674        NbSubDomain=size(Data.Civ2_SubRange,3);
     675        % get the guess from patch2
     676        for isub=1:NbSubDomain% for each sub-domain of Patch2
     677            nbvec_sub=Data.Civ2_NbCentres(isub);% nbre of Civ2 vectors in the subdomain
     678            ind_sel=find(par_civ3.Grid(:,1)>=Data.Civ2_SubRange(1,1,isub) & par_civ3.Grid(:,1)<=Data.Civ2_SubRange(1,2,isub) &...
     679                par_civ3.Grid(:,2)>=Data.Civ2_SubRange(2,1,isub) & par_civ3.Grid(:,2)<=Data.Civ2_SubRange(2,2,isub));
     680            epoints = par_civ3.Grid(ind_sel,:);% coordinates of interpolation sites
     681            ctrs=Data.Civ2_Coord_tps(1:nbvec_sub,:,isub) ;%(=initial points) ctrs
     682            nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
     683            EM = tps_eval(epoints,ctrs);
     684            Shiftx(ind_sel)=Shiftx(ind_sel)+EM*Data.Civ2_U_tps(1:nbvec_sub+3,isub);
     685            Shifty(ind_sel)=Shifty(ind_sel)+EM*Data.Civ2_V_tps(1:nbvec_sub+3,isub);
     686            if par_civ3.CheckDeformation
     687                [EMDX,EMDY] = tps_eval_dxy(epoints,ctrs);%2D matrix of distances between extrapolation points epoints and spline centres (=site points) ctrs
     688                DUDX(ind_sel)=DUDX(ind_sel)+EMDX*Data.Civ2_U_tps(1:nbvec_sub+3,isub);
     689                DUDY(ind_sel)=DUDY(ind_sel)+EMDY*Data.Civ2_U_tps(1:nbvec_sub+3,isub);
     690                DVDX(ind_sel)=DVDX(ind_sel)+EMDX*Data.Civ2_V_tps(1:nbvec_sub+3,isub);
     691                DVDY(ind_sel)=DVDY(ind_sel)+EMDY*Data.Civ2_V_tps(1:nbvec_sub+3,isub);
     692            end
     693        end
     694        mask='';
     695        if par_civ3.CheckMask&&~isempty(par_civ3.Mask)&& ~strcmp(maskname,par_civ3.Mask)% mask exist, not already read in Civ2
     696            mask=imread(par_civ3.Mask);
     697        end
     698        ibx2=ceil(par_civ3.CorrBoxSize(1)/2);
     699        iby2=ceil(par_civ3.CorrBoxSize(2)/2);
     700%         par_civ3.SearchBoxSize(1)=2*ibx2+9;% search ara +-4 pixels around the guess
     701%         par_civ3.SearchBoxSize(2)=2*iby2+9;
     702        par_civ3.SearchBoxShift=[Shiftx(nbval>=1)./nbval(nbval>=1) Shifty(nbval>=1)./nbval(nbval>=1)];
     703        par_civ3.Grid=[par_civ3.Grid(nbval>=1,1)-par_civ3.SearchBoxShift(:,1)/2 par_civ3.Grid(nbval>=1,2)-par_civ3.SearchBoxShift(:,2)/2];% grid taken at the extrapolated origin of the displacement vectors
     704        if par_civ3.CheckDeformation
     705            par_civ3.DUDX=DUDX./nbval;
     706            par_civ3.DUDY=DUDY./nbval;
     707            par_civ3.DVDX=DVDX./nbval;
     708            par_civ3.DVDY=DVDY./nbval;
     709        end
     710        % calculate velocity data (y and v in indices, reverse to y component)
     711        [xtable, ytable, utable, vtable, ctable, F] = civ (par_civ3);
     712        list_param=(fieldnames(Param.ActionInput.Civ3))';
     713        Civ3_param=regexprep(list_param,'^.+','Civ3_$0');% insert 'Civ3_' before  each string in list_param
     714        Civ3_param=[{'Civ3_ImageA','Civ3_ImageB','Civ3_Time','Civ3_Dt'} Civ3_param]; %insert the names of the two input images
     715        %indicate the values of all the global attributes in the output data
     716        Data.Civ3_ImageA=ImageName_A;
     717        Data.Civ3_ImageB=ImageName_B;
     718        Data.Civ3_Time=(time(i2+1,j2+1)+time(i1+1,j1+1))/2;
     719        Data.Civ3_Dt=0;
     720        for ilist=1:length(list_param)
     721            Data.(Civ3_param{4+ilist})=Param.ActionInput.Civ3.(list_param{ilist});
     722        end
     723        Data.ListGlobalAttribute=[Data.ListGlobalAttribute Civ3_param];
     724       
     725        nbvar=numel(Data.ListVarName);
     726        Data.ListVarName=[Data.ListVarName {'Civ3_X','Civ3_Y','Civ3_U','Civ3_V','Civ3_F','Civ3_C','Xphys','Yphys','Zphys','Civ3_E'}];%  cell array containing the names of the fields to record
     727        Data.VarDimName=[Data.VarDimName {'nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3'}];
     728        Data.VarAttribute{nbvar+1}.Role='coord_x';
     729        Data.VarAttribute{nbvar+2}.Role='coord_y';
     730        Data.VarAttribute{nbvar+3}.Role='vector_x';
     731        Data.VarAttribute{nbvar+4}.Role='vector_y';
     732        Data.VarAttribute{nbvar+5}.Role='warnflag';
     733        Data.Civ3_X=reshape(xtable,[],1);
     734        Data.Civ3_Y=reshape(size(par_civ3.ImageA,1)-ytable+1,[],1);
     735        Data.Civ3_U=reshape(utable,[],1);
     736        Data.Civ3_V=reshape(-vtable,[],1);
     737        Data.Civ3_C=reshape(ctable,[],1);
     738        Data.Civ3_F=reshape(F,[],1);
     739        Data.CivStage=Data.CivStage+1;
     740end
     741   
     742    %% Fix3
     743    if isfield (Param.ActionInput,'Fix3')
     744        ListFixParam=fieldnames(Param.ActionInput.Fix3);
     745        for ilist=1:length(ListFixParam)
     746            ParamName=ListFixParam{ilist};
     747            ListName=['Fix3_' ParamName];
     748            eval(['Data.ListGlobalAttribute=[Data.ListGlobalAttribute ''' ParamName '''];'])
     749            eval(['Data.' ListName '=Param.ActionInput.Fix3.' ParamName ';'])
     750        end
     751        if check_civx
     752            if ~isfield(Data,'fix3')
     753                Data.ListGlobalAttribute=[Data.ListGlobalAttribute 'fix3'];
     754                Data.fix3=1;
     755                Data.ListVarName=[Data.ListVarName {'vec3_FixFlag'}];
     756                Data.VarDimName=[Data.VarDimName {'nb_vectors3'}];
     757            end
     758            Data.vec_FixFlag=fix(Param.Fix3,Data.vec3_F,Data.vec3_C,Data.vec3_U,Data.vec3_V,Data.vec3_X,Data.vec3_Y);
     759        else
     760            Data.ListVarName=[Data.ListVarName {'Civ3_FF'}];
     761            Data.VarDimName=[Data.VarDimName {'nb_vec_3'}];
     762            nbvar=length(Data.ListVarName);
     763            Data.VarAttribute{nbvar}.Role='errorflag';
     764            Data.Civ3_FF=double(fix(Param.ActionInput.Fix3,Data.Civ3_F,Data.Civ3_C,Data.Civ3_U,Data.Civ3_V));
     765            Data.CivStage=Data.CivStage+1;
     766        end
     767       
     768    end
     769
     770   
     771     %% Patch3
     772    if isfield (Param.ActionInput,'Patch3')
     773        Data.ListGlobalAttribute=[Data.ListGlobalAttribute {'Patch3_Rho','Patch3_Threshold','Patch3_SubDomain'}];
     774        Data.Patch3_FieldSmooth=Param.ActionInput.Patch3.FieldSmooth;
     775        Data.Patch3_MaxDiff=Param.ActionInput.Patch3.MaxDiff;
     776        Data.Patch3_SubDomainSize=Param.ActionInput.Patch3.SubDomainSize;
     777        nbvar=length(Data.ListVarName);
     778        Data.ListVarName=[Data.ListVarName {'Civ3_U_smooth','Civ3_V_smooth','Civ3_SubRange','Civ3_NbCentres','Civ3_Coord_tps','Civ3_U_tps','Civ3_V_tps','Xmid','Ymid','Uphys','Vphys'}];
     779        Data.VarDimName=[Data.VarDimName {'nb_vec_3','nb_vec_3',{'nb_coord','nb_bounds','nb_subdomain_3'},{'nb_subdomain_3'},...
     780            {'nb_tps_3','nb_coord','nb_subdomain_3'},{'nb_tps_3','nb_subdomain_3'},{'nb_tps_3','nb_subdomain_3'},'nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3'}];
     781       
     782        Data.VarAttribute{nbvar+1}.Role='vector_x';
     783        Data.VarAttribute{nbvar+2}.Role='vector_y';
     784        Data.VarAttribute{nbvar+5}.Role='coord_tps';
     785        Data.VarAttribute{nbvar+6}.Role='vector_x';
     786        Data.VarAttribute{nbvar+7}.Role='vector_y';
     787        Data.Civ3_U_smooth=zeros(size(Data.Civ3_X));
     788        Data.Civ3_V_smooth=zeros(size(Data.Civ3_X));
     789        if isfield(Data,'Civ3_FF')
     790            ind_good=find(Data.Civ3_FF==0);
     791        else
     792            ind_good=1:numel(Data.Civ3_X);
     793        end
     794        [Data.Civ3_SubRange,Data.Civ3_NbCentres,Data.Civ3_Coord_tps,Data.Civ3_U_tps,Data.Civ3_V_tps,tild,Ures, Vres,tild,FFres]=...
     795            filter_tps([Data.Civ3_X(ind_good) Data.Civ3_Y(ind_good)],Data.Civ3_U(ind_good),Data.Civ3_V(ind_good),[],Data.Patch3_SubDomainSize,Data.Patch3_FieldSmooth,Data.Patch3_MaxDiff);
     796        Data.Civ3_U_smooth(ind_good)=Ures;
     797        Data.Civ3_V_smooth(ind_good)=Vres;
     798        Data.Civ3_FF(ind_good)=FFres;
     799        Data.CivStage=Data.CivStage+1;
     800       
     801           
    600802         % get z from u and v (displacements)
    601803       
    602         Data.Xphys=Rangx(1)+(Rangx(2)-Rangx(1))*(Data.Civ2_X-0.5)/(Npx-1);
    603         Data.Yphys=Rangy(1)+(Rangy(2)-Rangy(1))*(Data.Civ2_Y-0.5)/(Npy-1);
    604         U=Data.Civ2_U_smooth*(Rangx(2)-Rangx(1))/(Npx-1);
    605         V=Data.Civ2_V_smooth*(Rangy(2)-Rangy(1))/(Npy-1);
    606         [Data.Zphys,Data.Civ2_E]=shift2z(Data.Xphys,Data.Yphys,U,V,XmlData);
     804        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)
     805        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)
     806        Data.Uphys=Data.Civ3_U_smooth*(Rangx(2)-Rangx(1))/(Npx-1);
     807        Data.Vphys=Data.Civ3_V_smooth*(Rangy(2)-Rangy(1))/(Npy-1);
     808        [Data.Zphys,Data.Xphys,Data.Yphys,Data.Civ3_E]=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)
    607809        if ~isempty(errormsg)
    608810            disp_uvmat('ERROR',errormsg,checkrun)
     
    612814    end
    613815   
     816   
     817   
     818   
    614819    %% write result in a netcdf file if requested
    615     if exist('ncfile','var')
    616         errormsg=struct2nc(ncfile,Data);
    617         if isempty(errormsg)
    618             disp([ncfile ' written'])
    619         else
    620             disp(errormsg)
    621         end
    622     end
    623 end
     820    if LSM ~= 1 % store all data
     821        if exist('ncfile','var')
     822            errormsg=struct2nc(ncfile,Data);
     823            if isempty(errormsg)
     824                disp([ncfile ' written'])
     825            else
     826                disp(errormsg)
     827            end
     828        end
     829    else
     830       % store only phys data
     831        Data_light.ListVarName={'Xphys','Yphys','Zphys','Civ3_C','Xmid','Ymid','Uphys','Vphys'};
     832        Data_light.VarDimName={'nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3'};
     833        temp=find(Data.Civ3_FF==0);
     834        Data_light.Zphys=Data.Zphys(temp);
     835        Data_light.Yphys=Data.Yphys(temp);
     836        Data_light.Xphys=Data.Xphys(temp);
     837        Data_light.Civ3_C=Data.Civ3_C(temp);
     838        Data_light.Xmid=Data.Xmid(temp);
     839        Data_light.Ymid=Data.Ymid(temp);
     840        Data_light.Uphys=Data.Uphys(temp);
     841        Data_light.Vphys=Data.Vphys(temp);
     842       
     843       if exist('ncfile2','var')
     844            errormsg=struct2nc(ncfile2,Data_light);
     845            if isempty(errormsg)
     846                disp([ncfile2 ' written'])
     847            else
     848                disp(errormsg)
     849            end
     850       end
     851       
     852    end
     853   end
     854
     855  toc
     856
     857
    624858
    625859% 'civ': function piv.m adapted from PIVlab http://pivlab.blogspot.com/
     
    7691003end
    7701004% vector=[0 0];%default
     1005
    7711006for ivec=1:nbvec
    7721007    iref=round(par_civ.Grid(ivec,1)+0.5);% xindex on the image A for the middle of the correlation box
    7731008    jref=round(par_civ.ImageHeight-par_civ.Grid(ivec,2)+0.5);% yindex on the image B for the middle of the correlation box
     1009   
    7741010    %if ~(checkmask && par_civ.Mask(jref,iref)<=20) %velocity not set to zero by the black mask
    7751011    %         if jref-iby2<1 || jref+iby2>par_civ.ImageHeight|| iref-ibx2<1 || iref+ibx2>par_civ.ImageWidth||...
     
    8201056            image2_crop=interp2(image2_crop,xi,yi);
    8211057        end
    822         sum_square=(sum(sum(image1_crop.*image1_crop))+sum(sum(image2_crop.*image2_crop)))/2;
     1058        sum_square=(sum(sum(image1_crop.*image1_crop)));%+sum(sum(image2_crop.*image2_crop)))/2;
    8231059        %reference: Oliver Pust, PIV: Direct Cross-Correlation
    8241060        result_conv= conv2(image2_crop,flipdim(flipdim(image1_crop,2),1),'valid');
     
    8341070                    [vector,F(ivec)] = SUBPIX2DGAUSS (result_conv,x,y);
    8351071                end
    836 
    837 %                 if isfield(par_civ,'CheckDeformation')
    838 %                  utable(ivec)=shiftx(ivec);
    839 %                  vtable(ivec)=shifty(ivec);
    840 %                 else
    841                  utable(ivec)=vector(1)*mesh+shiftx(ivec);
     1072               
     1073
     1074                utable(ivec)=vector(1)*mesh+shiftx(ivec);
    8421075                vtable(ivec)=vector(2)*mesh+shifty(ivec);
    843 %                 end
    844 %                 xtable(ivec)=iref+utable(ivec)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2)
    845 %                 ytable(ivec)=jref+vtable(ivec)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge)
     1076
     1077               
     1078                xtable(ivec)=iref+utable(ivec)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2)
     1079                ytable(ivec)=jref+vtable(ivec)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge)
     1080               
    8461081                iref=round(xtable(ivec));% image index for the middle of the vector
    8471082                jref=round(ytable(ivec));
     
    10511286% ymid+ v/2: set of apparent phys y coordinates in the ref plane, image B
    10521287% XmlData: content of the xml files containing geometric calibration parameters
    1053 function [z,error]=shift2z(xmid, ymid, u, v,XmlData)
     1288function [z,Xphy,Yphy,error]=shift2z(xmid, ymid, u, v,XmlData)
    10541289z=0;
    10551290error=0;
     
    10591294Calib_A=XmlData{1}.GeometryCalib;
    10601295R=(Calib_A.R)';
    1061 x_a=xmid;
    1062 y_a=ymid;
     1296x_a=xmid- u/2;
     1297y_a=ymid- u/2;
    10631298z_a=R(7)*x_a+R(8)*y_a+Calib_A.Tx_Ty_Tz(1,3);
    10641299Xa=(R(1)*x_a+R(2)*y_a+Calib_A.Tx_Ty_Tz(1,1))./z_a;
     
    10781313Calib_B=XmlData{2}.GeometryCalib;
    10791314R=(Calib_B.R)';
    1080 x_b=xmid+ u;
    1081 y_b=ymid+ v;
     1315x_b=xmid+ u/2;
     1316y_b=ymid+ v/2;
    10821317z_b=R(7)*x_b+R(8)*y_b+Calib_B.Tx_Ty_Tz(1,3);
    10831318Xb=(R(1)*x_b+R(2)*y_b+Calib_B.Tx_Ty_Tz(1,1))./z_b;
     
    10961331Den=(Dxb-Dxa).*(Dxb-Dxa)+(Dyb-Dya).*(Dyb-Dya);
    10971332error=((Dyb-Dya).*u-(Dxb-Dxa).*v)./Den;
     1333
     1334xnew(1,:)=Dxa.*z+x_a;
     1335xnew(2,:)=Dxb.*z+x_b;
     1336ynew(1,:)=Dya.*z+y_a;
     1337ynew(2,:)=Dyb.*z+y_b;
     1338
     1339Xphy=mean(xnew,1); %on moyenne les 2 valeurs
     1340Yphy=mean(ynew,1);
    10981341z=((Dxb-Dxa).*u-(Dyb-Dya).*v)./Den;
    10991342
     1343
     1344
     1345
Note: See TracChangeset for help on using the changeset viewer.