Changeset 246 for trunk/src

Timestamp:

Apr 28, 2011, 10:52:31 AM (13 years ago)

Author:

sommeria

Message:

thin plate shell (patch) introduced

Location:

trunk/src

Files:

• calc_field.m (modified) (1 diff)
• civ.m (modified) (4 diffs)
• civ_uvmat.m (modified) (12 diffs)
• msgbox_uvmat.m (modified) (3 diffs)
• pivlab.m (modified) (2 diffs)
• read_civdata.m (modified) (7 diffs)
• series.m (modified) (2 diffs)
• struct2nc.m (modified) (1 diff)
• tps_coeff.m (added)
• tps_eval.m (added)
• tps_eval_dxy.m (added)
• uvmat.m (modified) (8 diffs)

trunk/src/calc_field.m

@@ -60 +60 @@
         nbcoord=2;
     end
-    DataOut=DataIn; %reproduce global attribute
     ListVarName={};
     ValueList={};
     RoleList={};
     units_cell={};
-    for ilist=1:length(FieldName)
-        if ~isempty(FieldName{ilist})
-            [VarName,Value,Role,units]=feval(FieldName{ilist},DataIn);%calculate field with appropriate function named FieldName{ilist}
-            ListVarName=[ListVarName VarName];
-            ValueList=[ValueList Value];
-            RoleList=[RoleList Role];
-            units_cell=[units_cell units];
-        end
-    end
-    %erase previous data (except coordinates)
-    for ivar=nbcoord+1:length(DataOut.ListVarName)
-        VarName=DataOut.ListVarName{ivar};
-        DataOut=rmfield(DataOut,VarName);
-    end
-    DataOut.ListVarName=DataOut.ListVarName(1:nbcoord);
-    if isfield(DataOut,'VarDimName')
-        DataOut.VarDimName=DataOut.VarDimName(1:nbcoord);
-    else
-        errormsg='element .VarDimName missing in input data';
-        return
-    end
-    DataOut.VarAttribute=DataOut.VarAttribute(1:nbcoord);
-    %append new data
-    DataOut.ListVarName=[DataOut.ListVarName ListVarName];
-    for ivar=1:length(ListVarName)
-        DataOut.VarDimName{nbcoord+ivar}=DataOut.VarDimName{1};
-        DataOut.VarAttribute{nbcoord+ivar}.Role=RoleList{ivar};
-        DataOut.VarAttribute{nbcoord+ivar}.units=units_cell{ivar};
-        eval(['DataOut.' ListVarName{ivar} '=ValueList{ivar};'])
-    end
-end
+    % new civ data
+    if isfield(DataIn,'X_SubRange')
+        XMax=max(max(DataIn.X_SubRange));
+        YMax=max(max(DataIn.Y_SubRange));
+        XMin=min(min(DataIn.X_SubRange));
+        YMin=min(min(DataIn.Y_SubRange));
+        Mesh=sqrt((YMax-YMin)*(XMax-XMin)/numel(DataIn.X_tps));%2D
+        xI=XMin:Mesh:XMax;
+        yI=YMin:Mesh:YMax;
+        [XI,YI]=meshgrid(xI,yI);
+        XI=reshape(XI,[],1);
+        YI=reshape(YI,[],1);
+        DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute; %reproduce global attribute
+        for ilist=1:numel(DataOut.ListGlobalAttribute)
+            eval(['DataOut.' DataOut.ListGlobalAttribute{ilist} '=DataIn.' DataIn.ListGlobalAttribute{ilist} ';'])
+        end
+        DataOut.ListVarName={'coord_y','coord_x','FF'};
+        DataOut.VarDimName{1}='coord_y';
+        DataOut.VarDimName{2}='coord_x';
+        DataOut.coord_y=[yI(1) yI(end)];
+        DataOut.coord_x=[xI(1) xI(end)];
+        DataOut.U=zeros(size(XI));
+        DataOut.V=zeros(size(XI));
+        DataOut.vort=zeros(size(XI));
+        DataOut.div=zeros(size(XI));
+        nbval=zeros(size(XI));
+        [NbSubDomain,xx]=size(DataIn.X_SubRange);
+        for isub=1:NbSubDomain
+            nbvec_sub=DataIn.NbSites(isub);         
+                ind_sel=find(XI>=DataIn.X_SubRange(isub,1) & XI<=DataIn.X_SubRange(isub,2) & YI>=DataIn.Y_SubRange(isub,1) & YI<=DataIn.Y_SubRange(isub,2));
+                %rho smoothing parameter
+                epoints = [XI(ind_sel) YI(ind_sel)];% coordinates of interpolation sites
+                ctrs=[DataIn.X_tps(1:nbvec_sub,isub) DataIn.Y_tps(1:nbvec_sub,isub)];%(=initial points) ctrs
+                nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
+%                 PM = [ones(size(epoints,1),1) epoints];
+                switch FieldName{1}
+                    case {'velocity','u','v'}
+                       % DM_eval = DistanceMatrix(epoints,ctrs);%2D matrix of distances between extrapolation points epoints and spline centres (=site points) ctrs
+                       % EM = tps(1,DM_eval);%values of thin plate
+                        EM = tps_eval(epoints,ctrs);
+                    case{'vort','div'}
+                        EMDXY = tps_eval_dxy(epoints,ctrs);%2D matrix of distances between extrapolation points epoints and spline centres (=site points) ctrs
+     
+                end
+                switch FieldName{1}
+                    case 'velocity'
+                        ListFields={'U', 'V'};
+                        VarAttributes{1}.Role='vector_x';
+                        VarAttributes{2}.Role='vector_y';
+                        DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
+                        DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
+                    case 'u'
+                        ListFields={'U'};
+                        VarAttributes{1}.Role='scalar';
+                        DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
+                    case 'v'
+                        ListFields={'V'};
+                        VarAttributes{1}.Role='scalar';
+                        DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
+                    case 'vort'
+                        ListFields={'vort'};
+                        VarAttributes{1}.Role='scalar';
+%                         EMX = [DMXY(:,:,1) PM];
+%                         EMY = [DMXY(:,:,2) PM];
+%                         DataIn.U_tps(nbvec_sub+1,isub)=0;%constant value suppressed by spatial derivatives
+%                         DataIn.V_tps(nbvec_sub+1,isub)=0;%constant value suppressed by spatial derivatives
+%                         DataIn.V_tps(nbvec_sub+2,isub)=0;% X coefficient suppressed for x wise derivatives
+%                         DataIn.U_tps(nbvec_sub+3,isub)=0;% Y coefficient suppressed for x wise derivatives
+                        DataOut.vort(ind_sel)=DataOut.vort(ind_sel)+EMDXY(:,:,2) *DataIn.U_tps(1:nbvec_sub+3,isub)-EMDXY(:,:,1) *DataIn.V_tps(1:nbvec_sub+3,isub);
+                    case 'div'
+                        ListFields={'div'};
+                        VarAttributes{1}.Role='scalar';
+%                         EMX = [DMXY(:,:,1) PM];
+%                         EMY = [DMXY(:,:,2) PM];
+%                         DataIn.U_tps(nbvec_sub+1,isub)=0;%constant value suppressed by spatial derivatives
+%                         DataIn.V_tps(nbvec_sub+1,isub)=0;%constant value suppressed by spatial derivatives
+%                         DataIn.V_tps(nbvec_sub+2,isub)=0;% X coefficient suppressed for x wise derivatives
+%                         DataIn.U_tps(nbvec_sub+3,isub)=0;% Y coefficient suppressed for x wise derivatives
+                        DataOut.div(ind_sel)=DataOut.div(ind_sel)+EMDXY(:,:,1)*DataIn.U_tps(1:nbvec_sub+3,isub)+EMDXY(:,:,2) *DataIn.V_tps(1:nbvec_sub+3,isub);
+                end
+        end
+        DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang  
+        
+        DataOut.FF=reshape(DataOut.FF,numel(yI),numel(xI));
+        nbval(nbval==0)=1;
+        DataOut.U=DataOut.U ./nbval;
+        DataOut.V=DataOut.V ./nbval;
+        DataOut.U=reshape(DataOut.U,numel(yI),numel(xI));
+        DataOut.V=reshape(DataOut.V,numel(yI),numel(xI));
+        DataOut.vort=reshape(DataOut.vort,numel(yI),numel(xI));
+        DataOut.div=reshape(DataOut.div,numel(yI),numel(xI));
+        DataOut.ListVarName=[DataOut.ListVarName ListFields];
+        for ilist=3:numel(DataOut.ListVarName)
+            DataOut.VarDimName{ilist}={'coord_y','coord_x'};
+        end
+        DataOut.VarAttribute={[],[]};
+        DataOut.VarAttribute{3}.Role='errorflag';
+        DataOut.VarAttribute=[DataOut.VarAttribute VarAttributes];      
+    else   
+        DataOut=DataIn;
+        for ilist=1:length(FieldName)
+            if ~isempty(FieldName{ilist})
+                [VarName,Value,Role,units]=feval(FieldName{ilist},DataIn);%calculate field with appropriate function named FieldName{ilist}
+                ListVarName=[ListVarName VarName];
+                ValueList=[ValueList Value];
+                RoleList=[RoleList Role];
+                units_cell=[units_cell units];
+            end
+        end
+        %erase previous data (except coordinates)
+        for ivar=nbcoord+1:length(DataOut.ListVarName)
+            VarName=DataOut.ListVarName{ivar};
+            DataOut=rmfield(DataOut,VarName);
+        end
+        DataOut.ListVarName=DataOut.ListVarName(1:nbcoord);
+        if isfield(DataOut,'VarDimName')
+            DataOut.VarDimName=DataOut.VarDimName(1:nbcoord);
+        else
+            errormsg='element .VarDimName missing in input data';
+            return
+        end
+        DataOut.VarAttribute=DataOut.VarAttribute(1:nbcoord);
+        %append new data
+        DataOut.ListVarName=[DataOut.ListVarName ListVarName];
+        for ivar=1:length(ListVarName)
+            DataOut.VarDimName{nbcoord+ivar}=DataOut.VarDimName{1};
+            DataOut.VarAttribute{nbcoord+ivar}.Role=RoleList{ivar};
+            DataOut.VarAttribute{nbcoord+ivar}.units=units_cell{ivar};
+            eval(['DataOut.' ListVarName{ivar} '=ValueList{ivar};'])
+        end
+    end
+end
+
 
 %%%%%%%%%%%%% velocity fieldn%%%%%%%%%%%%%%%%%%%%

