source: trunk/src/RUN_FIX.m @ 67

Last change on this file since 67 was 67, checked in by sommeria, 11 years ago

civ: RUN civ lounched out of the Matlab work space. RUN and BATCH now runned by a unique sub-function lounch.m.
FiLE PARAM.xml modified to provide different paths for Batch and Run !!!!!
RUN_FIX: minor error message modif
geometry_calib: calib point editing by the mouse
improvement of the interactions between the different GUIs. Close function...

File size: 7.7 KB
Line 
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 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)
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
23function error=RUN_FIX(filename,field,flagindex,iter,thresh_vecC,flag_mask,maskname,thresh_vel,inf_sup,fileref,fieldref)
24error=[]; %default
25vel_type{1}=field.vel_type;
26%check writing access
27[errorread,message]=fileattrib(filename);
28if ischar(message)
29    msgbox_uvmat('ERROR',[filename ':' message]);
30    return
31end
32if ~isequal(message.UserWrite,1)
33     msgbox_uvmat('ERROR',['no writting access to ' filename ' (RUN_FIX.m)']);
34    return
35end
36Field=read_civxdata(filename,'ima_cor',field.vel_type);
37if isfield(Field,'Txt')
38    error=Field.Txt; %error in reading
39    return
40end
41
42if ~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
45end
46if ~isfield(Field,'C')
47    Field.C=ones(size(Field.X));%correlation=1 by default
48end
49if ~isfield(Field,'F')
50    Field.F=ones(size(Field.X));%warning flag=1 by default
51end
52if ~isfield(Field,'FF')
53    Field.FF=zeros(size(Field.X));%fixflag=0 by default
54end
55
56vec_f_unit=abs(Field.F)-10*double(uint16(abs(Field.F)/10)); %unityterm of vec_F in abs value
57vec_f_sign=sign(Field.F).*vec_f_unit;% gives the unity digit of vec_f with correct sign
58flag1=(flagindex(1)==1)&(vec_f_sign==-2);%removed vectors vec_f=-2
59flag2=(flagindex(2)==1)&(vec_f_sign==3);%removed vectors vec_f=3
60if iter==1
61        flag3=(flagindex(3)==1)&(vec_f_sign==2); % Hart vectors
62elseif iter==2
63        flag3=(flagindex(3)==1)&(vec_f_sign==4); %
64end
65flag4=(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
68delta_u=Field.U;%default without ref file
69delta_v=Field.V;
70if 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.U-vec_U_ref;%take the difference with the interpolated ref field
95    delta_v=Field.V-vec_V_ref;
96end
97thresh_vel_x=thresh_vel;
98thresh_vel_y=thresh_vel;
99if 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);
101elseif 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);
103end
104
105            % flag7 introduce a grey mask, matrix M
106if 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.X-Field.U/2;%beginning points
113   vec_x2=Field.X+Field.U/2;%end points of vectors
114   vec_y1=Field.Y-Field.V/2;%beginning points
115   vec_y2=Field.Y+Field.V/2;%end points of vectors
116   indx=1+round((nxy(2)-1)*(vec_x1-rangx0(1))/(rangx0(2)-rangx0(1)));% image index x at abcissa vec_x
117   indy=1+round((nxy(1)-1)*(vec_y1-rangy0(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+(1-test_in); %replace indx by 1 out of the mask range
120   indy=indy.*test_in+(1-test_in); %replace indy by 1 out of the mask range
121   ICOMB=((indx-1)*nxy(1)+(nxy(1)+1-indy));%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_x2-rangx0(1))/(rangx0(2)-rangx0(1)));% image index x at abcissa Field.X
125   indy=1+round((nxy(1)-1)*(vec_y2-rangy0(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+(1-test_in); %replace indx by 1 out of the mask range
128   indy=indy.*test_in+(1-test_in); %replace indy by 1 out of the mask range
129   ICOMB=((indx-1)*nxy(1)+(nxy(1)+1-indy));%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=(flag7b|flag7e)';
133else
134   flag7=0;
135end   
136flagmagenta=flag1|flag2|flag3|flag4|flag5|flag7;
137fixflag_unit=Field.FF-10*floor(Field.FF/10); %unity term of fix_flag
138
139%write fix flags in the netcdf file
140hhh=which('netcdf.open');% look for built-in matlab netcdf library
141if ~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)
158else %old netcdf library
159    netcdf_toolbox(filename,Field,flagmagenta,iter,field)
160end
161
162function netcdf_toolbox(filename,Field,flagmagenta,iter,field)
163nc=netcdf(filename,'write'); %open netcdf file for writing
164result=redef(nc);
165if isempty(result), msgbox_uvmat('ERROR','##Bad redef operation.'),end 
166if iter==1
167    nc.fix=1;
168elseif iter==2
169    nc.fix2=1;
170end
171nc{field.fixflag}=ncfloat(field.nb);
172fixflag_unit=Field.FF-10*floor(Field.FF/10); %unity term of fix_flag
173nc{field.fixflag}(:)=fixflag_unit+10*flagmagenta;
174close(nc);
Note: See TracBrowser for help on using the repository browser.