Changeset 246


Ignore:
Timestamp:
Apr 28, 2011, 10:52:31 AM (10 years ago)
Author:
sommeria
Message:

thin plate shell (patch) introduced

Location:
trunk/src
Files:
3 added
9 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/calc_field.m

    r236 r246  
    6060        nbcoord=2;
    6161    end
    62     DataOut=DataIn; %reproduce global attribute
    6362    ListVarName={};
    6463    ValueList={};
    6564    RoleList={};
    6665    units_cell={};
    67     for ilist=1:length(FieldName)
    68         if ~isempty(FieldName{ilist})
    69             [VarName,Value,Role,units]=feval(FieldName{ilist},DataIn);%calculate field with appropriate function named FieldName{ilist}
    70             ListVarName=[ListVarName VarName];
    71             ValueList=[ValueList Value];
    72             RoleList=[RoleList Role];
    73             units_cell=[units_cell units];
    74         end
    75     end
    76     %erase previous data (except coordinates)
    77     for ivar=nbcoord+1:length(DataOut.ListVarName)
    78         VarName=DataOut.ListVarName{ivar};
    79         DataOut=rmfield(DataOut,VarName);
    80     end
    81     DataOut.ListVarName=DataOut.ListVarName(1:nbcoord);
    82     if isfield(DataOut,'VarDimName')
    83         DataOut.VarDimName=DataOut.VarDimName(1:nbcoord);
    84     else
    85         errormsg='element .VarDimName missing in input data';
    86         return
    87     end
    88     DataOut.VarAttribute=DataOut.VarAttribute(1:nbcoord);
    89     %append new data
    90     DataOut.ListVarName=[DataOut.ListVarName ListVarName];
    91     for ivar=1:length(ListVarName)
    92         DataOut.VarDimName{nbcoord+ivar}=DataOut.VarDimName{1};
    93         DataOut.VarAttribute{nbcoord+ivar}.Role=RoleList{ivar};
    94         DataOut.VarAttribute{nbcoord+ivar}.units=units_cell{ivar};
    95         eval(['DataOut.' ListVarName{ivar} '=ValueList{ivar};'])
    96     end
    97 end
     66    % new civ data
     67    if isfield(DataIn,'X_SubRange')
     68        XMax=max(max(DataIn.X_SubRange));
     69        YMax=max(max(DataIn.Y_SubRange));
     70        XMin=min(min(DataIn.X_SubRange));
     71        YMin=min(min(DataIn.Y_SubRange));
     72        Mesh=sqrt((YMax-YMin)*(XMax-XMin)/numel(DataIn.X_tps));%2D
     73        xI=XMin:Mesh:XMax;
     74        yI=YMin:Mesh:YMax;
     75        [XI,YI]=meshgrid(xI,yI);
     76        XI=reshape(XI,[],1);
     77        YI=reshape(YI,[],1);
     78        DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute; %reproduce global attribute
     79        for ilist=1:numel(DataOut.ListGlobalAttribute)
     80            eval(['DataOut.' DataOut.ListGlobalAttribute{ilist} '=DataIn.' DataIn.ListGlobalAttribute{ilist} ';'])
     81        end
     82        DataOut.ListVarName={'coord_y','coord_x','FF'};
     83        DataOut.VarDimName{1}='coord_y';
     84        DataOut.VarDimName{2}='coord_x';
     85        DataOut.coord_y=[yI(1) yI(end)];
     86        DataOut.coord_x=[xI(1) xI(end)];
     87        DataOut.U=zeros(size(XI));
     88        DataOut.V=zeros(size(XI));
     89        DataOut.vort=zeros(size(XI));
     90        DataOut.div=zeros(size(XI));
     91        nbval=zeros(size(XI));
     92        [NbSubDomain,xx]=size(DataIn.X_SubRange);
     93        for isub=1:NbSubDomain
     94            nbvec_sub=DataIn.NbSites(isub);         
     95                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));
     96                %rho smoothing parameter
     97                epoints = [XI(ind_sel) YI(ind_sel)];% coordinates of interpolation sites
     98                ctrs=[DataIn.X_tps(1:nbvec_sub,isub) DataIn.Y_tps(1:nbvec_sub,isub)];%(=initial points) ctrs
     99                nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
     100%                 PM = [ones(size(epoints,1),1) epoints];
     101                switch FieldName{1}
     102                    case {'velocity','u','v'}
     103                       % DM_eval = DistanceMatrix(epoints,ctrs);%2D matrix of distances between extrapolation points epoints and spline centres (=site points) ctrs
     104                       % EM = tps(1,DM_eval);%values of thin plate
     105                        EM = tps_eval(epoints,ctrs);
     106                    case{'vort','div'}
     107                        EMDXY = tps_eval_dxy(epoints,ctrs);%2D matrix of distances between extrapolation points epoints and spline centres (=site points) ctrs
     108     
     109                end
     110                switch FieldName{1}
     111                    case 'velocity'
     112                        ListFields={'U', 'V'};
     113                        VarAttributes{1}.Role='vector_x';
     114                        VarAttributes{2}.Role='vector_y';
     115                        DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
     116                        DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
     117                    case 'u'
     118                        ListFields={'U'};
     119                        VarAttributes{1}.Role='scalar';
     120                        DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
     121                    case 'v'
     122                        ListFields={'V'};
     123                        VarAttributes{1}.Role='scalar';
     124                        DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
     125                    case 'vort'
     126                        ListFields={'vort'};
     127                        VarAttributes{1}.Role='scalar';
     128%                         EMX = [DMXY(:,:,1) PM];
     129%                         EMY = [DMXY(:,:,2) PM];
     130%                         DataIn.U_tps(nbvec_sub+1,isub)=0;%constant value suppressed by spatial derivatives
     131%                         DataIn.V_tps(nbvec_sub+1,isub)=0;%constant value suppressed by spatial derivatives
     132%                         DataIn.V_tps(nbvec_sub+2,isub)=0;% X coefficient suppressed for x wise derivatives
     133%                         DataIn.U_tps(nbvec_sub+3,isub)=0;% Y coefficient suppressed for x wise derivatives
     134                        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);
     135                    case 'div'
     136                        ListFields={'div'};
     137                        VarAttributes{1}.Role='scalar';
     138%                         EMX = [DMXY(:,:,1) PM];
     139%                         EMY = [DMXY(:,:,2) PM];
     140%                         DataIn.U_tps(nbvec_sub+1,isub)=0;%constant value suppressed by spatial derivatives
     141%                         DataIn.V_tps(nbvec_sub+1,isub)=0;%constant value suppressed by spatial derivatives
     142%                         DataIn.V_tps(nbvec_sub+2,isub)=0;% X coefficient suppressed for x wise derivatives
     143%                         DataIn.U_tps(nbvec_sub+3,isub)=0;% Y coefficient suppressed for x wise derivatives
     144                        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);
     145                end
     146        end
     147        DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang 
     148       
     149        DataOut.FF=reshape(DataOut.FF,numel(yI),numel(xI));
     150        nbval(nbval==0)=1;
     151        DataOut.U=DataOut.U ./nbval;
     152        DataOut.V=DataOut.V ./nbval;
     153        DataOut.U=reshape(DataOut.U,numel(yI),numel(xI));
     154        DataOut.V=reshape(DataOut.V,numel(yI),numel(xI));
     155        DataOut.vort=reshape(DataOut.vort,numel(yI),numel(xI));
     156        DataOut.div=reshape(DataOut.div,numel(yI),numel(xI));
     157        DataOut.ListVarName=[DataOut.ListVarName ListFields];
     158        for ilist=3:numel(DataOut.ListVarName)
     159            DataOut.VarDimName{ilist}={'coord_y','coord_x'};
     160        end
     161        DataOut.VarAttribute={[],[]};
     162        DataOut.VarAttribute{3}.Role='errorflag';
     163        DataOut.VarAttribute=[DataOut.VarAttribute VarAttributes];     
     164    else   
     165        DataOut=DataIn;
     166        for ilist=1:length(FieldName)
     167            if ~isempty(FieldName{ilist})
     168                [VarName,Value,Role,units]=feval(FieldName{ilist},DataIn);%calculate field with appropriate function named FieldName{ilist}
     169                ListVarName=[ListVarName VarName];
     170                ValueList=[ValueList Value];
     171                RoleList=[RoleList Role];
     172                units_cell=[units_cell units];
     173            end
     174        end
     175        %erase previous data (except coordinates)
     176        for ivar=nbcoord+1:length(DataOut.ListVarName)
     177            VarName=DataOut.ListVarName{ivar};
     178            DataOut=rmfield(DataOut,VarName);
     179        end
     180        DataOut.ListVarName=DataOut.ListVarName(1:nbcoord);
     181        if isfield(DataOut,'VarDimName')
     182            DataOut.VarDimName=DataOut.VarDimName(1:nbcoord);
     183        else
     184            errormsg='element .VarDimName missing in input data';
     185            return
     186        end
     187        DataOut.VarAttribute=DataOut.VarAttribute(1:nbcoord);
     188        %append new data
     189        DataOut.ListVarName=[DataOut.ListVarName ListVarName];
     190        for ivar=1:length(ListVarName)
     191            DataOut.VarDimName{nbcoord+ivar}=DataOut.VarDimName{1};
     192            DataOut.VarAttribute{nbcoord+ivar}.Role=RoleList{ivar};
     193            DataOut.VarAttribute{nbcoord+ivar}.units=units_cell{ivar};
     194            eval(['DataOut.' ListVarName{ivar} '=ValueList{ivar};'])
     195        end
     196    end
     197end
     198
    98199
    99200%%%%%%%%%%%%% velocity fieldn%%%%%%%%%%%%%%%%%%%%
  • trunk/src/civ.m

    r244 r246  
    6262% Update handles structure
    6363guidata(hObject, handles);
    64 
     64set(hObject,'WindowButtonUpFcn',{'mouse_up_GUI',handles}) %set mouse action (zoom on uicontrols)
    6565%default initial parameters
    6666filebase=''; % root file name ('filebase'.civ)
     
    334334    browse.nom_type_nc=nom_type;
    335335    ind_opening=2;% propose 'fix' as the default option
    336     Data=nc2struct(fileinput,[]);
    337     if isfield(Data,'absolut_time_T0')%test for civx files
     336    Data=nc2struct(fileinput,'ListGlobalAttribute','CivStage','absolut_time_T0','fix','patch','civ2','fix2');
     337    if ~isempty(Data.CivStage)%test for civ files
     338        ind_opening=Data.CivStage;
     339        set(handles.CivAll,'Value',3)
     340    end
     341    if ~isempty(Data.absolut_time_T0)%test for civx files
     342        set(handles.CivAll,'Value',1)
    338343        if isfield(Data,'fix') && isequal(Data.fix,1)
    339344            ind_opening=3;
    340345        end
    341         if isfield(Data,'patch') && isequal(Data.patch,1)
     346        if isequal(Data.patch,1)
    342347            ind_opening=4;
    343348        end
    344         if isfield(Data,'civ2') && isequal(Data.civ2,1)
     349        if isequal(Data.civ2,1)
    345350            ind_opening=5;
    346351        end
    347         if isfield(Data,'fix2') && isequal(Data.fix2,1)
     352        if isequal(Data.fix2,1)
    348353            ind_opening=6;
    349354        end
     
    46834688    par_civ1.T0=0;
    46844689    par_civ1.Dt=1;
    4685     Data=civ_uvmat(par_civ1);
     4690    Param.Civ1=par_civ1;
     4691    Data=civ_uvmat(Param);
    46864692%     stepx=str2num(par_civ1.dx);
    46874693%     stepy=str2num(par_civ1.dy);
     
    48774883function open_view_field(hObject, eventdata)
    48784884%-------------------------------------------------------------------
    4879      list=get(gcbo,'String');
    4880      index=get(gcbo,'Value');
    4881      rootroot=get(gcbo,'UserData');
     4885     list=get(hObject,'String');
     4886     index=get(hObject,'Value');
     4887     rootroot=get(hObject,'UserData');
    48824888     filename=list{index};
    48834889     ind_dot=findstr(filename,'...');
    48844890     filename=filename(1:ind_dot-1);
    48854891      filename=fullfile(rootroot,filename);
     4892      delete(get(hObject,'parent'))%delete the display figure to stop the check process
    48864893      if exist(filename,'file')%visualise the vel field if it exists
    4887         %[Field,VelTypeOut]=read_civxdata(filename);
    4888         %view_field(Field)
    48894894        uvmat(filename)
    48904895        set(gcbo,'Value',1)
  • trunk/src/civ_uvmat.m

    r243 r246  
    11% To develop....
    22function [Data,errormsg]= civ_uvmat(Param,ncfile)
     3errormsg='';
    34Data.ListGlobalAttribute={'Conventions','Program','CivStage'};
    45Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes
    56Data.Program='civ_uvmat';
    67Data.CivStage=0;%default
     8ListVarCiv1={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F'};
     9ListVarFix1={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F','Civ1_FF'};
    710
    811%% Civ1
     
    2023    shiftx=str2num(par_civ1.shiftx);
    2124    shifty=str2num(par_civ1.shifty);
    22     miniy=max(1+isy2-shifty,1+iby2);
     25    miniy=max(1+isy2+shifty,1+iby2);
    2326    minix=max(1+isx2-shiftx,1+ibx2);
    24     maxiy=min(size(image1,1)-isy2-shifty,size(image1,1)-iby2);
     27    maxiy=min(size(image1,1)-isy2+shifty,size(image1,1)-iby2);
    2528    maxix=min(size(image1,2)-isx2-shiftx,size(image1,2)-ibx2);
    2629    [GridX,GridY]=meshgrid(minix:stepx:maxix,miniy:stepy:maxiy);
     
    2932   
    3033    % caluclate velocity data (y and v in indices, reverse to y component)
    31     [xtable ytable utable vtable ctable F] = pivlab (image1,image2,ibx2,iby2,isx2,isy2,shiftx,shifty,PointCoord,str2num(par_civ1.rho), []);
     34    [xtable ytable utable vtable ctable F] = pivlab (image1,image2,ibx2,iby2,isx2,isy2,shiftx,-shifty,PointCoord,str2num(par_civ1.rho), []);
    3235    list_param=(fieldnames(par_civ1))';
    3336    list_remove={'pxcmx','pxcmy','npx','npy','gridflag','maskflag','term_a','term_b','T0'};
     
    6770    Data.CivStage=1;
    6871else
    69     Data=nc2struct(ncfile)%read existing netcdf file
     72    if isfield(Param,'Fix1')
     73        Data=nc2struct(ncfile,ListVarCiv1);%read civ1 data in the existing netcdf file
     74    else
     75        Data=nc2struct(ncfile,ListVarFix1);%read civ1 and fix1 data in the existing netcdf file
     76    end
     77    if isfield(Data,'Txt')
     78        msgbox_uvmat('ERROR',Data.Txt)
     79        return
     80    end
     81    % read Civx data
    7082    if isfield(Data,'absolut_time_T0')%read civx file
    7183        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'};
     
    95107        ParamName=ListFixParam{ilist};
    96108        ListName=['Fix1_' ParamName];
    97         ['Data.ListGlobalAttribute=[Data.ListGlobalAttribute ''' ParamName '''];']
    98109        eval(['Data.ListGlobalAttribute=[Data.ListGlobalAttribute ''' ParamName '''];'])
    99110        eval(['Data.' ListName '=Param.Fix1.' ParamName ';'])
     
    107118    Data.VarDimName=[Data.VarDimName {'nbvec1'}];
    108119    nbvar=length(Data.ListVarName);
    109     Data.VarAttribute{nbvar}.Role='errorflag';
    110     [Data.Civ1_FF]=fix_uvmat(Param.Fix1,Data.Civ1_F,Data.Civ1_C,Data.Civ1_U,Data.Civ1_V);
     120    Data.VarAttribute{nbvar}.Role='errorflag';   
     121    Data.Civ1_FF=fix_uvmat(Param.Fix1,Data.Civ1_F,Data.Civ1_C,Data.Civ1_U,Data.Civ1_V);
    111122    Data.CivStage=2;                               
    112123end   
     
    117128    Data.Patch1_Threshold=str2double(Param.Patch1.Threshold);
    118129    Data.Patch1_SubDomain=str2double(Param.Patch1.SubDomain);
    119     Data.ListVarName=[Data.ListVarName {'Patch1_U','Patch1_V'}];
    120     Data.VarDimName=[Data.VarDimName {'nbvec1','nbvec1'}];
     130    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'}];
     131    Data.VarDimName=[Data.VarDimName {'NbVec1','NbVec1',{'NbSubDomain1','Two'},{'NbSubDomain1','Two'},'NbSubDomain1',...
     132             {'NbVec1Sub','NbSubDomain1'},{'NbVec1Sub','NbSubDomain1'},{'Nbtps1','NbSubDomain1'},{'Nbtps1','NbSubDomain1'}}];
    121133    nbvar=length(Data.ListVarName);
    122134    Data.VarAttribute{nbvar-1}.Role='vector_x';
    123135    Data.VarAttribute{nbvar}.Role='vector_y';
    124     Data.Patch1_U=zeros(size(Data.Civ1_X));
    125     Data.Patch1_V=zeros(size(Data.Civ1_X));
     136    Data.Civ1_U_Diff=zeros(size(Data.Civ1_X));
     137    Data.Civ1_V_Diff=zeros(size(Data.Civ1_X));
    126138    if isfield(Data,'Civ1_FF')
    127139        ind_good=find(Data.Civ1_FF==0);
     
    129141        ind_good=1:numel(Data.Civ1_X);
    130142    end
    131     Data.Civ1_X
    132     [Ures, Vres,FFres]=...
     143    [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]=...
    133144                            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);
    134                         size(Ures)
    135                         size(Vres)
    136                         size(FFres)
    137                         size(ind_good)
    138       Data.Patch1_U(ind_good)=Ures;
    139       Data.Patch1_V(ind_good)=Vres;
    140       Data.Civ1_FF(ind_good)=FFres
    141       Data.CivStage=3;                               
     145      Data.Civ1_U_Diff(ind_good)=Data.Civ1_U(ind_good)-Ures;
     146      Data.Civ1_V_Diff(ind_good)=Data.Civ1_V(ind_good)-Vres;
     147      Data.Civ1_FF(ind_good)=FFres;
     148      Data.CivStage=3;                             
    142149end   
    143150%% write result
     151% 'TESTcalc'
     152% [DataOut,errormsg]=calc_field('velocity',Data)
     153if exist('ncfile','var')
    144154errormsg=struct2nc(ncfile,Data);
    145 
     155end
    146156
    147157
     
    192202    FF=FF==1 | (U.*U+V.*V)>thresh;
    193203end
    194  
     204FF=double(FF);
    195205%
    196206% % criterium on velocity values
     
    270280%------------------------------------------------------------------------
    271281% patch function
    272 function [U_patch,V_patch,FF,SubRangx,SubRangy,X_ctrs,Y_ctrs] =patch_uvmat(X,Y,U,V,Rho,Threshold,SubDomain)
     282function [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)
    273283%subdomain decomposition
    274284warning off
     
    278288Y=reshape(Y,[],1);
    279289nbvec=numel(X);
    280 NbSubDomain=ceil(nbvec/SubDomain)
    281 MinX=min(X)
    282 MinY=min(Y)
    283 MaxX=max(X)
    284 MaxY=max(Y)
     290NbSubDomain=ceil(nbvec/SubDomain);
     291MinX=min(X);
     292MinY=min(Y);
     293MaxX=max(X);
     294MaxY=max(Y);
    285295RangX=MaxX-MinX;
    286296RangY=MaxY-MinY;
    287 AspectRatio=RangY/RangX
    288 NbSubDomainX=ceil(sqrt(NbSubDomain/AspectRatio))
    289 NbSubDomainY=ceil(sqrt(NbSubDomain*AspectRatio))
    290 
     297AspectRatio=RangY/RangX;
     298NbSubDomainX=max(floor(sqrt(NbSubDomain/AspectRatio)),1);
     299NbSubDomainY=max(floor(sqrt(NbSubDomain*AspectRatio)),1);
     300NbSubDomain=NbSubDomainX*NbSubDomainY;
    291301SizX=RangX/NbSubDomainX;%width of subdomains
    292302SizY=RangY/NbSubDomainY;%height of subdomains
    293303CentreX=linspace(MinX+SizX/2,MaxX-SizX/2,NbSubDomainX);
    294304CentreY=linspace(MinY+SizY/2,MaxY-SizY/2,NbSubDomainY);
    295 SubIndexX=ceil((X-MinX)/SizX);%subdomain index of vectors
    296 SubIndexY=ceil((Y-MinY)/SizY);
    297 rho=RangX*RangY*Rho;%optimum rho increase as teh area of the subdomain
    298 U_tps=zeros(length(X),NbSubDomainY,NbSubDomainX);%default spline
    299 V_tps=zeros(length(X),NbSubDomainY,NbSubDomainX);%default spline
    300 X_dist=zeros(length(X),NbSubDomainY,NbSubDomainX);%default spline
    301 Y_dist=zeros(length(X),NbSubDomainY,NbSubDomainX);%default spline
    302 U_smooth=zeros(length(X),NbSubDomainY,NbSubDomainX);
    303 V_smooth=zeros(length(X),NbSubDomainY,NbSubDomainX);
    304 dist_ctre=zeros(length(X),NbSubDomainY,NbSubDomainX);
     305[CentreX,CentreY]=meshgrid(CentreX,CentreY);
     306CentreY=reshape(CentreY,1,[]);
     307CentreX=reshape(CentreX,1,[]);
     308rho=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)
     309U_tps_sub=zeros(length(X),NbSubDomain);%default spline
     310V_tps_sub=zeros(length(X),NbSubDomain);%default spline
     311U_smooth=zeros(length(X),1);
     312V_smooth=zeros(length(X),1);
     313
     314nb_select=zeros(length(X),1);
    305315FF=zeros(length(X),1);
    306 for isubx=1:NbSubDomainX
    307     for isuby=1:NbSubDomainY
    308         SubRangx(isuby,isubx,:)=[CentreX(isubx)-SizX/2 CentreX(isubx)+SizX/2];
    309         SubRangy(isuby,isubx,:)=[CentreY(isubx)-SizY/2 CentreY(isubx)+SizY/2];
    310         for iter=1:3 %increase the subdomain during three iterations at most
    311             ind_sel=find(X>SubRangx(isuby,isubx,1) & X<SubRangx(isuby,isubx,2) & Y>SubRangy(isuby,isubx,1) & Y<SubRangy(isuby,isubx,2)); 
    312             size(ind_sel)
    313             if numel(ind_sel)<SubDomain/4;% too few selected vectors, increase the subrange for next iteration
    314                 SubRangx(isuby,isubx,1)=SubRangx(isuby,isubx,1)-SizX/4;
    315                 SubRangx(isuby,isubx,2)=SubRangx(isuby,isubx,2)+SizX/4;
    316                 SubRangy(isuby,isubx,1)=SubRangy(isuby,isubx,1)-SizY/4;
    317                 SubRangy(isuby,isubx,2)=SubRangy(isuby,isubx,2)+SizY/4;
     316test_empty=zeros(1,NbSubDomain);
     317for isub=1:NbSubDomain
     318    SubRangx(isub,:)=[CentreX(isub)-SizX/2 CentreX(isub)+SizX/2];
     319    SubRangy(isub,:)=[CentreY(isub)-SizY/2 CentreY(isub)+SizY/2];
     320    ind_sel_previous=[];
     321    ind_sel=0;
     322    while numel(ind_sel)>numel(ind_sel_previous) %increase the subdomain during four iterations at most
     323        ind_sel_previous=ind_sel;
     324        ind_sel=find(X>SubRangx(isub,1) & X<SubRangx(isub,2) & Y>SubRangy(isub,1) & Y<SubRangy(isub,2));
     325        % if no vector in the subdomain, skip the subdomain
     326        if isempty(ind_sel)
     327            test_empty(isub)=1;   
     328            U_tps(1,isub)=0;%define U_tps and V_tps by default
     329            V_tps(1,isub)=0;
     330            break
     331            % if too few selected vectors, increase the subrange for next iteration
     332        elseif numel(ind_sel)<SubDomain/4 && ~isequal( ind_sel,ind_sel_previous);
     333            SubRangx(isub,1)=SubRangx(isub,1)-SizX/4;
     334            SubRangx(isub,2)=SubRangx(isub,2)+SizX/4;
     335            SubRangy(isub,1)=SubRangy(isub,1)-SizY/4;
     336            SubRangy(isub,2)=SubRangy(isub,2)+SizY/4;
     337        else
     338            [U_smooth_sub,U_tps_sub]=tps_coeff(X(ind_sel),Y(ind_sel),U(ind_sel),rho);
     339            [V_smooth_sub,V_tps_sub]=tps_coeff(X(ind_sel),Y(ind_sel),V(ind_sel),rho);
     340            UDiff=U_smooth_sub-U(ind_sel);
     341            VDiff=V_smooth_sub-V(ind_sel);
     342            NormDiff=UDiff.*UDiff+VDiff.*VDiff;
     343            FF(ind_sel)=20*(NormDiff>Threshold);%put FF value to 20 to identify the criterium of elimmination
     344            ind_ind_sel=find(FF(ind_sel)==0); % select the indices of ind_sel corresponding to the remaining vectors
     345            % no value exceeds threshold, the result is recorded
     346            if isequal(numel(ind_ind_sel),numel(ind_sel))
     347                U_smooth(ind_sel)=U_smooth(ind_sel)+U_smooth_sub;
     348                V_smooth(ind_sel)=V_smooth(ind_sel)+V_smooth_sub;
     349                nbpoints(isub)=numel(ind_sel);
     350                X_tps(1:nbpoints(isub),isub)=X(ind_sel);
     351                Y_tps(1:nbpoints(isub),isub)=Y(ind_sel);
     352                U_tps(1:nbpoints(isub)+3,isub)=U_tps_sub;
     353                V_tps(1:nbpoints(isub)+3,isub)=V_tps_sub;         
     354                nb_select(ind_sel)=nb_select(ind_sel)+1;
     355                 display('good')
     356                break
     357                % too few selected vectors, increase the subrange for next iteration
     358            elseif numel(ind_ind_sel)<SubDomain/4 && ~isequal( ind_sel,ind_sel_previous);
     359                SubRangx(isub,1)=SubRangx(isub,1)-SizX/4;
     360                SubRangx(isub,2)=SubRangx(isub,2)+SizX/4;
     361                SubRangy(isub,1)=SubRangy(isub,1)-SizY/4;
     362                SubRangy(isub,2)=SubRangy(isub,2)+SizY/4;
     363%                 display('fewsmooth')
     364                % interpolation-smoothing is done again with the selected vectors
    318365            else
    319                 [U_smooth_sub,U_tps_sub]=tps_uvmat(X(ind_sel),Y(ind_sel),U(ind_sel),rho);
    320                 [V_smooth_sub,V_tps_sub]=tps_uvmat(X(ind_sel),Y(ind_sel),V(ind_sel),rho);
    321                 size(U_smooth_sub)
    322                 size(U(ind_sel))
    323                 UDiff=U_smooth_sub-U(ind_sel);
    324                 VDiff=V_smooth_sub-V(ind_sel);
    325                 NormDiff=UDiff.*UDiff+VDiff.*VDiff;
    326                 FF(ind_sel)=20*(NormDiff>Threshold);%put FF value to 20 to identify the criterium of elimmination
    327                 ind_ind_sel=find(FF(ind_sel)==0);
    328                 if isequal(numel(ind_ind_sel),numel(ind_sel))
    329                     U_smooth(ind_sel,isuby,isubx)=U_smooth_sub;
    330                     V_smooth(ind_sel,isuby,isubx)=V_smooth_sub;
    331                     break
    332                 elseif numel(ind_ind_sel)<SubDomain/4;% too few selected vectors, increase the subrange for next iteration
    333                     SubRangx(isuby,isubx,1)=SubRangx(isuby,isubx,1)-SizX/4;
    334                     SubRangx(isuby,isubx,2)=SubRangx(isuby,isubx,2)+SizX/4;
    335                     SubRangy(isuby,isubx,1)=SubRangy(isuby,isubx,1)-SizY/4;
    336                     SubRangy(isuby,isubx,2)=SubRangy(isuby,isubx,2)+SizY/4;
    337                 else
    338                     [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);
    339                     [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);
    340                      break                 
    341                 end
    342             end       
    343         end       
    344         X_ctrs(isuby,isubx)=mean(X(ind_sel));%gravity centre of selected points
    345         Y_ctrs(isuby,isubx)=mean(Y(ind_sel));%positions of tps sources for the subdomain i,     
    346         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)))));
    347     end
    348 end
    349 U_patch=sum(sum(U_smooth.*dist_ctre,3),2)./sum(sum(dist_ctre,3),2);
    350 V_patch=sum(sum(V_smooth.*dist_ctre,3),2)./sum(sum(dist_ctre,3),2);
    351 
    352 %------------------------------------------------------------------------
    353 %fasshauer@iit.edu MATH 590 ? Chapter 19 32
    354 % X,Y initial coordiantes
    355 % XI vector, YI column vector for the grid of interpolation points
    356 function [U_smooth,U_tps,U_tps3]=tps_uvmat(X,Y,U,rho)
    357 %------------------------------------------------------------------------
    358 %rho smoothing parameter
    359 ep = 1;
    360 X=reshape(X,[],1);
    361 Y=reshape(Y,[],1);
    362 rhs = reshape(U,[],1);
    363 % if exist('FF','var')
    364 % test_false=isnan(rhs)|FF~=0;
    365 % else
    366 %     test_false=isnan(rhs);
    367 % end
    368 % X(test_false)=[];
    369 % Y(test_false)=[];
    370 % rhs(test_false)=[];
    371 %randn('state',3); rhs = rhs + 0.03*randn(size(rhs));
    372 rhs = [rhs; zeros(3,1)];
    373 dsites = [X Y];% coordinates of measurement sites
    374 ctrs = dsites;%radial base functions are located at the measurement sites
    375 DM_data = DistanceMatrix(dsites,ctrs);%2D matrix of distances between spline centres (=initial points) ctrs
    376 % if size(XI,1)==1 && size(YI,2)==1 % XI vector, YI column vector
    377 %      [XI,YI]=meshgrid(XI,YI);
    378 % end
    379 % [npy,npx]=size(XI);
    380 % epoints = [reshape(XI,[],1) reshape(YI,[],1)];
    381 IM_sites = tps(ep,DM_data);%values of thin plate at site points
    382 IM = IM_sites + rho*eye(size(IM_sites));%  rho=1/(2*omega) , omega given by fasshauer;
    383 PM=[ones(size(dsites,1),1) dsites];
    384 IM=[IM PM; [PM' zeros(3,3)]];
    385 %fprintf('Condition number estimate: %e\n',condest(IM))
    386 %DM_eval = DistanceMatrix(epoints,ctrs);%2D matrix of distances between extrapolation points epoints and spline centres (=site points) ctrs
    387 %EM = tps(ep,DM_eval);%values of thin plate
    388 %PM = [ones(size(epoints,1),1) epoints];
    389 %EM = [EM PM];
    390 U_tps=(IM\rhs);
    391 PM = [ones(size(dsites,1),1) dsites];
    392 EM = [IM_sites PM];
    393 U_smooth=EM *U_tps;
    394 U_tps3=U_tps(end-2:end);
    395 U_tps=U_tps(1:end-3);
     366                [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);
     367                [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);
     368                U_smooth(ind_sel(ind_ind_sel))=U_smooth(ind_sel(ind_ind_sel))+U_smooth_sub;
     369                V_smooth(ind_sel(ind_ind_sel))=V_smooth(ind_sel(ind_ind_sel))+V_smooth_sub;
     370                nbpoints(isub)=numel(ind_ind_sel);
     371                X_tps(1:nbpoints(isub),isub)=X(ind_sel(ind_ind_sel));
     372                Y_tps(1:nbpoints(isub),isub)=Y(ind_sel(ind_ind_sel));
     373                U_tps(1:nbpoints(isub)+3,isub)=U_tps_sub;
     374                V_tps(1:nbpoints(isub)+3,isub)=V_tps_sub;
     375                nb_select(ind_sel(ind_ind_sel))=nb_select(ind_sel(ind_ind_sel))+1;
     376                display('good2')
     377                break
     378            end
     379        end
     380    end
     381end
     382ind_empty=find(test_empty);
     383%remove empty subdomains
     384if ~isempty(ind_empty)
     385    SubRangx(ind_empty,:)=[];
     386    SubRangy(ind_empty,:)=[];
     387    X_tps(:,ind_empty)=[];
     388    Y_tps(:,ind_empty)=[];
     389    U_tps(:,ind_empty)=[];
     390    V_tps(:,ind_empty)=[];
     391end
     392nb_select(nb_select==0)=1;%ones(size(find(nb_select==0)));
     393U_smooth=U_smooth./nb_select;
     394V_smooth=V_smooth./nb_select;
     395
     396
     397
     398
    396399
    397400
     
    421424
    422425
    423 %   DM:     MxN matrix whose i,j position contains the Euclidean
    424 %              distance between the i-th data site and j-th center
    425   function DM = DistanceMatrix(dsites,ctrs)
    426   [M,s] = size(dsites); [N,s] = size(ctrs);
    427   DM = zeros(M,N);
    428   % Accumulate sum of squares of coordinate differences
    429   % The ndgrid command produces two MxN matrices:
    430   %   dr, consisting of N identical columns (each containing
    431   %       the d-th coordinate of the M data sites)
    432   %   cc, consisting of M identical rows (each containing
    433   %       the d-th coordinate of the N centers)
    434   for d=1:s
    435      [dr,cc] = ndgrid(dsites(:,d),ctrs(:,d));
    436      DM = DM + (dr-cc).^2;
    437   end
    438   DM = sqrt(DM);
    439 
    440 
    441   % rbf = tps(e,r)
    442 % Defines thin plate spline RBF
    443 function rbf = tps(e,r)
    444 rbf = zeros(size(r));
    445 nz = find(r~=0);   % to deal with singularity at origin
    446 rbf(nz) = (e*r(nz)).^2.*log(e*r(nz));
    447 
     426
     427
  • trunk/src/msgbox_uvmat.m

    r244 r246  
    6363guidata(hObject, handles);
    6464testNo=0;
    65 testCancel=1;
     65testCancel=0;
    6666testinputstring=0;
    6767icontype='quest';%default question icon (text input asked)
     
    7171        case {'CONFIRMATION'}
    7272            icontype='';
    73             testCancel=0; %no cancel button
    7473        case 'ERROR'
    7574            icontype='error';
    76             testCancel=0; %no cancel button
    7775        case 'WARNING'
    7876            icontype='warn';
    79             testCancel=0; %no cancel button
    8077        case 'INPUT_Y-N'
    8178            icontype='quest';
     79            testCancel=1; %no cancel button
    8280            testNo=1; % button No activated
    8381        case {'RULER'}
    8482            icontype='';
    85             testCancel=0; %no cancel button
    8683            testinputstring=1;
    8784        case 'INPUT_TXT'
    8885            testinputstring=1;
     86            testCancel=1; %no cancel button
    8987        otherwise
    9088          %  testinputstring=1;
     
    189187
    190188% Get default command line output from handles structure
    191 if isequal(handles.output,'Cancel')
    192     varargout{1}='Cancel';
    193 elseif isequal(handles.output,'No')
    194     varargout{1}='No';
    195 else
    196     varargout{1}=get(handles.edit_box,'String');
    197     if isempty(varargout{1}) || isequal(varargout{1},'')
    198         varargout{1}='Yes';
    199     end
    200 end
    201 % The figure can be deleted now
    202 delete(handles.figure1);
    203 
     189if isfield(handles,'output')
     190    if isequal(handles.output,'Cancel')
     191        varargout{1}='Cancel';
     192    elseif isequal(handles.output,'No')
     193        varargout{1}='No';
     194    else
     195        varargout{1}=get(handles.edit_box,'String');
     196        if isempty(varargout{1}) || isequal(varargout{1},'')
     197            varargout{1}='Yes';
     198        end
     199    end
     200    % The figure can be deleted now
     201end
     202 delete(handles.figure1);
     203 
    204204% --- Executes on button press in OK.
    205205function OK_Callback(hObject, eventdata, handles)
  • trunk/src/pivlab.m

    r244 r246  
    1414% image1:first image (matrix)
    1515% image2: second image (matrix)
    16 % ibx,iby: size of the correlation box along x and y (in px)
     16% ibx2,iby2: half size of the correlation box along x and y, in px (size=(2*iby2+1,2*ibx2+1)
     17% isx2,isy2: half size of the search box along x and y, in px (size=(2*isy2+1,2*isx2+1)
     18% shiftx, shifty: shift of the search box (in pixel index, yshift reversed)
    1719% step: mesh of the measurement points (in px)
    1820% subpixfinder=1 or 2 controls the curve fitting of the image correlation
     
    119121    xtable(ivec)=iref+vector(1)/2;% convec flow (velocity taken at the point middle from imgae1 and 2)
    120122    ytable(ivec)=jref+vector(2)/2;
    121     utable(ivec)=vector(1);
    122     vtable(ivec)=vector(2);
     123    utable(ivec)=vector(1)+shiftx;
     124    vtable(ivec)=vector(2)+shifty;
    123125end
    124126result_conv=result_conv*corrmax/(255*sum_square);% keep the last correlation matrix for output
  • trunk/src/read_civdata.m

    r237 r246  
    4747% 'nc2struct': reads a netcdf file
    4848
    49 function [Field,VelTypeOut,errormsg]=read_civdata(filename,FieldNames,VelType)
     49function [Field,VelTypeOut,errormsg]=read_civdata(filename,FieldNames,VelType,XI,YI)
    5050errormsg='';
    5151VelTypeOut=VelType;%default
    5252%% default input
    5353if ~exist('VelType','var')
    54     VelType=[];
     54    VelType='';
    5555end
    5656if isequal(VelType,'*')
    57     VelType=[];
     57    VelType='';
     58end
     59if isempty(VelType)
     60    VelType='';
    5861end
    5962if ~exist('FieldNames','var')
     
    6265
    6366%% reading data
    64 %[var,role,units,vel_type_out_cell]=varcivx_generator(FieldNames,VelType);%determine the names of constants and variables to read
    65 % if isempty(VelType)
    66 %     VelType='civ1';
    67 % end
    68 % if isequal(VelType,'civ1')
    69 %     varlist=[{'X','Y','U','V','C','F','FF'};... %output names
    70 %         {'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F','Civ1_FF'}];%names in file
    71 %     role={'coord_x','coord_y','vector_x','vector_y','ancillary','warnflag','errorflag'};
    72 %     units={'pixel','pixel','pixel','pixel','','',''};
    73 % elseif isequal(VelType,'interp1')|| isequal(VelType,'filter1')% read diff between spline and initial data (change name interp1 -> splinediff
    74 %     varlist=[{'X','Y','U','V'};... %output names
    75 %         {'Civ1_X','Civ1_Y','Patch1_U','Patch1_V'}];
    76 %     role = {'coord_x','coord_y','vector_x','vector_y'};
    77 %     units={'pixel','pixel','pixel','pixel'};
    78 % end
    79 [varlist,role,units,vel_type_out_cell]=varcivx_generator(FieldNames,VelType);
     67[varlist,role,units,vel_type_out_cell]=varcivx_generator(FieldNames,VelType);
    8068[Field,vardetect,ichoice]=nc2struct(filename,varlist);%read the variables in the netcdf file
    8169if isfield(Field,'Txt')
     
    9381    Field.VarAttribute{ivar}.Mesh=0.1;%typical mesh for histograms O.1 pixel
    9482end
    95 
    9683if ~isempty(ichoice)
    9784    VelTypeOut=vel_type_out_cell{ichoice};
    9885end
    99 %  Field.CivStage=3;%for patch, to generalise
    100 % return
     86
    10187
    10288%% renaming for standard conventions
     
    162148
    163149
     150
    164151%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    165152% TAKEN FROM read_civxdata NOT USED
     
    178165
    179166%% default input values
    180 if ~exist('vel_type','var'),vel_type=[];end;
     167if ~exist('vel_type','var'),vel_type='';end;
    181168if iscell(vel_type),vel_type=vel_type{1}; end;%transform cell to string if needed
    182169if ~exist('FieldNames','var'),FieldNames={'ima_cor'};end;%default scalar
     
    185172%% select the priority order for automatic vel_type selection
    186173testder=0;
     174testpatch=0;
    187175for ilist=1:length(FieldNames)
    188176    if ~isempty(FieldNames{ilist})
    189177    switch FieldNames{ilist}
     178        case{'u','v'}
     179            testpatch=1;
    190180        case {'vort','div','strain'}
    191181            testder=1;
    192182    end
    193183    end
    194 end     
     184end   
     185switch vel_type
     186    case 'civ1'
     187        var={'X','Y','Z','U','V','W','C','F','FF';...
     188              'Civ1_X','Civ1_Y','Civ1_Z','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_F','Civ1_FF'};
     189        role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag'};
     190        units={'pixel','pixel','pixel','pixel','pixel','pixel',[],[],[]};
     191    case 'interp1'
     192         var={'X','Y','Z','U','V','W','FF';...
     193               'Civ1_X','Civ1_Y','','Civ1_U_Diff','Civ1_V_Diff','','Civ1_FF'};
     194         role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','errorflag'}; 
     195         units={'pixel','pixel','pixel','pixel','pixel','pixel',[]};
     196    case 'filter1'
     197        var={'X_tps','Y_tps','Z_tps','U_tps','V_tps','W_tps','X_SubRange','Y_SubRange','NbSites';...
     198            'Civ1_X_tps','Civ1_Y_tps','','Civ1_U_tps','Civ1_V_tps','','Civ1_X_SubRange','Civ1_Y_SubRange','Civ1_NbSites'};
     199         role={'','','','','','','','',''}; 
     200         units={'pixel','pixel','pixel','pixel','pixel','pixel','pixel','pixel',''};
     201    otherwise % if VelType=[]
     202        if testpatch
     203           var={'X_tps','Y_tps','Z_tps','U_tps','V_tps','W_tps','X_SubRange','Y_SubRange','NbSites';...
     204               'Civ2_X_tps','Civ2_Y_tps','','Civ2_U_tps','Civ2_V_tps','','Civ2_X_SubRange','Civ2_Y_SubRange','Civ2_NbSites';...
     205               'Civ1_X_tps','Civ1_Y_tps','','Civ1_U_tps','Civ1_V_tps','','Civ1_X_SubRange','Civ1_Y_SubRange','Civ1_NbSites'};
     206            role={'','','','','','','',''}; 
     207            units={'pixel','pixel','pixel','pixel','pixel','pixel','pixel','pixel'};
     208        else
     209             var={'X','Y','Z','U','V','W','C','F','FF';...
     210              'Civ2_X','Civ2_Y','Civ2_Z','Civ2_U','Civ2_V','Civ2_W','Civ2_C','Civ2_F','Civ2_FF';...
     211              'Civ1_X','Civ1_Y','Civ1_Z','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_F','Civ1_FF'};
     212            role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag'};
     213            units={'pixel','pixel','pixel','pixel','pixel','pixel',[],[],[]};
     214        end   
     215end
    195216if isempty(vel_type) || isequal(vel_type,'*') %undefined velocity type (civ1,civ2...)
    196217    if testder
     
    219240vel_type_out=vel_type_out';
    220241
    221 %% determine names of netcdf variables to read
    222 var={'X','Y','Z','U','V','W','C','F','FF'};
    223 role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag'};
    224 units={'pixel','pixel','pixel','pixel','pixel','pixel',[],[],[]};
    225 if testder
    226     var=[var {'DjUi(:,1,1)','DjUi(:,1,2)','DjUi(:,2,1)','DjUi(:,2,2)'}];
    227     role=[role {'tensor','tensor','tensor','tensor'}];
    228     units=[units {'pixel','pixel','pixel','pixel'}];
    229 end
    230 for ilist=1:length(vel_type_out)
    231     var=[var;varname1(vel_type_out{ilist},FieldNames)];
    232 end
    233 
    234 %------------------------------------------------------------------------ 
    235 % TAKEN FROM read_civxdata NOT USED
    236 %--- determine  var names to read
    237 
    238 function varin=varname1(vel_type,FieldNames)
    239 %------------------------------------------------------------------------
    240 testder=0;
    241 C1='';
    242 C2='';
    243 for ilist=1:length(FieldNames)
    244     if ~isempty(FieldNames{ilist})
    245     switch FieldNames{ilist}
    246         case 'ima_cor' %image correlation corresponding to a vel vector
    247             C1='vec_C';
    248             C2='vec2_C';
    249         case 'error'
    250             C1='vec_E';
    251             C2='vec2_E';
    252         case {'vort','div','strain'}
    253             testder=1;
    254     end
    255     end
    256 end     
    257 switch vel_type
    258     case 'civ1'
    259         varin={'Civ1_X','Civ1_Y','Civ1_Z','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_F','Civ1_FF'};
    260     case 'interp1'
    261         varin={'','','','','','','','',''};
    262     case 'filter1'
    263         varin={'Civ1_X','Civ1_Y','','Civ1_U_Spline','Civ1_V_Spline','','','',''};
    264     case 'civ2'
    265         varin={'vec2_X','vec2_Y','vec2_Z','vec2_U','vec2_V','vec2_W',C2,'vec2_F','vec2_FixFlag'};
    266     case 'interp2'
    267         varin={'vec2_patch_X','vec2_patch_Y','vec2_patch_Z','vec2_patch0_U','vec2_patch0_V','vec2_patch0_W','','',''};
    268     case 'filter2'
    269         varin={'vec2_patch_X','vec2_patch_Y','vec2_patch_Z','vec2_patch_U','vec2_patch_V','vec2_patch0_W','','',''};
    270 end
    271 if testder
    272      switch vel_type
    273         case 'filter1'
    274             varin=[varin {'vec_patch_DUDX','vec_patch_DVDX','vec_patch_DUDY','vec_patch_DVDY'}];
    275         case 'filter2'
    276             varin=[varin {'vec2_patch_DUDX','vec2_patch_DVDX','vec2_patch_DUDY','vec2_patch_DVDY'}];
    277     end   
    278 end
     242
     243
     244
     245
     246
  • trunk/src/series.m

    r244 r246  
    129129%     set(handles.VelTypeMenu,'String',param.VelTypeMenu)
    130130% end 
    131 set(hObject,'WindowButtonUpFcn',{@mouse_up_gui,handles})
     131set(hObject,'WindowButtonUpFcn',{'mouse_up_gui',handles})
    132132NomType_Callback(hObject, eventdata, handles)
    133133
     
    18581858end   
    18591859
    1860 %-----------------------------
    1861 function mouse_up_gui(hObject,eventdata,handles)
    1862 if isequal(get(hObject,'SelectionType'),'alt')
    1863     set(hObject,'Units','pixels')
    1864     series_pos=get(hObject,'Position');%position of the GUI series (in pixels)
    1865     set(hObject,'Units','normalized')
    1866     xy_fig=get(hObject,'CurrentPoint');% current point of the current figure (gcbo)
    1867     hchild=get(hObject,'Children');%handles of all objects in the current figure
    1868     %% loop on all the objects in the current figure (selected by the last mouse click)
    1869 
    1870     for ichild=1:length(hchild)
    1871         obj_pos=get(hchild(ichild),'Position');%position of the object       
    1872         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);         
    1873             htype=get(hchild(ichild),'Type');%type of object child of the current figure
    1874             %if the mouse is over an axis, look at the data
    1875             if isequal(htype,'uicontrol') && isequal(get(hchild(ichild),'Visible'),'on')
    1876                 msg_pos(1:2)=series_pos(1:2)+obj_pos(1:2).*series_pos(3:4);
    1877                 msgbox_uvmat(['uicontrol: ' get(hchild(ichild),'Tag')],'',get(hchild(ichild),'String'),msg_pos)
    1878                 break
    1879             end
    1880         end
    1881     end
    1882     set(hObject,'Units','pixels')
    1883 end
    18841860
    18851861%%%%%%%%%%%%%
  • trunk/src/struct2nc.m

    r189 r246  
    3535function errormsg=struct2nc(flname,Data)
    3636if ~ischar(flname)
    37     errormsg='no name input for the netcf file';
     37    errormsg='invalid input for the netcf file name';
    3838    return
    3939end
  • trunk/src/uvmat.m

    r243 r246  
    21502150        errormsg=['error in reading ' filename ': ' errormsg];
    21512151        return
    2152     end
    2153        
     2152    end       
    21542153    if isfield(ParamOut,'Npx')&& isfield(ParamOut,'Npy')
    21552154        set(handles.npx,'String',num2str(ParamOut.Npx));% display image size on the interface
     
    22962295    index_menu=strcmp(ParamOut.VelType,menu);
    22972296    set(handles.VelType,'Value',find(index_menu,1))
     2297    if ~get(handles.SubField,'value')
    22982298    set(handles.VelType,'String',menu)
    2299 %     set(handles.VelType_1,'Value',1)
    2300 %     set(handles.VelType_1,'String',[{''};menu])
     2299     set(handles.VelType_1,'Value',1)
     2300     set(handles.VelType_1,'String',[{''};menu])
     2301    end
    23012302end
    23022303field_index=strcmp(ParamOut.FieldName,ParamOut.FieldList);
     
    24172418    Field{1}=calc_field([{ParamOut.FieldName} {ParamOut.ColorVar}],Field{1});
    24182419end
    2419 if length(Field)==2 && ~test_keepdata_1 && isequal(FileType_1,'netcdf') && ~isequal(ParamOut_1.FieldName,'get_field...')%&&~isempty(FieldName_1)
     2420if numel(Field)==2 && ~test_keepdata_1 && isequal(FileType_1,'netcdf') && ~isequal(ParamOut_1.FieldName,'get_field...')%&&~isempty(FieldName_1)
    24202421    Field{2}=calc_field([{ParamOut_1.FieldName} {ParamOut_1.ColorVar}],Field{2});
    24212422end
     
    25102511            nbpoints_z=UvData.Field.DimValue(DimIndex(1));
    25112512            DZ=(ZMax-ZMin)/(nbpoints_z-1);
    2512             UvData.Field.Mesh=sqrt(DX*DY*DZ);
     2513            UvData.Field.Mesh=(DX*DY*DZ)^(1/3);
    25132514            UvData.Field.ZMax=ZMax;
    25142515            UvData.Field.ZMin=ZMin;
     
    31263127end
    31273128FileExt_1=get(handles.FileExt_1,'String');
    3128 if isequal(FileExt_1,'"'),FileExt_1=get(handles.FileExt,'String'); end;
     3129if isequal(get(handles.FileExt_1,'Visible'),'off') || isequal(FileExt_1,'"')
     3130    FileExt_1=get(handles.FileExt,'String');%read FileExt by default
     3131end
    31293132FileName_1=[FileName_1 FileIndices_1 FileExt_1];
    31303133
     
    32723275
    32733276%read the rootfile input display
    3274 [FileName,RootPath,FileBase,FileIndices,FileExt_prev]=read_file_boxes_1(handles);
    3275 [P,F,str1,str2,str_a,str_b,E,NomType]=name2display(['xxx' get(handles.FileIndex,'String') FileExt_prev]);
    3276 if isempty(FileExt_prev)|| strcmp(FileExt_prev,'')
    3277     FileExt_1=get(handles.FileExt,'String');
    3278 else
    3279     FileExt_1=FileExt_prev;
    3280 end
     3277[FileName,RootPath,FileBase,FileIndices,FileExt_1]=read_file_boxes_1(handles);
     3278[P,F,str1,str2,str_a,str_b,E,NomType]=name2display(['xxx' get(handles.FileIndex,'String') FileExt_1]);
     3279% if isempty(FileExt_prev)|| strcmp(FileExt_prev,'')
     3280%     FileExt_1=get(handles.FileExt,'String');
     3281% else
     3282%     FileExt_1=FileExt_prev;
     3283% end
    32813284NomType_1=get(handles.FileIndex_1,'UserData');
    32823285if isempty(NomType_1)|| strcmp(NomType_1,'')
     
    33503353else
    33513354    set(handles.SubDir_1,'Visible','on')
    3352     if ~isequal(FileExt_prev,'.nc') %find the new NomType if the previous display was not already a netcdf file
     3355    if ~isequal(FileExt_1,'.nc') %find the new NomType if the previous display was not already a netcdf file
    33533356%         veltype_handles=[handles.VelType_1 handles.interp1_1 handles.filter1_1 handles.civ2_1 handles.interp2_1 handles.filter2_1];
    33543357%         set_veltype_display(veltype_handles,6); % make all civ buttons visible
     
    34963499%---------------------------------------------
    34973500 
    3498 set(handles.FixVelType,'Value',1)
     3501set(handles.FixVelType,'Value',1)% the velocity type is now imposed by the GUI (not automatic)
    34993502%refresh field with a second filename=first fiel name
    35003503set(handles.run0,'BackgroundColor',[1 1 0])%paint the command button in yellow
    3501 drawnow
     3504drawnow   
    35023505filename=read_file_boxes(handles);
    3503 if get(handles.SubField,'Value')
    3504     filename_1=read_file_boxes_1(handles);
    3505 else
    3506     index=get(handles.VelType_1,'Value');
    3507     if index==1
     3506
     3507index=get(handles.VelType_1,'Value');
     3508if index==1
    35083509        filename_1='';% we plot the current field without the second field
    3509     else
    3510         filename_1=filename;
    3511     end
    3512 end
     3510        set(handles.SubField,'Value',0)
     3511        SubField_Callback(hObject, eventdata, handles)
     3512elseif get(handles.SubField,'Value')% if subfield is already 'on'
     3513    filename_1=read_file_boxes_1(handles); %read the current second field
     3514else
     3515     filename_1=filename;% we compare two fields in the same file
     3516     set(handles.SubField,'Value',1)
     3517end
     3518
    35133519num_i1=stra2num(get(handles.i1,'String'));
    35143520num_i2=stra2num(get(handles.i2,'String'));
Note: See TracChangeset for help on using the changeset viewer.