trunk/src/civ.m

@@ -62 +62 @@
 % Update handles structure
 guidata(hObject, handles);
-
+set(hObject,'WindowButtonUpFcn',{'mouse_up_GUI',handles}) %set mouse action (zoom on uicontrols)
 %default initial parameters
 filebase=''; % root file name ('filebase'.civ)
@@ -334 +334 @@
     browse.nom_type_nc=nom_type;
     ind_opening=2;% propose 'fix' as the default option
-    Data=nc2struct(fileinput,[]);
-    if isfield(Data,'absolut_time_T0')%test for civx files
+    Data=nc2struct(fileinput,'ListGlobalAttribute','CivStage','absolut_time_T0','fix','patch','civ2','fix2');
+    if ~isempty(Data.CivStage)%test for civ files
+        ind_opening=Data.CivStage;
+        set(handles.CivAll,'Value',3)
+    end
+    if ~isempty(Data.absolut_time_T0)%test for civx files
+        set(handles.CivAll,'Value',1)
         if isfield(Data,'fix') && isequal(Data.fix,1)
             ind_opening=3;
         end
-        if isfield(Data,'patch') && isequal(Data.patch,1)
+        if isequal(Data.patch,1)
             ind_opening=4;
         end
-        if isfield(Data,'civ2') && isequal(Data.civ2,1)
+        if isequal(Data.civ2,1)
             ind_opening=5;
         end
-        if isfield(Data,'fix2') && isequal(Data.fix2,1)
+        if  isequal(Data.fix2,1)
             ind_opening=6;
         end
@@ -4683 +4688 @@
     par_civ1.T0=0;
     par_civ1.Dt=1;
-    Data=civ_uvmat(par_civ1);
+    Param.Civ1=par_civ1;
+    Data=civ_uvmat(Param);
 %     stepx=str2num(par_civ1.dx);
 %     stepy=str2num(par_civ1.dy);
@@ -4877 +4883 @@
 function open_view_field(hObject, eventdata)
 %-------------------------------------------------------------------
-     list=get(gcbo,'String');
-     index=get(gcbo,'Value');
-     rootroot=get(gcbo,'UserData');
+     list=get(hObject,'String');
+     index=get(hObject,'Value');
+     rootroot=get(hObject,'UserData');
      filename=list{index};
      ind_dot=findstr(filename,'...');
      filename=filename(1:ind_dot-1);
       filename=fullfile(rootroot,filename);
+      delete(get(hObject,'parent'))%delete the display figure to stop the check process
       if exist(filename,'file')%visualise the vel field if it exists
-        %[Field,VelTypeOut]=read_civxdata(filename);
-        %view_field(Field)
         uvmat(filename)
         set(gcbo,'Value',1)

trunk/src/civ_uvmat.m

@@ -1 +1 @@
 % To develop....
 function [Data,errormsg]= civ_uvmat(Param,ncfile)
+errormsg='';
 Data.ListGlobalAttribute={'Conventions','Program','CivStage'};
 Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes
 Data.Program='civ_uvmat';
 Data.CivStage=0;%default
