Home > . > RUN_FIX.m

# RUN_FIX

## PURPOSE

'RUN_FIX': function for fixing velocity fields:

## SYNOPSIS

function error=RUN_FIX(filename,field,flagindex,iter,thresh_vecC,flag_mask,maskname,thresh_vel,inf_sup,fileref,fieldref)

## DESCRIPTION

    'RUN_FIX': function for fixing velocity fields:
    -----------------------------------------------
     RUN_FIX(filename,field,flagindex,thresh_vecC,thresh_vel,iter,flag_mask,maskname,fileref,fieldref)

    filename: name of the netcdf file (used as input and output)
    field: structure specifying the field to fix (field.vel_type='civ1' or 'civ2')
    flagindex: flag specifying which values of vec_f are removed:
     if flagindex(1)=1: vec_f=-2 vectors are removed
     if flagindex(2)=1: vec_f=3 vectors are removed
     if flagindex(3)=1: vec_f=2 vectors are removed (if iter=1) or vec_f=4 vectors are removed (if iter=2)
    iter=1 for civ1 fields and iter=2 for civ2 fields
    thresh_vecC: threshold in the image correlation vec_C
    flag_mask: =1 mask used to remove vectors (0 else)
    maskname: name of the mask image file for fix
    thresh_vel: threshold on velocity, or on the difference with the reference file fileref if exists
    inf_sup=1: remove values smaller than threshold thresh_vel, =2, larger than threshold
    fileref: .nc file name for a reference velocity (='': refrence 0 used)
    fieldref: 'civ1','filter1'...feld used in fileref

## CROSS-REFERENCE INFORMATION

This function calls:

- griddata_uvmat 'griddata_uvmat': function griddata_uvmat(vec2_X,vec2_Y,vec2_U,vec_X,vec_Y,'linear')
- read_civxdata 'read_civxdata': reads civx data from netcdf files
- warndlg_uvmat 'warndlg_uvmat': display warning message (error, warning, confirmation) in a given figure

This function is called by:

- civ 'civ': function associated with the interface 'civ.fig' for PIV, spline interpolation and stereo PIV (patch)

