[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)
|
---|
| 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.U-vec_U_ref;%take the difference with the interpolated ref field
|
---|
| 82 | delta_v=Field.V-vec_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.X-Field.U/2;%beginning points
|
---|
| 101 | vec_x2=Field.X+Field.U/2;%end points of vectors
|
---|
| 102 | vec_y1=Field.Y-Field.V/2;%beginning points
|
---|
| 103 | vec_y2=Field.Y+Field.V/2;%end points of vectors
|
---|
| 104 | indx=1+round((nxy(2)-1)*(vec_x1-rangx0(1))/(rangx0(2)-rangx0(1)));% image index x at abcissa vec_x
|
---|
| 105 | indy=1+round((nxy(1)-1)*(vec_y1-rangy0(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+(1-test_in); %replace indx by 1 out of the mask range
|
---|
| 108 | indy=indy.*test_in+(1-test_in); %replace indy by 1 out of the mask range
|
---|
| 109 | ICOMB=((indx-1)*nxy(1)+(nxy(1)+1-indy));%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_x2-rangx0(1))/(rangx0(2)-rangx0(1)));% image index x at abcissa Field.X
|
---|
| 113 | indy=1+round((nxy(1)-1)*(vec_y2-rangy0(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+(1-test_in); %replace indx by 1 out of the mask range
|
---|
| 116 | indy=indy.*test_in+(1-test_in); %replace indy by 1 out of the mask range
|
---|
| 117 | ICOMB=((indx-1)*nxy(1)+(nxy(1)+1-indy));%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=(flag7b|flag7e)';
|
---|
| 121 | else
|
---|
| 122 | flag7=0;
|
---|
| 123 | end
|
---|
| 124 | flagmagenta=flag1|flag2|flag3|flag4|flag5|flag7;
|
---|
| 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.FF-10*floor(Field.FF/10); %unity term of fix_flag
|
---|
| 144 | nc{field.fixflag}(:)=fixflag_unit+10*flagmagenta;
|
---|
| 145 | close(nc);
|
---|