+ListVarCiv1={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F'};
+ListVarFix1={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F','Civ1_FF'};
 
 %% Civ1
@@ -20 +23 @@
     shiftx=str2num(par_civ1.shiftx);
     shifty=str2num(par_civ1.shifty);
-    miniy=max(1+isy2-shifty,1+iby2);
+    miniy=max(1+isy2+shifty,1+iby2);
     minix=max(1+isx2-shiftx,1+ibx2);
-    maxiy=min(size(image1,1)-isy2-shifty,size(image1,1)-iby2);
+    maxiy=min(size(image1,1)-isy2+shifty,size(image1,1)-iby2);
     maxix=min(size(image1,2)-isx2-shiftx,size(image1,2)-ibx2);
     [GridX,GridY]=meshgrid(minix:stepx:maxix,miniy:stepy:maxiy);
@@ -29 +32 @@
     
     % caluclate velocity data (y and v in indices, reverse to y component)
-    [xtable ytable utable vtable ctable F] = pivlab (image1,image2,ibx2,iby2,isx2,isy2,shiftx,shifty,PointCoord,str2num(par_civ1.rho), []);
+    [xtable ytable utable vtable ctable F] = pivlab (image1,image2,ibx2,iby2,isx2,isy2,shiftx,-shifty,PointCoord,str2num(par_civ1.rho), []);
     list_param=(fieldnames(par_civ1))';
     list_remove={'pxcmx','pxcmy','npx','npy','gridflag','maskflag','term_a','term_b','T0'};
@@ -67 +70 @@
     Data.CivStage=1;
 else
-    Data=nc2struct(ncfile)%read existing netcdf file
+    if isfield(Param,'Fix1')
+        Data=nc2struct(ncfile,ListVarCiv1);%read civ1 data in the existing netcdf file
+    else
+        Data=nc2struct(ncfile,ListVarFix1);%read civ1 and fix1 data in the existing netcdf file
+    end
+    if isfield(Data,'Txt')
+        msgbox_uvmat('ERROR',Data.Txt)
+        return
+    end
+    % read Civx data
     if isfield(Data,'absolut_time_T0')%read civx file
         var={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F','Civ1_FF';'vec_X','vec_Y','vec_U','vec_V','vec_C','vec_F','vec_FixFlag'};
@@ -95 +107 @@
         ParamName=ListFixParam{ilist};
         ListName=['Fix1_' ParamName];
-        ['Data.ListGlobalAttribute=[Data.ListGlobalAttribute ''' ParamName '''];']
         eval(['Data.ListGlobalAttribute=[Data.ListGlobalAttribute ''' ParamName '''];'])
         eval(['Data.' ListName '=Param.Fix1.' ParamName ';'])
@@ -107 +118 @@
     Data.VarDimName=[Data.VarDimName {'nbvec1'}];
     nbvar=length(Data.ListVarName);
-    Data.VarAttribute{nbvar}.Role='errorflag';
-    [Data.Civ1_FF]=fix_uvmat(Param.Fix1,Data.Civ1_F,Data.Civ1_C,Data.Civ1_U,Data.Civ1_V);
+    Data.VarAttribute{nbvar}.Role='errorflag';    
+    Data.Civ1_FF=fix_uvmat(Param.Fix1,Data.Civ1_F,Data.Civ1_C,Data.Civ1_U,Data.Civ1_V);
     Data.CivStage=2;                               
 end   
@@ -117 +128 @@
     Data.Patch1_Threshold=str2double(Param.Patch1.Threshold);
     Data.Patch1_SubDomain=str2double(Param.Patch1.SubDomain);
-    Data.ListVarName=[Data.ListVarName {'Patch1_U','Patch1_V'}];
-    Data.VarDimName=[Data.VarDimName {'nbvec1','nbvec1'}];
+    Data.ListVarName=[Data.ListVarName {'Civ1_U_Diff','Civ1_V_Diff','Civ1_X_SubRange','Civ1_Y_SubRange','Civ1_NbSites','Civ1_X_tps','Civ1_Y_tps','Civ1_U_tps','Civ1_V_tps'}];
+    Data.VarDimName=[Data.VarDimName {'NbVec1','NbVec1',{'NbSubDomain1','Two'},{'NbSubDomain1','Two'},'NbSubDomain1',...
+             {'NbVec1Sub','NbSubDomain1'},{'NbVec1Sub','NbSubDomain1'},{'Nbtps1','NbSubDomain1'},{'Nbtps1','NbSubDomain1'}}];
     nbvar=length(Data.ListVarName);
     Data.VarAttribute{nbvar-1}.Role='vector_x';
     Data.VarAttribute{nbvar}.Role='vector_y';
-    Data.Patch1_U=zeros(size(Data.Civ1_X));
-    Data.Patch1_V=zeros(size(Data.Civ1_X));
+    Data.Civ1_U_Diff=zeros(size(Data.Civ1_X));
+    Data.Civ1_V_Diff=zeros(size(Data.Civ1_X));
     if isfield(Data,'Civ1_FF')
         ind_good=find(Data.Civ1_FF==0);
@@ -129 +141 @@
         ind_good=1:numel(Data.Civ1_X);
     end
-    Data.Civ1_X
-    [Ures, Vres,FFres]=...
+    [Data.Civ1_X_SubRange,Data.Civ1_Y_SubRange,Data.Civ1_NbSites,FFres,Ures, Vres,Data.Civ1_X_tps,Data.Civ1_Y_tps,Data.Civ1_U_tps,Data.Civ1_V_tps]=...
                             patch_uvmat(Data.Civ1_X(ind_good)',Data.Civ1_Y(ind_good)',Data.Civ1_U(ind_good)',Data.Civ1_V(ind_good)',Data.Patch1_Rho,Data.Patch1_Threshold,Data.Patch1_SubDomain); 
-                        size(Ures)
-                        size(Vres)
-                        size(FFres)
-                        size(ind_good)
-      Data.Patch1_U(ind_good)=Ures;
-      Data.Patch1_V(ind_good)=Vres;
-      Data.Civ1_FF(ind_good)=FFres
-      Data.CivStage=3;                               
+      Data.Civ1_U_Diff(ind_good)=Data.Civ1_U(ind_good)-Ures;
+      Data.Civ1_V_Diff(ind_good)=Data.Civ1_V(ind_good)-Vres;
+      Data.Civ1_FF(ind_good)=FFres;
+      Data.CivStage=3;                             
 end   
 %% write result
+% 'TESTcalc'
+% [DataOut,errormsg]=calc_field('velocity',Data)
+if exist('ncfile','var')
 errormsg=struct2nc(ncfile,Data);
-
+end
 
 
@@ -192 +202 @@
     FF=FF==1 | (U.*U+V.*V)>thresh;
 end
- 
+FF=double(FF);
 % 
 % % criterium on velocity values
@@ -270 +280 @@
 %------------------------------------------------------------------------
 % patch function
-function [U_patch,V_patch,FF,SubRangx,SubRangy,X_ctrs,Y_ctrs] =patch_uvmat(X,Y,U,V,Rho,Threshold,SubDomain)
+function [SubRangx,SubRangy,nbpoints,FF,U_smooth,V_smooth,X_tps,Y_tps,U_tps,V_tps] =patch_uvmat(X,Y,U,V,Rho,Threshold,SubDomain)
 %subdomain decomposition
 warning off
@@ -278 +288 @@
 Y=reshape(Y,[],1);
 nbvec=numel(X);
-NbSubDomain=ceil(nbvec/SubDomain)
-MinX=min(X)
-MinY=min(Y)
-MaxX=max(X)
-MaxY=max(Y)
+NbSubDomain=ceil(nbvec/SubDomain);
+MinX=min(X);
+MinY=min(Y);
+MaxX=max(X);
+MaxY=max(Y);
 RangX=MaxX-MinX;
 RangY=MaxY-MinY;
-AspectRatio=RangY/RangX
-NbSubDomainX=ceil(sqrt(NbSubDomain/AspectRatio))
-NbSubDomainY=ceil(sqrt(NbSubDomain*AspectRatio))
-
+AspectRatio=RangY/RangX;
+NbSubDomainX=max(floor(sqrt(NbSubDomain/AspectRatio)),1);
+NbSubDomainY=max(floor(sqrt(NbSubDomain*AspectRatio)),1);
+NbSubDomain=NbSubDomainX*NbSubDomainY;
 SizX=RangX/NbSubDomainX;%width of subdomains
 SizY=RangY/NbSubDomainY;%height of subdomains
 CentreX=linspace(MinX+SizX/2,MaxX-SizX/2,NbSubDomainX);
 CentreY=linspace(MinY+SizY/2,MaxY-SizY/2,NbSubDomainY);
-SubIndexX=ceil((X-MinX)/SizX);%subdomain index of vectors
-SubIndexY=ceil((Y-MinY)/SizY);
-rho=RangX*RangY*Rho;%optimum rho increase as teh area of the subdomain
-U_tps=zeros(length(X),NbSubDomainY,NbSubDomainX);%default spline
-V_tps=zeros(length(X),NbSubDomainY,NbSubDomainX);%default spline
-X_dist=zeros(length(X),NbSubDomainY,NbSubDomainX);%default spline
-Y_dist=zeros(length(X),NbSubDomainY,NbSubDomainX);%default spline
-U_smooth=zeros(length(X),NbSubDomainY,NbSubDomainX);
-V_smooth=zeros(length(X),NbSubDomainY,NbSubDomainX);
-dist_ctre=zeros(length(X),NbSubDomainY,NbSubDomainX);
+[CentreX,CentreY]=meshgrid(CentreX,CentreY);
+CentreY=reshape(CentreY,1,[]);
+CentreX=reshape(CentreX,1,[]);
+rho=SizX*SizY*Rho/1000000;%optimum rho increase as the area of the subdomain (division by 10^6 to reach good values with the default GUI input)
+U_tps_sub=zeros(length(X),NbSubDomain);%default spline
+V_tps_sub=zeros(length(X),NbSubDomain);%default spline
+U_smooth=zeros(length(X),1);
+V_smooth=zeros(length(X),1);
+
+nb_select=zeros(length(X),1);
 FF=zeros(length(X),1);
-for isubx=1:NbSubDomainX
-    for isuby=1:NbSubDomainY
-        SubRangx(isuby,isubx,:)=[CentreX(isubx)-SizX/2 CentreX(isubx)+SizX/2];
-        SubRangy(isuby,isubx,:)=[CentreY(isubx)-SizY/2 CentreY(isubx)+SizY/2];
-        for iter=1:3 %increase the subdomain during three iterations at most
-            ind_sel=find(X>SubRangx(isuby,isubx,1) & X<SubRangx(isuby,isubx,2) & Y>SubRangy(isuby,isubx,1) & Y<SubRangy(isuby,isubx,2));  
-            size(ind_sel)
-            if numel(ind_sel)<SubDomain/4;% too few selected vectors, increase the subrange for next iteration
-                SubRangx(isuby,isubx,1)=SubRangx(isuby,isubx,1)-SizX/4;
-                SubRangx(isuby,isubx,2)=SubRangx(isuby,isubx,2)+SizX/4;
-                SubRangy(isuby,isubx,1)=SubRangy(isuby,isubx,1)-SizY/4;
-                SubRangy(isuby,isubx,2)=SubRangy(isuby,isubx,2)+SizY/4;
+test_empty=zeros(1,NbSubDomain);
+for isub=1:NbSubDomain
+    SubRangx(isub,:)=[CentreX(isub)-SizX/2 CentreX(isub)+SizX/2];
+    SubRangy(isub,:)=[CentreY(isub)-SizY/2 CentreY(isub)+SizY/2];
+    ind_sel_previous=[];
+    ind_sel=0;
+    while numel(ind_sel)>numel(ind_sel_previous) %increase the subdomain during four iterations at most
+        ind_sel_previous=ind_sel;
+        ind_sel=find(X>SubRangx(isub,1) & X<SubRangx(isub,2) & Y>SubRangy(isub,1) & Y<SubRangy(isub,2));
+        % if no vector in the subdomain, skip the subdomain
+        if isempty(ind_sel)
+            test_empty(isub)=1;    
+            U_tps(1,isub)=0;%define U_tps and V_tps by default
+            V_tps(1,isub)=0;
+            break
+            % if too few selected vectors, increase the subrange for next iteration
+        elseif numel(ind_sel)<SubDomain/4 && ~isequal( ind_sel,ind_sel_previous);
+            SubRangx(isub,1)=SubRangx(isub,1)-SizX/4;
+            SubRangx(isub,2)=SubRangx(isub,2)+SizX/4;
+            SubRangy(isub,1)=SubRangy(isub,1)-SizY/4;
+            SubRangy(isub,2)=SubRangy(isub,2)+SizY/4;
+        else
+            [U_smooth_sub,U_tps_sub]=tps_coeff(X(ind_sel),Y(ind_sel),U(ind_sel),rho);
+            [V_smooth_sub,V_tps_sub]=tps_coeff(X(ind_sel),Y(ind_sel),V(ind_sel),rho);
+            UDiff=U_smooth_sub-U(ind_sel);
+            VDiff=V_smooth_sub-V(ind_sel);
+            NormDiff=UDiff.*UDiff+VDiff.*VDiff;
+            FF(ind_sel)=20*(NormDiff>Threshold);%put FF value to 20 to identify the criterium of elimmination
+            ind_ind_sel=find(FF(ind_sel)==0); % select the indices of ind_sel corresponding to the remaining vectors
+            % no value exceeds threshold, the result is recorded
+            if isequal(numel(ind_ind_sel),numel(ind_sel))
+                U_smooth(ind_sel)=U_smooth(ind_sel)+U_smooth_sub;
+                V_smooth(ind_sel)=V_smooth(ind_sel)+V_smooth_sub;
+                nbpoints(isub)=numel(ind_sel);
+                X_tps(1:nbpoints(isub),isub)=X(ind_sel);
+                Y_tps(1:nbpoints(isub),isub)=Y(ind_sel);
+                U_tps(1:nbpoints(isub)+3,isub)=U_tps_sub;
+                V_tps(1:nbpoints(isub)+3,isub)=V_tps_sub;         
+                nb_select(ind_sel)=nb_select(ind_sel)+1;
+                 display('good')
+                break
+                % too few selected vectors, increase the subrange for next iteration
+            elseif numel(ind_ind_sel)<SubDomain/4 && ~isequal( ind_sel,ind_sel_previous);
+                SubRangx(isub,1)=SubRangx(isub,1)-SizX/4;
+                SubRangx(isub,2)=SubRangx(isub,2)+SizX/4;
+                SubRangy(isub,1)=SubRangy(isub,1)-SizY/4;
+                SubRangy(isub,2)=SubRangy(isub,2)+SizY/4;
+%                 display('fewsmooth')
+                % interpolation-smoothing is done again with the selected vectors
             else
-                [U_smooth_sub,U_tps_sub]=tps_uvmat(X(ind_sel),Y(ind_sel),U(ind_sel),rho);
-                [V_smooth_sub,V_tps_sub]=tps_uvmat(X(ind_sel),Y(ind_sel),V(ind_sel),rho);
-                size(U_smooth_sub)
-                size(U(ind_sel))
-                UDiff=U_smooth_sub-U(ind_sel);
-                VDiff=V_smooth_sub-V(ind_sel);
-                NormDiff=UDiff.*UDiff+VDiff.*VDiff;
-                FF(ind_sel)=20*(NormDiff>Threshold);%put FF value to 20 to identify the criterium of elimmination
-                ind_ind_sel=find(FF(ind_sel)==0); 
-                if isequal(numel(ind_ind_sel),numel(ind_sel))
-                    U_smooth(ind_sel,isuby,isubx)=U_smooth_sub;
-                    V_smooth(ind_sel,isuby,isubx)=V_smooth_sub;
-                    break 
-                elseif numel(ind_ind_sel)<SubDomain/4;% too few selected vectors, increase the subrange for next iteration
-                    SubRangx(isuby,isubx,1)=SubRangx(isuby,isubx,1)-SizX/4;
-                    SubRangx(isuby,isubx,2)=SubRangx(isuby,isubx,2)+SizX/4;
-                    SubRangy(isuby,isubx,1)=SubRangy(isuby,isubx,1)-SizY/4;
-                    SubRangy(isuby,isubx,2)=SubRangy(isuby,isubx,2)+SizY/4;
-                else
-                    [U_smooth(ind_sel(ind_ind_sel),isuby,isubx),U_tps(ind_sel(ind_ind_sel),isuby,isubx),U_tps3(isuby,isubx)]=tps_uvmat(X(ind_sel(ind_ind_sel)),Y(ind_sel(ind_ind_sel)),U(ind_sel(ind_ind_sel)),rho);
-                    [V_smooth(ind_sel(ind_ind_sel),isuby,isubx),V_tps(ind_sel(ind_ind_sel),isuby,isubx),V_tps3(isuby,isubx)]=tps_uvmat(X(ind_sel(ind_ind_sel)),Y(ind_sel(ind_ind_sel)),V(ind_sel(ind_ind_sel)),rho);
-                     break                 
-                end
-            end        
-        end       
-        X_ctrs(isuby,isubx)=mean(X(ind_sel));%gravity centre of selected points
-        Y_ctrs(isuby,isubx)=mean(Y(ind_sel));%positions of tps sources for the subdomain i,      
-        dist_ctre(ind_sel,isuby,isubx)=sqrt(abs(((X(ind_sel)-X_ctrs(isuby,isubx)).*(X(ind_sel)-X_ctrs(isuby,isubx)))+((Y(ind_sel)-Y_ctrs(isuby,isubx)).*(Y(ind_sel)-Y_ctrs(isuby,isubx)))));
-    end
-end
-U_patch=sum(sum(U_smooth.*dist_ctre,3),2)./sum(sum(dist_ctre,3),2);
-V_patch=sum(sum(V_smooth.*dist_ctre,3),2)./sum(sum(dist_ctre,3),2);
-
-%------------------------------------------------------------------------
-%fasshauer@iit.edu MATH 590 ? Chapter 19 32
-% X,Y initial coordiantes
-% XI vector, YI column vector for the grid of interpolation points
-function [U_smooth,U_tps,U_tps3]=tps_uvmat(X,Y,U,rho)
-%------------------------------------------------------------------------
-%rho smoothing parameter
-ep = 1; 
-X=reshape(X,[],1);
-Y=reshape(Y,[],1);
-rhs = reshape(U,[],1);
-% if exist('FF','var')
-% test_false=isnan(rhs)|FF~=0;
-% else
-%     test_false=isnan(rhs);
-% end
-% X(test_false)=[];
-% Y(test_false)=[];
-% rhs(test_false)=[];
-%randn('state',3); rhs = rhs + 0.03*randn(size(rhs));
-rhs = [rhs; zeros(3,1)];
-dsites = [X Y];% coordinates of measurement sites
-ctrs = dsites;%radial base functions are located at the measurement sites
-DM_data = DistanceMatrix(dsites,ctrs);%2D matrix of distances between spline centres (=initial points) ctrs
-% if size(XI,1)==1 && size(YI,2)==1 % XI vector, YI column vector
-%      [XI,YI]=meshgrid(XI,YI);
-% end
-% [npy,npx]=size(XI);
-% epoints = [reshape(XI,[],1) reshape(YI,[],1)];
-IM_sites = tps(ep,DM_data);%values of thin plate at site points
-IM = IM_sites + rho*eye(size(IM_sites));%  rho=1/(2*omega) , omega given by fasshauer;
-PM=[ones(size(dsites,1),1) dsites];
-IM=[IM PM; [PM' zeros(3,3)]];
-%fprintf('Condition number estimate: %e\n',condest(IM))
-%DM_eval = DistanceMatrix(epoints,ctrs);%2D matrix of distances between extrapolation points epoints and spline centres (=site points) ctrs
-%EM = tps(ep,DM_eval);%values of thin plate 
-%PM = [ones(size(epoints,1),1) epoints]; 
-%EM = [EM PM];
-U_tps=(IM\rhs);
-PM = [ones(size(dsites,1),1) dsites]; 
-EM = [IM_sites PM];
-U_smooth=EM *U_tps;
-U_tps3=U_tps(end-2:end);
-U_tps=U_tps(1:end-3);
+                [U_smooth_sub,U_tps_sub]=tps_coeff(X(ind_sel(ind_ind_sel)),Y(ind_sel(ind_ind_sel)),U(ind_sel(ind_ind_sel)),rho);
+                [V_smooth_sub,V_tps_sub]=tps_coeff(X(ind_sel(ind_ind_sel)),Y(ind_sel(ind_ind_sel)),V(ind_sel(ind_ind_sel)),rho);
+                U_smooth(ind_sel(ind_ind_sel))=U_smooth(ind_sel(ind_ind_sel))+U_smooth_sub;
+                V_smooth(ind_sel(ind_ind_sel))=V_smooth(ind_sel(ind_ind_sel))+V_smooth_sub;
+                nbpoints(isub)=numel(ind_ind_sel);
+                X_tps(1:nbpoints(isub),isub)=X(ind_sel(ind_ind_sel));
+                Y_tps(1:nbpoints(isub),isub)=Y(ind_sel(ind_ind_sel));
+                U_tps(1:nbpoints(isub)+3,isub)=U_tps_sub;
+                V_tps(1:nbpoints(isub)+3,isub)=V_tps_sub;
+                nb_select(ind_sel(ind_ind_sel))=nb_select(ind_sel(ind_ind_sel))+1;
+                display('good2')
+                break
+            end
+        end
+    end
+end
+ind_empty=find(test_empty);
+%remove empty subdomains
+if ~isempty(ind_empty)
+    SubRangx(ind_empty,:)=[];
+    SubRangy(ind_empty,:)=[];
+    X_tps(:,ind_empty)=[];
+    Y_tps(:,ind_empty)=[];
+    U_tps(:,ind_empty)=[];
+    V_tps(:,ind_empty)=[];
+end
+nb_select(nb_select==0)=1;%ones(size(find(nb_select==0)));
+U_smooth=U_smooth./nb_select;
+V_smooth=V_smooth./nb_select;
+
+
+
+
 
 
@@ -421 +424 @@
 
 
-%   DM:     MxN matrix whose i,j position contains the Euclidean
-%              distance between the i-th data site and j-th center
-  function DM = DistanceMatrix(dsites,ctrs)
-  [M,s] = size(dsites); [N,s] = size(ctrs);
-  DM = zeros(M,N);
-  % Accumulate sum of squares of coordinate differences
-  % The ndgrid command produces two MxN matrices:
-  %   dr, consisting of N identical columns (each containing
-  %       the d-th coordinate of the M data sites)
-  %   cc, consisting of M identical rows (each containing
-  %       the d-th coordinate of the N centers)
-  for d=1:s
-     [dr,cc] = ndgrid(dsites(:,d),ctrs(:,d));
-     DM = DM + (dr-cc).^2;
-  end
-  DM = sqrt(DM);
-
-
-  % rbf = tps(e,r)
-% Defines thin plate spline RBF
-function rbf = tps(e,r) 
-rbf = zeros(size(r));
-nz = find(r~=0);   % to deal with singularity at origin
-rbf(nz) = (e*r(nz)).^2.*log(e*r(nz));
-
+
+

trunk/src/msgbox_uvmat.m

@@ -63 +63 @@
 guidata(hObject, handles);
 testNo=0;
-testCancel=1;
+testCancel=0;
 testinputstring=0;
 icontype='quest';%default question icon (text input asked)
@@ -71 +71 @@
         case {'CONFIRMATION'}
             icontype='';
-            testCancel=0; %no cancel button
         case 'ERROR'
             icontype='error';
-            testCancel=0; %no cancel button
         case 'WARNING'
             icontype='warn';
-            testCancel=0; %no cancel button
         case 'INPUT_Y-N'
             icontype='quest';
+            testCancel=1; %no cancel button
             testNo=1; % button No activated
         case {'RULER'}
             icontype='';
-            testCancel=0; %no cancel button
             testinputstring=1;
         case 'INPUT_TXT'
             testinputstring=1;
+            testCancel=1; %no cancel button
         otherwise
           %  testinputstring=1;
@@ -189 +187 @@
 
 % Get default command line output from handles structure
-if isequal(handles.output,'Cancel')
-    varargout{1}='Cancel';
-elseif isequal(handles.output,'No')
-    varargout{1}='No';
-else
-    varargout{1}=get(handles.edit_box,'String');
-    if isempty(varargout{1}) || isequal(varargout{1},'')
-        varargout{1}='Yes';
-    end
-end
-% The figure can be deleted now
-delete(handles.figure1);
-
+if isfield(handles,'output')
+    if isequal(handles.output,'Cancel')
+        varargout{1}='Cancel';
+    elseif isequal(handles.output,'No')
+        varargout{1}='No';
+    else
+        varargout{1}=get(handles.edit_box,'String');
+        if isempty(varargout{1}) || isequal(varargout{1},'')
+            varargout{1}='Yes';
+        end
+    end
+    % The figure can be deleted now
+end
+ delete(handles.figure1);
+ 
 % --- Executes on button press in OK.
 function OK_Callback(hObject, eventdata, handles)

trunk/src/pivlab.m

@@ -14 +14 @@
 % image1:first image (matrix)
 % image2: second image (matrix)
-% ibx,iby: size of the correlation box along x and y (in px)
+% ibx2,iby2: half size of the correlation box along x and y, in px (size=(2*iby2+1,2*ibx2+1)
+% isx2,isy2: half size of the search box along x and y, in px (size=(2*isy2+1,2*isx2+1)
+% shiftx, shifty: shift of the search box (in pixel index, yshift reversed)
 % step: mesh of the measurement points (in px)
 % subpixfinder=1 or 2 controls the curve fitting of the image correlation
@@ -119 +121 @@
     xtable(ivec)=iref+vector(1)/2;% convec flow (velocity taken at the point middle from imgae1 and 2)
     ytable(ivec)=jref+vector(2)/2;
-    utable(ivec)=vector(1);
-    vtable(ivec)=vector(2);
+    utable(ivec)=vector(1)+shiftx;
+    vtable(ivec)=vector(2)+shifty;
 end
 result_conv=result_conv*corrmax/(255*sum_square);% keep the last correlation matrix for output

trunk/src/read_civdata.m

@@ -47 +47 @@
 % 'nc2struct': reads a netcdf file 
 
-function [Field,VelTypeOut,errormsg]=read_civdata(filename,FieldNames,VelType)
+function [Field,VelTypeOut,errormsg]=read_civdata(filename,FieldNames,VelType,XI,YI)
 errormsg='';
 VelTypeOut=VelType;%default
 %% default input
 if ~exist('VelType','var')
-    VelType=[];
+    VelType='';
 end
 if isequal(VelType,'*')
-    VelType=[];
+    VelType='';
+end
+if isempty(VelType)
+    VelType='';
 end
 if ~exist('FieldNames','var') 
@@ -62 +65 @@
 
 %% reading data
-%[var,role,units,vel_type_out_cell]=varcivx_generator(FieldNames,VelType);%determine the names of constants and variables to read
-% if isempty(VelType)
-%     VelType='civ1';
-% end
-% if isequal(VelType,'civ1')
-%     varlist=[{'X','Y','U','V','C','F','FF'};... %output names
-%         {'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F','Civ1_FF'}];%names in file
-%     role={'coord_x','coord_y','vector_x','vector_y','ancillary','warnflag','errorflag'};
-%     units={'pixel','pixel','pixel','pixel','','',''};
-% elseif isequal(VelType,'interp1')|| isequal(VelType,'filter1')% read diff between spline and initial data (change name interp1 -> splinediff
-%     varlist=[{'X','Y','U','V'};... %output names
-%         {'Civ1_X','Civ1_Y','Patch1_U','Patch1_V'}];
-%     role = {'coord_x','coord_y','vector_x','vector_y'};
-%     units={'pixel','pixel','pixel','pixel'};
-% end
-[varlist,role,units,vel_type_out_cell]=varcivx_generator(FieldNames,VelType); 
+[varlist,role,units,vel_type_out_cell]=varcivx_generator(FieldNames,VelType);
 [Field,vardetect,ichoice]=nc2struct(filename,varlist);%read the variables in the netcdf file
 if isfield(Field,'Txt')
@@ -93 +81 @@
     Field.VarAttribute{ivar}.Mesh=0.1;%typical mesh for histograms O.1 pixel
 end
-
 if ~isempty(ichoice)
     VelTypeOut=vel_type_out_cell{ichoice};
 end
-%  Field.CivStage=3;%for patch, to generalise
-% return
+
 
 %% renaming for standard conventions
@@ -162 +148 @@
 
 
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % TAKEN FROM read_civxdata NOT USED
@@ -178 +165 @@
 
 %% default input values
-if ~exist('vel_type','var'),vel_type=[];end;
+if ~exist('vel_type','var'),vel_type='';end;
 if iscell(vel_type),vel_type=vel_type{1}; end;%transform cell to string if needed
 if ~exist('FieldNames','var'),FieldNames={'ima_cor'};end;%default scalar 
@@ -185 +172 @@
 %% select the priority order for automatic vel_type selection
 testder=0;
+testpatch=0;
 for ilist=1:length(FieldNames)
     if ~isempty(FieldNames{ilist})
     switch FieldNames{ilist}
+        case{'u','v'}
+            testpatch=1;
         case {'vort','div','strain'}
             testder=1;
     end
     end
-end      
+end    
+switch vel_type
+    case 'civ1'
+        var={'X','Y','Z','U','V','W','C','F','FF';...
+              'Civ1_X','Civ1_Y','Civ1_Z','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_F','Civ1_FF'};
+        role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag'};
+        units={'pixel','pixel','pixel','pixel','pixel','pixel',[],[],[]};
+    case 'interp1'
+         var={'X','Y','Z','U','V','W','FF';...
+               'Civ1_X','Civ1_Y','','Civ1_U_Diff','Civ1_V_Diff','','Civ1_FF'};
+         role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','errorflag'};  
+         units={'pixel','pixel','pixel','pixel','pixel','pixel',[]};
+    case 'filter1'
+        var={'X_tps','Y_tps','Z_tps','U_tps','V_tps','W_tps','X_SubRange','Y_SubRange','NbSites';...
+            'Civ1_X_tps','Civ1_Y_tps','','Civ1_U_tps','Civ1_V_tps','','Civ1_X_SubRange','Civ1_Y_SubRange','Civ1_NbSites'};
+         role={'','','','','','','','',''};  
+         units={'pixel','pixel','pixel','pixel','pixel','pixel','pixel','pixel',''};
+    otherwise % if VelType=[]
+        if testpatch
+           var={'X_tps','Y_tps','Z_tps','U_tps','V_tps','W_tps','X_SubRange','Y_SubRange','NbSites';...
+               'Civ2_X_tps','Civ2_Y_tps','','Civ2_U_tps','Civ2_V_tps','','Civ2_X_SubRange','Civ2_Y_SubRange','Civ2_NbSites';...
+               'Civ1_X_tps','Civ1_Y_tps','','Civ1_U_tps','Civ1_V_tps','','Civ1_X_SubRange','Civ1_Y_SubRange','Civ1_NbSites'};
+            role={'','','','','','','',''};  
+            units={'pixel','pixel','pixel','pixel','pixel','pixel','pixel','pixel'};
+        else
+             var={'X','Y','Z','U','V','W','C','F','FF';...
+              'Civ2_X','Civ2_Y','Civ2_Z','Civ2_U','Civ2_V','Civ2_W','Civ2_C','Civ2_F','Civ2_FF';...
+              'Civ1_X','Civ1_Y','Civ1_Z','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_F','Civ1_FF'};
+            role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag'};
+            units={'pixel','pixel','pixel','pixel','pixel','pixel',[],[],[]};
+        end   
+end
 if isempty(vel_type) || isequal(vel_type,'*') %undefined velocity type (civ1,civ2...)
     if testder
@@ -219 +240 @@
 vel_type_out=vel_type_out';
 
-%% determine names of netcdf variables to read
-var={'X','Y','Z','U','V','W','C','F','FF'};
-role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag'};
-units={'pixel','pixel','pixel','pixel','pixel','pixel',[],[],[]};
-if testder
-    var=[var {'DjUi(:,1,1)','DjUi(:,1,2)','DjUi(:,2,1)','DjUi(:,2,2)'}];
-    role=[role {'tensor','tensor','tensor','tensor'}];
-    units=[units {'pixel','pixel','pixel','pixel'}];
-end
-for ilist=1:length(vel_type_out)
-    var=[var;varname1(vel_type_out{ilist},FieldNames)];
-end
-
-%------------------------------------------------------------------------  
-% TAKEN FROM read_civxdata NOT USED
-%--- determine  var names to read
-
-function varin=varname1(vel_type,FieldNames)
-%------------------------------------------------------------------------
-testder=0;
-C1='';
-C2='';
-for ilist=1:length(FieldNames)
-    if ~isempty(FieldNames{ilist})
-    switch FieldNames{ilist}
-        case 'ima_cor' %image correlation corresponding to a vel vector
-            C1='vec_C';
-            C2='vec2_C';
-        case 'error'
-            C1='vec_E';
-            C2='vec2_E';
-        case {'vort','div','strain'}
-            testder=1;
-    end
-    end
-end      
-switch vel_type
-    case 'civ1'
-        varin={'Civ1_X','Civ1_Y','Civ1_Z','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_F','Civ1_FF'};
-    case 'interp1'
-        varin={'','','','','','','','',''};
-    case 'filter1'
-        varin={'Civ1_X','Civ1_Y','','Civ1_U_Spline','Civ1_V_Spline','','','',''};
-    case 'civ2'
-        varin={'vec2_X','vec2_Y','vec2_Z','vec2_U','vec2_V','vec2_W',C2,'vec2_F','vec2_FixFlag'};
-    case 'interp2'
-        varin={'vec2_patch_X','vec2_patch_Y','vec2_patch_Z','vec2_patch0_U','vec2_patch0_V','vec2_patch0_W','','',''};
-    case 'filter2'
-        varin={'vec2_patch_X','vec2_patch_Y','vec2_patch_Z','vec2_patch_U','vec2_patch_V','vec2_patch0_W','','',''};
-end
-if testder
-     switch vel_type
-        case 'filter1'
-            varin=[varin {'vec_patch_DUDX','vec_patch_DVDX','vec_patch_DUDY','vec_patch_DVDY'}];
-        case 'filter2'
-            varin=[varin {'vec2_patch_DUDX','vec2_patch_DVDX','vec2_patch_DUDY','vec2_patch_DVDY'}];
-    end   
-end
+
+
+
+
+

trunk/src/series.m

@@ -129 +129 @@
 %     set(handles.VelTypeMenu,'String',param.VelTypeMenu)
 % end  
-set(hObject,'WindowButtonUpFcn',{@mouse_up_gui,handles}) 
+set(hObject,'WindowButtonUpFcn',{'mouse_up_gui',handles}) 
 NomType_Callback(hObject, eventdata, handles)
 
@@ -1858 +1858 @@
 end   
 
-%-----------------------------
-function mouse_up_gui(hObject,eventdata,handles)
-if isequal(get(hObject,'SelectionType'),'alt')
-    set(hObject,'Units','pixels')
-    series_pos=get(hObject,'Position');%position of the GUI series (in pixels)
-    set(hObject,'Units','normalized')
-    xy_fig=get(hObject,'CurrentPoint');% current point of the current figure (gcbo)
-    hchild=get(hObject,'Children');%handles of all objects in the current figure
-    %% loop on all the objects in the current figure (selected by the last mouse click)
-
-    for ichild=1:length(hchild)
-        obj_pos=get(hchild(ichild),'Position');%position of the object        
-        if numel(obj_pos)>=4 && xy_fig(1) >=obj_pos(1) && xy_fig(2) >= obj_pos(2)&& xy_fig(1) <=obj_pos(1)+obj_pos(3) && xy_fig(2) <= obj_pos(2)+obj_pos(4);         
-            htype=get(hchild(ichild),'Type');%type of object child of the current figure
-            %if the mouse is over an axis, look at the data
-            if isequal(htype,'uicontrol') && isequal(get(hchild(ichild),'Visible'),'on')
-                msg_pos(1:2)=series_pos(1:2)+obj_pos(1:2).*series_pos(3:4);
-                msgbox_uvmat(['uicontrol: ' get(hchild(ichild),'Tag')],'',get(hchild(ichild),'String'),msg_pos)
-                break
-            end
-        end
-    end
-    set(hObject,'Units','pixels')
-end
 
 %%%%%%%%%%%%%

trunk/src/struct2nc.m

@@ -35 +35 @@
 function errormsg=struct2nc(flname,Data)
 if ~ischar(flname)
-    errormsg='no name input for the netcf file';
+    errormsg='invalid input for the netcf file name';
     return
 end

trunk/src/uvmat.m

@@ -2150 +2150 @@
         errormsg=['error in reading ' filename ': ' errormsg];
         return
-    end
-        
+    end        
     if isfield(ParamOut,'Npx')&& isfield(ParamOut,'Npy')
         set(handles.npx,'String',num2str(ParamOut.Npx));% display image size on the interface
@@ -2296 +2295 @@
     index_menu=strcmp(ParamOut.VelType,menu);
     set(handles.VelType,'Value',find(index_menu,1))
+    if ~get(handles.SubField,'value')
     set(handles.VelType,'String',menu)
-%     set(handles.VelType_1,'Value',1)
-%     set(handles.VelType_1,'String',[{''};menu])
+     set(handles.VelType_1,'Value',1)
+     set(handles.VelType_1,'String',[{''};menu])
+    end
 end
 field_index=strcmp(ParamOut.FieldName,ParamOut.FieldList);
@@ -2417 +2418 @@
     Field{1}=calc_field([{ParamOut.FieldName} {ParamOut.ColorVar}],Field{1});
 end
-if length(Field)==2 && ~test_keepdata_1 && isequal(FileType_1,'netcdf') && ~isequal(ParamOut_1.FieldName,'get_field...')%&&~isempty(FieldName_1)
+if numel(Field)==2 && ~test_keepdata_1 && isequal(FileType_1,'netcdf') && ~isequal(ParamOut_1.FieldName,'get_field...')%&&~isempty(FieldName_1)
     Field{2}=calc_field([{ParamOut_1.FieldName} {ParamOut_1.ColorVar}],Field{2});
 end
@@ -2510 +2511 @@
             nbpoints_z=UvData.Field.DimValue(DimIndex(1));
             DZ=(ZMax-ZMin)/(nbpoints_z-1);
-            UvData.Field.Mesh=sqrt(DX*DY*DZ);
+            UvData.Field.Mesh=(DX*DY*DZ)^(1/3);
             UvData.Field.ZMax=ZMax;
             UvData.Field.ZMin=ZMin;
@@ -3126 +3127 @@
 end
 FileExt_1=get(handles.FileExt_1,'String');
-if isequal(FileExt_1,'"'),FileExt_1=get(handles.FileExt,'String'); end;
+if isequal(get(handles.FileExt_1,'Visible'),'off') || isequal(FileExt_1,'"')
+    FileExt_1=get(handles.FileExt,'String');%read FileExt by default
+end
 FileName_1=[FileName_1 FileIndices_1 FileExt_1];
 
@@ -3272 +3275 @@
 
 %read the rootfile input display
-[FileName,RootPath,FileBase,FileIndices,FileExt_prev]=read_file_boxes_1(handles);
-[P,F,str1,str2,str_a,str_b,E,NomType]=name2display(['xxx' get(handles.FileIndex,'String') FileExt_prev]);
-if isempty(FileExt_prev)|| strcmp(FileExt_prev,'')
-    FileExt_1=get(handles.FileExt,'String');
-else
-    FileExt_1=FileExt_prev;
-end
+[FileName,RootPath,FileBase,FileIndices,FileExt_1]=read_file_boxes_1(handles);
+[P,F,str1,str2,str_a,str_b,E,NomType]=name2display(['xxx' get(handles.FileIndex,'String') FileExt_1]);
+% if isempty(FileExt_prev)|| strcmp(FileExt_prev,'')
+%     FileExt_1=get(handles.FileExt,'String');
+% else
+%     FileExt_1=FileExt_prev;
+% end
 NomType_1=get(handles.FileIndex_1,'UserData');
 if isempty(NomType_1)|| strcmp(NomType_1,'')
@@ -3350 +3353 @@
 else
     set(handles.SubDir_1,'Visible','on')
-    if ~isequal(FileExt_prev,'.nc') %find the new NomType if the previous display was not already a netcdf file
+    if ~isequal(FileExt_1,'.nc') %find the new NomType if the previous display was not already a netcdf file
 %         veltype_handles=[handles.VelType_1 handles.interp1_1 handles.filter1_1 handles.civ2_1 handles.interp2_1 handles.filter2_1];
 %         set_veltype_display(veltype_handles,6); % make all civ buttons visible
@@ -3496 +3499 @@
 %---------------------------------------------
   
-set(handles.FixVelType,'Value',1)
+set(handles.FixVelType,'Value',1)% the velocity type is now imposed by the GUI (not automatic)
 %refresh field with a second filename=first fiel name
 set(handles.run0,'BackgroundColor',[1 1 0])%paint the command button in yellow
-drawnow
+drawnow    
 filename=read_file_boxes(handles);
-if get(handles.SubField,'Value')
-    filename_1=read_file_boxes_1(handles);
-else
-    index=get(handles.VelType_1,'Value');
-    if index==1
+
+index=get(handles.VelType_1,'Value');
+if index==1
         filename_1='';% we plot the current field without the second field
-    else
-        filename_1=filename;
-    end
-end
+        set(handles.SubField,'Value',0)
+        SubField_Callback(hObject, eventdata, handles)
+elseif get(handles.SubField,'Value')% if subfield is already 'on'
+    filename_1=read_file_boxes_1(handles); %read the current second field
+else
+     filename_1=filename;% we compare two fields in the same file
+     set(handles.SubField,'Value',1)
+end
+
 num_i1=stra2num(get(handles.i1,'String'));
 num_i2=stra2num(get(handles.i2,'String'));