Changeset 363 for trunk/src/civ_matlab.m


Ignore:
Timestamp:
Jan 10, 2012, 9:13:31 AM (9 years ago)
Author:
sommeria
Message:

civ2 introduced in civ_matlab (still bugs)
bugs corrected in civ and fill_GUI, read_GUI
comments in get_field

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/civ_matlab.m

    r356 r363  
    7777    for ilist=1:length(list_param)
    7878        Civ1_param{ilist}=['Civ1_' list_param{ilist}];
    79         %         eval(['Data.Civ1_' list_param{ilist} '=Param.Civ1.' list_param{ilist} ';'])
    8079        Data.(['Civ1_' list_param{ilist}])=Param.Civ1.(list_param{ilist});
    8180    end
    8281    Data.ListGlobalAttribute=[Data.ListGlobalAttribute Civ1_param];% {'Civ1_Time','Civ1_Dt'}];
    83     %     if exist('ncfile','var')% TEST for use interactively with mouse_motion (no file created)
    8482    Data.ListVarName={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F'};%  cell array containing the names of the fields to record
    8583    Data.VarDimName={'NbVec1','NbVec1','NbVec1','NbVec1','NbVec1','NbVec1'};
     
    103101        CivFile=Param.Patch1.CivFile;
    104102    end
    105     Data=nc2struct(CivFile,'ListGlobalAttribute','absolut_time_T0') %look for the constant 'absolut_time_T0' to detect old civx data format
     103    Data=nc2struct(CivFile,'ListGlobalAttribute','absolut_time_T0'); %look for the constant 'absolut_time_T0' to detect old civx data format
    106104    if isfield(Data,'Txt')
    107105        errormsg=Data.Txt;
     
    112110        [Data,vardetect,ichoice]=nc2struct(CivFile);%read the variables in the netcdf file
    113111    else
    114         if isfield(Param,'Fix1')
    115             Data=nc2struct(CivFile,ListVarCiv1);%read civ1 data in the existing netcdf file
    116         else
    117             Data=nc2struct(CivFile,ListVarFix1);%read civ1 and fix1 data in the existing netcdf file
    118         end
     112        Data=nc2struct(CivFile);%read civ1 and fix1 data in the existing netcdf file
     113%         if isfield(Param,'Fix1')
     114%             Data=nc2struct(CivFile,ListVarCiv1);%read civ1 data in the existing netcdf file
     115%         else
     116%             Data=nc2struct(CivFile,ListVarFix1);%read civ1 and fix1 data in the existing netcdf file
     117%         end
    119118    end
    120119end
     
    153152    end
    154153    check_patch1=1;
    155     Param.Patch1
    156154    Data.ListGlobalAttribute=[Data.ListGlobalAttribute {'Patch1_Rho','Patch1_Threshold','Patch1_SubDomain'}];
    157155    Data.Patch1_Rho=Param.Patch1.SmoothingParam;
     
    168166    if isfield(Data,'Civ1_FF')
    169167        ind_good=find(Data.Civ1_FF==0);
    170     else
    171        
     168    else
    172169        ind_good=1:numel(Data.Civ1_X);
    173170    end
     
    183180if isfield (Param,'Civ2')
    184181    par_civ2=Param.Civ2;
    185     if ~check_civ1 || ~strcmp(par_civ1.filename_ima_a,par_civ2.filename_ima_a)
    186             par_civ2.ImageA=imread(par_civ2.filename_ima_a);%read first image if not already done for civ1
    187     end
    188     if ~check_civ1|| ~strcmp(par_civ1.filename_ima_b,par_civ2.filename_ima_b)
    189             par_civ2.ImageB=imread(par_civ2.filename_ima_b);%read second image if not already done for civ1
    190     end
    191 %     stepx=str2double(par_civ2.dx);
    192 %     stepy=str2double(par_civ2.dy);
    193     ibx2=ceil(str2double(par_civ2.ibx)/2);
    194     iby2=ceil(str2double(par_civ2.iby)/2);
     182    if ~isfield (Param,'Civ1') || ~strcmp(Param.Civ1.ImageA,par_civ2.ImageA)
     183        par_civ2.ImageA=imread(par_civ2.ImageA);%read first image if not already done for civ1
     184    else
     185        par_civ2.ImageA=Param.Civ1.ImageA;
     186    end
     187    if ~isfield (Param,'Civ1') || ~strcmp(Param.Civ1.ImageB,par_civ2.ImageB)
     188        par_civ2.ImageB=imread(par_civ2.ImageB);%read second image if not already done for civ1
     189         else
     190        par_civ2.ImageB=Param.Civ1.ImageB;
     191    end
     192    ibx2=ceil(par_civ2.Bx/2);
     193    iby2=ceil(par_civ2.By/2);
    195194    isx2=ibx2+2;
    196195    isy2=iby2+2;
    197 
    198 
    199 %         Data.Civ1_U_Diff=zeros(size(Data.Civ1_X));
    200 %         Data.Civ1_V_Diff=zeros(size(Data.Civ1_X));
    201 %         if isfield(Data,'Civ1_FF')
    202 %             ind_good=find(Data.Civ1_FF==0);
    203 %         else
    204 %             ind_good=1:numel(Data.Civ1_X);
    205 %         end
    206 %             [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]=...
    207 %                                 patch(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);
    208 %         end
    209 %     shiftx=str2num(par_civ1.shiftx);
    210 %     shifty=str2num(par_civ1.shifty);
    211 % TO GET shift from par_civ2.filename_nc1
    212     % shiftx=velocity interpolated at position
    213     miniy=max(1+isy2+shifty,1+iby2);
    214     minix=max(1+isx2-shiftx,1+ibx2);
    215     maxiy=min(size(par_civ2.ImageA,1)-isy2+shifty,size(par_civ2.ImageA,1)-iby2);
    216     maxix=min(size(par_civ2.ImageA,2)-isx2-shiftx,size(par_civ2.ImageA,2)-ibx2);
     196    % shift from par_civ2.filename_nc1
     197    % shiftx=velocity interpolated at position
     198    miniy=max(1+isy2,1+iby2);
     199    minix=max(1+isx2,1+ibx2);
     200    maxiy=min(size(par_civ2.ImageA,1)-isy2,size(par_civ2.ImageA,1)-iby2);
     201    maxix=min(size(par_civ2.ImageA,2)-isx2,size(par_civ2.ImageA,2)-ibx2);
    217202    [GridX,GridY]=meshgrid(minix:par_civ2.Dx:maxix,miniy:par_civ2.Dy:maxiy);
    218     PointCoord(:,1)=reshape(GridX,[],1);
    219     PointCoord(:,2)=reshape(GridY,[],1);
    220 
    221     Guess = tps_eval(PointCoord,[Data.Civ1_X_tps Data.Civ1_Y_tps])
    222     Shiftx=Guess*Data.Civ1_U_tps;
    223     Shifty=Guess*Data.Civ1_V_tps;
    224    
    225     if ~isempty(par_civ2.maskname)&& ~strcmp(maskname,par_civ2.maskname)% mask exist, not already read in civ1
     203    GridX=reshape(GridX,[],1);
     204    GridY=reshape(GridY,[],1);
     205    Shiftx=zeros(size(GridX));% shift expected from civ1 data
     206    Shifty=zeros(size(GridX));
     207    nbval=zeros(size(GridX));
     208    if par_civ2.CheckDeformation
     209        DUDX=zeros(size(GridX));
     210        DUDY=zeros(size(GridX));
     211        DVDX=zeros(size(GridX));
     212        DVDY=zeros(size(GridX));
     213    end
     214    [NbSubDomain,xx]=size(Data.Civ1_X_SubRange);
     215    % get the guess from patch1
     216    for isub=1:NbSubDomain
     217        nbvec_sub=Data.Civ1_NbSites(isub);
     218        ind_sel=find(GridX>=Data.Civ1_X_SubRange(isub,1) & GridX<=Data.Civ1_X_SubRange(isub,2) & GridY>=Data.Civ1_Y_SubRange(isub,1) & GridY<=Data.Civ1_Y_SubRange(isub,2));
     219        epoints = [GridX(ind_sel) GridY(ind_sel)];% coordinates of interpolation sites
     220        ctrs=[Data.Civ1_X_tps(1:nbvec_sub,isub) Data.Civ1_Y_tps(1:nbvec_sub,isub)];%(=initial points) ctrs
     221        nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
     222        EM = tps_eval(epoints,ctrs);
     223        Shiftx(ind_sel)=Shiftx(ind_sel)+EM*Data.Civ1_U_tps(1:nbvec_sub+3,isub);
     224        Shifty(ind_sel)=Shifty(ind_sel)+EM*Data.Civ1_V_tps(1:nbvec_sub+3,isub);
     225        if par_civ2.CheckDeformation
     226            [EMDX,EMDY] = tps_eval_dxy(epoints,ctrs);%2D matrix of distances between extrapolation points epoints and spline centres (=site points) ctrs
     227            DUDX(ind_sel)=DUDX(ind_sel)+EMDX*Data.Civ1_U_tps(1:nbvec_sub+3,isub);
     228            DUDY(ind_sel)=DUDY(ind_sel)+EMDY*Data.Civ1_U_tps(1:nbvec_sub+3,isub);
     229            DVDX(ind_sel)=DVDX(ind_sel)+EMDX*Data.Civ1_V_tps(1:nbvec_sub+3,isub);
     230            DVDY(ind_sel)=DVDY(ind_sel)+EMDY*Data.Civ1_V_tps(1:nbvec_sub+3,isub);
     231        end
     232    end
     233    mask='';
     234    if par_civ2.CheckMask&&~isempty(par_civ2.maskname)&& ~strcmp(maskname,par_civ2.maskname)% mask exist, not already read in civ1
    226235        mask=imread(par_civ2.maskname);
    227236    end
     237    par_civ2.Searchx=2*isx2+1;
     238    par_civ2.Searchy=2*isy2+1;
     239    par_civ2.Shiftx=Shiftx./nbval;
     240    par_civ2.Shifty=Shifty./nbval;
     241    par_civ2.Grid=[GridX GridY];   
     242    if par_civ2.CheckDeformation
     243        DUDX=DUDX./nbval;
     244        DUDY=DUDY./nbval;
     245        DVDX=DVDX./nbval;
     246        DVDY=DVDY./nbval;
     247    end
    228248    % caluclate velocity data (y and v in indices, reverse to y component)
    229     [xtable ytable utable vtable ctable F] = civ (par_civ2.ImageA,par_civ1.ImageB,ibx2,iby2,isx2,isy2,shiftx,-shifty,PointCoord,str2num(par_civ1.rho),mask);
    230     list_param=(fieldnames(par_civ1))';
     249    [xtable ytable utable vtable ctable F] = civ (par_civ2);
     250    list_param=(fieldnames(Param.Civ2))';
    231251    list_remove={'pxcmx','pxcmy','npx','npy','gridflag','maskflag','term_a','term_b','T0'};
    232     index=zeros(size(list_param));
    233252    for ilist=1:length(list_remove)
    234253        index=strcmp(list_remove{ilist},list_param);
     
    238257    end
    239258    for ilist=1:length(list_param)
    240         Civ1_param{ilist}=['Civ1_' list_param{ilist}];
    241         eval(['Data.Civ1_' list_param{ilist} '=Param.Civ1.' list_param{ilist} ';'])
    242     end
    243     if isfield(Data,'Civ1_gridname') && strcmp(Data.Civ1_gridname(1:6),'noFile')
     259        Civ2_param{ilist}=['Civ2_' list_param{ilist}];
     260        eval(['Data.Civ2_' list_param{ilist} '=Param.Civ2.' list_param{ilist} ';'])
     261    end
     262    if isfield(Data,'Civ2_gridname') && strcmp(Data.Civ1_gridname(1:6),'noFile')
    244263        Data.Civ1_gridname='';
    245264    end
    246     if isfield(Data,'Civ1_maskname') && strcmp(Data.Civ1_maskname(1:6),'noFile')
    247         Data.Civ1_maskname='';
    248     end
    249     Data.ListGlobalAttribute=[Data.ListGlobalAttribute Civ1_param {'Civ1_Time','Civ1_Dt'}];
    250     Data.Civ1_Time=str2double(par_civ1.T0);
    251     Data.Civ1_Dt=str2double(par_civ1.Dt);
    252     Data.ListVarName={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F'};%  cell array containing the names of the fields to record
    253     Data.VarDimName={'NbVec1','NbVec1','NbVec1','NbVec1','NbVec1','NbVec1'};
     265    if isfield(Data,'Civ2_maskname') && strcmp(Data.Civ1_maskname(1:6),'noFile')
     266        Data.Civ2_maskname='';
     267    end
     268    Data.ListGlobalAttribute=[Data.ListGlobalAttribute Civ2_param {'Civ2_Time','Civ2_Dt'}];
     269    Data.Civ2_Time=str2double(par_civ2.Time);
     270    Data.Civ2_Dt=str2double(par_civ2.Dt);
     271    Data.ListVarName=[Data.ListVarName {'Civ2_X','Civ2_Y','Civ2_U','Civ2_V','Civ2_C','Civ2_F'}];%  cell array containing the names of the fields to record
     272    Data.VarDimName=[Data.VarDimName {'NbVec2','NbVec2','NbVec2','NbVec2','NbVec2','NbVec2'}];
    254273    Data.VarAttribute{1}.Role='coord_x';
    255274    Data.VarAttribute{2}.Role='coord_y';
     
    257276    Data.VarAttribute{4}.Role='vector_y';
    258277    Data.VarAttribute{5}.Role='warnflag';
    259     Data.Civ1_X=reshape(xtable,[],1);
    260     Data.Civ1_Y=reshape(size(par_civ2.ImageA,1)-ytable+1,[],1);
    261     Data.Civ1_U=reshape(utable,[],1);
    262     Data.Civ1_V=reshape(-vtable,[],1);
    263     Data.Civ1_C=reshape(ctable,[],1);
    264     Data.Civ1_F=reshape(F,[],1);
     278    Data.Civ2_X=reshape(xtable,[],1);
     279    Data.Civ2_Y=reshape(size(par_civ2.ImageA,1)-ytable+1,[],1);
     280    Data.Civ2_U=reshape(utable,[],1);
     281    Data.Civ2_V=reshape(-vtable,[],1);
     282    Data.Civ2_C=reshape(ctable,[],1);
     283    Data.Civ2_F=reshape(F,[],1);
    265284    Data.CivStage=Data.CivStage+1;
    266285end
     
    361380if isfield(par_civ,'Grid')
    362381    if ischar(par_civ.Grid)
    363         par_civ.Grid;
    364382        par_civ.Grid=dlmread(par_civ.Grid);
    365383        par_civ.Grid(1,:)=[];%the first line must be removed (heading in the grid file)
     
    370388    isx2=ceil(par_civ.Searchx/2);
    371389    isy2=ceil(par_civ.Searchy/2);
    372     shiftx=par_civ.Shiftx;
    373     shifty=-par_civ.Shifty;
    374390    miniy=max(1+isy2+shifty,1+iby2);
    375391    minix=max(1+isx2-shiftx,1+ibx2);
     
    380396    par_civ.Grid(:,2)=reshape(GridY,[],1);
    381397end
    382 
     398nbvec=size(par_civ.Grid,1);
     399if numel(shiftx)==1
     400    shiftx=round(shiftx)*ones(nbvec,1);
     401    shifty=round(shifty)*ones(nbvec,1);
     402end
    383403%% Default output
    384 nbvec=size(par_civ.Grid,1);
    385404xtable=par_civ.Grid(:,1);
    386405ytable=par_civ.Grid(:,2);
     
    462481    %     ytable(ivec)=jref;%default
    463482    if ~(checkmask && par_civ.Mask(jref,iref)<=20) %velocity not set to zero by the black mask
    464         if jref-iby2<1 || jref+iby2>par_civ.ImageHeight|| iref-ibx2<1 || iref+ibx2>par_civ.ImageWidth% we are outside the image
     483        if jref-iby2<1 || jref+iby2>par_civ.ImageHeight|| iref-ibx2<1 || iref+ibx2>par_civ.ImageWidth||...
     484              jref+shifty(ivec)-isy2<1||jref+shifty(ivec)+isy2>par_civ.ImageHeight|| iref+shiftx(ivec)-isx2<1 || iref+shiftx(ivec)+isx2>par_civ.ImageWidth  % we are outside the image
    465485            F(ivec)=3;
    466486        else
    467487            image1_crop=par_civ.ImageA(jref-iby2:jref+iby2,iref-ibx2:iref+ibx2);%extract a subimage (correlation box) from image A
    468             image2_crop=par_civ.ImageB(jref+shifty-isy2:jref+shifty+isy2,iref+shiftx-isx2:iref+shiftx+isx2);%extract a larger subimage (search box) from image B
     488            image2_crop=par_civ.ImageB(jref+shifty(ivec)-isy2:jref+shifty(ivec)+isy2,iref+shiftx(ivec)-isx2:iref+shiftx(ivec)+isx2);%extract a larger subimage (search box) from image B
    469489            image1_mean=mean(mean(image1_crop));
    470490            image2_mean=mean(mean(image2_crop));
     
    496516                        [vector,F(ivec)] = SUBPIX2DGAUSS (result_conv,x,y);
    497517                    end
    498                     utable(ivec)=vector(1)+shiftx;
    499                     vtable(ivec)=vector(2)+shifty;
     518                    utable(ivec)=vector(1)+shiftx(ivec);
     519                    vtable(ivec)=vector(2)+shifty(ivec);
    500520                    xtable(ivec)=iref+utable(ivec)/2;% convec flow (velocity taken at the point middle from imgae1 and 2)
    501521                    ytable(ivec)=jref+vtable(ivec)/2;
     
    509529                    ctable(ivec)=corrmax/sum_square;% correlation value
    510530                catch ME
    511                     %                     vector=[0 0]; %if something goes wrong with cross correlation.....
    512531                    F(ivec)=3;
    513532                end
     
    619638function FF=fix(Param,F,C,U,V,X,Y)
    620639FF=zeros(size(F));%default
    621 Param
    622640
    623641%criterium on warn flags
Note: See TracChangeset for help on using the changeset viewer.