Changeset 1195 for trunk/src


Ignore:
Timestamp:
Feb 26, 2026, 4:16:03 PM (5 days ago)
Author:
sommeria
Message:

format of PIV data made more compact/4

Location:
trunk/src
Files:
14 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/civ.m

    r1179 r1195  
    44%
    55% OUTPUT:
    6 % xtable: set of x coordinates
    7 % ytable: set of y coordinates
    8 % utable: set of u displacements (along x)
     6% xtable: set of x positions of the first image in integer image coordinates , the measurement position is xtable-0.5+utable/2
     7% ytable: set of y coordinates of the first image in integer image coordinates (starting in the image bottom, image index= npy-ytable+1 , the measurement position is ytable-0.5+vtable/2
     8% utable: set of u displacements (along x), in pixels
    99% vtable: set of v displacements (along y)
    1010% ctable: max image correlation for each vector
    … …  
    4444%% prepare measurement grid if not given as input
    4545if ~isfield(par_civ,'Grid')% grid points defining central positions of the sub-images in image A
    46     nbinterv_x=floor((npx_ima-1)/par_civ.Dx);
     46    nbinterv_x=floor((npx_ima-1)/par_civ.Dx);%expected number of intervals Dx
    4747    gridlength_x=nbinterv_x*par_civ.Dx;
    4848    minix=ceil((npx_ima-gridlength_x)/2);
    49     nbinterv_y=floor((npy_ima-1)/par_civ.Dy);
     49    nbinterv_y=floor((npy_ima-1)/par_civ.Dy);%expected number of intervals Dy
    5050    gridlength_y=nbinterv_y*par_civ.Dy;
    5151    miniy=ceil((npy_ima-gridlength_y)/2);
    52     [GridX,GridY]=meshgrid(minix:par_civ.Dx:npx_ima-1,miniy:par_civ.Dy:npy_ima-1);
     52    [GridX,GridY]=meshgrid(minix+0.5:par_civ.Dx:npx_ima-0.5,miniy+0.5:par_civ.Dy:npy_ima-0.5);% grid of initial positions in pixel coordinates (y upward)
    5353    par_civ.Grid(:,1)=reshape(GridX,[],1);
    5454    par_civ.Grid(:,2)=reshape(GridY,[],1);% increases with array index
    … …  
    6666
    6767shiftx=par_civ.SearchBoxShift(:,1);%use the input shift estimate, rounded to the next integer value
    68 shifty=-par_civ.SearchBoxShift(:,2);% sign minus because image j index increases when y decreases
     68shifty=par_civ.SearchBoxShift(:,2);%
    6969if numel(shiftx)==1% case of a unique shift for the whole field( civ1)
    7070    shiftx=shiftx*ones(nbvec,1);
    … …  
    7373
    7474%% shift the grid by half the expected displacement to get the velocity closer to the initial grid
    75 par_civ.Grid(:,1)=par_civ.Grid(:,1)-shiftx/2;
    76 par_civ.Grid(:,2)=par_civ.Grid(:,2)+shifty/2;
     75par_civ.Grid(:,1)=par_civ.Grid(:,1)-shiftx/2;% initial positions in image coordinates (y upward)
     76par_civ.Grid(:,2)=par_civ.Grid(:,2)-shifty/2;
    7777
    7878%% Array initialisation and default output  if par_civ.CorrSmooth=0 (just the grid calculated, no civ computation)
    79 xtable=round(par_civ.Grid(:,1)+0.5)-0.5;
    80 ytable=round(npy_ima-par_civ.Grid(:,2)+0.5)-0.5;% y index corresponding to the position in image coordinates
     79xtable=round(par_civ.Grid(:,1)+0.5);% tables of image indices corresponding to the input grid
     80ytable=round(npy_ima-par_civ.Grid(:,2)+0.5);% y image indices corresponding to the position in image coordinates
    8181shiftx=round(shiftx);
    8282shifty=round(shifty);
    83 utable=shiftx;%zeros(nbvec,1);
    84 vtable=shifty;%zeros(nbvec,1);
     83utable=shiftx;%zeros(nbvec,1);% expected displacements in pixel indices (y downward)
     84vtable=-shifty;%zeros(nbvec,1);
    8585ctable=zeros(nbvec,1);
    8686FF=zeros(nbvec,1);
    … …  
    128128if par_civ.CorrSmooth~=0 % par_civ.CorrSmooth=0 implies no civ computation (just input image and grid points given)
    129129    for ivec=1:nbvec
    130         iref=round(par_civ.Grid(ivec,1)+0.5);% xindex on the image A for the middle of the correlation box
    131         jref=round(npy_ima-par_civ.Grid(ivec,2)+0.5);%  j index  for the middle of the correlation box in the image A
     130%         iref=round(par_civ.Grid(ivec,1));% xindex on the image A for the middle of the correlation box
     131%         jref=round(npy_ima-par_civ.Grid(ivec,2));%  j index  for the middle of the correlation box in the image A
     132         iref=xtable(ivec);% xindex on the image A for the middle of the correlation box
     133         jref=ytable(ivec);%  j index  for the middle of the correlation box in the image A
     134
    132135        FF(ivec)=0;
    133136        ibx2=floor(CorrBoxSizeX(ivec)/2);
    … …  
    138141        subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage
    139142        subrange2_x=iref+shiftx(ivec)-isx2:iref+shiftx(ivec)+isx2;%x indices defining the second subimage
    140         subrange2_y=jref+shifty(ivec)-isy2:jref+shifty(ivec)+isy2;%y indices defining the second subimage
     143        subrange2_y=jref-shifty(ivec)-isy2:jref-shifty(ivec)+isy2;%y indices defining the second subimage
    141144        image1_crop=MinA*ones(numel(subrange1_y),numel(subrange1_x));% default value=min of image A
    142145        image2_crop=MinA*ones(numel(subrange2_y),numel(subrange2_x));% default value=min of image A
    … …  
    224227                        end
    225228                        utable(ivec)=vector(1)*mesh+shiftx(ivec);
    226                         vtable(ivec)=-(vector(2)*mesh+shifty(ivec));
    227                         xtable(ivec)=iref+utable(ivec)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2)
    228                         ytable(ivec)=jref+vtable(ivec)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge)
    229                         iref=round(xtable(ivec)+0.5);% nearest image index for the middle of the vector
    230                         jref=round(ytable(ivec)+0.5);
     229                        vtable(ivec)=-vector(2)*mesh+shifty(ivec);% vtable and shifty in image coordinates (opposite to pixel shift)
     230                      %  xtable(ivec)=iref+utable(ivec)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2)
     231                       % ytable(ivec)=jref+vtable(ivec)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge)
     232%                         iref=round(xtable(ivec)+0.5);% nearest image index for the middle of the vector
     233%                         jref=round(ytable(ivec)+0.5);
    231234                        % eliminate vectors located in the mask
    232                         if  checkmask && (iref<1 || jref<1 ||iref>npx_ima || jref>npy_ima ||( par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100))
    233                             utable(ivec)=0;
    234                             vtable(ivec)=0;
    235                             FF(ivec)=1;
    236                         end
     235%                         if  checkmask && (iref<1 || jref<1 ||iref>npx_ima || jref>npy_ima ||( par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100))
     236%                             utable(ivec)=0;
     237%                             vtable(ivec)=0;
     238%                             FF(ivec)=1;
     239%                         end
    237240                        ctable(ivec)=corrmax/sum_square;% correlation value
    238241                    catch ME

  • trunk/src/filter_tps.m

    r1163 r1195  
    170170    end
    171171end
     172NbCentre=uint16(NbCentre);
    172173
    173174%% remove empty subdomains
    … …  
    188189%eliminate the vectors with diff>threshold not yet eliminated
    189190if exist('Threshold','var')&&~isempty(Threshold)
    190 UDiff=U_smooth-U;% difference between interpolated U component and initial value
    191 VDiff=V_smooth-V;% difference between interpolated V component and initial value
    192 NormDiff=UDiff.*UDiff+VDiff.*VDiff;% Square of difference norm
    193 FF(NormDiff>Threshold)=true;%put FF value to 1 to identify the criterium of elimmination
     191    UDiff=U_smooth-U;% difference between interpolated U component and initial value
     192    VDiff=V_smooth-V;% difference between interpolated V component and initial value
     193    NormDiff=UDiff.*UDiff+VDiff.*VDiff;% Square of difference norm
     194    FF(NormDiff>Threshold)=true;%put FF value to 1 to identify the criterium of elimmination
    194195end
    195196

  • trunk/src/get_field.m

    r1194 r1195  
    210210if isfield(ParamIn,'SeriesInput') && ParamIn.SeriesInput% case of call by series
    211211    set(handles.FieldOption,'value',1)
    212     if isfield(Field,'Conventions')&& strcmp(Field.Conventions,'uvmat/civdata')
     212    if isfield(Field,'Conventions')&& ~isempty(regexp(Field.Conventions,'^uvmat/civdata', 'once'))
    213213    set(handles.FieldOption,'String',{'scalar';'vectors';'civdata...'})
    214214    else
    … …  
    217217    checkseries=1;
    218218    set(handles.scalar,'Max',2)
    219 elseif isfield(Field,'Conventions')&& strcmp(Field.Conventions,'uvmat/civdata')
     219elseif isfield(Field,'Conventions')&& ~isempty(regexp(Field.Conventions,'^uvmat/civdata', 'once'))
    220220    set(handles.FieldOption,'String',{'1D plot';'scalar';'vectors';'civdata...'})% provides the possibility to come back to civdata
    221221    set(handles.scalar,'Max',1)

  • trunk/src/get_file_info.m

    r1191 r1195  
    151151            return
    152152        end
    153         XmlFile=fullfile(RootPath,[SubDir '.xml']);
    154         CheckWriteImaDoc=true;
     153%         XmlFile=fullfile(RootPath,[SubDir '.xml']);
     154%         CheckWriteImaDoc=true;
    155155%         if exist(XmlFile,'file')
    156156%             [XmlData,~,errormsg]=xml2struct(XmlFile);
    … …  
    207207                        [Data,tild,tild,errormsg]=nc2struct(fileinput,[]);
    208208                        if isempty(errormsg)
    209                             if isfield(Data,'Conventions') && strcmp(Data.Conventions,'uvmat/civdata')
     209                            if isfield(Data,'Conventions') && ~isempty(find(strcmp(Data.Conventions,{'uvmat/civdata','uvmat/civdata/compress'}), 1))
    210210                                FileInfo.FileType='civdata'; % test for civ velocity fields
    211211                                FileInfo.CivStage=Data.CivStage;

  • trunk/src/mouse_motion.m

    r1169 r1195  
    167167                                if isfield (CellInfo{icell},'XName')
    168168                                    XName=CellInfo{icell}.XName;
    169                                     text_displ_2=[XName '=' num2str(Field.(XName)(ivec),4)];
     169                                    text_displ_2=[XName '=' num2str(Field.(XName)(ivec),5)];
    170170                                end
    171171                                if isfield (CellInfo{icell},'YName')
    172172                                    YName=CellInfo{icell}.YName;
    173                                     text_displ_2=[text_displ_2 ',' YName '=' num2str(Field.(YName)(ivec),4)];
     173                                    text_displ_2=[text_displ_2 ',' YName '=' num2str(Field.(YName)(ivec),5)];
    174174                                end
    175175                            for ivar=1:numel(CellInfo{icell}.VarIndex)

  • trunk/src/nc2struct.m

    r1194 r1195  
    199199    end
    200200    ListVarName=ListVarName(:,logical(check_keep));
    201     if size(ListVarName,1)>1 %multiple choice of variable ranked by order of priority
    202         for iline=1:size(ListVarName,1)
    203             search_index=find(strcmp(ListVarName{iline,1},ListVarNameNetcdf),1);%look for the first variable name in the list of NetCDF variables
    204             if ~isempty(search_index)
    205                 break % go to the next line
    206             end
    207         end
    208         ichoice=iline-1;%selected line number in the list of input names of variables
    209     else
    210         iline=1;
    211     end
     201    iline=size(ListVarName,1);   
    212202    %ListVarName=ListVarName(iline,:);% select the appropriate option for input variable (lin ein the input name matrix)
    213203    if CheckTimeVar
    … …  
    274264        VarName=Data.ListVarName{ivar};
    275265        VarName=regexprep(VarName,'-','_'); %suppress '-' if it exists in the NetCDF variable name (leads to errors in matlab)
    276         %             CheckSub=0;
    277266        if input_index==4% if a dimension is selected as time
    278267            ind_vec=zeros(1,numel(var_dim{ivar}));% vector with zeros corresponding to al the dimensions of the variable VarName
    … …  
    290279                end
    291280            end
    292 %            tstart = tic;
    293           % disp(VarName)
    294281            Data.(VarName)=netcdf.getVar(nc,var_index(ivar)-1,ind_vec,ind_size); %read the variable data
    295 %              telapsed = toc(tstart)
    296282            Data.(VarName)=squeeze(Data.(VarName));%remove singeton dimension
    297 
    298283        else
     284                       disp(VarName)
     285           xtype(var_index(ivar))
    299286            Data.(VarName)=netcdf.getVar(nc,var_index(ivar)-1); %read the whole variable data
    300         end
    301         if xtype(var_index(ivar))==5% indicate single precision variable
    302             Data.(VarName)=double(Data.(VarName)); %transform to double for single pecision
     287        end       
     288        Data.(VarName)=double(Data.(VarName)); %transform all variables to double  pecision
     289        if isfield(Data,'VarAttribute') && numel(Data.VarAttribute)>=ivar && isfield(Data.VarAttribute{ivar},'scale_factor')
     290            Data.(VarName)=Data.VarAttribute{ivar}.scale_factor *Data.(VarName);
    303291        end
    304292    end

  • trunk/src/plot_field.m

    r1187 r1195  
    10051005                    A=uint8(255*double(A)/double(MaxA));
    10061006                end
    1007             otherwise
    1008                 colormap(ColorMap);
     1007%             otherwise
     1008%                 colormap(ColorMap);
    10091009        end
    10101010 

  • trunk/src/read_civdata.m

    r1169 r1195  
    9393
    9494%% reading data
    95 [Data,tild,tild,errormsg]=nc2struct(FileName,'ListGlobalAttribute','Conventions','CivStage');% read the global attributes to get Data.CivStage
     95[Data,~,~,errormsg]=nc2struct(FileName,'ListGlobalAttribute','Conventions','CivStage');% read the global attributes to get Data.CivStage
    9696if ~isempty(errormsg)
    9797     errormsg=['read_civdata: ' errormsg];
    9898    return
    9999end
     100
    100101% set the list of variables to read and their role
    101 [varlist,role,VelTypeOut]=varcivx_generator(ProjModeRequest,VelType,Data.CivStage);
     102CheckCompress=strcmp(Data.Conventions,'uvmat/civdata/compress');
     103[varlist,role,VelTypeOut]=varcivx_generator(ProjModeRequest,VelType,Data.CivStage,CheckCompress);
    102104if isempty(varlist)
    103105    errormsg=['read_civdata: unknow velocity type ' VelType];
    … …  
    105107else
    106108    if strcmp(Data.Conventions,'uvmat/civdata_3D')
    107         [Field,vardetect]=nc2struct(FileName,'TimeDimName','npz',frame_index,varlist);%read the variables in the netcdf file
    108     else   
    109     [Field,vardetect]=nc2struct(FileName,varlist);%read the variables in the netcdf file
    110     end
    111 end
    112 if isfield(Field,'Txt')
    113     errormsg=Field.Txt;
     109        [Field,vardetect,~,errormsg]=nc2struct(FileName,'TimeDimName','npz',frame_index,varlist);%read the variables in the netcdf file
     110    else
     111        [Field,vardetect,~,errormsg]=nc2struct(FileName,varlist);%read the variables in the netcdf file
     112    end
     113end
     114if ~isempty(errormsg)
     115     errormsg=['read_civdata: ' errormsg];
    114116    return
    115117end
    … …  
    117119     errormsg=[ 'requested field not available in ' FileName '/' VelType ': need to run patch'];
    118120     return
     121end
     122if strcmp(Data.Conventions,'uvmat/civdata/compress')
     123    Field.X=Field.X-0.5+Field.U/2;% shift to the convected position
     124    Field.Y=Field.Y-0.5+Field.V/2;
    119125end
    120126switch VelTypeOut
    … …  
    192198%            if vel_type=[] or'*', a  priority choice is done, civ2 considered better than civ1 )
    193199
    194 function [var,role,vel_type_out,errormsg]=varcivx_generator(ProjModeRequest,vel_type,CivStage)
     200function [var,role,vel_type_out,errormsg]=varcivx_generator(ProjModeRequest,vel_type,CivStage,CheckCompress)
    195201
    196202%% default input values
    … …  
    223229end
    224230
    225 var={};
    226 switch vel_type
    227     case{'civ1','civ2'}
    228         varout={'X','Y','Z','U','V','W','C','FF'};
    229         varin= {'Civ1_X','Civ1_Y','Civ1_Z','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_FF'};
    230          role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','errorflag'}; 
    231     case{'filter1','filter2'} 
    232         varout={'X','Y','Z','U','V','W','C','FF','Coord_tps','U_tps','V_tps','W_tps','SubRange','NbCentre','NbCentre','NbCentre'};
    233         varin={'Civ1_X','Civ1_Y','Civ1_Z','Civ1_U_smooth','Civ1_V_smooth','Civ1_W','Civ1_C','Civ1_FF',...
    234             'Civ1_Coord_tps','Civ1_U_tps','Civ1_V_tps','Civ1_W_tps','Civ1_SubRange','Civ1_NbCentre','Civ1_NbCentres','Civ1_NbSites'};
    235         role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','errorflag','coord_tps','vector_x_tps',...
    236             'vector_y_tps','vector_z_tps','ancillary','ancillary','ancillary','ancillary'};
    237           %  rmq: NbCentres and NbSites obsolete replaced by NbCentre, kept for consistency with previous data
    238 end
    239 switch vel_type
    240     case {'civ2','filter2'}
    241         varin=regexprep(varin,'1','2');
    242     % case {'civ3','filter3'}
    243     %     varin=regexprep(varin,'1','3');
     231if CheckCompress
     232    switch vel_type
     233        case{'civ1','civ2'}
     234            varout={'X','Y','Z','U','V','W','C','FF'};
     235            varin= {'X','Y','Z','U','V','W','C','FF'};
     236            role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','errorflag'};
     237        case{'filter1','filter2'}
     238            varout={'X','Y','Z','U','V','W','C','FF'};
     239            varin={'X','Y','Z','U_smooth','V_smooth','W_smooth','C','FF'};
     240            role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','errorflag'};
     241            %  rmq: NbCentres and NbSites obsolete replaced by NbCentre, kept for consistency with previous data
     242    end
     243else
     244    switch vel_type
     245        case{'civ1','civ2'}
     246            varout={'X','Y','Z','U','V','W','C','FF'};
     247            varin= {'Civ1_X','Civ1_Y','Civ1_Z','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_FF'};
     248            role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','errorflag'};
     249        case{'filter1','filter2'}
     250            varout={'X','Y','Z','U','V','W','C','FF','Coord_tps','U_tps','V_tps','W_tps','SubRange','NbCentre','NbCentre','NbCentre'};
     251            varin={'Civ1_X','Civ1_Y','Civ1_Z','Civ1_U_smooth','Civ1_V_smooth','Civ1_W','Civ1_C','Civ1_FF',...
     252                'Civ1_Coord_tps','Civ1_U_tps','Civ1_V_tps','Civ1_W_tps','Civ1_SubRange','Civ1_NbCentre','Civ1_NbCentres','Civ1_NbSites'};
     253            role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','errorflag','coord_tps','vector_x_tps',...
     254                'vector_y_tps','vector_z_tps','ancillary','ancillary','ancillary','ancillary'};
     255            %  rmq: NbCentres and NbSites obsolete replaced by NbCentre, kept for consistency with previous data
     256    end
     257    switch vel_type
     258        case {'civ2','filter2'}
     259            varin=regexprep(varin,'1','2');
     260    end
    244261end
    245262var=[varout;varin];

  • trunk/src/read_civxdata.m

    r1127 r1195  
    1 %'read_civxdata': reads civx data from netcdf files
     1%'read_civxdata': reads civx data from netcdf files (OBSOLETE)
    22%------------------------------------------------------------------
    33%

  • trunk/src/read_field.m

    r1189 r1195  
    8181%% distingush different input file types
    8282switch FileType
    83     case {'civdata','civdata_3D'}% new format for civ results
     83    case {'civdata','civdata_3D'}% format for civ results
    8484        [Field,ParamOut.VelType,errormsg]=read_civdata(FileName,InputField,ParamIn.VelType,frame_index);
    8585        if ~isempty(errormsg),errormsg=['read_civdata / ' errormsg];return,end

  • trunk/src/series.m

    r1194 r1195  
    12521252
    12531253%% display the set of existing files as an image with black bands for gaps showing gaps in the series
    1254 set(handles.FileStatus,'Units','pixels')
    1255 Position=get(handles.FileStatus,'Position');
    1256 set(handles.FileStatus,'Units','normalized')
     1254
     1255%set(handles.FileStatus,'Units','normalized')
    12571256nbview=numel(SeriesData.i1_series);
    12581257i_max=cell(1,nbview);
    … …  
    12671266        MaxIndex_i(iline)=find(i_max{iline}, 1, 'last' )-1; % max ref index i
    12681267        MinIndex_i(iline)=find(i_max{iline}, 1 )-1; % min ref index i
    1269          missing_indices{iline}= find(i_max{iline}(2:end)==0);         
     1268         exist_indices{iline}= find(i_max{iline}(2:end)~=0);
     1269         index_series=i_max{iline}( exist_indices{iline});
    12701270    end
    12711271end
    … …  
    12731273MaxIndex_i=max(MaxIndex_i);
    12741274range_index=MaxIndex_i-MinIndex_i+1;
     1275
     1276set(handles.FileStatus,'Units','pixels')
     1277Position=get(handles.FileStatus,'Position');
    12751278range_y=max(1,floor(Position(4)/nbview));
    12761279npx=floor(Position(3));%length of the bar image FileStatus in pixels
    … …  
    12781281%file_indices=MinIndex_i+floor(((0.5:npx-0.5)/npx)*range_index)+1;
    12791282CData=ones(nbview*range_y,npx); % initiate the image representing the existing files
    1280 LineData=ones(1,npx);
    1281 for iline=1:nbview
    1282     ind_y=1+(iline-1)*range_y:iline*range_y;
    1283     missing_pixels=floor((missing_indices{iline}-MinIndex_i+1)*npx/range_index)+1;
    1284     LineData(missing_pixels)=0;
    1285 %     LineData=zeros(size(file_indices));
    1286 %     file_select=file_indices(file_indices<=numel(i_max{iline}));
    1287 %     ind_select=file_indices<=numel(i_max{iline});
    1288 %     LineData(ind_select)=i_max{iline}(file_select)~=0;
    1289     CData(ind_y,:)=ones(numel(ind_y),1)*LineData;%create an image band with width numel(ind_y)
     1283if MaxIndex_i>MinIndex_i
     1284    LineData=ones(1,npx);
     1285    for iline=1:nbview
     1286        ind_y=1+(iline-1)*range_y:iline*range_y;
     1287%         missing_pixels=floor((missing_indices{iline}-MinIndex_i+1)*npx/range_index)+1;
     1288%         LineData(missing_pixels)=0;
     1289        %     LineData=zeros(size(file_indices));
     1290        %     file_select=file_indices(file_indices<=numel(i_max{iline}));
     1291        %     ind_select=file_indices<=numel(i_max{iline});
     1292        %     LineData(ind_select)=i_max{iline}(file_select)~=0;
     1293        CData(ind_y,:)=ones(numel(ind_y),1)*LineData;%create an image band with width numel(ind_y)
     1294    end
    12901295end
    12911296CData=cat(3,zeros(size(CData)),CData,zeros(size(CData))); % make color images r=0,g,b=0
    12921297set(handles.FileStatus,'CData',CData);
    1293 
     1298set(handles.FileStatus,'Units','normalized')
    12941299%-----------------------------------------------------------guide -------------
    12951300%------------------------------------------------------------------------

  • trunk/src/series/civ_series.m

    r1188 r1195  
    9696    return
    9797end
    98 %iview_A=0;%default values
    99 NbField=1;
    100 RUNHandle=[];
    101 % CheckInputFile=isfield(Param,'InputTable');%= 1 in test use for TestCiv (no nc file involved)
    102 % CheckOutputFile=isfield(Param,'OutputSubDir');%= 1 in test use for TestPatch (no nc file produced)
    103 
    104 %% input files and indexing
     98
     99inv_scale_factor=100; % scale factor of displacements for uin16 records in netcdf files (dx expressed in pixels)
     100%inv_scale_factor=[];% no scale factor, displacements written as single precision real
     101
     102%% input files and indexing
    105103hseries=findobj(allchild(0),'Tag','series');
    106104RUNHandle=findobj(hseries,'Tag','RUN');%handle of RUN button in GUI series
    … …  
    188186            j2_series_Civ1=ones(size(i1_series_Civ1));
    189187        end
     188        Check_j_Civ1=false;
    190189    else % movie for each burst or volume (index j)
    191190        FrameIndex_A_Civ1=j1_series_Civ1;
    192191        FrameIndex_B_Civ1=j2_series_Civ1;
     192        Check_j_Civ1=true;
    193193    end
    194194    if isempty(PairCiv2)
    … …  
    203203                j2_series_Civ2=ones(size(i1_series_Civ2));
    204204            end
     205            Check_j_Civ2=false;
    205206        else
    206207            FrameIndex_A_Civ2=j1_series_Civ2;
    207208            FrameIndex_B_Civ2=j2_series_Civ2;
     209            Check_j_Civ2=true;
    208210        end
    209211    end
    … …  
    226228OutputDir=[Param.OutputSubDir Param.OutputDirExt];
    227229ListGlobalAttribute={'Conventions','Program','CivStage'};
    228 Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes
     230Data.Conventions='uvmat/civdata/compress';% states the conventions used for the description of field variables and attributes
    229231Data.Program='civ_series';
    230232if isfield(Param,'UvmatRevision')
    … …  
    282284    OutputPath=fullfile(Param.OutputPath,Param.Experiment,Param.Device);
    283285    if CheckRelabel
    284          RootFileOut=index2filename(XmlData.FileSeries,1,1,MaxIndex_j);
     286        RootFileOut=index2filename(XmlData.FileSeries,1,1,MaxIndex_j);
    285287    else
    286288        RootFileOut=RootFile_A;
    … …  
    322324    ImageName_A='';ImageName_B='';%default
    323325    VideoObject_A=[];VideoObject_B=[];
     326    Civ_FF=[];%default
    324327   
    325328    %% Civ1
    … …  
    388391            end
    389392            [par_civ1.ImageB,VideoObject_B] = read_image(ImageName_B,FileType_B,VideoObject_B,FrameIndex_B);
    390 
     393           
    391394        catch ME % display errors in reading input images
    392395            if ~isempty(ME.message)
    … …  
    395398            end
    396399        end
    397 
    398  % case of background image to subtract
     400       
     401        % case of background image to subtract
    399402        if par_civ1.CheckBackground &&~isempty(par_civ1.Background)
    400403            [RootPath_background,SubDir_background,RootFile_background,~,~,~,~,Ext_background]=fileparts_uvmat(Param.ActionInput.Civ1.Background);
    401             j1=1;
    402             if ~isempty(j1_series_Civ1)
    403                 j1=j1_series_Civ1(ifield);
    404             end
    405404            if ~isempty(i2_series_Civ1)% case of volume,backgrounds act on different j levels
    406                 backgroundname=fullfile_uvmat(RootPath_background,SubDir_background,RootFile_background,Ext_background,'_1',j1);
     405                backgroundname=fullfile_uvmat(RootPath_background,SubDir_background,RootFile_background,Ext_background,'_1',j1_series_Civ1(ifield));
    407406            elseif isfield(par_civ1,'NbSlice')
    408407                i1_background=mod(i1-1,par_civ1.NbSlice)+1;
    … …  
    436435            par_civ1.ImageB=par_civ1.ImageB-par_civ1.Background;
    437436        end
    438 
    439  % case of image luminosity rescaling
    440       if par_civ1.CheckRescale &&~isempty(par_civ1.Maxtanh)
    441                par_civ1.ImageA =par_civ1.Maxtanh*tanh(double(par_civ1.ImageA)/par_civ1.Maxtanh);
    442                par_civ1.ImageB=par_civ1.Maxtanh*tanh(double(par_civ1.ImageB)/par_civ1.Maxtanh);
     437       
     438        % case of image luminosity rescaling
     439        if par_civ1.CheckRescale &&~isempty(par_civ1.Maxtanh)
     440            par_civ1.ImageA =par_civ1.Maxtanh*tanh(double(par_civ1.ImageA)/par_civ1.Maxtanh);
     441            par_civ1.ImageB=par_civ1.Maxtanh*tanh(double(par_civ1.ImageB)/par_civ1.Maxtanh);
    443442        end
    444443       
    … …  
    457456            i2=i2_series_Civ1(ifield);
    458457        end
    459         j1=1;
    460         if ~isempty(j1_series_Civ1)
    461             j1=j1_series_Civ1(ifield);
    462         end
     458       
     459        j1=j1_series_Civ1(ifield);
     460       
    463461        j2=j1;
    464         if ~isempty(j2_series_Civ1)
     462        if Check_j_Civ1
    465463            j2=j2_series_Civ1(ifield);
    466464        end
    … …  
    479477       
    480478        % set the list of variables
    481         Data.ListVarName={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_FF'};%  cell array containing the names of the fields to record
    482         Data.VarDimName={'nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1'};
    483         Data.VarAttribute{1}.Role='coord_x';
    484         Data.VarAttribute{2}.Role='coord_y';
    485         Data.VarAttribute{3}.Role='vector_x';
    486         Data.VarAttribute{4}.Role='vector_y';
    487         Data.VarAttribute{5}.Role='ancillary';
    488         Data.VarAttribute{6}.Role='errorflag';
     479        %         Data.ListVarName={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_FF'};%  cell array containing the names of the fields to record
     480        %         Data.VarDimName={'nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1'};
     481        %         Data.VarAttribute{1}.Role='coord_x';
     482        %         Data.VarAttribute{2}.Role='coord_y';
     483        %         Data.VarAttribute{3}.Role='vector_x';
     484        %         Data.VarAttribute{4}.Role='vector_y';
     485        %         Data.VarAttribute{5}.Role='ancillary';
     486        %         Data.VarAttribute{6}.Role='errorflag';
    489487       
    490488        % case of mask
    491489        if par_civ1.CheckMask&&~isempty(par_civ1.Mask)
    492490            [RootPath_mask,SubDir_mask,RootFile_mask,~,~,~,~,Ext_mask]=fileparts_uvmat(Param.ActionInput.Civ1.Mask);
    493             j1=1;
    494             if ~isempty(j1_series_Civ1)
     491            if Check_j_Civ1
    495492                j1=j1_series_Civ1(ifield);
     493            else
     494                j1=[];
    496495            end
    497496            if ~isempty(i2_series_Civ1)&& ~isequal(i1_series_Civ1,i2_series_Civ1)% case of volume,masks act on different j levels
    … …  
    528527            end
    529528        end
    530 
     529       
    531530        % case of input grid
    532531        if par_civ1.CheckGrid &&~isempty(par_civ1.Grid)
    … …  
    537536       
    538537        % caluclate velocity data
    539         [Data.Civ1_X,Data.Civ1_Y,Data.Civ1_U,Data.Civ1_V,Data.Civ1_C,Data.Civ1_FF, result_conv, errormsg] = civ (par_civ1);
     538        %   [Data.Civ1_X,Data.Civ1_Y,Data.Civ1_U,Data.Civ1_V,Data.Civ1_C,Data.Civ1_FF, result_conv, errormsg] = civ (par_civ1);
     539        [Civ_X,Civ_Y,Civ_U,Civ_V,Civ_C,Civ_FF, result_conv, errormsg] = civ (par_civ1);
     540        Civ_X_shifted=Civ_X-0.5+Civ_U/2;% get the exact positions
     541        Civ_Y_shifted=Civ_Y-0.5+Civ_V/2;
    540542        if ~isempty(errormsg)
    541543            disp_uvmat('ERROR',errormsg,checkrun)
    … …  
    562564        end
    563565        Data.ListGlobalAttribute=[Data.ListGlobalAttribute Fix1_param];
    564         Data.Civ1_FF=uint8(detect_false(Param.ActionInput.Fix1,Data.Civ1_C,Data.Civ1_U,Data.Civ1_V,Data.Civ1_FF));
     566        % Data.Civ1_FF=uint8(detect_false(Param.ActionInput.Fix1,Data.Civ1_C,Data.Civ1_U,Data.Civ1_V,Data.Civ1_FF));
     567        Civ_FF=uint8(detect_false(Param.ActionInput.Fix1,Civ_C,Civ_U,Civ_V,Civ_FF));
    565568        Data.CivStage=2;
    566569    end
    … …  
    581584       
    582585        % list the variables to record
    583         nbvar=length(Data.ListVarName);
    584         Data.ListVarName=[Data.ListVarName {'Civ1_U_smooth','Civ1_V_smooth','Civ1_SubRange','Civ1_NbCentres','Civ1_Coord_tps','Civ1_U_tps','Civ1_V_tps'}];
    585         Data.VarDimName=[Data.VarDimName {'nb_vec_1','nb_vec_1',{'nb_coord','nb_bounds','nb_subdomain_1'},'nb_subdomain_1',...
    586             {'nb_tps_1','nb_coord','nb_subdomain_1'},{'nb_tps_1','nb_subdomain_1'},{'nb_tps_1','nb_subdomain_1'}}];
    587         Data.VarAttribute{nbvar+1}.Role='vector_x';
    588         Data.VarAttribute{nbvar+2}.Role='vector_y';
    589         Data.VarAttribute{nbvar+5}.Role='coord_tps';
    590         Data.VarAttribute{nbvar+6}.Role='vector_x';
    591         Data.VarAttribute{nbvar+7}.Role='vector_y';
    592         Data.Civ1_U_smooth=Data.Civ1_U; % zeros(size(Data.Civ1_X));
    593         Data.Civ1_V_smooth=Data.Civ1_V; %zeros(size(Data.Civ1_X));
    594         if isfield(Data,'Civ1_FF')
    595             ind_good=find(Data.Civ1_FF==0);
     586        %         nbvar=length(Data.ListVarName);
     587        %         Data.ListVarName=[Data.ListVarName {'Civ1_U_smooth','Civ1_V_smooth','Civ1_SubRange','Civ1_NbCentres','Civ1_Coord_tps','Civ1_U_tps','Civ1_V_tps'}];
     588        %         Data.VarDimName=[Data.VarDimName {'nb_vec_1','nb_vec_1',{'nb_coord','nb_bounds','nb_subdomain_1'},'nb_subdomain_1',...
     589        %             {'nb_tps_1','nb_coord','nb_subdomain_1'},{'nb_tps_1','nb_subdomain_1'},{'nb_tps_1','nb_subdomain_1'}}];
     590        %         Data.VarAttribute{nbvar+1}.Role='vector_x';
     591        %         Data.VarAttribute{nbvar+2}.Role='vector_y';
     592        %         Data.VarAttribute{nbvar+5}.Role='coord_tps';
     593        %         Data.VarAttribute{nbvar+6}.Role='vector_x';
     594        %         Data.VarAttribute{nbvar+7}.Role='vector_y';
     595        %Data.Civ1_U_smooth=Data.Civ1_U; % zeros(size(Data.Civ1_X));
     596        %Data.Civ1_V_smooth=Data.Civ1_V; %zeros(size(Data.Civ1_X));
     597        %         if isfield(Data,'Civ1_FF')
     598        if isempty(Civ_FF)
     599            ind_good=1:numel(Civ_X);
    596600        else
    597             ind_good=1:numel(Data.Civ1_X);
     601            ind_good=find(Civ_FF==0);
    598602        end
    599603        if isempty(ind_good)
    … …  
    603607       
    604608        % perform Patch calculation using the UVMAT fct 'filter_tps'
    605         [Data.Civ1_SubRange,Data.Civ1_NbCentres,Data.Civ1_Coord_tps,Data.Civ1_U_tps,Data.Civ1_V_tps,~,Ures, Vres,~,FFres]=...
    606             filter_tps([Data.Civ1_X(ind_good) Data.Civ1_Y(ind_good)],Data.Civ1_U(ind_good),Data.Civ1_V(ind_good),[],Data.Patch1_SubDomainSize,Data.Patch1_FieldSmooth,Data.Patch1_MaxDiff);
    607         Data.Civ1_U_smooth(ind_good)=Ures;% take the interpolated (smoothed) velocity values for good vectors, keep civ1 data for the other
    608         Data.Civ1_V_smooth(ind_good)=Vres;
    609         Data.Civ1_FF(ind_good)=uint8(4*FFres);%set FF to value =4 for vectors eliminated by filter_tps
     609        [SubRange,NbCentres,Coord_tps,U_tps,V_tps,~,Ures, Vres,~,FFres]=...
     610            filter_tps([Civ_X_shifted(ind_good),Civ_Y_shifted(ind_good)],Civ_U(ind_good),Civ_V(ind_good),[],Data.Patch1_SubDomainSize,Data.Patch1_FieldSmooth,Data.Patch1_MaxDiff);
     611        Civ_U_smooth=Civ_U;% false vectors kept unchanged
     612        Civ_V_smooth=Civ_V;
     613        Civ_U_smooth(ind_good)=Ures;% take the interpolated (smoothed) velocity values for good vectors, keep civ1 data for the other
     614        Civ_V_smooth(ind_good)=Vres;
     615        Civ_FF(ind_good)=uint8(4*FFres);%set FF to value =4 for vectors eliminated by filter_tps
    610616        time_patch1=toc(tstart_patch1);
    611617        disp('patch1 performed')
    … …  
    666672        end
    667673       
    668 %                 %% user defined image transform
    669 %         if ~isempty(transform_fct)
    670 %                par_civ2 =transform_fct(par_civ2,Param);
    671 %         end
    672          %% case of image luminosity rescaling
    673       if par_civ2.CheckRescale &&~isempty(par_civ2.Maxtanh)
    674                par_civ2.ImageA =par_civ2.Maxtanh*tanh(double(par_civ2.ImageA)/par_civ2.Maxtanh);
    675                par_civ2.ImageB=par_civ2.Maxtanh*tanh(double(par_civ2.ImageB)/par_civ2.Maxtanh);
     674        %                 %% user defined image transform
     675        %         if ~isempty(transform_fct)
     676        %                par_civ2 =transform_fct(par_civ2,Param);
     677        %         end
     678        %% case of image luminosity rescaling
     679        if par_civ2.CheckRescale &&~isempty(par_civ2.Maxtanh)
     680            par_civ2.ImageA =par_civ2.Maxtanh*tanh(double(par_civ2.ImageA)/par_civ2.Maxtanh);
     681            par_civ2.ImageB=par_civ2.Maxtanh*tanh(double(par_civ2.ImageB)/par_civ2.Maxtanh);
    676682        end
    677683       
    … …  
    688694            U_tps=DataIn.Civ2_U_tps;
    689695            V_tps=DataIn.Civ2_V_tps;
    690             %CivStage=DataIn.CivStage;%store the current CivStage
    691696            Civ1_Dt=DataIn.Civ2_Dt;
    692697            Data=[];%reinitialise the result structure Data
    … …  
    694699            Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes
    695700            Data.Program='civ_series';
    696            % Data.CivStage=CivStage+1;%update the current civStage after reinitialisation of Data
    697701            Data.ListVarName={};
    698702            Data.VarDimName={};
    699703        else % get the guess from patch1
    700             SubRange= Data.Civ1_SubRange;
    701             NbCentres=Data.Civ1_NbCentres;
    702             Coord_tps=Data.Civ1_Coord_tps;
    703             U_tps=Data.Civ1_U_tps;
    704             V_tps=Data.Civ1_V_tps;
     704            %             SubRange= Data.Civ_SubRange;
     705            %             NbCentres=Data.Civ_NbCentres;
     706            %             Coord_tps=Data.Civ_Coord_tps;
     707            %             U_tps=Data.Civ_U_tps;
     708            %             V_tps=Data.Civ_V_tps;
    705709            Civ1_Dt=Data.Civ1_Dt;
    706 %             Data.CivStage=4;
    707         end
    708          Data.CivStage=4;
    709         %             SubRange= par_civ2.Civ1_SubRange;
    710         %             NbCentres=par_civ2.Civ1_NbCentres;
    711         %             Coord_tps=par_civ2.Civ1_Coord_tps;
    712         %             U_tps=par_civ2.Civ1_U_tps;
    713         %             V_tps=par_civ2.Civ1_V_tps;
    714         %             Civ1_Dt=par_civ2.Civ1_Dt;
    715         %             Civ2_Dt=par_civ2.Civ1_Dt;
    716         %             Data.ListVarName={};
    717         %             Data.VarDimName={};
    718         %         end
     710        end
     711        Data.CivStage=4;
    719712        Shiftx=zeros(size(par_civ2.Grid,1),1);% initialise the shift expected from civ1 data
    720713        Shifty=zeros(size(par_civ2.Grid,1),1);
    … …  
    761754            [RootPath_mask,SubDir_mask,RootFile_mask,~,~,~,~,Ext_mask]=fileparts_uvmat(Param.ActionInput.Civ2.Mask);
    762755            if ~isempty(i2_series_Civ2) && ~isequal(i1_series_Civ2,i2_series_Civ2) % we do PIV among indices i,  at given indices j (volume scan), mask depends on position j
    763                 j1=1;
    764                 if ~isempty(j1_series_Civ2)
    765                     j1=j1_series_Civ1(ifield);
     756                if Check_j_Civ2
     757                    j1=j1_series_Civ2(ifield);
     758                else
     759                    j1=[];
    766760                end
    767761                maskname=fullfile_uvmat(RootPath_mask,SubDir_mask,RootFile_mask,Ext_mask,'_1',j1);
    … …  
    859853        % calculate velocity data (y and v in image indices, reverse to y component)
    860854       
    861         [Data.Civ2_X,Data.Civ2_Y,Data.Civ2_U,Data.Civ2_V,Data.Civ2_C,Data.Civ2_FF,~, errormsg] = civ (par_civ2);
    862        
     855        [Civ_X,Civ_Y,Civ_U,Civ_V,Civ_C,Civ_FF,~, errormsg] = civ (par_civ2);
     856        Civ_X_shifted=Civ_X-0.5+Civ_U/2;% get the exact positions
     857        Civ_Y_shifted=Civ_Y-0.5+Civ_V/2;
    863858        list_param=(fieldnames(Param.ActionInput.Civ2))';
    864859        list_param(strcmp('TestCiv2',list_param))=[];% remove the parameter TestCiv2 from the list
    … …  
    882877        Data.ListGlobalAttribute=[Data.ListGlobalAttribute Civ2_param];
    883878       
    884         nbvar=numel(Data.ListVarName);
     879        %        nbvar=numel(Data.ListVarName);
    885880        % define the Civ2 variable (if Civ2 data are not replaced from previous calculation)
    886         if isempty(find(strcmp('Civ2_X',Data.ListVarName),1))
    887             Data.ListVarName=[Data.ListVarName {'Civ2_X','Civ2_Y','Civ2_U','Civ2_V','Civ2_C','Civ2_FF'}];%  cell array containing the names of the fields to record
    888             Data.VarDimName=[Data.VarDimName {'nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2'}];
    889             Data.VarAttribute{nbvar+1}.Role='coord_x';
    890             Data.VarAttribute{nbvar+2}.Role='coord_y';
    891             Data.VarAttribute{nbvar+3}.Role='vector_x';
    892             Data.VarAttribute{nbvar+4}.Role='vector_y';
    893             Data.VarAttribute{nbvar+5}.Role='ancillary';
    894             Data.VarAttribute{nbvar+6}.Role='errorflag';
    895         end
     881        %         if isempty(find(strcmp('Civ2_X',Data.ListVarName),1))
     882        %             Data.ListVarName=[Data.ListVarName {'Civ2_X','Civ2_Y','Civ2_U','Civ2_V','Civ2_C','Civ2_FF'}];%  cell array containing the names of the fields to record
     883        %             Data.VarDimName=[Data.VarDimName {'nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2'}];
     884        %             Data.VarAttribute{nbvar+1}.Role='coord_x';
     885        %             Data.VarAttribute{nbvar+2}.Role='coord_y';
     886        %             Data.VarAttribute{nbvar+3}.Role='vector_x';
     887        %             Data.VarAttribute{nbvar+4}.Role='vector_y';
     888        %             Data.VarAttribute{nbvar+5}.Role='ancillary';
     889        %             Data.VarAttribute{nbvar+6}.Role='errorflag';
     890        %         end
    896891        disp('civ2 performed')
    897892        time_civ2=toc(tstart_civ2);
    898     elseif ~isfield(Data,'ListVarName') % we start there, using existing Civ2 data
     893    elseif Param.ActionInput.CheckFix2 && isfield (Param.ActionInput,'Fix2') % we start there, using existing Civ2 data
    899894        if exist('ncfile','var')
    900895            CivFile=ncfile;
    … …  
    908903   
    909904    %% Fix2
    910     if Param.ActionInput.CheckFix2 && isfield (Param.ActionInput,'Fix2')% if Fix2 computation is requested
     905    if (Param.ActionInput.CheckFix2 && isfield (Param.ActionInput,'Fix2'))||(Param.ActionInput.CheckPatch2 && isfield (Param.ActionInput,'Patch2'))%% if Fix2 computation is requested
    911906        disp('detect_false2 started')
    912907        list_param=fieldnames(Param.ActionInput.Fix2)';
    … …  
    917912        end
    918913        Data.ListGlobalAttribute=[Data.ListGlobalAttribute Fix2_param];
    919         Data.Civ2_FF=double(detect_false(Param.ActionInput.Fix2,Data.Civ2_C,Data.Civ2_U,Data.Civ2_V,Data.Civ2_FF));
     914        Civ_FF=uint8(detect_false(Param.ActionInput.Fix2,Civ_C,Civ_U,Civ_V,Civ_FF));
    920915        Data.CivStage=Data.CivStage+1;
    921916    end
    … …  
    935930        Data.ListGlobalAttribute=[Data.ListGlobalAttribute Patch2_param];
    936931       
    937         nbvar=length(Data.ListVarName);
    938         Data.ListVarName=[Data.ListVarName {'Civ2_U_smooth','Civ2_V_smooth','Civ2_SubRange','Civ2_NbCentres','Civ2_Coord_tps','Civ2_U_tps','Civ2_V_tps'}];
    939         Data.VarDimName=[Data.VarDimName {'nb_vec_2','nb_vec_2',{'nb_coord','nb_bounds','nb_subdomain_2'},{'nb_subdomain_2'},...
    940             {'nb_tps_2','nb_coord','nb_subdomain_2'},{'nb_tps_2','nb_subdomain_2'},{'nb_tps_2','nb_subdomain_2'}}];
    941        
    942         Data.VarAttribute{nbvar+1}.Role='vector_x';
    943         Data.VarAttribute{nbvar+2}.Role='vector_y';
    944         Data.VarAttribute{nbvar+5}.Role='coord_tps';
    945         Data.VarAttribute{nbvar+6}.Role='vector_x';
    946         Data.VarAttribute{nbvar+7}.Role='vector_y';
    947         Data.Civ2_U_smooth=Data.Civ2_U;
    948         Data.Civ2_V_smooth=Data.Civ2_V;
    949         if isfield(Data,'Civ2_FF')
    950             ind_good=find(Data.Civ2_FF==0);
     932        %  nbvar=length(Data.ListVarName);
     933        %         Data.ListVarName=[Data.ListVarName {'Civ2_U_smooth','Civ2_V_smooth','Civ2_SubRange','Civ2_NbCentres','Civ2_Coord_tps','Civ2_U_tps','Civ2_V_tps'}];
     934        %         Data.VarDimName=[Data.VarDimName {'nb_vec_2','nb_vec_2',{'nb_coord','nb_bounds','nb_subdomain_2'},{'nb_subdomain_2'},...
     935        %             {'nb_tps_2','nb_coord','nb_subdomain_2'},{'nb_tps_2','nb_subdomain_2'},{'nb_tps_2','nb_subdomain_2'}}];
     936       
     937        %         Data.VarAttribute{nbvar+1}.Role='vector_x';
     938        %         Data.VarAttribute{nbvar+2}.Role='vector_y';
     939        %         Data.VarAttribute{nbvar+5}.Role='coord_tps';
     940        %         Data.VarAttribute{nbvar+6}.Role='vector_x';
     941        %         Data.VarAttribute{nbvar+7}.Role='vector_y';
     942       
     943        if isempty(Civ_FF)
     944            ind_good=1:numel(Civ_X);
    951945        else
    952             ind_good=1:numel(Data.Civ2_X);
     946            ind_good=find(Civ_FF==0);
    953947        end
    954948        if isempty(ind_good)
    … …  
    957951        end
    958952       
    959         [Data.Civ2_SubRange,Data.Civ2_NbCentres,Data.Civ2_Coord_tps,Data.Civ2_U_tps,Data.Civ2_V_tps,tild,Ures,Vres,tild,FFres]=...
    960             filter_tps([Data.Civ2_X(ind_good) Data.Civ2_Y(ind_good)],Data.Civ2_U(ind_good),Data.Civ2_V(ind_good),[],Data.Patch2_SubDomainSize,Data.Patch2_FieldSmooth,Data.Patch2_MaxDiff);
    961         Data.Civ2_U_smooth(ind_good)=Ures;
    962         Data.Civ2_V_smooth(ind_good)=Vres;
    963         Data.Civ2_FF(ind_good)=uint8(4*FFres);
     953        [Civ_SubRange,Civ_NbCentres,Civ_Coord_tps,Civ_U_tps,Civ_V_tps,~,Ures,Vres,~,FFres]=...
     954            filter_tps([Civ_X_shifted(ind_good) Civ_Y_shifted(ind_good)],Civ_U(ind_good),Civ_V(ind_good),[],Data.Patch2_SubDomainSize,Data.Patch2_FieldSmooth,Data.Patch2_MaxDiff);
     955        Civ_U_smooth=Civ_U;% keep the false vectors unchanged
     956        Civ_V_smooth=Civ_V;
     957        Civ_U_smooth(ind_good)=Ures;
     958        Civ_V_smooth(ind_good)=Vres;
     959        Civ_FF(ind_good)=uint8(4*FFres);
    964960        Data.CivStage=Data.CivStage+1;
    965961        time_patch2=toc(tstart_patch2);
    … …  
    968964   
    969965    %% write result in a netcdf file
     966   
     967    Data.ListVarName={'X','Y','U','V','C','FF'};%  cell array containing the names of the fields to record
     968    Data.VarDimName={'nb_vec','nb_vec','nb_vec','nb_vec','nb_vec','nb_vec'};
     969    Data.VarAttribute{1}.Role='coord_x';
     970    Data.VarAttribute{2}.Role='coord_y';
     971    Data.VarAttribute{3}.Role='vector_x';
     972    Data.VarAttribute{3}.scale_factor=1/inv_scale_factor;
     973    Data.VarAttribute{4}.Role='vector_y';
     974    Data.VarAttribute{4}.scale_factor=1/inv_scale_factor;
     975    Data.VarAttribute{5}.Role='ancillary';
     976    Data.VarAttribute{5}.scale_factor=1/inv_scale_factor;
     977    Data.VarAttribute{6}.Role='errorflag';
     978    Data.X=uint16(Civ_X);
     979    Data.Y=uint16(Civ_Y);
     980    Data.U=int16(inv_scale_factor*Civ_U);
     981    Data.V=int16(inv_scale_factor*Civ_V);
     982    Data.C=uint8(inv_scale_factor*Civ_C);
     983    Data.FF=uint8(Civ_FF);
     984    if (Param.ActionInput.CheckPatch1 && ~Param.ActionInput.CheckCiv2) ||Param.ActionInput.CheckPatch2
     985        nbvar=6;
     986        %     Data.ListVarName=[Data.ListVarName {'U_smooth','V_smooth','SubRange','NbCentres','Coord_tps','U_tps','V_tps'}];
     987        %         Data.VarDimName=[Data.VarDimName {'nb_vec','nb_vec',{'nb_coord','nb_bounds','nb_subdomain'},{'nb_subdomain'},...
     988        %             {'nb_tps','nb_coord','nb_subdomain'},{'nb_tps','nb_subdomain'},{'nb_tps','nb_subdomain'}}];
     989        Data.ListVarName=[Data.ListVarName {'U_smooth','V_smooth'}];
     990        Data.VarDimName=[Data.VarDimName {'nb_vec','nb_vec'}];
     991        Data.VarAttribute{nbvar+1}.Role='vector_x';
     992        Data.VarAttribute{nbvar+1}.scale_factor=1/inv_scale_factor;
     993        Data.VarAttribute{nbvar+2}.Role='vector_y';
     994        Data.VarAttribute{nbvar+2}.scale_factor=1/inv_scale_factor;
     995        Data.U_smooth=int16(inv_scale_factor*Civ_U_smooth);
     996        Data.V_smooth=int16(inv_scale_factor*Civ_V_smooth);
     997        %         Data.VarAttribute{nbvar+5}.Role='coord_tps';
     998        %         Data.VarAttribute{nbvar+6}.Role='vector_x';
     999        %         Data.VarAttribute{nbvar+7}.Role='vector_y';
     1000        %         Data.U_smooth=int16(inv_scale_factor*U_smooth);
     1001        %         Data.V_smooth=int16(inv_scale_factor*V_smooth);
     1002        %         Data.SubRange=SubRange;
     1003        %         Data.NbCentres=NbCentres;
     1004        %         Data.Coord_tps=Coord_tps;
     1005        %         Data.U_tps=U_tps;
     1006        %         Data.V_tps=V_tps;
     1007    end
     1008    %
     1009    %      if ~isempty(inv_scale_factor)
     1010    %              Data=compress_data(Data,inv_scale_factor);% compress the data using integers instead of (single precision)floating reals
     1011    %      end
    9701012    errormsg=struct2nc(ncfile_out,Data);
    9711013    if isempty(errormsg)
    … …  
    10741116end
    10751117
    1076 
    1077 
     1118% % --- compress the data using integers instead of (single precision)floating reals
     1119% function DataOut=compress_data(DataIn,inv_scale_factor)
     1120% DataOut=DataIn; %default
     1121% DataOut.Conventions= 'uvmat/civdata/compress';
     1122% ListVarName=DataIn.ListVarName;
     1123% for ilist=1:numel(ListVarName)
     1124%     switch ListVarName{ilist}
     1125%         case {'Civ1_C','Civ2_C'}
     1126%             DataOut.(ListVarName{ilist})=uint8(inv_scale_factor*DataIn.(ListVarName{ilist}));
     1127%             DataOut.VarAttribute{ilist}.scale_factor=1/inv_scale_factor;
     1128%         case {'Civ1_U','Civ1_V','Civ2_U','Civ2_V','Civ1_U_smooth','Civ1_V_smooth','Civ2_U_smooth','Civ2_V_smooth'}
     1129%             DataOut.(ListVarName{ilist})=int16(inv_scale_factor*DataIn.(ListVarName{ilist}));
     1130%             DataOut.VarAttribute{ilist}.scale_factor=1/inv_scale_factor;
     1131%         case 'Civ1_X'
     1132%             DataOut.Civ1_X=uint16(DataIn.Civ1_X);
     1133%         case 'Civ1_Y'
     1134%             DataOut.Civ1_Y=uint16(DataIn.Civ1_Y);
     1135%        case 'Civ2_X'
     1136%             DataOut.Civ2_X=uint16(DataIn.Civ2_X);
     1137%         case 'Civ2_Y'
     1138%             DataOut.Civ2_Y=uint16(DataIn.Civ2_Y);
     1139%     end
     1140% end

  • trunk/src/series/test_filter_tps.m

    r1161 r1195  
    191191FieldSmooth=(Param.ActionInput.FieldSmooth)*[0.1 0.2 0.5 1 2 5 10];%scan the smoothing param from 1/10 to 10 current value
    192192NbSmooth=numel(FieldSmooth);
    193 % for irho=1:NbSmooth
    194 %     str=num2str(FieldSmooth(irho));
    195 %     str=regexprep(str,'\.','p');
    196 %     Ustr{irho}=['U_' str];
    197 %     Vstr{irho}=['V_' str];
    198 %     Xstr{irho}=['X_' str];
    199 %     Ystr{irho}=['Y_' str];
    200 %     Dimstr{irho}='NbVec';
    201 %     str_i{irho}=str;
    202 % end
     193
    203194
    204195%% Prepare the structure of output netcdf file

  • trunk/src/struct2nc.m

    r1194 r1195  
    114114        VarType='';
    115115        switch VarClass{ivar}
    116             case {'single','double'}
     116            case {'single','double','int32','uint32','int64','uint64'}
    117117                VarType='nc_float'; % store all floating reals as single
    118             case {'int8','uint8','int16','uint16','int32','uint32','int64','uint64'}
    119                 VarType='nc_int';
     118            case {'int16','uint16'}
     119                  VarType='nc_short';
     120            case {'int8','uint8'}
     121                VarType='nc_byte';
    120122            case 'logical'
    121                 VarType='nc_int';
     123                VarType='nc_byte';
    122124                Data.(ListVarName{ivar})=uint8(Data.(ListVarName{ivar}));
    123125        end
Note: See TracChangeset for help on using the changeset viewer.