Index: trunk/src/geometry_calib.m =================================================================== --- trunk/src/geometry_calib.m (revision 1147) +++ trunk/src/geometry_calib.m (revision 1148) @@ -982,4 +982,12 @@ end + + +%% read the current image, displayed in the GUI uvmat +huvmat=findobj(allchild(0),'Name','uvmat'); +UvData=get(huvmat,'UserData'); +A=UvData.Field.A;%currently displayed image +npxy=size(A); + %% initiate the grid in phys coordinates CalibData=get(handles.geometry_calib,'UserData');%get information stored on the GUI geometry_calib @@ -987,4 +995,10 @@ if isfield(CalibData,'grid') grid_input=CalibData.grid;%retrieve the previously used grid +else + %S=skewness(double(reshape(A,1,[]))); + A=double(A); + A=A-mean(mean(A)); + S=mean(mean(A.*A.*A))/(mean(mean(A.*A)))^1.5 + grid_input.CheckWhite=sign(S);%propose white markers if image skewness>0, black markers otherwise end [T,CalibData.grid,CalibData.grid.CheckWhite,CalibData.grid.FilterWindow]=create_grid(grid_input,'detect_grid');%display the GUI create_grid, read the set of phys coordinates T @@ -992,10 +1006,4 @@ X=[CalibData.grid.x_0 CalibData.grid.x_1 CalibData.grid.x_0 CalibData.grid.x_1]';%corner absissa in the phys coordinates (cm) Y=[CalibData.grid.y_0 CalibData.grid.y_0 CalibData.grid.y_1 CalibData.grid.y_1]';%corner ordinates in the phys coordinates (cm) - -%% read the current image, displayed in the GUI uvmat -huvmat=findobj(allchild(0),'Name','uvmat'); -UvData=get(huvmat,'UserData'); -A=UvData.Field.A;%currently displayed image -npxy=size(A); %% calculate transform matrices for plane projection: rectangle assumed to be viewed in perspective @@ -1034,4 +1042,7 @@ Rangx=DataOut.Coord_x;% x coordinates of first and last pixel centres in phys Rangy=DataOut.Coord_y;% y coordinates of first and last pixel centres in phys + + +%% inverse the image in the case of black lines if CalibData.grid.CheckWhite Amod=double(Amod);%case of white grid markers: will look for image maxima Index: trunk/src/series.m =================================================================== --- trunk/src/series.m (revision 1147) +++ trunk/src/series.m (revision 1148) @@ -2321,5 +2321,5 @@ % remove old Action options in the menu (keeping a menu length NbBuiltinAction+5; % nb_builtin_ACTION=nbre of functions always remaining in the initial menu + if length(ActionList)>NbBuiltinAction+5% nb_builtin_ACTION=nbre of functions always remaining in the initial menu nbremove=length(ActionList)-NbBuiltinAction-5; ActionList(NbBuiltinAction+1:end-5)=[]; @@ -2354,5 +2354,5 @@ %% check the current ActionPath to the selected function ActionPath=ActionPathList{ActionIndex}; % current recorded path -set(handles.ActionPath,'String',ActionPath); % show the path to the senlected function +set(handles.ActionPath,'String',ActionPath); % show the path to the selected function %% reinitialise the waitbar @@ -2361,5 +2361,5 @@ %% Put the first line of the selected Action fct as tooltip help try - [fid,errormsg] =fopen([ActionName '.m']); + [fid,errormsg] =fopen([fullfile(ActionPath,ActionName) '.m']); InputText=textscan(fid,'%s',1,'delimiter','\n'); fclose(fid); @@ -2563,21 +2563,32 @@ %% Impose the whole input file index range if requested -if isfield(ParamOut,'WholeIndexRange')&&isequal(ParamOut.WholeIndexRange,'on') +if ~isfield(ParamOut,'WholeIndexRange') + ParamOut.WholeIndexRange='off'; +end +if ~isfield(ParamOut,'WholeIndexRange_j') + ParamOut.WholeIndexRange_j='off'; +end + +if strcmp(ParamOut.WholeIndexRange,'on') MinIndex_i=get(handles.MinIndex_i,'Data'); - MinIndex_j=get(handles.MinIndex_j,'Data'); MaxIndex_i=get(handles.MaxIndex_i,'Data'); - MaxIndex_j=get(handles.MaxIndex_j,'Data'); set(handles.num_first_i,'String',num2str(MinIndex_i(1)))% set first as the min index (for the first line) set(handles.num_last_i,'String',num2str(MaxIndex_i(1)))% set last as the max index (for the first line) set(handles.num_incr_i,'String','1') +else + first_i=1;last_i=1; + if isfield(Param.IndexRange,'first_i') + first_i=Param.IndexRange.first_i; + last_i=Param.IndexRange.last_i; + end +end +if strcmp(ParamOut.WholeIndexRange,'on')||strcmp(ParamOut.WholeIndexRange_j,'on') + MinIndex_j=get(handles.MinIndex_j,'Data'); + MaxIndex_j=get(handles.MaxIndex_j,'Data'); set(handles.num_first_j,'String',num2str(MinIndex_j(1)))% set first as the min index (for the first line) set(handles.num_last_j,'String',num2str(MaxIndex_j(1)))% set last as the max index (for the first line) set(handles.num_incr_j,'String','1') else % check index ranges - first_i=1;last_i=1;first_j=1;last_j=1; - if isfield(Param.IndexRange,'first_i') - first_i=Param.IndexRange.first_i; - last_i=Param.IndexRange.last_i; - end + first_j=1;last_j=1; if isfield(Param.IndexRange,'first_j') first_j=Param.IndexRange.first_j; Index: trunk/src/series/aver_stat.m =================================================================== --- trunk/src/series/aver_stat.m (revision 1147) +++ trunk/src/series/aver_stat.m (revision 1148) @@ -136,8 +136,10 @@ Param=xml2struct(Param);% read Param as input file (batch case) checkrun=0; -end + WaitbarHandle=[]; +else hseries=findobj(allchild(0),'Tag','series'); RUNHandle=findobj(hseries,'Tag','RUN');%handle of RUN button in GUI series WaitbarHandle=findobj(hseries,'Tag','Waitbar');%handle of waitbar in GUI series +end %% define the directory for result file (with path=RootPath{1}) Index: trunk/src/series/civ_3D.m =================================================================== --- trunk/src/series/civ_3D.m (revision 1147) +++ trunk/src/series/civ_3D.m (revision 1148) @@ -1,9 +1,6 @@ -%'civ_3D': PIV function activated by the general GUI series -% --- call the sub-functions: -% civ: PIV function itself -% fix: removes false vectors after detection by various criteria -% filter_tps: make interpolation-smoothing +%'civ_3D': 3D PIV from image scan in a volume + %------------------------------------------------------------------------ -% function [Data,errormsg,result_conv]= civ_series(Param) +% function [Data,errormsg,result_conv]= civ_3D(Param) % %OUTPUT @@ -21,6 +18,6 @@ % if absent no file is produced, result in the output structure Data (test mode) % Param.ActionInput: substructure with the parameters provided by the GUI civ_input -% .Civ1: parameters for civ1 -% .Fix1: parameters for fix1 +% .Civ1: parameters for civ1cc +% .Fix1: parameters for detect_false1 % .Patch1: % .Civ2: for civ2 @@ -53,13 +50,10 @@ path_series=fileparts(which('series')); addpath(fullfile(path_series,'series')) - Data=civ_input_3D(Param);% introduce the civ parameters using the GUI civ_input - %Data=civ_input_App(Param);% introduce the civ parameters using the GUI - %civ_input %%%% TO MOVE - if isempty(Data) - Data=Param;% if civ_input has been cancelled, keep previous parameters - end + + Data=civ_input(Param) + Data.Program=mfilename;%gives the name of the current function Data.AllowInputSort='off';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default) - Data.WholeIndexRange='off';% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default) + Data.WholeIndexRange_j='on';% prescribes the file index ranges j from min to max (options 'off'/'on', 'off' by default) Data.NbSlice='off'; %nbre of slices ('off' by default) Data.VelType='off';% menu for selecting the velocity type (options 'off'/'one'/'two', 'off' by default) @@ -68,9 +62,9 @@ Data.ProjObject='off';%can use projection object(option 'off'/'on', Data.Mask='off';%can use mask option (option 'off'/'on', 'off' by default) - Data.OutputDirExt='.civ';%set the output dir extension + Data.OutputDirExt='.civ_3D';%set the output dir extension Data.OutputSubDirMode='last'; %select the last subDir in the input table as root of the output subdir name (option 'all'/'first'/'last', 'all' by default) Data.OutputFileMode='NbInput_i';% one output file expected per value of i index (used for waitbar) Data.CheckOverwriteVisible='on'; % manage the overwrite of existing files (default=1) - if isfield(Data.ActionInput,'PairIndices') && strcmp(Data.ActionInput.PairIndices.ListPairMode,'pair j1-j2') + if isfield(Data,'ActionInput') && isfield(Data.ActionInput,'PairIndices') && strcmp(Data.ActionInput.PairIndices.ListPairMode,'pair j1-j2') if isfield(Data.ActionInput.PairIndices,'ListPairCiv2') str_civ=Data.ActionInput.PairIndices.ListPairCiv2; @@ -90,51 +84,27 @@ end end - % estimated CPUTime - CPUtime_unit=0.01;%estimated time for a multiplication (in microsecond) - if isfield(Param.SeriesData,'FileInfo')&&isfield(Param.SeriesData.FileInfo{1},'Height')&&isfield(Param.SeriesData.FileInfo{1},'Width') - pixnbre=Param.SeriesData.FileInfo{1}.Height*Param.SeriesData.FileInfo{1}.Width; % total number of pxels for input images - CPUtime=0; - if isfield(Data.ActionInput,'Civ1') - %BoxSize=Data.ActionInput.Civ1.CorrBoxSize(1)*Data.ActionInput.Civ1.CorrBoxSize(2); - tic - testboxa=rand(Data.ActionInput.Civ1.CorrBoxSize(1),Data.ActionInput.Civ1.CorrBoxSize(2)); - testboxb=rand(Data.ActionInput.Civ1.SearchBoxSize(1),Data.ActionInput.Civ1.SearchBoxSize(2)); - anss=conv2(testboxa,testboxb); - CPUtime_unit=toc; - nb_box=pixnbre/(Data.ActionInput.Civ1.Dx*Data.ActionInput.Civ1.Dy); - %nbpos=Data.ActionInput.Civ1.SearchBoxSize-Data.ActionInput.Civ1.CorrBoxSize; - CPUtime=2*CPUtime_unit*nb_box%*BoxSize*nbpos(1)*nbpos(2);% adjustement factor 2 used - end - if isfield(Data.ActionInput,'Patch1') - CPUtime=2*CPUtime; - end - if isfield(Data.ActionInput,'Civ2') - tic - testboxa=rand(Data.ActionInput.Civ2.CorrBoxSize(1),Data.ActionInput.Civ2.CorrBoxSize(2)); - testboxb=rand(Data.ActionInput.Civ2.SearchBoxSize(1),Data.ActionInput.Civ2.SearchBoxSize(2)); - anss=conv2(testboxa,testboxb); - CPUtime_unit=toc; - nb_box=pixnbre/(Data.ActionInput.Civ2.Dx*Data.ActionInput.Civ2.Dy); - %BoxSize=Data.ActionInput.Civ2.CorrBoxSize(1)*Data.ActionInput.Civ2.CorrBoxSize(2); - %nbpos=Data.ActionInput.Civ2.SearchBoxSize-Data.ActionInput.Civ2.CorrBoxSize; - CPUtime=CPUtime+2*CPUtime_unit*nb_box;%*BoxSize*nbpos(1)*nbpos(2); - end - if isfield(Data.ActionInput,'Patch2') - CPUtime=(4/3)*CPUtime; - end - Data.CPUTime=ceil(CPUtime/6); % estimated CPU time per field pair in minute - Data.CPUTime=Data.CPUTime/10; % displqy CPU time with 1 digit beyond dot - end + + return end - -%% read input parameters from an xml file if input is a file name (batch mode) -checkrun=1; +%% END OF ENTERING INPUT PARAMETER MODE + +%% RUN MODE: read input parameters from an xml file if input is a file name (batch mode) if ischar(Param) Param=xml2struct(Param);% read Param as input file (batch case) checkrun=0; -end - -%% test input + RUNHandle=[]; +else + hseries=findobj(allchild(0),'Tag','series'); + RUNHandle=findobj(hseries,'Tag','RUN');%handle of RUN button in GUI series + WaitbarHandle=findobj(hseries,'Tag','Waitbar');%handle of waitbar in GUI series + checkrun=1; +end + +%% test input Param +if ~isfield(Param,'InputTable') + disp('ERROR: no input file entered') + return +end if ~isfield(Param,'ActionInput') disp_uvmat('ERROR','no parameter set for PIV',checkrun) @@ -142,145 +112,57 @@ end iview_A=0;%default values -NbField=1; -RUNHandle=[]; -CheckInputFile=isfield(Param,'InputTable');%= 1 in test use for TestCiv (no nc file involved) -CheckOutputFile=isfield(Param,'OutputSubDir');%= 1 in test use for TestPatch (no nc file produced) - -%% input files and indexing (skipped in Test mode) -if CheckInputFile - hseries=findobj(allchild(0),'Tag','series'); - RUNHandle=findobj(hseries,'Tag','RUN');%handle of RUN button in GUI series - WaitbarHandle=findobj(hseries,'Tag','Waitbar');%handle of waitbar in GUI series - MaxIndex_i=Param.IndexRange.MaxIndex_i; - MinIndex_i=Param.IndexRange.MinIndex_i; - MaxIndex_j=ones(size(MaxIndex_i));MinIndex_j=ones(size(MinIndex_i)); - if isfield(Param.IndexRange,'MaxIndex_j')&& isfield(Param.IndexRange,'MinIndex_j') - MaxIndex_j=Param.IndexRange.MaxIndex_j; - MinIndex_j=Param.IndexRange.MinIndex_j; - end - if isfield(Param,'InputTable') - [tild,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param); - iview_A=0;% series index (iview) for the first image series - iview_B=0;% series index (iview) for the second image series (only non zero for option 'shift' comparing two image series ) - if Param.ActionInput.CheckCiv1 - iview_A=1;% usual PIV, the image series is on the first line of the table - elseif Param.ActionInput.CheckCiv2 % civ2 is performed without Civ1, a netcdf file series is needed in the first table line - iview_A=2;% the second line is used for the input images of Civ2 - end -% if strcmp(Param.ActionInput.ListCompareMode,'shift') -% iview_B=iview_A+1; % the second image series is on the next line of the input table -% end - if iview_A~=0 - RootPath_A=Param.InputTable{iview_A,1}; - RootFile_A=Param.InputTable{iview_A,3}; - SubDir_A=Param.InputTable{iview_A,2}; - NomType_A=Param.InputTable{iview_A,4}; - FileExt_A=Param.InputTable{iview_A,5}; - if iview_B==0 - iview_B=iview_A;% the second image series is the same as the first - end - RootPath_B=Param.InputTable{iview_B,1}; - RootFile_B=Param.InputTable{iview_B,3}; - SubDir_B=Param.InputTable{iview_B,2}; - NomType_B=Param.InputTable{iview_B,4}; - FileExt_B=Param.InputTable{iview_B,5}; - end - - PairCiv2=''; - switch Param.ActionInput.ListCompareMode - case 'PIV' - PairCiv1=Param.ActionInput.PairIndices.ListPairCiv1; - if isfield(Param.ActionInput.PairIndices,'ListPairCiv2') - PairCiv2=Param.ActionInput.PairIndices.ListPairCiv2;%string which determines the civ2 pair - end - if iview_A==1% if Civ1 is performed - [i1_series_Civ1,i2_series_Civ1,j1_series_Civ1,j2_series_Civ1,check_bounds,NomTypeNc]=... + +if isfield(Param,'OutputSubDir')&& isfield(Param,'OutputDirExt') + OutputDir=[Param.OutputSubDir Param.OutputDirExt]; + OutputPath=fullfile(Param.OutputPath,num2str(Param.Experiment),num2str(Param.Device)); +else + disp_uvmat('ERROR','no output folder defined',checkrun) + return +end + +%% input files and indexing +MaxIndex_i=Param.IndexRange.MaxIndex_i; +MinIndex_i=Param.IndexRange.MinIndex_i; +MaxIndex_j=ones(size(MaxIndex_i));MinIndex_j=ones(size(MinIndex_i)); +if isfield(Param.IndexRange,'MaxIndex_j')&& isfield(Param.IndexRange,'MinIndex_j') + MaxIndex_j=Param.IndexRange.MaxIndex_j; + MinIndex_j=Param.IndexRange.MinIndex_j; +end + +[tild,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param); +iview_A=0;% series index (iview) for the first image series +iview_B=0;% series index (iview) for the second image series (only non zero for option 'shift' comparing two image series ) +if Param.ActionInput.CheckCiv1 + iview_A=1;% usual PIV, the image series is on the first line of the table +elseif Param.ActionInput.CheckCiv2 % civ2 is performed without Civ1, a netcdf file series is needed in the first table line + iview_A=2;% the second line is used for the input images of Civ2 +end +if iview_A~=0 + RootPath_A=Param.InputTable{iview_A,1}; + RootFile_A=Param.InputTable{iview_A,3}; + SubDir_A=Param.InputTable{iview_A,2}; + NomType_A=Param.InputTable{iview_A,4}; + FileExt_A=Param.InputTable{iview_A,5}; + if iview_B==0 + iview_B=iview_A;% the second image series is the same as the first + end + RootPath_B=Param.InputTable{iview_B,1}; + RootFile_B=Param.InputTable{iview_B,3}; + SubDir_B=Param.InputTable{iview_B,2}; + NomType_B=Param.InputTable{iview_B,4}; + FileExt_B=Param.InputTable{iview_B,5}; +end + +PairCiv1=Param.ActionInput.PairIndices.ListPairCiv1; + +[i1_series_Civ1,i2_series_Civ1,j1_series_Civ1,j2_series_Civ1,check_bounds,NomTypeNc]=... find_pair_indices(PairCiv1,i1_series{1},j1_series{1},MinIndex_i,MaxIndex_i,MinIndex_j,MaxIndex_j); - if ~isempty(PairCiv2) - [i1_series_Civ2,i2_series_Civ2,j1_series_Civ2,j2_series_Civ2,check_bounds_Civ2]=... - find_pair_indices(PairCiv2,i1_series{1},j1_series{1},MinIndex_i(1),MaxIndex_i(1),MinIndex_j(1),MaxIndex_j(1)); - check_bounds=check_bounds | check_bounds_Civ2; - end - else% we start from an existing Civ1 file - i1_series_Civ1=i1_series{1}; - i2_series_Civ1=i2_series{1}; - j1_series_Civ1=j1_series{1}; - j2_series_Civ1=j2_series{1}; - NomTypeNc=Param.InputTable{1,4}; - if ~isempty(PairCiv2) - [i1_series_Civ2,i2_series_Civ2,j1_series_Civ2,j2_series_Civ2,check_bounds,NomTypeNc]=... - find_pair_indices(PairCiv2,i1_series{2},j1_series{2},MinIndex_i(2),MaxIndex_i(2),MinIndex_j(2),MaxIndex_j(2)); - end - end - case 'displacement' - if isfield(Param.ActionInput,'OriginIndex') - i1_series_Civ1=Param.ActionInput.OriginIndex*ones(size(i1_series{1})); - else - i1_series_Civ1=ones(size(i1_series{1})); - end - i1_series_Civ2=i1_series_Civ1; - i2_series_Civ1=i1_series{1}; - i2_series_Civ2=i1_series{1}; - j1_series_Civ1=[];% no j index variation for the ref image - j1_series_Civ2=[]; - if isempty(j1_series{1}) - j2_series_Civ1=ones(size(i1_series_Civ1)); - else - j2_series_Civ1=j1_series{1};% if j index exist - end - j2_series_Civ2=j2_series_Civ1; - NomTypeNc='_1'; - case 'PIV volume' - % TODO, TODO - end - %determine frame indices for input with movie or other multiframe input file - if isempty(j1_series_Civ1)% simple movie with index i - FrameIndex_A_Civ1=i1_series_Civ1; - FrameIndex_B_Civ1=i2_series_Civ1; - j1_series_Civ1=ones(size(i1_series_Civ1)); - if strcmp(Param.ActionInput.ListCompareMode,'PIV') - j2_series_Civ1=ones(size(i1_series_Civ1)); - end - else % movie for each burst or volume (index j) - FrameIndex_A_Civ1=j1_series_Civ1; - FrameIndex_B_Civ1=j2_series_Civ1; - end - if isempty(PairCiv2) - FrameIndex_A_Civ2=FrameIndex_A_Civ1; - FrameIndex_B_Civ2=FrameIndex_B_Civ1; - else - if isempty(j1_series_Civ2) - FrameIndex_A_Civ2=i1_series_Civ2; - FrameIndex_B_Civ2=i2_series_Civ2; - j1_series_Civ2=ones(size(i1_series_Civ2)); - if strcmp(Param.ActionInput.ListCompareMode,'PIV') - j2_series_Civ2=ones(size(i1_series_Civ2)); - end - else - FrameIndex_A_Civ2=j1_series_Civ2; - FrameIndex_B_Civ2=j2_series_Civ2; - end - end - if isempty(i1_series_Civ1)||(~isempty(PairCiv2) && isempty(i1_series_Civ2)) - disp_uvmat('ERROR','no image pair for civ in the input file index range',checkrun) - return - end - end - - %% check the first image pair - if Param.ActionInput.CheckCiv1% Civ1 is performed - NbField=numel(i1_series_Civ1); - elseif Param.ActionInput.CheckCiv2 % Civ2 is performed without Civ1 - NbField=numel(i1_series_Civ2); - else - NbField=numel(i1_series_Civ1);% no image used (only fix or patch) TO CHECK - end - - %% Output directory - OutputDir=''; - if CheckOutputFile - OutputDir=[Param.OutputSubDir Param.OutputDirExt]; - end -end + +if isempty(i1_series_Civ1) + disp_uvmat('ERROR','no image pair for civ in the input file index range',checkrun) + return +end +NbField_i=size(i1_series_Civ1,2); +NbSlice=size(i1_series_Civ1,1); %% prepare output Data @@ -289,5 +171,19 @@ Data.Program='civ_series'; Data.CivStage=0;%default -check_civx=0;%default +list_param=(fieldnames(Param.ActionInput.Civ1))'; + list_param(strcmp('TestCiv1',list_param))=[];% remove the parameter TestCiv1 from the list + Civ1_param=regexprep(list_param,'^.+','Civ1_$0');% insert 'Civ1_' before each string in list_param + Civ1_param=[{'Civ1_ImageA','Civ1_ImageB','Civ1_Time','Civ1_Dt'} Civ1_param]; %insert the names of the two input images +Data.ListGlobalAttribute=[ListGlobalAttribute Civ1_param]; + % set the list of variables + 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 + Data.VarDimName={'nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1'}; + Data.VarAttribute{1}.Role='coord_x'; + Data.VarAttribute{2}.Role='coord_y'; + Data.VarAttribute{3}.Role='vector_x'; + Data.VarAttribute{4}.Role='vector_y'; + Data.VarAttribute{5}.Role='vector_z'; + Data.VarAttribute{6}.Role='ancillary'; + Data.VarAttribute{7}.Role='errorflag'; %% get timing from the ImaDoc file or input video @@ -314,5 +210,4 @@ end -%%%%% MAIN LOOP %%%%%% maskoldname='';% initiate the mask name FileType_A=''; @@ -322,8 +217,26 @@ CheckOverwrite=Param.CheckOverwrite; end -for ifield=1:NbField - tic - if ~isempty(RUNHandle)% update the waitbar in interactive mode with GUI series (checkrun=1) - update_waitbar(WaitbarHandle,ifield/NbField) + Data.Civ1_ImageA=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1),[],j1_series_Civ1(1,1)); + Data.Civ1_ImageB=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i1_series_Civ1(1),[],j2_series_Civ1(1,1)); + FileInfo=get_file_info(Data.Civ1_ImageA); + par_civ1=Param.ActionInput.Civ1;% parameters for civ1 + par_civ1.ImageHeight=FileInfo.Height;npy=FileInfo.Height; + par_civ1.ImageWidth=FileInfo.Width;npx=FileInfo.Width; +SearchRange_z=floor(Param.ActionInput.Civ1.SearchBoxSize(3)/2); + par_civ1.Dz=Param.ActionInput.Civ1.Dz; + par_civ1.ImageA=zeros(2*SearchRange_z+1,npy,npx); +par_civ1.ImageB=zeros(2*SearchRange_z+1,npy,npx); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%% MAIN LOOP %%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +for ifield=1:NbField_i + tstart=tic; + time_civ1=0; + time_patch1=0; + time_civ2=0; + time_patch2=0; + if checkrun% update the waitbar in interactive mode with GUI series (checkrun=1) + update_waitbar(WaitbarHandle,ifield/NbField_i) if checkrun && ~strcmp(get(RUNHandle,'BusyAction'),'queue') disp('program stopped by user') @@ -331,155 +244,72 @@ end end - if CheckInputFile - if iview_A==0 % no nc file has been entered - ncfile=fullfile_uvmat(Param.InputTable{1,1},Param.InputTable{1,2},Param.InputTable{1,3},Param.InputTable{1,5},... - NomTypeNc,i1_series_Civ1(ifield),i2_series_Civ1(ifield),j1_series_Civ1(ifield),j2_series_Civ1(ifield)); - else% an existing nc file has been entered - if iview_A==1% if Civ1 is performed - Civ1Dir=OutputDir; - else - Civ1Dir=Param.InputTable{1,2}; - end - if strcmp(Param.ActionInput.ListCompareMode,'PIV') - ncfile=fullfile_uvmat(RootPath_A,Civ1Dir,RootFile_A,'.nc',NomTypeNc,i1_series_Civ1(ifield),i2_series_Civ1(ifield),... - j1_series_Civ1(ifield),j2_series_Civ1(ifield)); - else - ncfile=fullfile_uvmat(RootPath_A,Civ1Dir,RootFile_A,'.nc',NomTypeNc,i2_series_Civ1(ifield),[],... - j1_series_Civ1(ifield),j2_series_Civ1(ifield)); - end - end - ncfile_out=ncfile;% by default - if isfield (Param.ActionInput,'Civ2') - i1_civ2=i1_series_Civ2(ifield); - i2_civ2=i1_civ2; - if ~isempty(i2_series_Civ2) - i2_civ2=i2_series_Civ2(ifield); - end - j1_civ2=1; - if ~isempty(j1_series_Civ2) - j1_civ2=j1_series_Civ2(ifield); - end - j2_civ2=i1_civ2; - if ~isempty(j2_series_Civ2) - j2_civ2=j2_series_Civ2(ifield); - end - if strcmp(Param.ActionInput.ListCompareMode,'PIV') - ncfile_out=fullfile_uvmat(RootPath_A,OutputDir,RootFile_A,'.nc',NomTypeNc,i1_civ2,i2_civ2,j1_civ2,j2_civ2); - else % displacement - ncfile_out=fullfile_uvmat(RootPath_A,OutputDir,RootFile_A,'.nc',NomTypeNc,i2_civ2,[],j2_civ2); - end - end - if ~CheckOverwrite && exist(ncfile_out,'file') - disp(['existing output file ' ncfile_out ' already exists, skip to next field']) - continue% skip iteration if the mode overwrite is desactivated and the result file already exists - end - end + + %indicate the values of all the global attributes in the output data + i1=i1_series_Civ1(ifield); + i2=i1; + if ~isempty(i2_series_Civ1) + i2=i2_series_Civ1(ifield); + end + j1=1; + if ~isempty(j1_series_Civ1) + j1=j1_series_Civ1(ifield); + end + j2=j1; + if ~isempty(j2_series_Civ1) + j2=j2_series_Civ1(ifield); + end + + Data.Civ1_Time=(Time(i2+1,j2+1)+Time(i1+1,j1+1))/2;% the Time is the Time at the middle of the image pair + Data.Civ1_Dt=Time(i2+1,j2+1)-Time(i1+1,j1+1); + + for ilist=1:length(list_param) + Data.(Civ1_param{4+ilist})=Param.ActionInput.Civ1.(list_param{ilist}); + end + + Data.CivStage=1; + % output nc file + ncfile_out=fullfile_uvmat(OutputPath,Param.InputTable{1,2},Param.InputTable{1,3},Param.InputTable{1,5},... + NomTypeNc,i1_series_Civ1(ifield),i2_series_Civ1(ifield),j1_series_Civ1(ifield),j2_series_Civ1(ifield)); + + if ~CheckOverwrite && exist(ncfile_out,'file') + disp(['existing output file ' ncfile_out ' already exists, skip to next field']) + continue% skip iteration if the mode overwrite is desactivated and the result file already exists + end + + %% Civ1 % if Civ1 computation is requested if isfield (Param.ActionInput,'Civ1') - if CheckInputFile - disp('civ1 started') - end - par_civ1=Param.ActionInput.Civ1;% parameters for civ1 - if CheckInputFile % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1.ImageNameA and B) - try - if strcmp(Param.ActionInput.ListCompareMode,'displacement') - ImageName_A=Param.ActionInput.RefFile; - else - ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(ifield),[],j1_series_Civ1(ifield)); - end - if strcmp(FileExt_A,'.nc')% case of input images in format netcdf - FieldName_A=Param.InputFields.FieldName; - [DataIn,tild,tild,errormsg]=nc2struct(ImageName_A,{FieldName_A}); - par_civ1.ImageA=DataIn.(FieldName_A); - else % usual image formats for image A - if isempty(FileType_A)% open the image object if not already done in case of movie input - [FileInfo_A,VideoObject_A]=get_file_info(ImageName_A); - FileType_A=FileInfo_A.FileType; - if isempty(Time) && ~isempty(find(strcmp(FileType_A,{'mmreader','video','cine_phantom'})))% case of video input - Time=zeros(FileInfo_A.NumberOfFrames+1,2); - Time(:,2)=(0:1/FileInfo_A.FrameRate:(FileInfo_A.NumberOfFrames)/FileInfo_A.FrameRate)'; - end - if ~isempty(FileType_A) && isempty(Time)% Time = index i +0.001 index j by default - MaxIndex_i=max(i2_series_Civ1); - MaxIndex_j=max(j2_series_Civ1); - Time=(1:MaxIndex_i)'*ones(1,MaxIndex_j); - Time=Time+0.001*ones(MaxIndex_i,1)*(1:MaxIndex_j); - Time=[zeros(1,MaxIndex_j);Time];% insert a first line of zeros - Time=[zeros(MaxIndex_i+1,1) Time];% insert a first column of zeros - end - end - if ~exist(ImageName_A,'file') - disp([ImageName_A ' missing']) - continue - end - [par_civ1.ImageA,VideoObject_A] = read_image(ImageName_A,FileType_A,VideoObject_A,FrameIndex_A_Civ1(ifield)); - end - ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(ifield),[],j2_series_Civ1(ifield)); - if strcmp(FileExt_B,'.nc') % case of input images in format netcdf - FieldName_B=Param.InputFields.FieldName; - [DataIn,tild,tild,errormsg]=nc2struct(ImageName_B,{FieldName_B}); - par_civ1.ImageB=DataIn.(FieldName_B); - else % usual image formats for image B - if isempty(FileType_B) - [FileInfo_B,VideoObject_B]=get_file_info(ImageName_B); - FileType_B=FileInfo_B.FileType; - end - if ~exist(ImageName_B,'file') - disp([ImageName_B ' missing']) - continue - end - [par_civ1.ImageB,VideoObject_B] = read_image(ImageName_B,FileType_B,VideoObject_B,FrameIndex_B_Civ1(ifield)); - end - catch ME % display errors in reading input images - if ~isempty(ME.message) - disp_uvmat('ERROR', ['error reading input image: ' ME.message],checkrun) - continue - end - end - par_civ1.ImageWidth=size(par_civ1.ImageA,2); - par_civ1.ImageHeight=size(par_civ1.ImageA,1); - list_param=(fieldnames(Param.ActionInput.Civ1))'; - list_param(strcmp('TestCiv1',list_param))=[];% remove the parameter TestCiv1 from the list - Civ1_param=regexprep(list_param,'^.+','Civ1_$0');% insert 'Civ1_' before each string in list_param - Civ1_param=[{'Civ1_ImageA','Civ1_ImageB','Civ1_Time','Civ1_Dt'} Civ1_param]; %insert the names of the two input images - %indicate the values of all the global attributes in the output data - Data.Civ1_ImageA=ImageName_A; - Data.Civ1_ImageB=ImageName_B; - i1=i1_series_Civ1(ifield); - i2=i1; - if ~isempty(i2_series_Civ1) - i2=i2_series_Civ1(ifield); - end - j1=1; - if ~isempty(j1_series_Civ1) - j1=j1_series_Civ1(ifield); - end - j2=j1; - if ~isempty(j2_series_Civ1) - j2=j2_series_Civ1(ifield); - end - if strcmp(Param.ActionInput.ListCompareMode,'displacement') - Data.Civ1_Time=Time(i2+1,j2+1);% the Time is the Time of the secodn image - Data.Civ1_Dt=1;% Time interval is 1, to yield displacement instead of velocity=displacement/Dt at reading - else - Data.Civ1_Time=(Time(i2+1,j2+1)+Time(i1+1,j1+1))/2;% the Time is the Time at the middle of the image pair - Data.Civ1_Dt=Time(i2+1,j2+1)-Time(i1+1,j1+1); - end - for ilist=1:length(list_param) - Data.(Civ1_param{4+ilist})=Param.ActionInput.Civ1.(list_param{ilist}); - end - Data.ListGlobalAttribute=[ListGlobalAttribute Civ1_param]; - Data.CivStage=1; - end - % set the list of variables - Data.ListVarName={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_F','Civ1_C'};% cell array containing the names of the fields to record - Data.VarDimName={'nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1'}; - Data.VarAttribute{1}.Role='coord_x'; - Data.VarAttribute{2}.Role='coord_y'; - Data.VarAttribute{3}.Role='vector_x'; - Data.VarAttribute{4}.Role='vector_y'; - Data.VarAttribute{5}.Role='warnflag'; - % case of mask + + disp('civ1 started') + + + % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1.ImageNameA and B) + + par_civ1.ImageA=zeros(2*SearchRange_z+1,npy,npx); + par_civ1.ImageB=zeros(2*SearchRange_z+1,npy,npx); + + for islice=ceil(par_civ1.Dz/2):par_civ1.Dz:NbSlice + if par_civ1.Dz<2*SearchRange_z+1 + par_civ1.ImageA=circshift(par_civ1.ImageA,-par_civ1.Dz); + par_civ1.ImageB=circshift(par_civ1.ImageA,-par_civ1.Dz); + end + for iz=1:par_civ1.Dz + ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(ifield),[],j1_series_Civ1(ifield,islice)); + A= read_image(ImageName_A,FileType_A); + ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(ifield),[],j2_series_Civ1(ifield,islice)); + B= read_image(ImageName_B,FileType_B); + par_civ1.ImageA(2*SearchRange_z+2-iz,:,:) = A; + par_civ1.ImageB(2*SearchRange_z+2-iz,:,:) = B; + end + % caluclate velocity data (y and v in indices, reverse to y component) + [Data.Civ1_X(islice,:,:),Data.Civ1_Y(islice,:,:), utable, vtable,wtable, ctable, FF, result_conv, errormsg] = civ3D (par_civ1); + if ~isempty(errormsg) + disp_uvmat('ERROR',errormsg,checkrun) + return + end + + end + % case of mask TO ADAPT if par_civ1.CheckMask&&~isempty(par_civ1.Mask) if isfield(par_civ1,'NbSlice') @@ -493,5 +323,5 @@ par_civ1.Mask=mask; %use mask already opened else - if exist(maskname,'file') + if ~isempty(regexp(maskname,'(^http://)|(^https://)'))|| exist(maskname,'file') try par_civ1.Mask=imread(maskname);%update the mask, an store it for future use @@ -510,61 +340,22 @@ end end - if strcmp(Param.ActionInput.ListCompareMode, 'PIV volume') - Data.ListVarName=[Data.ListVarName 'Civ1_Z']; - Data.Civ1_X=[];Data.Civ1_Y=[];Data.Civ1_Z=[]; - Data.Civ1_U=[];Data.Civ1_V=[];Data.Civ1_C=[];Data.Civ1_F=[]; - for ivol=1:NbSlice - % caluclate velocity data (y and v in indices, reverse to y component) - [xtable, ytable, utable, vtable, ctable, F, result_conv, errormsg] = civ (par_civ1); - if ~isempty(errormsg) - disp_uvmat('ERROR',errormsg,checkrun) - return - end - Data.Civ1_X=[Data.Civ1_X reshape(xtable,[],1)]; - Data.Civ1_Y=[Data.Civ1_Y reshape(Param.Civ1.ImageHeight-ytable+1,[],1)]; - Data.Civ1_Z=[Data.Civ1_Z ivol*ones(numel(xtable),1)];% z=image index in image coordinates - Data.Civ1_U=[Data.Civ1_U reshape(utable,[],1)]; - Data.Civ1_V=[Data.Civ1_V reshape(-vtable,[],1)]; - Data.Civ1_C=[Data.Civ1_C reshape(ctable,[],1)]; - Data.Civ1_F=[Data.Civ1_C reshape(F,[],1)]; - end - else %usual PIV - % caluclate velocity data (y and v in indices, reverse to y component) - [xtable, ytable, utable, vtable, ctable, F, result_conv, errormsg] = civ (par_civ1); - if ~isempty(errormsg) - disp_uvmat('ERROR',errormsg,checkrun) - return - end - Data.Civ1_X=reshape(xtable,[],1); - Data.Civ1_Y=reshape(par_civ1.ImageHeight-ytable+1,[],1); - Data.Civ1_U=reshape(utable,[],1); - Data.Civ1_V=reshape(-vtable,[],1); - Data.Civ1_C=reshape(ctable,[],1); - Data.Civ1_F=reshape(F,[],1); - end - else% we use existing Civ1 data - if exist('ncfile','var') - CivFile=ncfile; - [Data,tild,tild,errormsg]=nc2struct(CivFile,'ListGlobalAttribute','absolut_time_T0'); %look for the constant 'absolut_time_T0' to detect old civx data format - if ~isempty(errormsg) - disp_uvmat('ERROR',errormsg,checkrun) - return - end - [Data,tild,tild,errormsg]=nc2struct(CivFile);%read civ1 and fix1 data in the existing netcdf file - elseif isfield(Param,'Civ1_X') - Data.ListGlobalAttribute={}; - Data.ListVarName={}; - Data.VarDimName={}; - Data.Civ1_X=Param.Civ1_X; - Data.Civ1_Y=Param.Civ1_Y; - Data.Civ1_U=Param.Civ1_U; - Data.Civ1_V=Param.Civ1_V; - Data.Civ1_FF=Param.Civ1_FF; - end - end - + % Data.ListVarName=[Data.ListVarName 'Civ1_Z']; + % Data.Civ1_X=[];Data.Civ1_Y=[];Data.Civ1_Z=[]; + % Data.Civ1_U=[];Data.Civ1_V=[];Data.Civ1_C=[]; + % + % + % Data.Civ1_X=[Data.Civ1_X reshape(xtable,[],1)]; + % Data.Civ1_Y=[Data.Civ1_Y reshape(Param.Civ1.ImageHeight-ytable+1,[],1)]; + % Data.Civ1_Z=[Data.Civ1_Z ivol*ones(numel(xtable),1)];% z=image index in image coordinates + % Data.Civ1_U=[Data.Civ1_U reshape(utable,[],1)]; + % Data.Civ1_V=[Data.Civ1_V reshape(-vtable,[],1)]; + % Data.Civ1_C=[Data.Civ1_C reshape(ctable,[],1)]; + % Data.Civ1_FF=[Data.Civ1_FF reshape(F,[],1)]; + + end + %% Fix1 if isfield (Param.ActionInput,'Fix1') - disp('fix1 started') + disp('detect_false1 started') if ~isfield (Param.ActionInput,'Civ1')% if we use existing Civ1, remove previous data beyond Civ1 Fix1_attr=find(strcmp('Fix1',Data.ListGlobalAttribute)); @@ -581,9 +372,5 @@ end Data.ListGlobalAttribute=[Data.ListGlobalAttribute Fix1_param]; - Data.ListVarName=[Data.ListVarName {'Civ1_FF'}]; - Data.VarDimName=[Data.VarDimName {'nb_vec_1'}]; - nbvar=length(Data.ListVarName); - Data.VarAttribute{nbvar}.Role='errorflag'; - Data.Civ1_FF=int8(fix(Param.ActionInput.Fix1,Data.Civ1_F,Data.Civ1_C,Data.Civ1_U,Data.Civ1_V)); + Data.Civ1_FF=uint8(detect_false(Param.ActionInput.Fix1,Data.Civ1_C,Data.Civ1_U,Data.Civ1_V,Data.Civ1_FF)); Data.CivStage=2; end @@ -591,10 +378,6 @@ if isfield (Param.ActionInput,'Patch1') disp('patch1 started') - if check_civx - errormsg='Civ Matlab input needed for patch'; - disp_uvmat('ERROR',errormsg,checkrun) - return - end - + tstart_patch1=tic; + % record the processing parameters of Patch1 as global attributes in the result nc file list_param=fieldnames(Param.ActionInput.Patch1)'; @@ -606,5 +389,5 @@ Data.CivStage=3;% record the new state of processing Data.ListGlobalAttribute=[Data.ListGlobalAttribute Patch1_param]; - + % list the variables to record nbvar=length(Data.ListVarName); @@ -624,92 +407,100 @@ ind_good=1:numel(Data.Civ1_X); end - + if isempty(ind_good) + disp_uvmat('ERROR','all vectors of civ1 are bad, check input parameters' ,checkrun) + return + end + % perform Patch calculation using the UVMAT fct 'filter_tps' + [Data.Civ1_SubRange,Data.Civ1_NbCentres,Data.Civ1_Coord_tps,Data.Civ1_U_tps,Data.Civ1_V_tps,tild,Ures, Vres,tild,FFres]=... filter_tps([Data.Civ1_X(ind_good) Data.Civ1_Y(ind_good)],Data.Civ1_U(ind_good),Data.Civ1_V(ind_good),[],Data.Patch1_SubDomainSize,Data.Patch1_FieldSmooth,Data.Patch1_MaxDiff); Data.Civ1_U_smooth(ind_good)=Ures;% take the interpolated (smoothed) velocity values for good vectors, keep civ1 data for the other Data.Civ1_V_smooth(ind_good)=Vres; - Data.Civ1_FF(ind_good)=int8(FFres); + Data.Civ1_FF(ind_good)=uint8(FFres); + time_patch1=toc(tstart_patch1); disp('patch1 performed') end - + %% Civ2 if isfield (Param.ActionInput,'Civ2') disp('civ2 started') + tstart_civ2=tic; par_civ2=Param.ActionInput.Civ2; - if CheckInputFile % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1.ImageNameA and B) - par_civ2.ImageA=[]; - par_civ2.ImageB=[]; - if strcmp(Param.ActionInput.ListCompareMode,'displacement') - ImageName_A_Civ2=Param.ActionInput.RefFile; - else + % read input images + par_civ2.ImageA=[]; + par_civ2.ImageB=[]; + if strcmp(Param.ActionInput.ListCompareMode,'displacement') + ImageName_A_Civ2=Param.ActionInput.RefFile; + else ImageName_A_Civ2=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_civ2,[],j1_civ2); - end - if strcmp(ImageName_A_Civ2,ImageName_A) && isequal(FrameIndex_A_Civ1(ifield),FrameIndex_A_Civ2(ifield)) - par_civ2.ImageA=par_civ1.ImageA; - else - [par_civ2.ImageA,VideoObject_A] = read_image(ImageName_A_Civ2,FileType_A,VideoObject_A,FrameIndex_A_Civ2(ifield)); - end - ImageName_B_Civ2=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_civ2,[],j2_civ2); - if strcmp(ImageName_B_Civ2,ImageName_B) && isequal(FrameIndex_B_Civ1(ifield),FrameIndex_B_Civ2) - par_civ2.ImageB=par_civ1.ImageB; - else - [par_civ2.ImageB,VideoObject_B] = read_image(ImageName_B_Civ2,FileType_B,VideoObject_B,FrameIndex_B_Civ2(ifield)); - end - par_civ2.ImageWidth=FileInfo_A.Width; - par_civ2.ImageHeight=FileInfo_A.Height; - if isfield(par_civ2,'Grid')% grid points set as input file - if ischar(par_civ2.Grid)%read the grid file if the input is a file name - par_civ2.Grid=dlmread(par_civ2.Grid); - par_civ2.Grid(1,:)=[];%the first line must be removed (heading in the grid file) - end - else% automatic grid - minix=floor(par_civ2.Dx/2)-0.5; - maxix=minix+par_civ2.Dx*floor((par_civ2.ImageWidth-1)/par_civ2.Dx); - miniy=floor(par_civ2.Dy/2)-0.5; - maxiy=minix+par_civ2.Dy*floor((par_civ2.ImageHeight-1)/par_civ2.Dy); - [GridX,GridY]=meshgrid(minix:par_civ2.Dx:maxix,miniy:par_civ2.Dy:maxiy); - par_civ2.Grid(:,1)=reshape(GridX,[],1); - par_civ2.Grid(:,2)=reshape(GridY,[],1); - end - end - + end + if strcmp(ImageName_A_Civ2,ImageName_A) && isequal(FrameIndex_A_Civ1(ifield),FrameIndex_A_Civ2(ifield)) + par_civ2.ImageA=par_civ1.ImageA; + else + [par_civ2.ImageA,VideoObject_A] = read_image(ImageName_A_Civ2,FileType_A,VideoObject_A,FrameIndex_A_Civ2(ifield)); + end + ImageName_B_Civ2=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_civ2,[],j2_civ2); + if strcmp(ImageName_B_Civ2,ImageName_B) && isequal(FrameIndex_B_Civ1(ifield),FrameIndex_B_Civ2) + par_civ2.ImageB=par_civ1.ImageB; + else + [par_civ2.ImageB,VideoObject_B] = read_image(ImageName_B_Civ2,FileType_B,VideoObject_B,FrameIndex_B_Civ2(ifield)); + end + [FileInfo_A,VideoObject_A]=get_file_info(ImageName_A_Civ2); + par_civ2.ImageWidth=FileInfo_A.Width; + par_civ2.ImageHeight=FileInfo_A.Height; + if isfield(par_civ2,'Grid')% grid points set as input file + if ischar(par_civ2.Grid)%read the grid file if the input is a file name + par_civ2.Grid=dlmread(par_civ2.Grid); + par_civ2.Grid(1,:)=[];%the first line must be removed (heading in the grid file) + end + else% automatic grid + minix=floor(par_civ2.Dx/2)-0.5; + maxix=minix+par_civ2.Dx*floor((par_civ2.ImageWidth-1)/par_civ2.Dx); + miniy=floor(par_civ2.Dy/2)-0.5; + maxiy=minix+par_civ2.Dy*floor((par_civ2.ImageHeight-1)/par_civ2.Dy); + [GridX,GridY]=meshgrid(minix:par_civ2.Dx:maxix,miniy:par_civ2.Dy:maxiy); + par_civ2.Grid(:,1)=reshape(GridX,[],1); + par_civ2.Grid(:,2)=reshape(GridY,[],1); + end + % end + % get the guess from patch1 or patch2 (case 'CheckCiv3') - if CheckInputFile % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1.ImageNameA and B) - if isfield (par_civ2,'CheckCiv3') && par_civ2.CheckCiv3 %get the guess from patch2 - SubRange= Data.Civ2_SubRange; - NbCentres=Data.Civ2_NbCentres; - Coord_tps=Data.Civ2_Coord_tps; - U_tps=Data.Civ2_U_tps; - V_tps=Data.Civ2_V_tps; - CivStage=Data.CivStage;%store the current CivStage - Civ1_Dt=Data.Civ2_Dt; - Data=[];%reinitialise the result structure Data - Data.ListGlobalAttribute={'Conventions','Program','CivStage'}; - Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes - Data.Program='civ_series'; - Data.CivStage=CivStage+1;%update the current civStage after reinitialisation of Data - Data.ListVarName={}; - Data.VarDimName={}; - else % get the guess from patch1 - SubRange= Data.Civ1_SubRange; - NbCentres=Data.Civ1_NbCentres; - Coord_tps=Data.Civ1_Coord_tps; - U_tps=Data.Civ1_U_tps; - V_tps=Data.Civ1_V_tps; - Civ1_Dt=Data.Civ1_Dt; - Data.CivStage=4; - end - else - SubRange= par_civ2.Civ1_SubRange; - NbCentres=par_civ2.Civ1_NbCentres; - Coord_tps=par_civ2.Civ1_Coord_tps; - U_tps=par_civ2.Civ1_U_tps; - V_tps=par_civ2.Civ1_V_tps; - Civ1_Dt=par_civ2.Civ1_Dt; - Civ2_Dt=par_civ2.Civ1_Dt; + + if isfield (par_civ2,'CheckCiv3') && par_civ2.CheckCiv3 %get the guess from patch2 + SubRange= Data.Civ2_SubRange; + NbCentres=Data.Civ2_NbCentres; + Coord_tps=Data.Civ2_Coord_tps; + U_tps=Data.Civ2_U_tps; + V_tps=Data.Civ2_V_tps; + CivStage=Data.CivStage;%store the current CivStage + Civ1_Dt=Data.Civ2_Dt; + Data=[];%reinitialise the result structure Data + Data.ListGlobalAttribute={'Conventions','Program','CivStage'}; + Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes + Data.Program='civ_series'; + Data.CivStage=CivStage+1;%update the current civStage after reinitialisation of Data Data.ListVarName={}; Data.VarDimName={}; - end + else % get the guess from patch1 + SubRange= Data.Civ1_SubRange; + NbCentres=Data.Civ1_NbCentres; + Coord_tps=Data.Civ1_Coord_tps; + U_tps=Data.Civ1_U_tps; + V_tps=Data.Civ1_V_tps; + Civ1_Dt=Data.Civ1_Dt; + Data.CivStage=4; + end + % else + % SubRange= par_civ2.Civ1_SubRange; + % NbCentres=par_civ2.Civ1_NbCentres; + % Coord_tps=par_civ2.Civ1_Coord_tps; + % U_tps=par_civ2.Civ1_U_tps; + % V_tps=par_civ2.Civ1_V_tps; + % Civ1_Dt=par_civ2.Civ1_Dt; + % Civ2_Dt=par_civ2.Civ1_Dt; + % Data.ListVarName={}; + % Data.VarDimName={}; + % end Shiftx=zeros(size(par_civ2.Grid,1),1);% shift expected from civ1 data Shifty=zeros(size(par_civ2.Grid,1),1); @@ -742,5 +533,5 @@ end end - if par_civ2.CheckMask&&~isempty(par_civ2.Mask) + if par_civ2.CheckMask&&~isempty(par_civ2.Mask) if isfield(par_civ2,'NbSlice') [RootPath_mask,SubDir_mask,RootFile_mask,i1_mask,i2_mask,j1_mask,j2_mask,Ext_mask]=fileparts_uvmat(Param.ActionInput.Civ2.Mask); @@ -754,13 +545,13 @@ else if exist(maskname,'file') - try - par_civ2.Mask=imread(maskname);%update the mask, an store it for future use - catch ME - if ~isempty(ME.message) - errormsg=['error reading input image: ' ME.message]; - disp_uvmat('ERROR',errormsg,checkrun) - return + try + par_civ2.Mask=imread(maskname);%update the mask, an store it for future use + catch ME + if ~isempty(ME.message) + errormsg=['error reading input image: ' ME.message]; + disp_uvmat('ERROR',errormsg,checkrun) + return + end end - end else par_civ2.Mask=[]; @@ -771,12 +562,12 @@ end - if CheckInputFile % else Dt given by par_civ2 - if strcmp(Param.ActionInput.ListCompareMode,'displacement') - Civ1_Dt=1; - Civ2_Dt=1; - else - Civ2_Dt=Time(i2_civ2+1,j2_civ2+1)-Time(i1_civ2+1,j1_civ2+1); - end - end + + if strcmp(Param.ActionInput.ListCompareMode,'displacement') + Civ1_Dt=1; + Civ2_Dt=1; + else + Civ2_Dt=Time(i2_civ2+1,j2_civ2+1)-Time(i1_civ2+1,j1_civ2+1); + end + par_civ2.SearchBoxShift=(Civ2_Dt/Civ1_Dt)*[Shiftx(nbval>=1)./nbval(nbval>=1) Shifty(nbval>=1)./nbval(nbval>=1)]; % shift the grid points by half the expected shift to provide the correlation box position in image A @@ -788,8 +579,9 @@ par_civ2.DVDY=DVDY(nbval>=1)./nbval(nbval>=1); end - + % calculate velocity data (y and v in image indices, reverse to y component) + [xtable, ytable, utable, vtable, ctable, F,result_conv,errormsg] = civ (par_civ2); - + list_param=(fieldnames(Param.ActionInput.Civ2))'; list_param(strcmp('TestCiv2',list_param))=[];% remove the parameter TestCiv2 from the list @@ -812,9 +604,9 @@ end Data.ListGlobalAttribute=[Data.ListGlobalAttribute Civ2_param]; - + nbvar=numel(Data.ListVarName); % define the Civ2 variable (if Civ2 data are not replaced from previous calculation) if isempty(find(strcmp('Civ2_X',Data.ListVarName),1)) - 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 + Data.ListVarName=[Data.ListVarName {'Civ2_X','Civ2_Y','Civ2_U','Civ2_V','Civ2_C','Civ2_FF'}];% cell array containing the names of the fields to record Data.VarDimName=[Data.VarDimName {'nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2'}]; Data.VarAttribute{nbvar+1}.Role='coord_x'; @@ -822,5 +614,6 @@ Data.VarAttribute{nbvar+3}.Role='vector_x'; Data.VarAttribute{nbvar+4}.Role='vector_y'; - Data.VarAttribute{nbvar+5}.Role='warnflag'; + Data.VarAttribute{nbvar+5}.Role='ancillary'; + Data.VarAttribute{nbvar+6}.Role='errorflag'; end Data.Civ2_X=reshape(xtable,[],1); @@ -829,29 +622,30 @@ Data.Civ2_V=reshape(-vtable,[],1); Data.Civ2_C=reshape(ctable,[],1); - Data.Civ2_F=reshape(F,[],1); + Data.Civ2_FF=reshape(F,[],1); disp('civ2 performed') + time_civ2=toc(tstart_civ2); elseif ~isfield(Data,'ListVarName') % we start there, using existing Civ2 data if exist('ncfile','var') CivFile=ncfile; - [Data,tild,tild,errormsg]=nc2struct(CivFile);%read civ1 and fix1 data in the existing netcdf file + [Data,tild,tild,errormsg]=nc2struct(CivFile);%read civ1 and detect_false1 data in the existing netcdf file if ~isempty(errormsg) disp_uvmat('ERROR',errormsg,checkrun) return end - elseif isfield(Param,'Civ2_X')% use Civ2 data as input in Param (test mode) - Data.ListGlobalAttribute={}; - Data.ListVarName={}; - Data.VarDimName={}; - Data.Civ2_X=Param.Civ2_X; - Data.Civ2_Y=Param.Civ2_Y; - Data.Civ2_U=Param.Civ2_U; - Data.Civ2_V=Param.Civ2_V; - Data.Civ2_FF=Param.Civ2_FF; - end - end - + % elseif isfield(Param,'Civ2_X')% use Civ2 data as input in Param (test mode) + % Data.ListGlobalAttribute={}; + % Data.ListVarName={}; + % Data.VarDimName={}; + % Data.Civ2_X=Param.Civ2_X; + % Data.Civ2_Y=Param.Civ2_Y; + % Data.Civ2_U=Param.Civ2_U; + % Data.Civ2_V=Param.Civ2_V; + % Data.Civ2_FF=Param.Civ2_FF; + end + end + %% Fix2 if isfield (Param.ActionInput,'Fix2') - disp('fix2 started') + disp('detect_false2 started') list_param=fieldnames(Param.ActionInput.Fix2)'; Fix2_param=regexprep(list_param,'^.+','Fix2_$0');% insert 'Fix1_' before each string in ListFixParam @@ -861,33 +655,12 @@ end Data.ListGlobalAttribute=[Data.ListGlobalAttribute Fix2_param]; - % - % ListFixParam=fieldnames(Param.ActionInput.Fix2); - % for ilist=1:length(ListFixParam) - % ParamName=ListFixParam{ilist}; - % ListName=['Fix2_' ParamName]; - % eval(['Data.ListGlobalAttribute=[Data.ListGlobalAttribute ''' ParamName '''];']) - % eval(['Data.' ListName '=Param.ActionInput.Fix2.' ParamName ';']) - % end - if check_civx - if ~isfield(Data,'fix2') - Data.ListGlobalAttribute=[Data.ListGlobalAttribute 'fix2']; - Data.fix2=1; - Data.ListVarName=[Data.ListVarName {'vec2_FixFlag'}]; - Data.VarDimName=[Data.VarDimName {'nb_vectors2'}]; - end - Data.vec_FixFlag=fix(Param.Fix2,Data.vec2_F,Data.vec2_C,Data.vec2_U,Data.vec2_V,Data.vec2_X,Data.vec2_Y); - else - Data.ListVarName=[Data.ListVarName {'Civ2_FF'}]; - Data.VarDimName=[Data.VarDimName {'nb_vec_2'}]; - nbvar=length(Data.ListVarName); - Data.VarAttribute{nbvar}.Role='errorflag'; - Data.Civ2_FF=double(fix(Param.ActionInput.Fix2,Data.Civ2_F,Data.Civ2_C,Data.Civ2_U,Data.Civ2_V)); - Data.CivStage=Data.CivStage+1; - end - end - + Data.Civ2_FF=double(detect_false(Param.ActionInput.Fix2,Data.Civ2_C,Data.Civ2_U,Data.Civ2_V,Data.Civ2_FF)); + Data.CivStage=Data.CivStage+1; + end + %% Patch2 if isfield (Param.ActionInput,'Patch2') disp('patch2 started') + tstart_patch2=tic; list_param=fieldnames(Param.ActionInput.Patch2)'; list_param(strcmp('TestPatch2',list_param))=[];% remove the parameter TestCiv1 from the list @@ -898,10 +671,10 @@ end Data.ListGlobalAttribute=[Data.ListGlobalAttribute Patch2_param]; - + nbvar=length(Data.ListVarName); Data.ListVarName=[Data.ListVarName {'Civ2_U_smooth','Civ2_V_smooth','Civ2_SubRange','Civ2_NbCentres','Civ2_Coord_tps','Civ2_U_tps','Civ2_V_tps'}]; Data.VarDimName=[Data.VarDimName {'nb_vec_2','nb_vec_2',{'nb_coord','nb_bounds','nb_subdomain_2'},{'nb_subdomain_2'},... {'nb_tps_2','nb_coord','nb_subdomain_2'},{'nb_tps_2','nb_subdomain_2'},{'nb_tps_2','nb_subdomain_2'}}]; - + Data.VarAttribute{nbvar+1}.Role='vector_x'; Data.VarAttribute{nbvar+2}.Role='vector_y'; @@ -916,4 +689,9 @@ ind_good=1:numel(Data.Civ2_X); end + if isempty(ind_good) + disp_uvmat('ERROR','all vectors of civ2 are bad, check input parameters' ,checkrun) + return + end + [Data.Civ2_SubRange,Data.Civ2_NbCentres,Data.Civ2_Coord_tps,Data.Civ2_U_tps,Data.Civ2_V_tps,tild,Ures,Vres,tild,FFres]=... filter_tps([Data.Civ2_X(ind_good) Data.Civ2_Y(ind_good)],Data.Civ2_U(ind_good),Data.Civ2_V(ind_good),[],Data.Patch2_SubDomainSize,Data.Patch2_FieldSmooth,Data.Patch2_MaxDiff); @@ -922,18 +700,26 @@ Data.Civ2_FF(ind_good)=FFres; Data.CivStage=Data.CivStage+1; + time_patch2=toc(tstart_patch2); disp('patch2 performed') end - + %% write result in a netcdf file if requested - if CheckOutputFile - errormsg=struct2nc(ncfile_out,Data); - if isempty(errormsg) - disp([ncfile_out ' written']) - %[success,msg] = fileattrib(ncfile_out ,'+w','g');% done in struct2nc - else - disp(errormsg) - end - disp(['ellapsed time ' num2str(toc/60,2) ' minutes']) - end + % if CheckOutputFile + errormsg=struct2nc(ncfile_out,Data); + if isempty(errormsg) + disp([ncfile_out ' written']) + %[success,msg] = fileattrib(ncfile_out ,'+w','g');% done in struct2nc + else + disp(errormsg) + end + time_total=toc(tstart); + disp(['ellapsed time ' num2str(time_total/60,2) ' minutes']) + disp(['time image reading ' num2str(time_input,2) ' s']) + disp(['time civ1 ' num2str(time_civ1,2) ' s']) + disp(['time patch1 ' num2str(time_patch1,2) ' s']) + disp(['time civ2 ' num2str(time_civ2,2) ' s']) + disp(['time patch2 ' num2str(time_patch2,2) ' s']) + disp(['time other ' num2str((time_total-time_input-time_civ1-time_patch1-time_civ2-time_patch2),2) ' s']) + % end end @@ -971,23 +757,25 @@ % .DVDY -function [xtable,ytable,utable,vtable,ctable,F,result_conv,errormsg] = civ (par_civ) +function [xtable,ytable,utable,vtable,wtable,ctable,FF,result_conv,errormsg] = civ3D (par_civ) %% prepare measurement grid -if isfield(par_civ,'Grid')% grid points set as input, central positions of the sub-images in image A - if ischar(par_civ.Grid)%read the grid file if the input is a file name (grid in x, y image coordinates) - par_civ.Grid=dlmread(par_civ.Grid); - par_civ.Grid(1,:)=[];%the first line must be removed (heading in the grid file) - end - % else par_civ.Grid is already an array, no action here -else% automatic grid in x, y image coordinates - minix=floor(par_civ.Dx/2)-0.5; - maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx); - miniy=floor(par_civ.Dy/2)-0.5;% first automatic grid point at half the mesh Dy - maxiy=minix+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy); - [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy); - par_civ.Grid(:,1)=reshape(GridX,[],1); - par_civ.Grid(:,2)=reshape(GridY,[],1);% increases with array index -end -nbvec=size(par_civ.Grid,1); + +minix=floor(par_civ.Dx/2)-0.5; +maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx); +miniy=floor(par_civ.Dy/2)-0.5;% first automatic grid point at half the mesh Dy +maxiy=miniy+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy); +[GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy); +par_civ.Grid(:,:,1)=GridX; +par_civ.Grid(:,:,2)=GridY;% increases with array index, +[nbvec_y,nbvec_x,~]=size(par_civ.Grid); +% +% +% minix=floor(par_civ.Dx/2)-0.5; +% maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx); +% miniy=floor(par_civ.Dy/2)-0.5;% first automatic grid point at half the mesh Dy +% maxiy=miniy+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy); +% [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy); +% par_civ.Grid(:,1)=reshape(GridX,[],1); +% par_civ.Grid(:,2)=reshape(GridY,[],1);% increases with array index %% prepare correlation and search boxes @@ -996,36 +784,30 @@ isx2=floor(par_civ.SearchBoxSize(1)/2); isy2=floor(par_civ.SearchBoxSize(2)/2); +isz2=floor(par_civ.SearchBoxSize(3)/2); +kref=isz2+1;%middle index of the z slice shiftx=round(par_civ.SearchBoxShift(:,1));%use the input shift estimate, rounded to the next integer value shifty=-round(par_civ.SearchBoxShift(:,2));% sign minus because image j index increases when y decreases if numel(shiftx)==1% case of a unique shift for the whole field( civ1) - shiftx=shiftx*ones(nbvec,1); - shifty=shifty*ones(nbvec,1); + shiftx=shiftx*ones(nbvec_y,nbvec_x,1); + shifty=shifty*ones(nbvec_y,nbvec_x,1); end %% Array initialisation and default output if par_civ.CorrSmooth=0 (just the grid calculated, no civ computation) -xtable=round(par_civ.Grid(:,1)+0.5)-0.5; -ytable=round(par_civ.ImageHeight-par_civ.Grid(:,2)+0.5)-0.5;% y index corresponding to the position in image coordiantes +xtable=round(par_civ.Grid(:,:,1)+0.5)-0.5; +ytable=round(par_civ.ImageHeight-par_civ.Grid(:,:,2)+0.5)-0.5;% y index corresponding to the position in image coordiantes utable=shiftx;%zeros(nbvec,1); vtable=shifty;%zeros(nbvec,1); -ctable=zeros(nbvec,1); -F=zeros(nbvec,1); +wtable=zeros(size(utable)); +ctable=zeros(nbvec_y,nbvec_x,1); +FF=zeros(nbvec_y,nbvec_x,1); result_conv=[]; errormsg=''; %% prepare mask -if isfield(par_civ,'Mask') && ~isempty(par_civ.Mask) - if strcmp(par_civ.Mask,'all') - return % get the grid only, no civ calculation - elseif ischar(par_civ.Mask) - par_civ.Mask=imread(par_civ.Mask);% read the mask if not allready done - end -end check_MinIma=isfield(par_civ,'MinIma');% test for image luminosity threshold check_MaxIma=isfield(par_civ,'MaxIma') && ~isempty(par_civ.MaxIma); -par_civ.ImageA=sum(double(par_civ.ImageA),3);%sum over rgb component for color images -par_civ.ImageB=sum(double(par_civ.ImageB),3); -[npy_ima npx_ima]=size(par_civ.ImageA); -if ~isequal(size(par_civ.ImageB),[npy_ima npx_ima]) +[npz,npy_ima,npx_ima]=size(par_civ.ImageA); +if ~isequal(size(par_civ.ImageB),[npz npy_ima npx_ima]) errormsg='image pair with unequal size'; return @@ -1033,5 +815,5 @@ %% Apply mask -% Convention for mask IDEAS TO IMPLEMENT ? +% Convention for mask, IDEAS NOT IMPLEMENTED % mask >200 : velocity calculated % 200 >=mask>150;velocity not calculated, interpolation allowed (bad spots) @@ -1040,5 +822,5 @@ % 20>=mask: velocity=0 checkmask=0; -MinA=min(min(par_civ.ImageA)); +MinA=min(min(min(par_civ.ImageA))); %MinB=min(min(par_civ.ImageB)); %check_undefined=false(size(par_civ.ImageA)); @@ -1050,6 +832,4 @@ end check_undefined=(par_civ.Mask<200 & par_civ.Mask>=20 ); - % par_civ.ImageA(check_undefined)=0;% put image A to zero (i.e. the min image value) in the undefined area - % par_civ.ImageB(check_undefined)=0;% put image B to zero (i.e. the min image value) in the undefined area end @@ -1065,112 +845,110 @@ if par_civ.CorrSmooth~=0 % par_civ.CorrSmooth=0 implies no civ computation (just input image and grid points given) - for ivec=1:nbvec - iref=round(par_civ.Grid(ivec,1)+0.5);% xindex on the image A for the middle of the correlation box - jref=round(par_civ.ImageHeight-par_civ.Grid(ivec,2)+0.5);% j index for the middle of the correlation box in the image A - F(ivec)=0; - subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage - subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage - subrange2_x=iref+shiftx(ivec)-isx2:iref+shiftx(ivec)+isx2;%x indices defining the second subimage - subrange2_y=jref+shifty(ivec)-isy2:jref+shifty(ivec)+isy2;%y indices defining the second subimage - image1_crop=MinA*ones(numel(subrange1_y),numel(subrange1_x));% default value=min of image A - image2_crop=MinA*ones(numel(subrange2_y),numel(subrange2_x));% default value=min of image A - check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1 - check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight; - check2_x=subrange2_x>=1 & subrange2_x<=par_civ.ImageWidth;% check which points in the subimage 2 are contained in the initial image 2 - check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight; - image1_crop(check1_y,check1_x)=par_civ.ImageA(subrange1_y(check1_y),subrange1_x(check1_x));%extract a subimage (correlation box) from image A - image2_crop(check2_y,check2_x)=par_civ.ImageB(subrange2_y(check2_y),subrange2_x(check2_x));%extract a larger subimage (search box) from image B - if checkmask - mask1_crop=ones(numel(subrange1_y),numel(subrange1_x));% default value=1 for mask - mask2_crop=ones(numel(subrange2_y),numel(subrange2_x));% default value=1 for mask - mask1_crop(check1_y,check1_x)=check_undefined(subrange1_y(check1_y),subrange1_x(check1_x));%extract a mask subimage (correlation box) from image A - mask2_crop(check2_y,check2_x)=check_undefined(subrange2_y(check2_y),subrange2_x(check2_x));%extract a mask subimage (search box) from image B - sizemask=sum(sum(mask1_crop))/(numel(subrange1_y)*numel(subrange1_x));%size of the masked part relative to the correlation sub-image - if sizemask > 1/2% eliminate point if more than half of the correlation box is masked - F(ivec)=3; % - utable(ivec)=0; - vtable(ivec)=0; + for ivec_x=1:nbvec_x + for ivec_y=1:nbvec_y + ivec_y + iref=round(par_civ.Grid(ivec_y,ivec_x,1)+0.5)% xindex on the image A for the middle of the correlation box + jref=round(par_civ.ImageHeight-par_civ.Grid(ivec_y,ivec_x,2)+0.5)% j index for the middle of the correlation box in the image A + subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage + subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage + subrange2_x=iref+shiftx(ivec_y,ivec_x)-isx2:iref+shiftx(ivec_y,ivec_x)+isx2;%x indices defining the second subimage + subrange2_y=jref+shifty(ivec_y,ivec_x)-isy2:jref+shifty(ivec_y,ivec_x)+isy2;%y indices defining the second subimage + image1_crop=MinA*ones(npz,numel(subrange1_y),numel(subrange1_x));% default value=min of image A + image2_crop=MinA*ones(npz,numel(subrange2_y),numel(subrange2_x));% default value=min of image A + check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1 + check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight; + check2_x=subrange2_x>=1 & subrange2_x<=par_civ.ImageWidth;% check which points in the subimage 2 are contained in the initial image 2 + check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight; + image1_crop(:,check1_y,check1_x)=par_civ.ImageA(:,subrange1_y(check1_y),subrange1_x(check1_x));%extract a subimage (correlation box) from image A + image2_crop(:,check2_y,check2_x)=par_civ.ImageB(:,subrange2_y(check2_y),subrange2_x(check2_x));%extract a larger subimage (search box) from image B + if checkmask + mask1_crop=ones(numel(subrange1_y),numel(subrange1_x));% default value=1 for mask + mask2_crop=ones(numel(subrange2_y),numel(subrange2_x));% default value=1 for mask + mask1_crop(check1_y,check1_x)=check_undefined(subrange1_y(check1_y),subrange1_x(check1_x));%extract a mask subimage (correlation box) from image A + mask2_crop(check2_y,check2_x)=check_undefined(subrange2_y(check2_y),subrange2_x(check2_x));%extract a mask subimage (search box) from image B + sizemask=sum(sum(mask1_crop))/(numel(subrange1_y)*numel(subrange1_x));%size of the masked part relative to the correlation sub-image + if sizemask > 1/2% eliminate point if more than half of the correlation box is masked + FF(ivec_y,ivec_x)=1; % + utable(ivec_y,ivec_x)=NaN; + vtable(ivec_y,ivec_x)=NaN; + else + FF(ivec_y,ivec_x)=0; + image1_crop=image1_crop.*~mask1_crop;% put to zero the masked pixels (mask1_crop='true'=1) + image2_crop=image2_crop.*~mask2_crop; + image1_mean=mean(mean(image1_crop))/(1-sizemask); + image2_mean=mean(mean(image2_crop))/(1-sizemask); + end else - image1_crop=image1_crop.*~mask1_crop;% put to zero the masked pixels (mask1_crop='true'=1) - image2_crop=image2_crop.*~mask2_crop; - image1_mean=mean(mean(image1_crop))/(1-sizemask); - image2_mean=mean(mean(image2_crop))/(1-sizemask); - end - else - image1_mean=mean(mean(image1_crop)); - image2_mean=mean(mean(image2_crop)); - end - %threshold on image minimum - if F(ivec)~=3 - if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma) - F(ivec)=3; - %threshold on image maximum - elseif check_MaxIma && (image1_mean > par_civ.MaxIma || image2_mean > par_civ.MaxIma) - F(ivec)=3; - end - if F(ivec)==3 - utable(ivec)=0; - vtable(ivec)=0; - else - %mask - if checkmask - image1_crop=(image1_crop-image1_mean).*~mask1_crop;%substract the mean, put to zero the masked parts - image2_crop=(image2_crop-image2_mean).*~mask2_crop; + image1_mean=mean(mean(image1_crop)); + image2_mean=mean(mean(image2_crop)); + end + %threshold on image minimum + if FF(ivec_y,ivec_x)==0 + if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma) + FF(ivec_y,ivec_x)=1; + %threshold on image maximum + elseif check_MaxIma && (image1_mean > par_civ.MaxIma || image2_mean > par_civ.MaxIma) + FF(ivec_y,ivec_x)=1; + end + if FF(ivec_y,ivec_x)==1 + utable(ivec_y,ivec_x)=NaN; + vtable(ivec_y,ivec_x)=NaN; else - image1_crop=(image1_crop-image1_mean); - image2_crop=(image2_crop-image2_mean); - end - %deformation - if CheckDeformation - xi=(1:mesh:size(image1_crop,2)); - yi=(1:mesh:size(image1_crop,1))'; - [XI,YI]=meshgrid(xi-ceil(size(image1_crop,2)/2),yi-ceil(size(image1_crop,1)/2)); - XIant=XI-par_civ.DUDX(ivec)*XI+par_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2); - YIant=YI+par_civ.DVDX(ivec)*XI-par_civ.DVDY(ivec)*YI+ceil(size(image1_crop,1)/2); - image1_crop=interp2(image1_crop,XIant,YIant); - image1_crop(isnan(image1_crop))=0; - xi=(1:mesh:size(image2_crop,2)); - yi=(1:mesh:size(image2_crop,1))'; - image2_crop=interp2(image2_crop,xi,yi,'*spline'); - image2_crop(isnan(image2_crop))=0; - end - sum_square=sum(sum(image1_crop.*image1_crop)); - %reference: Oliver Pust, PIV: Direct Cross-Correlation - result_conv= conv2(image2_crop,flipdim(flipdim(image1_crop,2),1),'valid'); - corrmax= max(max(result_conv)); - result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255 - %Find the correlation max, at 255 - [y,x] = find(result_conv==255,1); - subimage2_crop=image2_crop(y:y+2*iby2/mesh,x:x+2*ibx2/mesh);%subimage of image 2 corresponding to the optimum displacement of first image - sum_square=sum_square*sum(sum(subimage2_crop.*subimage2_crop));% product of variances of image 1 and 2 - sum_square=sqrt(sum_square);% srt of the variance product to normalise correlation - if ~isempty(y) && ~isempty(x) - try - if par_civ.CorrSmooth==1 - [vector,F(ivec)] = SUBPIXGAUSS (result_conv,x,y); - elseif par_civ.CorrSmooth==2 - [vector,F(ivec)] = SUBPIX2DGAUSS (result_conv,x,y); - else - [vector,F(ivec)] = quadr_fit(result_conv,x,y); - end - utable(ivec)=vector(1)*mesh+shiftx(ivec); - vtable(ivec)=vector(2)*mesh+shifty(ivec); - xtable(ivec)=iref+utable(ivec)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2) - ytable(ivec)=jref+vtable(ivec)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge) - iref=round(xtable(ivec)+0.5);% nearest image index for the middle of the vector - jref=round(ytable(ivec)+0.5); - % eliminate vectors located in the mask - if checkmask && (iref<1 || jref<1 ||iref>npx_ima || jref>npy_ima ||( par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100)) - utable(ivec)=0; - vtable(ivec)=0; - F(ivec)=3; - end - ctable(ivec)=corrmax/sum_square;% correlation value - catch ME - F(ivec)=3; - disp(ME.message) + %mask + if checkmask + image1_crop=(image1_crop-image1_mean).*~mask1_crop;%substract the mean, put to zero the masked parts + image2_crop=(image2_crop-image2_mean).*~mask2_crop; + else + image1_crop=(image1_crop-image1_mean); + image2_crop=(image2_crop-image2_mean); end - else - F(ivec)=3; + + + %reference: Oliver Pust, PIV: Direct Cross-Correlation + for kz=1:par_civ.SearchBoxSize(3) + subima2=squeeze(image2_crop(kz,:,:)); + subima1=squeeze(image1_crop(kref,:,:)); + correl_xy=conv2(subima2,flip(flip(subima1,2),1),'valid'); + result_conv(kz,:,:)= correl_xy; + max_xy(kz)=max(max(correl_xy)); + [xk(kz),yk(kz)]=find(correl_xy==max_xy(kz),1); + + end + [corrmax,z]=max(max_xy); + + x=xk(z); + y=yk(z); + result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255 + 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 + sum_square=sum(sum(squeeze(image1_crop(z,:,:).*image1_crop(z,:,:)))); + sum_square=sum_square*sum(sum(subimage2_crop.*subimage2_crop));% product of variances of image 1 and 2 + sum_square=sqrt(sum_square);% srt of the variance product to normalise correlation + if ~isempty(y) && ~isempty(x) + + if par_civ.CorrSmooth==1 + [vector,FF(ivec_y,ivec_x)] = SUBPIXGAUSS (result_conv(z,:,:),x,y);%TODO: improve by max optimisation along z + elseif par_civ.CorrSmooth==2 + [vector,FF(ivec_y,ivec_x)] = SUBPIX2DGAUSS (result_conv(z,:,:),x,y); + else + [vector,FF(ivec_y,ivec_x)] = quadr_fit(result_conv(z,:,:),x,y); + end + utable(ivec_y,ivec_x)=vector(1)*mesh+shiftx(ivec_y,ivec_x); + vtable(ivec_y,ivec_x)=vector(2)*mesh+shifty(ivec_y,ivec_x); + xtable(ivec_y,ivec_x)=iref+utable(ivec_y,ivec_x)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2) + ytable(ivec_y,ivec_x)=jref+vtable(ivec_y,ivec_x)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge) + iref=round(xtable(ivec_y,ivec_x)+0.5);% nearest image index for the middle of the vector + jref=round(ytable(ivec_y,ivec_x)+0.5); + wtable(ivec_y,ivec_x)=z-kref; + % eliminate vectors located in the mask + if checkmask && (iref<1 || jref<1 ||iref>npx_ima || jref>npy_ima ||( par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100)) + utable(ivec_y,ivec_x)=0; + vtable(ivec_y,ivec_x)=0; + FF(ivec_y,ivec_x)=1; + end + ctable(ivec_y,ivec_x)=corrmax/sum_square;% correlation value + + else + FF(ivec_y,ivec_x)=1; + end end end @@ -1204,5 +982,5 @@ peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2); else - F=-2; % warning flag for vector truncated by the limited search box + F=1; % warning flag for vector truncated by the limited search box end peakx=x; @@ -1213,5 +991,5 @@ peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2); else - F=-2; % warning flag for vector truncated by the limited search box + F=1; % warning flag for vector truncated by the limited search box end vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1]; @@ -1222,5 +1000,5 @@ %------------------------------------------------------------------------ % vector=[0 0]; %default -F=-2; +F=1; peaky=y; peakx=x; @@ -1265,5 +1043,5 @@ [npy,npx]=size(result_conv); if x<4 || y<4 || npx-x<4 ||npy-y <4 - F=-2; + F=1; vector=[x y]; else @@ -1298,51 +1076,26 @@ end -%'RUN_FIX': function for fixing velocity fields: -%----------------------------------------------- -% RUN_FIX(filename,field,flagindex,thresh_vecC,thresh_vel,iter,flag_mask,maskname,fileref,fieldref) -% -%filename: name of the netcdf file (used as input and output) -%field: structure specifying the names of the fields to fix (depending on civ1 or civ2) -%.vel_type='civ1' or 'civ2'; -%.nb=name of the dimension common to the field to fix ('nb_vectors' for civ1); -%.fixflag=name of fix flag variable ('vec_FixFlag' for civ1) -%flagindex: flag specifying which values of vec_f are removed: -% if flagindex(1)=1: vec_f=-2 vectors are removed -% if flagindex(2)=1: vec_f=3 vectors are removed -% if flagindex(3)=1: vec_f=2 vectors are removed (if iter=1) or vec_f=4 vectors are removed (if iter=2) -%iter=1 for civ1 fields and iter=2 for civ2 fields -%thresh_vecC: threshold in the image correlation vec_C -%flag_mask: =1 mask used to remove vectors (0 else) -%maskname: name of the mask image file for fix -%thresh_vel: threshold on velocity, or on the difference with the reference file fileref if exists -%inf_sup=1: remove values smaller than threshold thresh_vel, =2, larger than threshold -%fileref: .nc file name for a reference velocity (='': refrence 0 used) -%fieldref: 'civ1','filter1'...feld used in fileref - -function FF=fix(Param,F,C,U,V,X,Y) -FF=zeros(size(F));%default - -%criterium on warn flags -FlagName={'CheckFmin2','CheckF2','CheckF3','CheckF4'}; -FlagVal=[-2 2 3 4]; -for iflag=1:numel(FlagName) - if isfield(Param,FlagName{iflag}) && Param.(FlagName{iflag}) - FF=(FF==1| F==FlagVal(iflag)); - end -end -%criterium on correlation values + +function FF=detect_false(Param,C,U,V,FFIn) +FF=FFIn;%default, good vectors +% FF=1, for correlation max at edge, not set in this function +% FF=2, for too small correlation +% FF=3, for velocity outside bounds +% FF=4 for exclusion by difference with the smoothed field, not set in this function + if isfield (Param,'MinCorr') - FF=FF==1 | C(Param.MaxVel*Param.MaxVel); - end -end - + U2Max=Param.MaxVel*Param.MaxVel; + FF(Umod>U2Max & FFIn==0)=3; + end +end %------------------------------------------------------------------------ Index: trunk/src/series/civ_input.m =================================================================== --- trunk/src/series/civ_input.m (revision 1147) +++ trunk/src/series/civ_input.m (revision 1148) @@ -67,21 +67,4 @@ %% set visibility options depending on the calling function (Param.Action.ActionName): -if strcmp(Param.Action.ActionName,'civ_series')||strcmp(Param.Action.ActionName,'stereo_civ') - set(handles.num_MaxDiff,'Visible','on') - set(handles.num_MaxVel,'Visible','on') - set(handles.title_MaxVel,'Visible','on') - set(handles.title_MaxDiff,'Visible','on') - set(handles.num_Nx,'Visible','off') - set(handles.num_Ny,'Visible','off') - set(handles.title_Nx,'Visible','off') - set(handles.title_Ny,'Visible','off') - set(handles.num_CorrSmooth,'Style','popupmenu') - set(handles.num_CorrSmooth,'Value',1) - set(handles.num_CorrSmooth,'String',{'1';'2'}) - set(handles.CheckThreshold,'Visible','on') - set(handles.CheckDeformation,'Value',0)% desactivate -% set(handles.num_SubDomainSize(1),'String','250') -% set(handles.num_SubDomainSize(2),'String','500') -end switch Param.Action.ActionName case 'stereo_civ' @@ -91,4 +74,13 @@ set(handles.ListCompareMode,'Visible','on') set(handles.PairIndices,'Visible','on') + case 'civ_3D' + set(handles.ListCompareMode,'Visible','on') + set(handles.PairIndices,'Visible','on') + set(handles.title_z,'Visible','on') + set(handles.num_CorrBoxSize_3,'Visible','on') + set(handles.num_SearchBoxSize_3,'Visible','on') + set(handles.num_SearchBoxShift_3,'Visible','on') + set(handles.num_Dz,'Visible','on') + set(handles.title_Dz,'Visible','on') end @@ -157,12 +149,4 @@ end end - -%% reinitialise menus -set(handles.ListPairMode,'Value',1) -set(handles.ListPairMode,'String',{''}) -set(handles.ListPairCiv1,'Value',1) -set(handles.ListPairCiv1,'String',{''}) -set(handles.ListPairCiv2,'Value',1) -set(handles.ListPairCiv2,'String',{''}) %% prepare the GUI with input parameters @@ -321,6 +305,14 @@ end +%% reinitialise pair menus +set(handles.ListPairMode,'Value',1) +set(handles.ListPairMode,'String',{''}) +set(handles.ListPairCiv1,'Value',1) +set(handles.ListPairCiv1,'String',{''}) +set(handles.ListPairCiv2,'Value',1) +set(handles.ListPairCiv2,'String',{''}) + %% set the menu and default choice of civ pairs -if isequal(MaxIndex_j,MinIndex_j)% no possibility of j pairs +if isequal(MaxIndex_j,MinIndex_j)|| strcmp(Param.Action.ActionName,'civ_3D')% no possibility of j pairs PairMenu={'series(Di)'}; elseif MaxIndex_j-MinIndex_j==1 @@ -338,22 +330,26 @@ PairIndex=find(strcmp(Param.ActionInput.PairIndices.ListPairMode,PairMenu));%retrieve the previous option end -if isempty(PairIndex) - if ~isfield(Param.IndexRange,'first_j')||isequal(MaxIndex_j,MinIndex_j)% no possibility of j pairs - PairIndex=1; - elseif MaxIndex_i==1 && MaxIndex_j>1% simple series in j - if MaxIndex_j <= 10 - PairIndex=1;% advice 'pair j1-j2' except in MaxIndex_j is large - end - else - if strcmp(NomTypeNc,'_1-2_1') - PairIndex=3;% advise 'series(Di)' - elseif MaxIndex_j <= 10 - PairIndex=1;% advice 'pair j1-j2' except in MaxIndex_j is large +if strcmp(Param.Action.ActionName,'civ_3D') + PairIndex=1 +else + if isempty(PairIndex) + if ~isfield(Param.IndexRange,'first_j')||isequal(MaxIndex_j,MinIndex_j)% no possibility of j pairs + PairIndex=1; + elseif MaxIndex_i==1 && MaxIndex_j>1% simple series in j + if MaxIndex_j <= 10 + PairIndex=1;% advice 'pair j1-j2' except in MaxIndex_j is large + end else - PairIndex=2;% advice 'Dj' - end - end -end -set(handles.ListPairMode,'Value',PairIndex); + if strcmp(NomTypeNc,'_1-2_1') + PairIndex=3;% advise 'series(Di)' + elseif MaxIndex_j <= 10 + PairIndex=1;% advice 'pair j1-j2' except in MaxIndex_j is large + else + PairIndex=2;% advice 'Dj' + end + end + end +end +set(handles.ListPairMode,'Value',PairIndex); %% indicate the min and max indices i and j on the GUI @@ -407,6 +403,6 @@ %% Civ1 parameters %Param.CheckCiv1=1; -Param.Civ1.CorrBoxSize=[25 25]; -Param.Civ1.SearchBoxSize=[55 55]; +Param.Civ1.CorrBoxSize=[25 25 1]; +Param.Civ1.SearchBoxSize=[55 55 5]; Param.Civ1.SearchBoxShift=[0 0]; Param.Civ1.CorrSmooth=1; @@ -420,15 +416,11 @@ %% Fix1 parameters -%Param.CheckFix1=1; -Param.Fix1.CheckFmin2=1; -Param.Fix1.CheckF3=1; Param.Fix1.MinCorr=0.2000; %% Patch1 parameters %Param.CheckPatch1=1; -Param.Patch1.FieldSmooth=10; -Param.Patch1.MaxDiff=1.5000; -Param.Patch1.SubDomainSize=250; -Param.Patch1.TestPatch1=0; +Param.Patch1.FieldSmooth=20; +Param.Patch1.MaxDiff=2; +Param.Patch1.SubDomainSize=125; %% Civ2 parameters @@ -443,18 +435,12 @@ Param.Civ2.Mask=''; Param.Civ2.CheckThreshold=0; -Param.Civ2.TestCiv2=0; %% Fix2 parameters -%Param.CheckFix2=1; -Param.Fix2.CheckFmin2=1; -Param.Fix2.CheckF4=1; -Param.Fix2.CheckF3=1; Param.Fix2.MinCorr=0.2000; %% Patch2 parameters -%Param.CheckPatch2=1; -Param.Patch2.FieldSmooth=2; +Param.Patch2.FieldSmooth=5; Param.Patch2.MaxDiff=1.5000; -Param.Patch2.SubDomainSize=500; +Param.Patch2.SubDomainSize=250; Param.Patch2.TestPatch2=0; @@ -611,5 +597,5 @@ checkeven=(mod(ActionInput.Civ1.CorrBoxSize,2)==0); ActionInput.Civ1.CorrBoxSize(checkeven)=ActionInput.Civ1.CorrBoxSize(checkeven)+1;% set correlation box sizes to odd values - ActionInput.Civ1.SearchBoxSize=max(ActionInput.Civ1.SearchBoxSize,ActionInput.Civ1.CorrBoxSize+8);% insure that the search box size is large enough + ActionInput.Civ1.SearchBoxSize(1:2)=max(ActionInput.Civ1.SearchBoxSize(1:2),ActionInput.Civ1.CorrBoxSize(1:2)+8);% insure that the search box size is large enough checkeven=(mod(ActionInput.Civ1.SearchBoxSize,2)==0); ActionInput.Civ1.SearchBoxSize(checkeven)=ActionInput.Civ1.SearchBoxSize(checkeven)+1;% set search box sizes to odd values @@ -669,9 +655,9 @@ case 'PIV' PairIndices='on';% needs to define index pairs for PIV - case 'PIV volume' - PairIndices='on';% needs to define index pairs for PIV - set(handles.ListPairMode,'Value',1) - set(handles.ListPairMode,'String',{'series(Di)'}) - ListPairMode_Callback(hObject, eventdata, handles) + % case 'PIV volume' + % PairIndices='on';% needs to define index pairs for PIV + % set(handles.ListPairMode,'Value',1) + % set(handles.ListPairMode,'String',{'series(Di)'}) + % ListPairMode_Callback(hObject, eventdata, handles) case 'displacement' OriginIndex='on';%define a frame origin for displacement @@ -2363,2 +2349,18 @@ drawnow end + + +function num_CorrBoxSize_3_Callback(hObject, eventdata, handles) + + +function num_SearchBoxSize_3_Callback(hObject, eventdata, handles) + + +function MinIndex_j_Callback(hObject, eventdata, handles) + + +% --- Executes on selection change in field_ref2. +function field_ref2_Callback(hObject, eventdata, handles) + + + Index: trunk/src/series/civ_series.m =================================================================== --- trunk/src/series/civ_series.m (revision 1147) +++ trunk/src/series/civ_series.m (revision 1148) @@ -56,8 +56,10 @@ Data=civ_input(Param);% introduce the civ parameters using the GUI civ_input + % TODO: change from guide to App: modify the input procedure, adapt read_GUI function + %App=civ_input_App %Data=civ_input_App(Param);% introduce the civ parameters using the GUI civ_input - if isempty(Data) - Data=Param;% if civ_input has been cancelled, keep previous parameters - end + % if isempty(App) + % Data=Param;% if civ_input has been cancelled, keep previous parameters + % end Data.Program=mfilename;%gives the name of the current function Data.AllowInputSort='off';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default) Index: trunk/src/series/extract_rdvision.m =================================================================== --- trunk/src/series/extract_rdvision.m (revision 1147) +++ trunk/src/series/extract_rdvision.m (revision 1148) @@ -1,3 +1,3 @@ -%'extunningact_rdvision': relabel an image series with two indices, and correct errors from the RDvision transfer program +%'extract_rdvision': relabel an image series with two indices, and correct errors from the RDvision transfer program %------------------------------------------------------------------------ % function ParamOut=extract_rdvision(Param) @@ -152,5 +152,5 @@ end difftime=timestamp-timexml; - if max(difftime)>0.01 + %if max(difftime)>0.01 figure plot(timestamp,difftime) @@ -158,5 +158,5 @@ ylabel('time difference(s)') title('discrepency timestamps-timexml') - end + %end return end @@ -308,4 +308,12 @@ % m.Data=data; %%%%%%% + [XmlData,errormsg]=imadoc2struct(filexml);% check reading of the xml file + if ~isempty(errormsg) + disp(errormsg) + return + end + [nbfield1,nbfield2]=size(XmlData.Time); + nbfield1=nbfield1-1;nbfield2=nbfield2-1; + timestamp=zeros(1,numel(m.Data)); for ii=1: numel(m.Data) @@ -313,10 +321,7 @@ end if isequal(Param.IndexRange.first_i,1) - [nbfield1,nbfield2,msg]=copyfile_modif(filexml,timestamp,newxml); %copy the xml file in the upper folder - [XmlData,errormsg]=imadoc2struct(newxml);% check reading of the new xml file - if ~isempty(errormsg) - disp(errormsg) - return - end + %[nbfield1,nbfield2,msg]=copyfile_modif(filexml,timestamp,newxml); %copy the xml file in the upper folder + % [XmlData,errormsg]=imadoc2struct(newxml);% check reading of the new xml file + if numel(timestamp)~=nbfield1*nbfield2 disp('WARNING: total image number defined by the xml file differs from the number of frames ') @@ -329,6 +334,6 @@ end end - else - [nbfield1,nbfield2,msg]=copyfile_modif(filexml,timestamp,''); + % else + % [nbfield1,nbfield2,msg]=copyfile_modif(filexml,timestamp,''); end if nbfield2>1 @@ -468,9 +473,9 @@ fclose(fid) end - - - - -function [nbfield1,nbfield2,msg]=copyfile_modif(filexml,timestamp,newxml) +disp('END EXTRACT') + + +%OBSOLETE +function [nbfield1,nbfield2,msg]=copyfile_modif(filexml,timestamp,newxml)% msg=''; t=xmltree(filexml);