source: trunk/src/series/merge_proj.m @ 29

Last change on this file since 29 was 29, checked in by sommeria, 11 years ago

read_imadoc suppressed (obsolete, replaced by imadoc2struct)
update_imadoc: bug repaired: existing xml file was erased
various cleaning (deals with non existing input file for series and uvmat)

File size: 14.3 KB
Line 
1function GUI_input=merge_proj(num_i1,num_i2,num_j1,num_j2,Series);
2
3%requests for the visibility of input windows in the GUI series  (activated directly by the selection in the menu ACTION)
4if ~exist('num_i1','var')
5    GUI_input={'RootPath';'two';...%nbre of possible input series (options 'on'/'two'/'many', default:'one')
6        'SubDir';'on';... % subdirectory of derived files (PIV fields), ('on' by default)
7        'RootFile';'on';... %root input file name ('on' by default)
8        'FileExt';'on';... %input file extension ('on' by default)
9        'NomType';'on';...%type of file indexing ('on' by default)
10        'NbSlice';'on'; ...%nbre of slices ('off' by default)
11        'VelTypeMenu';'one';...% menu for selecting the velocity type (civ1,..) options 'off'/'one'/'two', 'off' by default)
12        'FieldMenu';'one';...% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
13        'CoordType';'on';...%can use a transform function 'off' by default
14        'GetObject';'on';...%can use projection object ,'off' by default
15        %'GetMask';'on'...%can use mask option   ,'off' by default
16        %'PARAMETER'; options: name of the user defined parameter',repeat a line for each parameter
17               ''};
18    return %exit the function
19end
20
21%-------------------------------------------------
22hseries=guidata(Series.hseries);%handles of the GUI series
23WaitbarPos=get(hseries.waitbar_frame,'Position'); %positiopn of waitbar frame
24%-------------------------------------------------
25
26%numbers of view fields (nbre of inputs in RootPath)
27testcell=iscell(Series.RootFile);
28if ~testcell
29    Series.RootPath={Series.RootPath};
30    Series.RootFile={Series.RootFile};
31    Series.SubDir={Series.SubDir};
32    Series.FileExt={Series.FileExt};
33    Series.NomType={Series.NomType};
34    num_i1={num_i1};
35    num_i2={num_i2};
36    num_j1={num_j1};
37    num_j2={num_j2};
38end
39nbview=length(Series.RootFile);%number of views (file series to merge)
40nbfield=size(num_i1{1},1)*size(num_i1{1},2);%number of fields in the time series
41transform=Series.CoordType; %  field transform function
42hhh=which('mmreader');
43for iview=1:nbview
44    test_movie(iview)=0;
45    if ~isequal(hhh,'')&& mmreader.isPlatformSupported()
46        if isequal(lower(FileExt{iview}),'.avi')
47            MovieObject{iview}=mmreader(fullfile(RootPath{iview},[RootFile{iview} FileExt{iview}]));
48            test_movie(iview)=1;
49        end
50    end
51end
52
53%Calibration data and timing: read the ImaDoc files
54mode=''; %default
55timecell={};
56itime=0;
57NbSlice_calib={}; %test for z index
58for iview=1:nbview%Loop on views
59    XmlData{iview}=[];%default
60    filebase{iview}=fullfile(Series.RootPath{iview},Series.RootFile{iview});
61    if exist([filebase{iview} '.xml'],'file')
62        [XmlData{iview},error]=imadoc2struct([filebase{iview} '.xml']);
63        if isfield(XmlData{iview},'Time')
64            itime=itime+1;
65            timecell{itime}=XmlData{iview}.Time;
66        end
67        if isfield(XmlData{iview},'GeometryCalib') && isfield(XmlData{iview}.GeometryCalib,'SliceCoord')
68            NbSlice_calib{iview}=size(XmlData{iview}.GeometryCalib.SliceCoord,1);
69            if ~isequal(NbSlice_calib{iview},NbSlice_calib{1})
70                msgbox_uvmat('WARNING','inconsistent number of Z indices for the two field series');
71            end
72        end   
73    elseif exist([filebase{iview} '.civ'],'file')
74        [error,time,TimeUnit,mode,npx,npy,pxcmx,pxcmy]=read_imatext([filebase{iview} '.civ']);
75        itime=itime+1;
76        timecell{itime}=time;
77        XmlData{iview}.Time=time;
78        GeometryCalib.R=[pxcmx 0 0; 0 pxcmy 0;0 0 0];
79        GeometryCalib.Tx=0;
80        GeometryCalib.Ty=0;
81        GeometryCalib.Tz=1;
82        GeometryCalib.dpx=1;
83        GeometryCalib.dpy=1;
84        GeometryCalib.sx=1;
85        GeometryCalib.Cx=0;
86        GeometryCalib.Cy=0;
87        GeometryCalib.f=1;
88        GeometryCalib.kappa1=0;
89        GeometryCalib.CoordUnit='cm';
90        XmlData{iview}.GeometryCalib=GeometryCalib;
91        if error==1
92            msgbox_uvmat('WARNING','inconsistent number of fields in the .civ file');
93        end
94    end
95end
96
97%check coincidence in time
98multitime=0;
99if length(timecell)==0
100    time=[];
101elseif length(timecell)==1
102    time=timecell{1};
103elseif length(timecell)>1
104    multitime=1;
105    for icell=1:length(timecell)
106        if ~isequal(size(timecell{icell}),size(timecell{1}))
107            msgbox_uvmat('WARNING','inconsistent time array dimensions in ImaDoc fields, the time for the first series is used')
108            time=timecell{1};
109            multitime=0;
110            break
111        end
112    end
113end
114if multitime
115    for icell=1:length(timecell)
116        time(icell,:,:)=timecell{icell};
117    end
118    diff_time=max(max(diff(time)));
119    if diff_time>0
120        msgbox_uvmat('WARNING',['times of series differ by more than ' num2str(diff_time)])
121    end   
122end
123if size(time,2) < num_i2{1}(end) || size(time,3) < num_j2{1}(end)% ime array absent or too short in ImaDoc xml file'
124    time=[];
125end
126
127% Field and velocity type (the same for all views)
128Field_str=get(hseries.FieldMenu,'String');
129val=get(hseries.FieldMenu,'Value');
130FieldName=Field_str(val);%the same set of fields for all views
131VelType_str=get(hseries.VelTypeMenu,'String');
132VelType_val=get(hseries.VelTypeMenu,'Value');
133VelType=VelType_str{VelType_val}; %the same for all views
134if isequal(FieldName,'get_field...')
135    hget_field=findobj(allchild(0),'Name','get_field');%find the get_field... GUI
136   % hhget_field=guidata(hget_field);%handles of GUI elements in get_field
137    SubField=get_field('read_get_field',hObject,eventdata,hget_field); %read the names of the variables to plot in the get_field GUI
138%     if isequal(get(hhget_field.menu_coord,'Visible'),'on')
139%         list_transform=get(hhget_field.menu_coord,'String');
140%         val_list=get(hhget_field.menu_coord,'Value');
141%         transform=list_transform{val_list};
142%     end
143end
144%detect whether all the files are 'images' or 'netcdf'
145testima=0;
146testvol=0;
147testcivx=0;
148testnc=0;
149FileExt=get(hseries.FileExt,'String');
150for iview=1:nbview
151     ext=FileExt{iview};
152     form=imformats(ext([2:end]));
153     if isequal(lower(ext),'.vol')
154         testvol=testvol+1;
155     elseif ~isempty(form)||isequal(lower(ext),'.avi')% if the extension corresponds to an image format recognized by Matlab
156         testima=testima+1;
157     elseif isequal(ext,'.nc')
158         testnc=testnc+1;
159     end
160end
161if testvol
162    msgbox_uvmat('ERROR','volume images not implemented yet')
163    return
164end
165if testnc~=nbview && testima~=nbview && testvol~=nbview
166    msgbox_uvmat('ERROR','need a set of images or a set of netcdf files with the same fields as input')
167    return
168end
169if ~isequal(FieldName,'get_field...')
170    testcivx=testnc;
171end
172%name of output files and directory:
173% res_subdir=fullfile(Series.RootPath{1},[Series.SubDir{1} '_STAT']);
174ProjectDir=fileparts(fileparts(Series.RootPath{1}));% preoject directory (GERK)
175prompt={['result directory (in' ProjectDir ')']};
176RootPath=get(hseries.RootPath,'String');
177SubDir=get(hseries.SubDir,'String');
178if isequal(length(RootPath),1)
179    fulldir=RootPath{1};
180    subdir='GRID';
181    res_subdir=fullfile(fulldir,subdir);
182else
183    def={fullfile(ProjectDir,'0_RESULTS')};
184    dlgTitle='result directory';
185    lineNo=1;
186    answer=msgbox_uvmat('INPUT_TXT',dlgTitle,def);
187    fulldir=answer{1};
188    subdir=[];
189    dirlist=sort(Series.RootFile);
190    for iview=1:nbview
191        if ~isempty(subdir)
192            subdir=[subdir '-'];
193        end
194        subdir=[subdir dirlist{iview}];
195    end 
196    res_subdir=fullfile(fulldir,subdir);
197end
198ext=FileExt{1};
199if ~exist(fulldir,'dir')
200    msgbox_uvmat('ERROR',['directory ' fulldir ' needs to be created'])
201    return
202end
203if ~exist(res_subdir,'dir')
204    dircur=pwd;
205    cd(fulldir)
206    error=mkdir(subdir);
207    cd(dircur)
208end
209filebasesub=fullfile(res_subdir,Series.RootFile{1});
210filebase_merge=fullfile(res_subdir,'merged');%root name for the merged files
211
212%projection object
213if isfield(Series,'sethandles')
214    if ishandle(Series.sethandles.set_object)
215        Series.ProjObject=read_set_object(Series.sethandles);
216        if ~isfield(Series.ProjObject,'Style')
217            msgbox_uvmat('ERROR','Undefined projection object style')
218            return
219        end
220        if ~isequal(Series.ProjObject.Style,'plane')
221            msgbox_uvmat('ERROR','The projection object must be a plane')
222            return
223        end
224    end
225end
226
227    %MAIN LOOP
228for ifile=1:nbfield               
229    stopstate=get(hseries.RUN,'BusyAction');
230    if isequal(stopstate,'queue')% enable STOP command from the 'series' interface
231         update_waitbar(hseries.waitbar,WaitbarPos,ifile/nbfield)
232         Amerge=0;
233         
234         %----------LOOP ON VIEWS----------------------
235        nbtime=0;
236        for iview=1:nbview
237            %name of the current file
238            filename=name_generator(filebase{iview},num_i1{iview}(ifile),num_j1{iview}(ifile),Series.FileExt{iview},Series.NomType{iview},1,num_i2{iview}(ifile),num_j2{iview}(ifile),SubDir{iview});
239            if ~exist(filename,'file')
240                msgbox_uvmat('ERROR',['missing input file' filename])
241                break
242            end
243
244            %reading the current file
245            if testima
246                if test_movie(iview)
247                    Field{iview}.A=read(MovieObject{iview},num_i1{iview}(ifile));
248                else
249                    Field{iview}.A=read_image(filename,Series.NomType{iview},num_i1{iview}(ifile));
250                end % TODO: introduce ListVarName
251                npxy=size(Field{iview}.A);
252                Field{iview}.AX=[0.5 npxy(2)-0.5]; % coordinates of the first and last pixel centers
253                Field{iview}.AY=[npxy(1)-0.5 0.5];
254                Field{iview}.CoordType='px';
255                Field{iview}.AName='image';
256            else
257                if testcivx
258                    [Field{iview},VelTypeOut]=read_civxdata(filename,FieldName,VelType);
259                else
260                    [Field{iview},var_detect]=nc2struct(filename,SubField.ListVarName); %read the corresponding input data               
261                    Field{iview}.VarAttribute=SubField.VarAttribute;
262                end
263                if isfield(Field{iview},'Time')
264                    timeread(iview)=Field{iview}.Time;
265                    nbtime=nbtime+1;
266                end
267            end
268            % coord transform
269            % z index
270            if ~isempty(NbSlice_calib)
271                Field{iview}.ZIndex=mod(num_i1{iview}(ifile)-1,NbSlice_calib{1})+1;
272            end
273            if ~isequal(transform,'')
274                Field{iview}=feval(Series.CoordType,Field{iview},XmlData{iview});%transform to phys if requested
275            end
276            if testcivx
277                    Field{iview}=calc_field(FieldName,Field{iview});
278            end
279
280            %projection on object (gridded plane)
281            if isfield(Series,'ProjObject')
282                Field{iview}=proj_field(Field{iview},Series.ProjObject);
283            end
284        end   
285       
286         %----------END LOOP ON VIEWS----------------------
287         
288        %merge the nbview fields
289        MergeData=merge_field(Field);
290        if isfield(MergeData,'Txt')
291            msgbox_uvmat('ERROR',MergeData.Txt)
292            return
293        end
294       
295        % generating the name of the merged field
296        mergename=name_generator(filebase_merge,num_i1{iview}(ifile),num_j1{iview}(ifile),Series.FileExt{iview},Series.NomType{iview},1,num_i2{iview}(ifile),num_j2{iview}(ifile));
297       
298        % time:
299        time_i=0;%default
300        if isempty(time)% time from ImaDoc prevails
301            time_i=sum(timeread)/nbtime;
302        else
303            time_i=(time(iview,num_i1{iview}(ifile),num_j1{iview}(ifile))+time(iview,num_i2{iview}(ifile),num_j2{iview}(ifile)))/2;
304        end
305       
306        % recording the merged field
307        if testima    %in case of input images an image is produced   
308            if isa(MergeData.A,'uint8')
309                bitdepth=8;
310            elseif isa(MergeData.A,'uint16')
311                bitdepth=16;
312            end
313            imwrite(MergeData.A,mergename,'BitDepth',bitdepth);
314            %write xml calibration file
315            siz=size(MergeData.A);
316            npy=siz(1);
317            npx=siz(2);
318            if isfield(MergeData,'VarAttribute')&&isfield(MergeData.VarAttribute{1},'Coord_2')&&isfield(MergeData.VarAttribute{1},'Coord_1')
319                Rangx=MergeData.VarAttribute{1}.Coord_2;
320                Rangy=MergeData.VarAttribute{1}.Coord_1;
321            elseif isfield(MergeData,'AX')&& isfield(MergeData,'AY')
322                Rangx=[MergeData.AX(1) MergeData.AX(end)];
323                Rangy=[MergeData.AY(1) MergeData.AY(end)];
324            else
325                Rangx=[0.5 npx-0.5];
326                Rangy=[npy-0.5 0.5];%default
327            end
328            pxcmx=(npx-1)/(Rangx(2)-Rangx(1));
329            pxcmy=(npy-1)/(Rangy(1)-Rangy(2));
330            T_x=-pxcmx*Rangx(1)+0.5;
331            T_y=-pxcmy*Rangy(2)+0.5;
332            GeometryCal.focal=1;
333            GeometryCal.R=[pxcmx,0,0;0,pxcmy,0;0,0,1];
334            GeometryCal.Tx_Ty_Tz=[T_x T_y 1];
335            ImaDoc.GeometryCalib=GeometryCal;
336            t=struct2xml(ImaDoc);
337            t=set(t,1,'name','ImaDoc');
338            save(t,[filebase_merge '.xml'])     
339            display([filebase_merge '.xml saved'])
340        else
341            MergeData.ListGlobalAttribute={'Project','InputFile_1','InputFile_end','nb_coord','nb_dim','dt','Time','civ'};       
342            MergeData.nb_coord=2;
343            MergeData.nb_dim=2;
344            MergeData.dt=1;
345            MergeData.Time=time_i;
346            error=struct2nc(mergename,MergeData); %save result file
347            if isempty(error)
348                display(['output file ' mergename ' written'])
349            else
350                display(error)
351            end
352        end
353    end
354end
355
356%--------------------------------------------------------------------------   
357function MergeData=merge_field(Data)
358% initiate Matlab  structure for physical field
359if isempty(Data)||~iscell(Data)
360    MergeData=[];
361    return
362end
363MergeData=Data{1};%default
364error=0;
365nbview=length(Data);
366if nbview==1
367    return
368end
369for iview=1:nbview
370    if ~isequal(MergeData.ListDimName,Data{iview}.ListDimName)
371        error=1;
372    end
373    if ~isequal(MergeData.ListVarName,Data{iview}.ListVarName)
374        error=1;
375    end
376%      if ~isequal(MergeData.VarDimIndex,Data{iview}.VarDimIndex)
377%         error=1;
378%      end
379end
380if error
381    MergeData.Txt='ERROR: attempt at merging fields of incompatible type';
382    return
383end
Note: See TracBrowser for help on using the repository browser.