[8]  1  %'RUN_FIX': function for fixing velocity fields:


 2  %


 3  % RUN_FIX(filename,field,flagindex,thresh_vecC,thresh_vel,iter,flag_mask,maskname,fileref,fieldref)


 4  %


 5  %filename: name of the netcdf file (used as input and output)


[47]  6  %field: structure specifying the names of the fields to fix (depending on civ1 or civ2)


 7  %.vel_type='civ1' or 'civ2';


 8  %.nb=name of the dimension common to the field to fix ('nb_vectors' for civ1);


 9  %.fixflag=name of fix flag variable ('vec_FixFlag' for civ1)


[8]  10  %flagindex: flag specifying which values of vec_f are removed:


 11  % if flagindex(1)=1: vec_f=2 vectors are removed


 12  % if flagindex(2)=1: vec_f=3 vectors are removed


 13  % if flagindex(3)=1: vec_f=2 vectors are removed (if iter=1) or vec_f=4 vectors are removed (if iter=2)


 14  %iter=1 for civ1 fields and iter=2 for civ2 fields


 15  %thresh_vecC: threshold in the image correlation vec_C


 16  %flag_mask: =1 mask used to remove vectors (0 else)


 17  %maskname: name of the mask image file for fix


 18  %thresh_vel: threshold on velocity, or on the difference with the reference file fileref if exists


 19  %inf_sup=1: remove values smaller than threshold thresh_vel, =2, larger than threshold


 20  %fileref: .nc file name for a reference velocity (='': refrence 0 used)


 21  %fieldref: 'civ1','filter1'...feld used in fileref


 22 


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


 24  error=[]; %default


 25  vel_type{1}=field.vel_type;


[47]  26  %check writing access


 27  [errorread,message]=fileattrib(filename);


[67]  28  if ischar(message)


 29  msgbox_uvmat('ERROR',[filename ':' message]);


 30  return


 31  end


 32  if ~isequal(message.UserWrite,1)


[47]  33  msgbox_uvmat('ERROR',['no writting access to ' filename ' (RUN_FIX.m)']);


 34  return


 35  end


 36  Field=read_civxdata(filename,'ima_cor',field.vel_type);


