 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)


 6  %field: structure specifying the field to fix (field.vel_type='civ1' or 'civ2')


 7  %flagindex: flag specifying which values of vec_f are removed:


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


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


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


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


 12  %thresh_vecC: threshold in the image correlation vec_C


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


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


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


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


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


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


 19 


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


 21  error=[]; %default


 22  vel_type{1}=field.vel_type;


 23  Field=read_civxdata(filename,[],field.vel_type);


 24  if isfield(Field,'Txt')


 25  error=Field.Txt; %error in reading


 26  return


 27  end


 28 


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


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


 31  return


 32  end


 33  if ~isfield(Field,'C')


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


 35  end


 36  if ~isfield(Field,'F')


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


 38  end


 39  if ~isfield(Field,'FF')


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


 41  end


 42 


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


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


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


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


 47  if iter==1


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


 49  elseif iter==2


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


 51  end


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


 53 


 54  % criterium on velocity values


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


 56  delta_v=Field.V;


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


 58  if ~exist(fileref,'file')


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


 60  display(error);


 61  return


 62  end


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


 64  if isfield(FieldRef,'FF')


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


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


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


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


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


 70  end


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


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


 73  return


 74  end


 75  if length(FieldRef.X)<=1


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


 77  return


 78  end


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


 80  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


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


 82  delta_v=Field.Vvec_V_ref;


 83  end


 84 


 85  thresh_vel_x=thresh_vel;


 86  thresh_vel_y=thresh_vel;


 87  if isequal(inf_sup,1)


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


 89  elseif isequal(inf_sup,2)


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


 91  end


 92 


 93  % flag7 introduce a grey mask, matrix M


 94  if isequal (flag_mask,1)


 95  M=imread(maskname);


 96  nxy=size(M);


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


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


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


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


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


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


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


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


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


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


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


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


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


 110  Mvalues=M(ICOMB);


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


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


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


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


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


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


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


 118  Mvalues=M(ICOMB);


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


 120  flag7=(flag7bflag7e)';


 121  else


 122  flag7=0;


 123  end


 124  flagmagenta=flag1flag2flag3flag4flag5flag7;


 125 


 126  %write fix flags


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


 128 


 129  if ~isempty(message) && ~isequal(message.UserWrite,1)


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


 131  return


 132  end


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


 134  result=redef(nc);


 135  if isempty(result), errordlg('##Bad redef operation.'),end


 136  if iter==1


 137  nc.fix=1;


 138  elseif iter==2


 139  nc.fix2=1;


 140  end


 141  %theDim=nc(field.nb) ;% get the number of velocity vectors


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


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


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


 145  close(nc);