## SOURCE CODE

    0001 %'RUN_FIX': function for fixing velocity fields:
    0002 %-----------------------------------------------
    0003 % RUN_FIX(filename,field,flagindex,thresh_vecC,thresh_vel,iter,flag_mask,maskname,fileref,fieldref)
    0004 %
    0005 %filename: name of the netcdf file (used as input and output)
    0006 %field: structure specifying the field to fix (field.vel_type='civ1' or 'civ2')
    0007 %flagindex: flag specifying which values of vec_f are removed:
    0008         % if flagindex(1)=1: vec_f=-2 vectors are removed
    0009         % if flagindex(2)=1: vec_f=3 vectors are removed
    0010         % if flagindex(3)=1: vec_f=2 vectors are removed (if iter=1) or vec_f=4 vectors are removed (if iter=2)
    0011 %iter=1 for civ1 fields and iter=2 for civ2 fields
    0012 %thresh_vecC: threshold in the image correlation vec_C
    0013 %flag_mask: =1 mask used to remove vectors (0 else)
    0014 %maskname: name of the mask image file for fix
    0015 %thresh_vel: threshold on velocity, or on the difference with the reference file fileref if exists
    0016 %inf_sup=1: remove values smaller than threshold thresh_vel, =2, larger than threshold
    0017 %fileref: .nc file name for a reference velocity (='': refrence 0 used)
    0018 %fieldref: 'civ1','filter1'...feld used in fileref
    0019
    0020 function error=RUN_FIX(filename,field,flagindex,iter,thresh_vecC,flag_mask,maskname,thresh_vel,inf_sup,fileref,fieldref)
    0021 error=[]; %default
    0022 vel_type{1}=field.vel_type;
    0023 %[nbc,nbd,Civ,CivStage,tt,Field]=read_ncfield(filename,field.vel_type,'vec_color','ima_cor');
    0024 Field=read_civxdata(filename,[],field.vel_type);
    0025 if isfield(Field,'Txt')
    0026     error=Field.Txt; %error in reading
    0027     return
    0028 end
    0029
    0030 if ~isfield(Field,'X') || ~isfield(Field,'Y') || ~isfield(Field,'U') || ~isfield(Field,'V')
    0031     error=['input file ' filename ' does not contain vectors in RUN_FIX.m']; %bad input file
    0032     return
    0033 end
    0034 if ~isfield(Field,'C')
    0035     Field.C=ones(size(Field.X));%correlation=1 by default
    0036 end
    0037 if ~isfield(Field,'F')
    0038     Field.F=ones(size(Field.X));%warning flag=1 by default
    0039 end
    0040 if ~isfield(Field,'FF')
    0041     Field.FF=zeros(size(Field.X));%fixflag=0 by default
    0042 end
    0043
    0044 vec_f_unit=abs(Field.F)-10*double(uint16(abs(Field.F)/10)); %unityterm of vec_F in abs value
    0045 vec_f_sign=sign(Field.F).*vec_f_unit;% gives the unity digit of vec_f with correct sign
    0046 flag1=(flagindex(1)==1)&(vec_f_sign==-2);%removed vectors vec_f=-2
    0047 flag2=(flagindex(2)==1)&(vec_f_sign==3);%removed vectors vec_f=3
    0048 if iter==1
    0049         flag3=(flagindex(3)==1)&(vec_f_sign==2); % Hart vectors
    0050 elseif iter==2
    0051         flag3=(flagindex(3)==1)&(vec_f_sign==4); %
    0052 end
    0053             %  flag4=(vec_Fmin2_2_Fmin2_2_c <
    0054             %  colvec1)&(vec_Fmin2_2_Fmin2_2_f_sign~=2)&(vec_Fmin2_2_Fmin2_
    0055             %  2_f~=0); % =1 for red vectors not flaged
    0056
    0057 flag4=(Field.C < thresh_vecC)&(flag1~=1)&(flag2~=1)&(flag3~=1); % =1 for low vec_C vectors not previously removed
    0058
    0059 % criterium on velocity values
    0060 delta_u=Field.U;%default without ref file
    0061 delta_v=Field.V;
    0062 if exist('fileref','var') && ~isempty(fileref)
    0063     if ~exist(fileref,'file')
    0064         error='reference file not found in RUN_FIX.m';
    0065         display(error);
    0066         return
    0067     end
    0068
    0069 %     [nbc,nbd,Civ,CivStage,tt,FieldRef]=read_ncfield(fileref,fieldref);
    0070     FieldRef=read_civxdata(fileref,[],fieldref);
    0071
    0072     if isfield(FieldRef,'FF')
    0073         index_true=find(FieldRef.FF==0);
    0074         FieldRef.X=FieldRef.X(index_true);
    0075         FieldRef.Y=FieldRef.Y(index_true);
    0076         FieldRef.U=FieldRef.U(index_true);
    0077         FieldRef.V=FieldRef.V(index_true);
    0078     end
    0079     if ~isfield(FieldRef,'X') || ~isfield(FieldRef,'Y') || ~isfield(FieldRef,'U') || ~isfield(FieldRef,'V')
    0080         error='reference file is not a velocity field in RUN_FIX.m '; %bad input file
    0081         return
    0082     end
    0083     if length(FieldRef.X)<=1
    0084         errordlg('reference field with one vector or less in RUN_FIX.m')
    0085         return
    0086     end
    0087     vec_U_ref=griddata_uvmat(FieldRef.X,FieldRef.Y,FieldRef.U,Field.X,Field.Y);  %interpolate vectors in the ref field
    0088     vec_V_ref=griddata_uvmat(FieldRef.X,FieldRef.Y,FieldRef.V,Field.X,Field.Y);  %interpolate vectors in the ref field to the positions  of the main field
    0089     delta_u=Field.U-vec_U_ref;%take the difference with the interpolated ref field
    0090     delta_v=Field.V-vec_V_ref;
    0091 end
    0092
    0093 thresh_vel_x=thresh_vel;
    0094 thresh_vel_y=thresh_vel;
    0095 if isequal(inf_sup,1)
    0096     flag5=abs(delta_u)<thresh_vel_x & abs(delta_v)<thresh_vel_y &(flag1~=1)&(flag2~=1)&(flag3~=1)&(flag4~=1);
    0097 elseif isequal(inf_sup,2)
    0098     flag5=(abs(delta_u)>thresh_vel_x | abs(delta_v)>thresh_vel_y) &(flag1~=1)&(flag2~=1)&(flag3~=1)&(flag4~=1);
    0099 end
    0100
    0101             % flag7 introduce a grey mask, matrix M
    0102 if isequal (flag_mask,1)
    0103    M=imread(maskname);
    0104    nxy=size(M);
    0105    M=reshape(M,1,nxy(1)*nxy(2));
    0106    rangx0=[0.5 nxy(2)-0.5];
    0107    rangy0=[0.5 nxy(1)-0.5];
    0108    vec_x1=Field.X-Field.U/2;%beginning points
    0109    vec_x2=Field.X+Field.U/2;%end points of vectors
    0110    vec_y1=Field.Y-Field.V/2;%beginning points
    0111    vec_y2=Field.Y+Field.V/2;%end points of vectors
    0112    indx=1+round((nxy(2)-1)*(vec_x1-rangx0(1))/(rangx0(2)-rangx0(1)));% image index x at abcissa vec_x
    0113    indy=1+round((nxy(1)-1)*(vec_y1-rangy0(1))/(rangy0(2)-rangy0(1)));% image index y at ordinate vec_y
    0114    test_in=~(indx < 1 |indy < 1 | indx > nxy(2) |indy > nxy(1)); %=0 out of the mask image, 1 inside
    0115    indx=indx.*test_in+(1-test_in); %replace indx by 1 out of the mask range
    0116    indy=indy.*test_in+(1-test_in); %replace indy by 1 out of the mask range
    0117    ICOMB=((indx-1)*nxy(1)+(nxy(1)+1-indy));%determine the indices in the image reshaped in a Matlab vector
    0118    Mvalues=M(ICOMB);
    0119 %    size(indx)
    0120 %     size(indy)
    0121 %    Ma=M(indy,indx);% relevant mask values
    0122 %    Mvalues=(diag(Ma,0))'; % mask values for vector indices, flag7=1 for masked vectors
    0123    flag7b=((20 < Mvalues) & (Mvalues < 200))| ~test_in';
    0124    indx=1+round((nxy(2)-1)*(vec_x2-rangx0(1))/(rangx0(2)-rangx0(1)));% image index x at abcissa Field.X
    0125    indy=1+round((nxy(1)-1)*(vec_y2-rangy0(1))/(rangy0(2)-rangy0(1)));% image index y at ordinate vec_y
    0126    test_in=~(indx < 1 |indy < 1 | indx > nxy(2) |indy > nxy(1)); %=0 out of the mask image, 1 inside
    0127    indx=indx.*test_in+(1-test_in); %replace indx by 1 out of the mask range
    0128    indy=indy.*test_in+(1-test_in); %replace indy by 1 out of the mask range
    0129    ICOMB=((indx-1)*nxy(1)+(nxy(1)+1-indy));%determine the indices in the image reshaped in a Matlab vector
    0130    Mvalues=M(ICOMB);
    0131 %    Ma=M(indy,indx);% relevant mask values
    0132 %    Mvalues=(diag(Ma,0))'; % mask values for vector indices, flag7=1 for masked vectors
    0133    flag7e=((Mvalues > 20) & (Mvalues < 200))| ~test_in';
    0134    flag7=(flag7b|flag7e)';
    0135 else
    0136    flag7=0;
    0137 end
    0138 % flagmagenta=flag1|flag2|flag3|flag4|flag7;
    0139 flagmagenta=flag1|flag2|flag3|flag4|flag5|flag7;
    0140 %write fix flags
    0141
    0142 [errorread,message]=fileattrib(filename);
    0143
    0144 if ~isempty(message) & ~isequal(message.UserWrite,1)
    0145      warndlg_uvmat(['no writting access to ' filename ' (RUN_FIX.m)'],'ERROR');
    0146     return
    0147 end
    0148 nc=netcdf(filename,'write'); %open netcdf file for writing
    0149 result=redef(nc);
    0150 if isempty(result), errordlg('##Bad redef operation.'),end
    0151 if iter==1
    0152     nc.fix=1;
    0153 elseif iter==2
    0154     nc.fix2=1;
    0155 end
    0156 theDim=nc(field.nb) ;% get the number of velocity vectors
    0157 % nb_vectors=size(theDim);
    0158 % nb_vectors=nb_vectors(1);
    0159 nc{field.fixflag}=ncfloat(field.nb);
    0160 fixflag_unit=Field.FF-10*floor(Field.FF/10); %unity term of fix_flag
    0161 nc{field.fixflag}(:)=fixflag_unit+10*flagmagenta;
    0162 close(nc);