Changeset 905

Timestamp:

May 30, 2015, 8:44:20 PM (10 years ago)

Author:

sommeria

Message:

projection with tps corrected + minor changes

Location:

trunk/src

Files:

• calc_field_tps.m (modified) (1 diff)
• mouse_up.m (modified) (2 diffs)
• proj_field.m (modified) (6 diffs)
• series.m (modified) (8 diffs)
• series/civ_series.m (modified) (1 diff)
• series/merge_proj.m (modified) (4 diffs)
• uvmat.m (modified) (4 diffs)

trunk/src/calc_field_tps.m

@@ -113 +113 @@
     nbvec_sub=NbCentre(isub);
     check_range=(Coord_interp >=ones(nb_sites,1)*SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*SubRange(:,2,isub)');
-    ind_sel=find(sum(check_range,2)==nb_coord);
+    ind_sel=find(sum(check_range,2)==nb_coord);% select points whose all coordinates are in the prescribed range
     nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
     if check_grid

trunk/src/mouse_up.m

@@ -141 +141 @@
         if  ~isempty(ObjectData)
             % plot the field projected on the object
-            ProjData= proj_field(UvData.Field,ObjectData);%project the current interface field on ObjectData
-            if ~isempty(ProjData)
+            [ProjData,errormsg]= proj_field(UvData.Field,ObjectData);%project the current interface field on ObjectData
+            if isempty(errormsg) && ~isempty(ProjData)
                 if strcmp(FigTag,'uvmat')% uvmat plot selected, projection plot seen in view_field
                     hview_field=findobj(allchild(0),'tag','view_field');
@@ -170 +170 @@
                 end
             end
+            if ~isempty(errormsg)
+                msgbox_uvmat('ERROR',errormsg)
+                return
+            end
             set(hhuvmat.CheckViewField,'Value',1);%
             set(hhuvmat.CheckEditObject,'Value',1);%   

trunk/src/proj_field.m

@@ -522 +522 @@
 function  [ProjData,errormsg] = proj_line(FieldData, ObjectData)
 %-----------------------------------------------------------------
+
+%% prepare heading for the projected field
 [ProjData,errormsg]=proj_heading(FieldData,ObjectData);%transfer global attributes
 if ~isempty(errormsg)
@@ -544 +546 @@
 ListIndex={};
 
+%% group the variables (fields of 'FieldData') in cells of variables with the same dimensions
+[CellInfo,NbDim,errormsg]=find_field_cells(FieldData);
+if ~isempty(errormsg)
+    errormsg=['error in proj_field/proj_line:' errormsg];
+    return
+end
+CellInfo=CellInfo(NbDim==2); %keep only the 2D cells
+%%%%%% TODO: treat 1D fields: project as identity so that P o P=P for projection operation
+cell_select=true(size(CellInfo));
+
+for icell=1:length(CellInfo)
+    if isfield(CellInfo{icell},'ProjModeRequest')
+        if ~strcmp(CellInfo{icell}.ProjModeRequest, ProjMode)
+            cell_select(icell)=0;
+        end
+        if strcmp(ProjMode,'interp_tps')&& ~strcmp(CellInfo{icell}.CoordType,'tps')
+            cell_select(icell)=0;
+        end
+    end
+end
+if isempty(find(cell_select))
+    errormsg=[' invalid projection mode ''' ProjMode ''': use ''interp_tps'' to interpolate spatial derivatives'];
+    return
+end
+CellInfo=CellInfo(cell_select);
 
 %% projection line: object types selected from  proj_field='line','polyline','polygon','rectangle','ellipse':
@@ -633 +660 @@
     end
 end
-
-
-%% group the variables (fields of 'FieldData') in cells of variables with the same dimensions
-[CellInfo,NbDim,errormsg]=find_field_cells(FieldData);
-if ~isempty(errormsg)
-    errormsg=['error in proj_field/proj_line:' errormsg];
-    return
-end
-CellInfo=CellInfo(NbDim==2); %keep only the 2D cells
-%%%%%% TODO: treat 1D fields: project as identity so that P o P=P for projection operation
 
 %% loop on variable cells with the same space dimension 2
@@ -1287 +1304 @@
                     Dist=Distx.*Distx+Disty.*Disty;
                     for ivar=1:numel(VarVal)
-                        VarVal{ivar}(Dist>4*ProjData.CoordMesh)=NaN;% put to NaN interpolated positions too far from initial data
+                        VarVal{ivar}(Dist>16*ProjData.CoordMesh)=NaN;% % put to NaN interpolated positions further than 4 meshes from initial data
                     end  
                     
@@ -1325 +1342 @@
                 
                 % set to NaN interpolation points which are too far from any initial data (more than 2 CoordMesh)
-                    if exist('scatteredInterpolant','file')%recent Matlab versions
-                        F=scatteredInterpolant(coord_X, coord_Y,coord_X,'nearest');
-                        G=scatteredInterpolant(coord_X, coord_Y,coord_Y,'nearest');
-                    else
-                        F=TriScatteredInterp([coord_X coord_Y],coord_X,'nearest');
-                        G=TriScatteredInterp([coord_X coord_Y],coord_Y,'nearest');
-                    end
-                    Distx=F(XI,YI)-XI;% diff of x coordinates with the nearest measurement point
-                    Disty=G(XI,YI)-YI;% diff of y coordinates with the nearest measurement point
-                    Dist=Distx.*Distx+Disty.*Disty;
-                    for ivar=1:numel(VarVal)
-                        VarVal{ivar}(Dist>4*ProjData.CoordMesh)=NaN;% put to NaN interpolated positions too far from initial data
-                    end  
-                    
-                
+                Coord=permute(Coord,[1 3 2]);
+                Coord=reshape(Coord,size(Coord,1)*size(Coord,2),2);
+                if exist('scatteredInterpolant','file')%recent Matlab versions
+                    F=scatteredInterpolant(Coord,Coord(:,1),'nearest');
+                    G=scatteredInterpolant(Coord,Coord(:,2),'nearest');
+                else
+                    F=TriScatteredInterp(Coord,Coord(:,1),'nearest');
+                    G=TriScatteredInterp(Coord,Coord(:,2),'nearest');
+                end
+                Distx=F(XI,YI)-XI;% diff of x coordinates with the nearest measurement point
+                Disty=G(XI,YI)-YI;% diff of y coordinates with the nearest measurement point
+                Dist=Distx.*Distx+Disty.*Disty;
                 ListVarName=(fieldnames(DataOut))';
                 VarDimName=cell(size(ListVarName));
@@ -1345 +1359 @@
                     VarName=ListVarName{ilist};
                     ProjData.(VarName)=DataOut.(VarName);
+                    ProjData.(VarName)(Dist>16*ProjData.CoordMesh)=NaN;% put to NaN interpolated positions further than 4 meshes from initial data
                     VarDimName{ilist}={'coord_y','coord_x'};
                 end

trunk/src/series.m

@@ -1607 +1607 @@
     detect=exist(fullfile(Param.InputTable{1,1},SubDirOutNew),'dir');% test if  the dir  already exist
     check_create=1; %need to create the result directory by default
+    CheckOverwrite=1;
+    if isfield(Param,'CheckOverwrite')
+        CheckOverwrite=Param.CheckOverwrite;
+    end
     while detect
-        if Param.CheckOverwrite
+        if CheckOverwrite
             comment=', possibly overwrite previous data';
         else
@@ -1719 +1723 @@
 end
 nbfield_j=numel(ref_j); % number of j indices
+BlockLength=numel(ref_i);%default
 if isempty(Param.IndexRange.NbSlice)
     NbProcess=NbCore;% choose one process per core by default if NbSlice is not imposed
@@ -1805 +1810 @@
     %% processing on a different session of the same computer (background) or cluster, create executable files
     batch_file_list=cell(NbProcess,1);% initiate the list of executable files
-    DirBat=fullfile(OutputDir,'0_EXE');
+    DirExe=fullfile(OutputDir,'0_EXE');%directory name for executable files
     switch computer
         case {'PCWIN','PCWIN64'} %Windows system
@@ -1813 +1818 @@
     end
     %create subdirectory for executable files
-    if ~exist(DirBat,'dir')
-        [tild,msg1]=mkdir(DirBat);
+    if ~exist(DirExe,'dir')
+        [tild,msg1]=mkdir(DirExe);
         if ~strcmp(msg1,'')
-            errormsg=['cannot create ' DirBat ': ' msg1];%error message for directory creation
+            errormsg=['cannot create ' DirExe ': ' msg1];%error message for directory creation
             return
         end
@@ -1852 +1857 @@
         
         %create the executable file
-         filebat=fullfile_uvmat(DirBat,'',Param.InputTable{1,3},ExeExt,OutputNomType,...
+         filebat=fullfile_uvmat(DirExe,'',Param.InputTable{1,3},ExeExt,OutputNomType,...
            Param.IndexRange.first_i,Param.IndexRange.last_i,first_j,last_j);
         batch_file_list{iprocess}=filebat;
@@ -1919 +1924 @@
     case 'cluster_oar' % option 'oar-parexec' used
         %create subdirectory for oar command and log files
-        DirOAR=fullfile(OutputDir,'0_LOG');
-        if exist(DirOAR,'dir')% delete the content of the dir 0_LOG to allow new input
-            curdir=pwd;
-            cd(DirOAR)
-            delete('*')
-            cd(curdir)
-        else
-            [tild,msg1]=mkdir(DirOAR);
-            if ~strcmp(msg1,'')
-                errormsg=['cannot create ' DirOAR ': ' msg1];%error message for directory creation
-                return
-            end
-        end
-        filename_joblist=fullfile(DirOAR,'0_job_list.txt');%create name of the global executable file
+        %DirOARLog=fullfile(OutputDir,'0_LOG');
+        %DirOARExe=fullfile(OutputDir,'0_EXE');
+%         if exist(DirOAR,'dir')% delete the content of the dir 0_LOG to allow new input
+%             curdir=pwd;
+%             cd(DirOAR)
+%             delete('*')
+%             cd(curdir)
+%         else
+%             [tild,msg1]=mkdir(DirOAR);
+%             if ~strcmp(msg1,'')
+%                 errormsg=['cannot create ' DirOAR ': ' msg1];%error message for directory creation
+%                 return
+%             end
+%         end
+        filename_joblist=fullfile(DirExe,'0_job_list.txt');%create name of the global executable file
+        filename_log=fullfile(DirLog,'0_job_list.stdout');%file for output messages of the master oar process
+        filename_errors=fullfile(DirLog,'0_job_list.stderr');%file for error messages of the master oar process
+        
         fid=fopen(filename_joblist,'w');
         for p=1:length(batch_file_list)
@@ -1958 +1967 @@
             '-l /core=' num2str(NbCore) ','...
             'walltime=' datestr(WallTimeOneJob/24,13) ' '...
-            '-E ' regexprep(filename_joblist,'\.txt\>','.stderr') ' '...
-            '-O ' regexprep(filename_joblist,'\.txt\>','.stdout') ' '...
+            '-E ' filename_errors ' '...
+            '-O ' filename_log ' '...
             extra_oar ' '...
             '"oar-parexec -s -f ' filename_joblist ' '...
             '-l ' filename_joblist '.log"\n'];
         
-        filename_oarcommand=fullfile(DirOAR,'0_oar_command');
+        filename_oarcommand=fullfile(DirExe,'0_oar_command');
         fid=fopen(filename_oarcommand,'w');
         fprintf(fid,oar_command);
@@ -3597 +3606 @@
 % --- Executes on button press in CheckOverwrite.
 function CheckOverwrite_Callback(hObject, eventdata, handles)
-% hObject    handle to CheckOverwrite (see GCBO)
-% eventdata  reserved - to be defined in a future version of MATLAB
-% handles    structure with handles and user data (see GUIDATA)
-
-% Hint: get(hObject,'Value') returns toggle state of CheckOverwrite
+
 
 

trunk/src/series/civ_series.m

@@ -333 +333 @@
     if ~isempty(RUNHandle)% update the waitbar in interactive mode with GUI series  (checkrun=1)
         update_waitbar(WaitbarHandle,ifield/NbField)
-        if  ~strcmp(get(RUNHandle,'BusyAction'),'queue')
+        if  checkrun && ~strcmp(get(RUNHandle,'BusyAction'),'queue')
             disp('program stopped by user')
             break

trunk/src/series/merge_proj.m

@@ -73 +73 @@
     ParamOut.OutputFileMode='NbInput';% '=NbInput': 1 output file per input file index, '=NbInput_i': 1 file per input file index i, '=NbSlice': 1 file per slice
       %check the input files
+    ParamOut.CheckOverwriteVisible='on'; % manage the overwrite of existing files (default=1)
     first_j=[];
     if isfield(Param.IndexRange,'first_j'); first_j=Param.IndexRange.first_j; end
@@ -90 +91 @@
 %%%%%%%%%%%% STANDARD PART (DO NOT EDIT) %%%%%%%%%%%%
 ParamOut=[]; %default output
+RUNHandle=[];
+WaitbarHandle=[];
 %% read input parameters from an xml file if input is a file name (batch mode)
 checkrun=1;
@@ -95 +98 @@
     Param=xml2struct(Param);% read Param as input file (batch case)
     checkrun=0;
-end
-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
+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})
@@ -407 +411 @@
     end
 end
-disp(['total ellapsed time ' num2str(toc(tstart))])
+ellapsed_time=toc(tstart);
+disp(['total ellapsed time ' num2str(ellapsed_time/60,2) ' minutes'])
+disp([ num2str(ellapsed_time/(60*NbField),3) ' minutes per iteration'])
 
 %'merge_field': concatene fields

trunk/src/uvmat.m

@@ -3735 +3735 @@
         [ObjectData,errormsg]=proj_field(UvData.Field,UvData.ProjObject{iobj});% project field on the object
         if ~isempty(errormsg)
+            errormsg=['projection on ' UvData.ProjObject{iobj}.Type ': ' errormsg ]; 
             return
         end
@@ -5378 +5379 @@
 
 %% update the projection plot on uvmat
-ProjData= proj_field(UvData.Field,ObjectData);%project the current input field on object ObjectData
+[ProjData,errormsg]= proj_field(UvData.Field,ObjectData);%project the current input field on object ObjectData
+    if ~isempty(errormsg)
+        msgbox_uvmat('ERROR',['projection on ' UvData.ProjObject{iobj}.Type ': ' errormsg ])
+        return
+    end
 plot_field(ProjData,handles.PlotAxes,read_GUI(handles.uvmat));% plot the projected field;
 UvData.PlotAxes=ProjData;% store the plotted field for further update
@@ -5446 +5451 @@
     Data=get(hview_field,'UserData');
     hhview_field=guidata(hview_field);
-    ProjData= proj_field(UvData.Field,ObjectData);%project the current interface field on ObjectData
+    [ProjData,errormsg]= proj_field(UvData.Field,ObjectData);%project the current interface field on ObjectData
+    if ~isempty(errormsg)
+        msgbox_uvmat('ERROR',['projection on ' UvData.ProjObject{iobj}.Type ': ' errormsg ])
+        return
+    end
     [PlotType,PlotParam]=plot_field(ProjData,hhview_field.PlotAxes,read_GUI(hview_field));%read plotting parameters on the uvmat interface
     haxes=findobj(hview_field,'tag','axes3');
@@ -5614 +5623 @@
     
     %% show the projection of the selected object on view_field
-    ProjData= proj_field(UvData.Field,UvData.ProjObject{IndexObj});%project the current field on ObjectData
+    [ProjData,errormsg]= proj_field(UvData.Field,UvData.ProjObject{IndexObj});%project the current field on ObjectData
+    if ~isempty(errormsg)
+        msgbox_uvmat('ERROR',['projection on ' UvData.ProjObject{IndexObj}.Type ': ' errormsg])
+        return
+    end
     hview_field=findobj(allchild(0),'tag','view_field');
     if isempty(hview_field)