[8]  37  if isfield(Field,'Txt')


 38  error=Field.Txt; %error in reading


 39  return


 40  end


 41 


 42  if ~isfield(Field,'X')  ~isfield(Field,'Y')  ~isfield(Field,'U')  ~isfield(Field,'V')


 43  error=['input file ' filename ' does not contain vectors in RUN_FIX.m']; %bad input file


 44  return


 45  end


 46  if ~isfield(Field,'C')


 47  Field.C=ones(size(Field.X));%correlation=1 by default


 48  end


 49  if ~isfield(Field,'F')


 50  Field.F=ones(size(Field.X));%warning flag=1 by default


 51  end


 52  if ~isfield(Field,'FF')


 53  Field.FF=zeros(size(Field.X));%fixflag=0 by default


 54  end


 55 


 56  vec_f_unit=abs(Field.F)10*double(uint16(abs(Field.F)/10)); %unityterm of vec_F in abs value


 57  vec_f_sign=sign(Field.F).*vec_f_unit;% gives the unity digit of vec_f with correct sign


 58  flag1=(flagindex(1)==1)&(vec_f_sign==2);%removed vectors vec_f=2


 59  flag2=(flagindex(2)==1)&(vec_f_sign==3);%removed vectors vec_f=3


 60  if iter==1


 61  flag3=(flagindex(3)==1)&(vec_f_sign==2); % Hart vectors


 62  elseif iter==2


 63  flag3=(flagindex(3)==1)&(vec_f_sign==4); %


 64  end


 65  flag4=(Field.C < thresh_vecC)&(flag1~=1)&(flag2~=1)&(flag3~=1); % =1 for low vec_C vectors not previously removed


 66 


 67  % criterium on velocity values


 68  delta_u=Field.U;%default without ref file


 69  delta_v=Field.V;


 70  if exist('fileref','var') && ~isempty(fileref)


 71  if ~exist(fileref,'file')


 72  error='reference file not found in RUN_FIX.m';


 73  display(error);


 74  return


 75  end


 76  FieldRef=read_civxdata(fileref,[],fieldref);


 77  if isfield(FieldRef,'FF')


 78  index_true=find(FieldRef.FF==0);


 79  FieldRef.X=FieldRef.X(index_true);


 80  FieldRef.Y=FieldRef.Y(index_true);


 81  FieldRef.U=FieldRef.U(index_true);


 82  FieldRef.V=FieldRef.V(index_true);


 83  end


 84  if ~isfield(FieldRef,'X')  ~isfield(FieldRef,'Y')  ~isfield(FieldRef,'U')  ~isfield(FieldRef,'V')


 85  error='reference file is not a velocity field in RUN_FIX.m '; %bad input file


 86  return


 87  end


 88  if length(FieldRef.X)<=1


 89  errordlg('reference field with one vector or less in RUN_FIX.m')


 90  return


 91  end


 92  vec_U_ref=griddata_uvmat(FieldRef.X,FieldRef.Y,FieldRef.U,Field.X,Field.Y); %interpolate vectors in the ref field


 93  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


 94  delta_u=Field.Uvec_U_ref;%take the difference with the interpolated ref field


 95  delta_v=Field.Vvec_V_ref;


 96  end


 97  thresh_vel_x=thresh_vel;


 98  thresh_vel_y=thresh_vel;


 99  if isequal(inf_sup,1)


 100  flag5=abs(delta_u)<thresh_vel_x & abs(delta_v)<thresh_vel_y &(flag1~=1)&(flag2~=1)&(flag3~=1)&(flag4~=1);


 101  elseif isequal(inf_sup,2)


 102  flag5=(abs(delta_u)>thresh_vel_x  abs(delta_v)>thresh_vel_y) &(flag1~=1)&(flag2~=1)&(flag3~=1)&(flag4~=1);


 103  end


 104 


 105  % flag7 introduce a grey mask, matrix M


 106  if isequal (flag_mask,1)


 107  M=imread(maskname);


 108  nxy=size(M);


 109  M=reshape(M,1,nxy(1)*nxy(2));


 110  rangx0=[0.5 nxy(2)0.5];


 111  rangy0=[0.5 nxy(1)0.5];


 112  vec_x1=Field.XField.U/2;%beginning points


 113  vec_x2=Field.X+Field.U/2;%end points of vectors


 114  vec_y1=Field.YField.V/2;%beginning points


 115  vec_y2=Field.Y+Field.V/2;%end points of vectors


 116  indx=1+round((nxy(2)1)*(vec_x1rangx0(1))/(rangx0(2)rangx0(1)));% image index x at abcissa vec_x


 117  indy=1+round((nxy(1)1)*(vec_y1rangy0(1))/(rangy0(2)rangy0(1)));% image index y at ordinate vec_y


 118  test_in=~(indx < 1 indy < 1  indx > nxy(2) indy > nxy(1)); %=0 out of the mask image, 1 inside


 119  indx=indx.*test_in+(1test_in); %replace indx by 1 out of the mask range


 120  indy=indy.*test_in+(1test_in); %replace indy by 1 out of the mask range


 121  ICOMB=((indx1)*nxy(1)+(nxy(1)+1indy));%determine the indices in the image reshaped in a Matlab vector


 122  Mvalues=M(ICOMB);


 123  flag7b=((20 < Mvalues) & (Mvalues < 200)) ~test_in';


 124  indx=1+round((nxy(2)1)*(vec_x2rangx0(1))/(rangx0(2)rangx0(1)));% image index x at abcissa Field.X


 125  indy=1+round((nxy(1)1)*(vec_y2rangy0(1))/(rangy0(2)rangy0(1)));% image index y at ordinate vec_y


 126  test_in=~(indx < 1 indy < 1  indx > nxy(2) indy > nxy(1)); %=0 out of the mask image, 1 inside


 127  indx=indx.*test_in+(1test_in); %replace indx by 1 out of the mask range


 128  indy=indy.*test_in+(1test_in); %replace indy by 1 out of the mask range


 129  ICOMB=((indx1)*nxy(1)+(nxy(1)+1indy));%determine the indices in the image reshaped in a Matlab vector


 130  Mvalues=M(ICOMB);


 131  flag7e=((Mvalues > 20) & (Mvalues < 200)) ~test_in';


 132  flag7=(flag7bflag7e)';


 133  else


 134  flag7=0;


 135  end


 136  flagmagenta=flag1flag2flag3flag4flag5flag7;


[47]  137  fixflag_unit=Field.FF10*floor(Field.FF/10); %unity term of fix_flag


[8]  138 


[47]  139  %write fix flags in the netcdf file


 140  hhh=which('netcdf.open');% look for builtin matlab netcdf library


 141  if ~isequal(hhh,'')% case of new builtin Matlab netcdf library


 142  nc=netcdf.open(filename,'NC_WRITE');


 143  netcdf.reDef(nc)


 144  if iter==1


 145  netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'fix1',1)


 146  elseif iter==2


 147  netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'fix2',1)


 148  end


 149  dimid = netcdf.inqDimID(nc,field.nb);


 150  try


 151  varid = netcdf.inqVarID(nc,field.fixflag);% look for already existing fixflag variable


 152  catch


 153  varid=netcdf.defVar(nc,field.fixflag,'double',dimid);%create fixflag variable if it does not exist


 154  end


 155  netcdf.endDef(nc)


 156  netcdf.putVar(nc,varid,fixflag_unit+10*flagmagenta);


 157  netcdf.close(nc)


 158  else %old netcdf library


[56]  159  netcdf_toolbox(filename,Field,flagmagenta,iter,field)


[8]  160  end


[56]  161 


 162  function netcdf_toolbox(filename,Field,flagmagenta,iter,field)


 163  nc=netcdf(filename,'write'); %open netcdf file for writing


 164  result=redef(nc);


 165  if isempty(result), msgbox_uvmat('ERROR','##Bad redef operation.'),end


 166  if iter==1


 167  nc.fix=1;


 168  elseif iter==2


 169  nc.fix2=1;


 170  end


 171  nc{field.fixflag}=ncfloat(field.nb);


 172  fixflag_unit=Field.FF10*floor(Field.FF/10); %unity term of fix_flag


 173  nc{field.fixflag}(:)=fixflag_unit+10*flagmagenta;


 174  close(nc);

