1 | %'proj_field': projects the field on a projection object |
---|
2 | %-------------------------------------------------------------------------- |
---|
3 | % function [ProjData,errormsg]=proj_field(FieldData,ObjectData,VarMesh) |
---|
4 | % |
---|
5 | % OUTPUT: |
---|
6 | % ProjData structure containing the fields of the input field FieldData, |
---|
7 | % transmitted or projected on the object, plus the additional fields |
---|
8 | % .UMax, .UMin, .VMax, .VMin: min and max of velocity components in a domain |
---|
9 | % .UMean,VMean: mean of the velocity components in a domain |
---|
10 | % .AMin, AMax: min and max of a scalar |
---|
11 | % .AMean: mean of a scalar in a domain |
---|
12 | % .NbPix; |
---|
13 | % .DimName= names of the matrix dimensions (matlab cell) |
---|
14 | % .VarName= names of the variables [ProjData.VarName {'A','AMean','AMin','AMax'}]; |
---|
15 | % .VarDimNameIndex= dimensions of the variables, indicated by indices in the list .DimName; |
---|
16 | % |
---|
17 | %INPUT |
---|
18 | % ObjectData: structure characterizing the projection object |
---|
19 | % .Type : type of projection object |
---|
20 | % .ProjMode=mode of projection ; |
---|
21 | % .CoordUnit: 'px', 'cm' units for the coordinates defining the object |
---|
22 | % .Angle ( angles of rotation (=[0 0 0] by default) |
---|
23 | % .ProjAngle=angle of projection; |
---|
24 | % .DX,.DY,.DZ=increments along each coordinate |
---|
25 | % .Coord(nbpoints,3): set of coordinates defining the object position; |
---|
26 | |
---|
27 | %FieldData: data of the field to be projected on the projection object, with optional fields |
---|
28 | % .Txt: error message, transmitted to the projection |
---|
29 | % .FieldList: cell array of strings representing the fields to calculate |
---|
30 | % .CoordMesh: typical distance between data points (used for mouse action or display), transmitted |
---|
31 | % .CoordUnit, .TimeUnit, .dt: transmitted |
---|
32 | % standardised description of fields, nc-formated Matlab structure with fields: |
---|
33 | % .ListGlobalAttribute: cell listing the names of the global attributes |
---|
34 | % .Att_1,Att_2... : values of the global attributes |
---|
35 | % .ListVarName: cell listing the names of the variables |
---|
36 | % .VarAttribute: cell of structures s containing names and values of variable attributes (s.name=value) for each variable of .ListVarName |
---|
37 | % .Var1, .Var2....: variables (Matlab arrays) with names listed in .ListVarName |
---|
38 | % The variables are grouped in 'fields', made of a set of variables with common dimensions (using the function find_field_cells) |
---|
39 | % The variable attribute 'Role' is used to define the role for plotting: |
---|
40 | % Role = 'scalar': (default) represents a scalar field |
---|
41 | % = 'coord': represents a set of unstructured coordinates, whose |
---|
42 | % space dimension is given by the last array dimension (called 'NbDim'). |
---|
43 | % = 'coord_x', 'coord_y', 'coord_z': represents a separate set of |
---|
44 | % unstructured coordinate x, y or z |
---|
45 | % = 'vector': represents a vector field whose number of components |
---|
46 | % is given by the last dimension (called 'NbDim') |
---|
47 | % = 'vector_x', 'vector_y', 'vector_z' :represents the x, y or z component of a vector |
---|
48 | % = 'warnflag' : provides a warning flag about the quality of data in a 'Field', default=0, no warning |
---|
49 | % = 'errorflag': provides an error flag marking false data, |
---|
50 | % default=0, no error. Different non zero values can represent different criteria of elimination. |
---|
51 | % |
---|
52 | % Default role of variables (by name) |
---|
53 | % vector field: |
---|
54 | % .X,.Y: position of the velocity vectors, projected on the object |
---|
55 | % .U, .V, .W: velocity components, projected on the object |
---|
56 | % .C, .CName: scalar associated to the vector |
---|
57 | % .F : equivalent to 'warnflag' |
---|
58 | % .FF: equivalent to 'errorflag' |
---|
59 | % scalar field or image: |
---|
60 | % .AName: name of a scalar (to be calculated from velocity fields after projection), transmitted |
---|
61 | % .A: scalar, projected on the object |
---|
62 | % .AX, .AY: positions for the scalar |
---|
63 | % case of a structured grid: A is a dim 2 matrix and .AX=[first last] (length 2 vector) represents the first and last abscissa of the grid |
---|
64 | % case of an unstructured scalar: A is a vector, AX and AY the corresponding coordinates |
---|
65 | |
---|
66 | %======================================================================= |
---|
67 | % Copyright 2008-2016, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France |
---|
68 | % http://www.legi.grenoble-inp.fr |
---|
69 | % Joel.Sommeria - Joel.Sommeria (A) legi.cnrs.fr |
---|
70 | % |
---|
71 | % This file is part of the toolbox UVMAT. |
---|
72 | % |
---|
73 | % UVMAT is free software; you can redistribute it and/or modify |
---|
74 | % it under the terms of the GNU General Public License as published |
---|
75 | % by the Free Software Foundation; either version 2 of the license, |
---|
76 | % or (at your option) any later version. |
---|
77 | % |
---|
78 | % UVMAT is distributed in the hope that it will be useful, |
---|
79 | % but WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
80 | % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
---|
81 | % GNU General Public License (see LICENSE.txt) for more details. |
---|
82 | %======================================================================= |
---|
83 | |
---|
84 | function [ProjData,errormsg]=proj_field(FieldData,ObjectData,VarMesh) |
---|
85 | errormsg='';%default |
---|
86 | ProjData=[]; |
---|
87 | |
---|
88 | %% check input projection object: type, projection mode and Coord: |
---|
89 | if ~isfield(ObjectData,'Type')||~isfield(ObjectData,'ProjMode') |
---|
90 | return |
---|
91 | end |
---|
92 | % check list of effective projection modes |
---|
93 | if ~ismember(ObjectData.ProjMode,{'projection','interp_lin','interp_tps','inside','outside'}) |
---|
94 | return |
---|
95 | end |
---|
96 | if ~isfield(ObjectData,'Coord')||isempty(ObjectData.Coord) |
---|
97 | if strcmp(ObjectData.Type,'plane') |
---|
98 | ObjectData.Coord=[0 0];%default |
---|
99 | else |
---|
100 | return |
---|
101 | end |
---|
102 | end |
---|
103 | |
---|
104 | %% apply projection depending on the object type |
---|
105 | switch ObjectData.Type |
---|
106 | case 'points' |
---|
107 | [ProjData,errormsg]=proj_points(FieldData,ObjectData); |
---|
108 | case {'line','polyline'} |
---|
109 | [ProjData,errormsg] = proj_line(FieldData,ObjectData); |
---|
110 | case {'polygon','rectangle','ellipse'} |
---|
111 | if isequal(ObjectData.ProjMode,'inside')||isequal(ObjectData.ProjMode,'outside') |
---|
112 | if ~exist('VarMesh','var') |
---|
113 | VarMesh=[]; |
---|
114 | end |
---|
115 | [ProjData,errormsg] = proj_patch(FieldData,ObjectData,VarMesh); |
---|
116 | else |
---|
117 | [ProjData,errormsg] = proj_line(FieldData,ObjectData); |
---|
118 | end |
---|
119 | case {'plane','plane_z'} |
---|
120 | [ProjData,errormsg] = proj_plane(FieldData,ObjectData); |
---|
121 | case 'volume' |
---|
122 | [ProjData,errormsg] = proj_volume(FieldData,ObjectData); |
---|
123 | end |
---|
124 | |
---|
125 | %----------------------------------------------------------------- |
---|
126 | %project on a set of points |
---|
127 | function [ProjData,errormsg]=proj_points(FieldData,ObjectData)%% |
---|
128 | %------------------------------------------------------------------- |
---|
129 | |
---|
130 | siz=size(ObjectData.Coord); |
---|
131 | width=0; |
---|
132 | if isfield(ObjectData,'Range') |
---|
133 | width=ObjectData.Range(1,2); |
---|
134 | end |
---|
135 | if isfield(ObjectData,'RangeX')&&~isempty(ObjectData.RangeX) |
---|
136 | width=max(ObjectData.RangeX); |
---|
137 | end |
---|
138 | if isfield(ObjectData,'RangeY')&&~isempty(ObjectData.RangeY) |
---|
139 | width=max(width,max(ObjectData.RangeY)); |
---|
140 | end |
---|
141 | if isfield(ObjectData,'RangeZ')&&~isempty(ObjectData.RangeZ) |
---|
142 | width=max(width,max(ObjectData.RangeZ)); |
---|
143 | end |
---|
144 | if isequal(ObjectData.ProjMode,'projection') |
---|
145 | if width==0 |
---|
146 | errormsg='projection range around points needed'; |
---|
147 | return |
---|
148 | end |
---|
149 | elseif ~isequal(ObjectData.ProjMode,'interp_lin') |
---|
150 | errormsg=(['ProjMode option ' ObjectData.ProjMode ' not available in proj_field']); |
---|
151 | return |
---|
152 | end |
---|
153 | [ProjData,errormsg]=proj_heading(FieldData,ObjectData); |
---|
154 | if ~isempty(errormsg) |
---|
155 | return |
---|
156 | end |
---|
157 | ProjData.NbDim=0; |
---|
158 | [CellInfo,NbDimArray,errormsg]=find_field_cells(FieldData); |
---|
159 | if ~isempty(errormsg) |
---|
160 | errormsg=['error in proj_field/proj_points:' errormsg]; |
---|
161 | return |
---|
162 | end |
---|
163 | %LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS |
---|
164 | for icell=1:length(CellInfo) |
---|
165 | if NbDimArray(icell)<=1 |
---|
166 | continue %projection only for multidimensional fields |
---|
167 | end |
---|
168 | VarIndex=CellInfo{icell}.VarIndex;% indices of the selected variables in the list FieldData.ListVarName |
---|
169 | ivar_X=CellInfo{icell}.CoordIndex(end); |
---|
170 | ivar_Y=CellInfo{icell}.CoordIndex(end-1); |
---|
171 | ivar_Z=[]; |
---|
172 | if NbDimArray(icell)==3 |
---|
173 | ivar_Z=CellInfo{icell}.CoordIndex(1); |
---|
174 | end |
---|
175 | ivar_FF=[]; |
---|
176 | if isfield(CellInfo{icell},'VarIndex_errorflag') |
---|
177 | ivar_FF=CellInfo{icell}.VarIndex_errorflag; |
---|
178 | if numel(ivar_FF)>1 |
---|
179 | errormsg='multiple error flag input'; |
---|
180 | return |
---|
181 | end |
---|
182 | end |
---|
183 | % select types of variables to be projected |
---|
184 | ListProj={'VarIndex_scalar','VarIndex_image','VarIndex_color','VarIndex_vector_x','VarIndex_vector_y'}; |
---|
185 | check_proj=false(size(FieldData.ListVarName)); |
---|
186 | for ilist=1:numel(ListProj) |
---|
187 | if isfield(CellInfo{icell},ListProj{ilist}) |
---|
188 | check_proj(CellInfo{icell}.(ListProj{ilist}))=1; |
---|
189 | end |
---|
190 | end |
---|
191 | VarIndex=find(check_proj); |
---|
192 | ProjData.ListVarName={'Y','X','NbVal'}; |
---|
193 | ProjData.VarDimName={'nb_points','nb_points','nb_points'}; |
---|
194 | ProjData.VarAttribute{1}.Role='ancillary'; |
---|
195 | ProjData.VarAttribute{2}.Role='ancillary'; |
---|
196 | ProjData.VarAttribute{3}.Role='ancillary'; |
---|
197 | for ivar=VarIndex |
---|
198 | VarName=FieldData.ListVarName{ivar}; |
---|
199 | ProjData.ListVarName=[ProjData.ListVarName {VarName}];% add the current variable to the list of projected variables |
---|
200 | ProjData.VarDimName=[ProjData.VarDimName {'nb_points'}]; % projected VarName has a single dimension called 'nb_points' (set of projection points) |
---|
201 | |
---|
202 | end |
---|
203 | if strcmp( CellInfo{icell}.CoordType,'scattered') |
---|
204 | coord_x=FieldData.(FieldData.ListVarName{ivar_X}); |
---|
205 | coord_y=FieldData.(FieldData.ListVarName{ivar_Y}); |
---|
206 | test3D=0;% TEST 3D CASE : NOT COMPLETED , 3D CASE : NOT COMPLETED |
---|
207 | if length(ivar_Z)==1 |
---|
208 | coord_z=FieldData.(FieldData.ListVarName{ivar_Z}); |
---|
209 | test3D=1; |
---|
210 | end |
---|
211 | |
---|
212 | for ipoint=1:siz(1) |
---|
213 | Xpoint=ObjectData.Coord(ipoint,:); |
---|
214 | distX=coord_x-Xpoint(1); |
---|
215 | distY=coord_y-Xpoint(2); |
---|
216 | dist=distX.*distX+distY.*distY; |
---|
217 | indsel=find(dist<width*width); |
---|
218 | ProjData.X(ipoint,1)=Xpoint(1); |
---|
219 | ProjData.Y(ipoint,1)=Xpoint(2); |
---|
220 | if isequal(length(ivar_FF),1) |
---|
221 | FFName=FieldData.ListVarName{ivar_FF}; |
---|
222 | FF=FieldData.(FFName)(indsel); |
---|
223 | indsel=indsel(~FF); |
---|
224 | end |
---|
225 | ProjData.NbVal(ipoint,1)=length(indsel); |
---|
226 | for ivar=VarIndex |
---|
227 | VarName=FieldData.ListVarName{ivar}; |
---|
228 | if isempty(indsel) |
---|
229 | ProjData.(VarName)(ipoint,1)=NaN; |
---|
230 | else |
---|
231 | Var=FieldData.(VarName)(indsel); |
---|
232 | ProjData.(VarName)(ipoint,1)=mean(Var); |
---|
233 | if isequal(ObjectData.ProjMode,'interp_lin') |
---|
234 | ProjData.(VarName)(ipoint,1)=griddata_uvmat(coord_x(indsel),coord_y(indsel),Var,Xpoint(1),Xpoint(2)); |
---|
235 | end |
---|
236 | end |
---|
237 | end |
---|
238 | end |
---|
239 | else %case of structured coordinates |
---|
240 | if strcmp( CellInfo{icell}.CoordType,'grid') |
---|
241 | AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)}; |
---|
242 | AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)}; |
---|
243 | eval(['AX=FieldData.' AXName ';']);% set of x positions |
---|
244 | eval(['AY=FieldData.' AYName ';']);% set of y positions |
---|
245 | AName=FieldData.ListVarName{VarIndex(1)};% a single variable assumed in the current cell |
---|
246 | eval(['A=FieldData.' AName ';']);% scalar |
---|
247 | npxy=size(A); |
---|
248 | %update VarDimName in case of components (non coordinate dimensions e;g. color components) |
---|
249 | if numel(npxy)>NbDimArray(icell) |
---|
250 | ProjData.VarDimName{end}={'nb_points','component'}; |
---|
251 | end |
---|
252 | for idim=1:NbDimArray(icell) %loop on space dimensions |
---|
253 | test_interp(idim)=0;%test for coordiate interpolation (non regular grid), =0 by default |
---|
254 | test_coord(idim)=0;%test for defined coordinates, =0 by default |
---|
255 | ivar=CellInfo{icell}.CoordIndex(idim); |
---|
256 | Coord{idim}=FieldData.(FieldData.ListVarName{ivar}); % position for the first index |
---|
257 | if numel(Coord{idim})==2 |
---|
258 | DCoord_min(idim)= (Coord{idim}(2)-Coord{idim}(1))/(npxy(idim)-1); |
---|
259 | test_direct(idim)=DCoord_min(idim)>0;% =1 for increasing values, 0 otherwise |
---|
260 | else |
---|
261 | DCoord=diff(Coord{idim}); |
---|
262 | DCoord_min(idim)=min(DCoord); |
---|
263 | DCoord_max=max(DCoord); |
---|
264 | test_direct(idim)=DCoord_max>0;% =1 for increasing values, 0 otherwise |
---|
265 | test_direct_min=DCoord_min(idim)>0;% =1 for increasing values, 0 otherwise |
---|
266 | if ~isequal(test_direct(idim),test_direct_min) |
---|
267 | errormsg=['non monotonic dimension variable # ' num2str(idim) ' in proj_field.m']; |
---|
268 | return |
---|
269 | end |
---|
270 | test_interp(idim)=(DCoord_max-DCoord_min(idim))> 0.0001*abs(DCoord_max);% test grid regularity |
---|
271 | test_coord(idim)=1; |
---|
272 | end |
---|
273 | end |
---|
274 | DX=DCoord_min(2); |
---|
275 | DY=DCoord_min(1); |
---|
276 | for ipoint=1:siz(1) |
---|
277 | xwidth=width/(abs(DX)); |
---|
278 | ywidth=width/(abs(DY)); |
---|
279 | i_min=round((ObjectData.Coord(ipoint,1)-Coord{2}(1))/DX+0.5-xwidth); %minimum index of the selected region |
---|
280 | i_min=max(1,i_min);%restrict to field limit |
---|
281 | i_plus=round((ObjectData.Coord(ipoint,1)-Coord{2}(1))/DX+0.5+xwidth); |
---|
282 | i_plus=min(npxy(2),i_plus); %restrict to field limit |
---|
283 | j_min=round((ObjectData.Coord(ipoint,2)-Coord{1}(1))/DY-ywidth+0.5); |
---|
284 | j_min=max(1,j_min); |
---|
285 | j_plus=round((ObjectData.Coord(ipoint,2)-Coord{1}(1))/DY+ywidth+0.5); |
---|
286 | j_plus=min(npxy(1),j_plus); |
---|
287 | ProjData.X(ipoint,1)=ObjectData.Coord(ipoint,1); |
---|
288 | ProjData.Y(ipoint,1)=ObjectData.Coord(ipoint,2); |
---|
289 | i_int=(i_min:i_plus); |
---|
290 | j_int=(j_min:j_plus); |
---|
291 | ProjData.NbVal(ipoint,1)=length(j_int)*length(i_int); |
---|
292 | if isempty(i_int) || isempty(j_int) |
---|
293 | for ivar=VarIndex |
---|
294 | eval(['ProjData.' FieldData.ListVarName{ivar} '(ipoint,:)=NaN;']); |
---|
295 | end |
---|
296 | errormsg=['no data points in the selected projection range ' num2str(width) ]; |
---|
297 | else |
---|
298 | %TODO: introduce circle in the selected subregion |
---|
299 | %[I,J]=meshgrid([1:j_int],[1:i_int]); |
---|
300 | for ivar=VarIndex |
---|
301 | Avalue=FieldData.(FieldData.ListVarName{ivar})(j_int,i_int,:); |
---|
302 | ProjData.(FieldData.ListVarName{ivar})(ipoint,:)=mean(mean(Avalue)); |
---|
303 | end |
---|
304 | end |
---|
305 | end |
---|
306 | end |
---|
307 | end |
---|
308 | end |
---|
309 | |
---|
310 | %----------------------------------------------------------------- |
---|
311 | %project in a patch |
---|
312 | function [ProjData,errormsg]=proj_patch(FieldData,ObjectData,VarMesh)%% |
---|
313 | %------------------------------------------------------------------- |
---|
314 | [ProjData,errormsg]=proj_heading(FieldData,ObjectData); |
---|
315 | if ~isempty(errormsg) |
---|
316 | return |
---|
317 | end |
---|
318 | %objectfield=fieldnames(ObjectData); |
---|
319 | widthx=0; |
---|
320 | widthy=0; |
---|
321 | if isfield(ObjectData,'RangeX') && ~isempty(ObjectData.RangeX) |
---|
322 | widthx=max(ObjectData.RangeX); |
---|
323 | end |
---|
324 | if isfield(ObjectData,'RangeY') && ~isempty(ObjectData.RangeY) |
---|
325 | widthy=max(ObjectData.RangeY); |
---|
326 | end |
---|
327 | |
---|
328 | %A REVOIR, GENERALISER: UTILISER proj_line |
---|
329 | ProjData.NbDim=1; |
---|
330 | ProjData.ListVarName={}; |
---|
331 | ProjData.VarDimName={}; |
---|
332 | ProjData.VarAttribute={}; |
---|
333 | |
---|
334 | CoordMesh=zeros(1,numel(FieldData.ListVarName)); |
---|
335 | if isfield (FieldData,'VarAttribute') |
---|
336 | for iattr=1:length(FieldData.VarAttribute)%initialization of variable attribute values |
---|
337 | if isfield(FieldData.VarAttribute{iattr},'Unit') |
---|
338 | unit{iattr}=FieldData.VarAttribute{iattr}.Unit; |
---|
339 | end |
---|
340 | if isfield(FieldData.VarAttribute{iattr},'CoordMesh') |
---|
341 | CoordMesh(iattr)=FieldData.VarAttribute{iattr}.CoordMesh; |
---|
342 | end |
---|
343 | end |
---|
344 | end |
---|
345 | |
---|
346 | %group the variables (fields of 'FieldData') in cells of variables with the same dimensions |
---|
347 | [CellInfo,NbDim,errormsg]=find_field_cells(FieldData); |
---|
348 | if ~isempty(errormsg) |
---|
349 | errormsg=['error in proj_field/proj_patch:' errormsg]; |
---|
350 | return |
---|
351 | end |
---|
352 | |
---|
353 | %LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS |
---|
354 | for icell=1:length(CellInfo) |
---|
355 | CoordType=CellInfo{icell}.CoordType; |
---|
356 | test_Amat=0; |
---|
357 | if NbDim(icell)~=2% proj_patch acts only on fields of space dimension 2 |
---|
358 | continue |
---|
359 | end |
---|
360 | ivar_FF=[]; |
---|
361 | testfalse=isfield(CellInfo{icell},'VarIndex_errorflag'); |
---|
362 | if testfalse |
---|
363 | ivar_FF=CellInfo{icell}.VarIndex_errorflag; |
---|
364 | FFName=FieldData.ListVarName{ivar_FF}; |
---|
365 | errorflag=FieldData.(FFName); |
---|
366 | end |
---|
367 | % select types of variables to be projected |
---|
368 | ListProj={'VarIndex_scalar','VarIndex_image','VarIndex_color','VarIndex_vector_x','VarIndex_vector_y'}; |
---|
369 | check_proj=false(size(FieldData.ListVarName)); |
---|
370 | for ilist=1:numel(ListProj) |
---|
371 | if isfield(CellInfo{icell},ListProj{ilist}) |
---|
372 | check_proj(CellInfo{icell}.(ListProj{ilist}))=1; |
---|
373 | end |
---|
374 | end |
---|
375 | VarIndex=find(check_proj); |
---|
376 | |
---|
377 | ivar_X=CellInfo{icell}.CoordIndex(end); |
---|
378 | ivar_Y=CellInfo{icell}.CoordIndex(end-1); |
---|
379 | ivar_Z=[]; |
---|
380 | if NbDim(icell)==3 |
---|
381 | ivar_Z=CellInfo{icell}.CoordIndex(1); |
---|
382 | end |
---|
383 | switch CellInfo{icell}.CoordType |
---|
384 | case 'scattered' %case of unstructured coordinates |
---|
385 | for ivar=[VarIndex ivar_X ivar_Y ivar_FF] |
---|
386 | VarName=FieldData.ListVarName{ivar}; |
---|
387 | FieldData.(VarName)=reshape(FieldData.(VarName),[],1); |
---|
388 | end |
---|
389 | XName=FieldData.ListVarName{ivar_X}; |
---|
390 | YName=FieldData.ListVarName{ivar_Y}; |
---|
391 | coord_x=FieldData.(FieldData.ListVarName{ivar_X}); |
---|
392 | coord_y=FieldData.(FieldData.ListVarName{ivar_Y}); |
---|
393 | % image or 2D matrix |
---|
394 | case 'grid' %case of structured coordinates |
---|
395 | test_Amat=1;% test for image or 2D matrix |
---|
396 | AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)}; |
---|
397 | AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)}; |
---|
398 | AX=FieldData.(AXName);% x coordinate |
---|
399 | AY=FieldData.(AYName);% y coordinate |
---|
400 | VarName=FieldData.ListVarName{VarIndex(1)}; |
---|
401 | DimValue=size(FieldData.(VarName)); |
---|
402 | if length(AX)==2 |
---|
403 | AX=linspace(AX(1),AX(end),DimValue(2)); |
---|
404 | end |
---|
405 | if length(AY)==2 |
---|
406 | AY=linspace(AY(1),AY(end),DimValue(1)); |
---|
407 | end |
---|
408 | if length(DimValue)==3 |
---|
409 | testcolor=1; |
---|
410 | npxy(3)=3; |
---|
411 | else |
---|
412 | testcolor=0; |
---|
413 | npxy(3)=1; |
---|
414 | end |
---|
415 | [Xi,Yi]=meshgrid(AX,AY); |
---|
416 | npxy(1)=length(AY); |
---|
417 | npxy(2)=length(AX); |
---|
418 | Xi=reshape(Xi,npxy(1)*npxy(2),1); |
---|
419 | Yi=reshape(Yi,npxy(1)*npxy(2),1); |
---|
420 | for ivar=1:length(VarIndex) |
---|
421 | VarName=FieldData.ListVarName{VarIndex(ivar)}; |
---|
422 | FieldData.(VarName)=reshape(FieldData.(VarName),npxy(1)*npxy(2),npxy(3)); % keep only non false vectors |
---|
423 | end |
---|
424 | end |
---|
425 | %select the indices in the range of action |
---|
426 | testin=[];%default |
---|
427 | if isequal(ObjectData.Type,'rectangle') |
---|
428 | if strcmp(CellInfo{icell}.CoordType,'scattered') |
---|
429 | distX=abs(coord_x-ObjectData.Coord(1,1)); |
---|
430 | distY=abs(coord_y-ObjectData.Coord(1,2)); |
---|
431 | testin=distX<widthx & distY<widthy; |
---|
432 | elseif test_Amat |
---|
433 | distX=abs(Xi-ObjectData.Coord(1,1)); |
---|
434 | distY=abs(Yi-ObjectData.Coord(1,2)); |
---|
435 | testin=distX<widthx & distY<widthy; |
---|
436 | end |
---|
437 | elseif isequal(ObjectData.Type,'polygon') |
---|
438 | if strcmp(CoordType,'scattered') |
---|
439 | testin=inpolygon(coord_x,coord_y,ObjectData.Coord(:,1),ObjectData.Coord(:,2)); |
---|
440 | elseif strcmp(CoordType,'grid') |
---|
441 | testin=inpolygon(Xi,Yi,ObjectData.Coord(:,1),ObjectData.Coord(:,2)); |
---|
442 | else%calculate the scalar |
---|
443 | testin=[]; %A REVOIR |
---|
444 | end |
---|
445 | elseif isequal(ObjectData.Type,'ellipse') |
---|
446 | X2Max=widthx*widthx; |
---|
447 | Y2Max=(widthy)*(widthy); |
---|
448 | if strcmp(CoordType,'scattered') |
---|
449 | distX=(coord_x-ObjectData.Coord(1,1)); |
---|
450 | distY=(coord_y-ObjectData.Coord(1,2)); |
---|
451 | testin=(distX.*distX/X2Max+distY.*distY/Y2Max)<1; |
---|
452 | elseif strcmp(CoordType,'grid') %case of usual 2x2 matrix |
---|
453 | distX=(Xi-ObjectData.Coord(1,1)); |
---|
454 | distY=(Yi-ObjectData.Coord(1,2)); |
---|
455 | testin=(distX.*distX/X2Max+distY.*distY/Y2Max)<1; |
---|
456 | end |
---|
457 | end |
---|
458 | %selected indices |
---|
459 | if isequal(ObjectData.ProjMode,'outside') |
---|
460 | testin=~testin; |
---|
461 | end |
---|
462 | if testfalse |
---|
463 | testin=testin & (errorflag==0); % keep only non false vectors |
---|
464 | end |
---|
465 | indsel=find(testin); |
---|
466 | nbvar=0; |
---|
467 | VarSize=zeros(size(VarIndex)); |
---|
468 | for ivar=VarIndex |
---|
469 | VarName=FieldData.ListVarName{ivar}; |
---|
470 | ProjData.([VarName 'Mean'])=mean(double(FieldData.(VarName)(indsel,:))); % take the mean in the selected region, for each co |
---|
471 | ProjData.([VarName 'Min'])=min(double(FieldData.(VarName)(indsel,:))); % take the min in the selected region , for each color component |
---|
472 | ProjData.([VarName 'Max'])=max(double(FieldData.(VarName)(indsel,:))); % take the max in the selected region , for each color co |
---|
473 | nbvar=nbvar+1; |
---|
474 | VarSize(nbvar)=mean((ProjData.([VarName 'Max'])-ProjData.([VarName 'Min']))/100); |
---|
475 | end |
---|
476 | if isempty(VarMesh)% || isnan(VarMesh) % mesh not specified as input, estimate from the bounds |
---|
477 | VarMesh=mean(VarSize); |
---|
478 | ord=10^(floor(log10(VarMesh)));%order of magnitude |
---|
479 | if VarMesh/ord >=5 |
---|
480 | VarMesh=5*ord; |
---|
481 | elseif VarMesh/ord >=2 |
---|
482 | VarMesh=2*ord; |
---|
483 | else |
---|
484 | VarMesh=ord; |
---|
485 | end |
---|
486 | end |
---|
487 | for ivar=VarIndex |
---|
488 | VarName=FieldData.ListVarName{ivar}; |
---|
489 | LowBound=VarMesh*ceil(ProjData.([VarName 'Min'])/VarMesh); |
---|
490 | UpperBound=VarMesh*floor(ProjData.([VarName 'Max'])/VarMesh); |
---|
491 | if numel(indsel)<=1 |
---|
492 | errormsg='only one data point or less for histogram'; |
---|
493 | return |
---|
494 | elseif isequal(LowBound,UpperBound) |
---|
495 | errormsg='attempt histogram of uniform field: low bound = high bound'; |
---|
496 | return |
---|
497 | end |
---|
498 | ProjData.(VarName)=LowBound:VarMesh:UpperBound; % list of bin values |
---|
499 | ProjData.([VarName 'Histo'])=hist(double(FieldData.(VarName)(indsel,:)),ProjData.(VarName)); % histogram at predefined bin positions |
---|
500 | ProjData.ListVarName=[ProjData.ListVarName {VarName} {[VarName 'Histo']} {[VarName 'Mean']} {[VarName 'Min']} {[VarName 'Max']}]; |
---|
501 | if test_Amat && testcolor |
---|
502 | ProjData.VarDimName=[ProjData.VarDimName {VarName} {{VarName,'rgb'}} {'rgb'} {'rgb'} {'rgb'}];%{{'nb_point','rgb'}}; |
---|
503 | else |
---|
504 | ProjData.VarDimName=[ProjData.VarDimName {VarName} {VarName} {'one'} {'one'} {'one'}]; |
---|
505 | end |
---|
506 | VarAttribute_var=[]; |
---|
507 | if isfield(FieldData,'VarAttribute')&& numel(FieldData.VarAttribute)>=ivar |
---|
508 | VarAttribute_var=FieldData.VarAttribute{ivar}; |
---|
509 | end |
---|
510 | % VarAttribute_var.Role='coord_x';% the variable is now used as an absissa |
---|
511 | VarAttribute_histo.Role='histo'; |
---|
512 | ProjData.VarAttribute=[ProjData.VarAttribute {VarAttribute_var} {VarAttribute_histo} {[]} {[]} {[]}]; |
---|
513 | end |
---|
514 | end |
---|
515 | |
---|
516 | %----------------------------------------------------------------- |
---|
517 | %project on a line |
---|
518 | % AJOUTER flux,circul,error |
---|
519 | % OUTPUT: |
---|
520 | % ProjData: projected field |
---|
521 | % |
---|
522 | function [ProjData,errormsg] = proj_line(FieldData, ObjectData) |
---|
523 | %----------------------------------------------------------------- |
---|
524 | |
---|
525 | %% prepare heading for the projected field |
---|
526 | [ProjData,errormsg]=proj_heading(FieldData,ObjectData);%transfer global attributes |
---|
527 | if ~isempty(errormsg) |
---|
528 | return |
---|
529 | end |
---|
530 | ProjData.NbDim=1; |
---|
531 | %initialisation of the input parameters and defaultoutput |
---|
532 | ProjMode=ObjectData.ProjMode; %rmq: ProjMode always defined from input={'projection','interp_lin','interp_tps'} |
---|
533 | % ProjAngle=90; %90 degrees projection by default |
---|
534 | width=0; |
---|
535 | if isfield(ObjectData,'RangeY') |
---|
536 | width=max(ObjectData.RangeY);%Rangey needed bfor mode 'projection' |
---|
537 | end |
---|
538 | % default output |
---|
539 | errormsg='';%default |
---|
540 | Xline=[]; |
---|
541 | flux=0; |
---|
542 | circul=0; |
---|
543 | liny=ObjectData.Coord(:,2); |
---|
544 | NbPoints=size(ObjectData.Coord,1); |
---|
545 | testfalse=0; |
---|
546 | ListIndex={}; |
---|
547 | |
---|
548 | %% group the variables (fields of 'FieldData') in cells of variables with the same dimensions |
---|
549 | [CellInfo,NbDim,errormsg]=find_field_cells(FieldData); |
---|
550 | if ~isempty(errormsg) |
---|
551 | errormsg=['error in proj_field/proj_line:' errormsg]; |
---|
552 | return |
---|
553 | end |
---|
554 | CellInfo=CellInfo(NbDim>=2); %keep only the 2D cells |
---|
555 | %%%%%% TODO: treat 1D fields: project as identity so that P o P=P for projection operation |
---|
556 | cell_select=true(size(CellInfo)); |
---|
557 | |
---|
558 | for icell=1:length(CellInfo) |
---|
559 | if isfield(CellInfo{icell},'ProjModeRequest') |
---|
560 | if ~strcmp(CellInfo{icell}.ProjModeRequest, ProjMode) |
---|
561 | cell_select(icell)=0; |
---|
562 | end |
---|
563 | if strcmp(ProjMode,'interp_tps')&& ~strcmp(CellInfo{icell}.CoordType,'tps') |
---|
564 | cell_select(icell)=0; |
---|
565 | end |
---|
566 | end |
---|
567 | end |
---|
568 | if isempty(find(cell_select)) |
---|
569 | errormsg=[' invalid projection mode ''' ProjMode ''': use ''interp_tps'' to interpolate spatial derivatives']; |
---|
570 | return |
---|
571 | end |
---|
572 | CellInfo=CellInfo(cell_select); |
---|
573 | |
---|
574 | %% projection line: object types selected from proj_field='line','polyline','polygon','rectangle','ellipse': |
---|
575 | LineCoord=ObjectData.Coord; |
---|
576 | switch ObjectData.Type |
---|
577 | case 'ellipse' |
---|
578 | LineLength=2*pi*ObjectData.RangeX*ObjectData.RangeY; |
---|
579 | NbSegment=0; |
---|
580 | case 'rectangle' |
---|
581 | LineCoord([1 4],1)=ObjectData.Coord(1,1)-ObjectData.RangeX; |
---|
582 | LineCoord([1 2],2)=ObjectData.Coord(1,2)-ObjectData.RangeY; |
---|
583 | LineCoord([2 3],1)=ObjectData.Coord(1,1)+ObjectData.RangeX; |
---|
584 | LineCoord([4 1],2)=ObjectData.Coord(1,2)+ObjectData.RangeY; |
---|
585 | case 'polygon' |
---|
586 | LineCoord(NbPoints+1)=LineCoord(1); |
---|
587 | end |
---|
588 | if ~strcmp(ObjectData.Type,'ellipse') |
---|
589 | if ~strcmp(ObjectData.Type,'rectangle') && NbPoints<2 |
---|
590 | return% line needs at least 2 points to be defined |
---|
591 | end |
---|
592 | dlinx=diff(LineCoord(:,1)); |
---|
593 | dliny=diff(LineCoord(:,2)); |
---|
594 | [theta,dlength]=cart2pol(dlinx,dliny);%angle and length of each segment |
---|
595 | LineLength=sum(dlength); |
---|
596 | NbSegment=numel(LineLength); |
---|
597 | end |
---|
598 | CheckClosedLine=~isempty(find(strcmp(ObjectData.Type,{'rectangle','ellipse','polygon'}))); |
---|
599 | |
---|
600 | % x = a \ \cosh \mu \ \cos \nu |
---|
601 | % |
---|
602 | % y = a \ \sinh \mu \ \sin \nu |
---|
603 | |
---|
604 | %% angles of the polyline and boundaries of action for mode 'projection' |
---|
605 | |
---|
606 | % determine a rectangles at +-width from the line (only used for the ProjMode='projection or 'interp_tps') |
---|
607 | xsup=zeros(1,NbPoints); xinf=zeros(1,NbPoints); ysup=zeros(1,NbPoints); yinf=zeros(1,NbPoints); |
---|
608 | if isequal(ProjMode,'projection') |
---|
609 | if strcmp(ObjectData.Type,'line') |
---|
610 | xsup=ObjectData.Coord(:,1)-width*sin(theta); |
---|
611 | xinf=ObjectData.Coord(:,1)+width*sin(theta); |
---|
612 | ysup=ObjectData.Coord(:,2)+width*cos(theta); |
---|
613 | yinf=ObjectData.Coord(:,2)-width*cos(theta); |
---|
614 | else |
---|
615 | errormsg='mode projection only available for simple line, use interpolation otherwise'; |
---|
616 | return |
---|
617 | end |
---|
618 | else % need to define the set of interpolation points |
---|
619 | if isfield(ObjectData,'DX') && ~isempty(ObjectData.DX) |
---|
620 | DX=abs(ObjectData.DX);%mesh of interpolation points along the line |
---|
621 | if CheckClosedLine |
---|
622 | NbPoint=ceil(LineLength/DX); |
---|
623 | DX=LineLength/NbPoint;%adjust DX to get an integer nbre of intervals in a closed line |
---|
624 | DX_edge=DX/2; |
---|
625 | else |
---|
626 | DX_edge=(LineLength-DX*floor(LineLength/DX))/2;%margin from the first point and first interpolation point, the same for the end point |
---|
627 | end |
---|
628 | XI=[]; |
---|
629 | YI=[]; |
---|
630 | ThetaI=[]; |
---|
631 | dlengthI=[]; |
---|
632 | if strcmp(ObjectData.Type,'ellipse') |
---|
633 | phi=(DX_edge:DX:LineLength)*2*pi/LineLength; |
---|
634 | XI=ObjectData.RangeX*cos(phi); |
---|
635 | YI=ObjectData.RangeY*sin(phi); |
---|
636 | dphi=2*pi*DX/LineLength; |
---|
637 | [ThetaI,dlengthI]=cart2pol(-ObjectData.RangeX*sin(phi)*dphi,ObjectData.RangeY*cos(phi)*dphi); |
---|
638 | else |
---|
639 | for isegment=1:NbSegment |
---|
640 | costheta=cos(theta(isegment)); |
---|
641 | sintheta=sin(theta(isegment)); |
---|
642 | % XIsegment=LineCoord(isegment,1)+DX_edge*costheta:DX*costheta:LineCoord(isegment+1,1)); |
---|
643 | % YIsegment=(LineCoord(isegment,2)+DX_edge*sintheta:DX*sintheta:LineCoord(isegment+1,2)); |
---|
644 | NbInterval=floor((dlength(isegment)-DX_edge)/DX); |
---|
645 | LastX=DX_edge+DX*NbInterval; |
---|
646 | NbPoint=NbInterval+1; |
---|
647 | XIsegment=linspace(LineCoord(isegment,1)+DX_edge*costheta,LineCoord(isegment,1)+LastX*costheta,NbPoint); |
---|
648 | YIsegment=linspace(LineCoord(isegment,2)+DX_edge*sintheta,LineCoord(isegment,2)+LastX*sintheta,NbPoint); |
---|
649 | XI=[XI XIsegment]; |
---|
650 | YI=[YI YIsegment]; |
---|
651 | ThetaI=[ThetaI theta(isegment)*ones(1,numel(XIsegment))]; |
---|
652 | dlengthI=[dlengthI DX*ones(1,numel(XIsegment))]; |
---|
653 | DX_edge=DX-(dlength(isegment)-LastX);%edge for the next segment set to keep DX=DX_end+DX_edge between two segments |
---|
654 | end |
---|
655 | end |
---|
656 | Xproj=cumsum(dlengthI); |
---|
657 | else |
---|
658 | errormsg='abscissa mesh along line DX needed for interpolation'; |
---|
659 | return |
---|
660 | end |
---|
661 | end |
---|
662 | |
---|
663 | %% loop on variable cells with the same space dimension 2 |
---|
664 | ProjData.ListVarName={}; |
---|
665 | ProjData.VarDimName={}; |
---|
666 | check_abscissa=0; |
---|
667 | for icell=1:length(CellInfo) |
---|
668 | % list of variable types to be projected |
---|
669 | ListProj={'VarIndex_scalar','VarIndex_image','VarIndex_color','VarIndex_vector_x','VarIndex_vector_y'}; |
---|
670 | check_proj=false(size(FieldData.ListVarName)); |
---|
671 | for ilist=1:numel(ListProj) |
---|
672 | if isfield(CellInfo{icell},ListProj{ilist}) |
---|
673 | check_proj(CellInfo{icell}.(ListProj{ilist}))=1; |
---|
674 | end |
---|
675 | end |
---|
676 | VarIndex=find(check_proj);% indices of the variables to be projected |
---|
677 | |
---|
678 | %identify error flag |
---|
679 | errorflag=0; %default, no error flag |
---|
680 | if isfield(CellInfo{icell},'VarIndex_errorflag');% test for error flag |
---|
681 | FFName=FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag}; |
---|
682 | errorflag=FieldData.(FFName); |
---|
683 | end |
---|
684 | VarName=FieldData.ListVarName(VarIndex);% cell array of the names of variables to pje |
---|
685 | ivar_U=[]; |
---|
686 | ivar_V=[]; |
---|
687 | %% check needed object properties for unstructured positions (position given by the variables with role coord_x, coord_y |
---|
688 | |
---|
689 | % circul=0; |
---|
690 | % flux=0; |
---|
691 | %%%%%%% % A FAIRE CALCULER MEAN DES QUANTITES %%%%%% |
---|
692 | switch CellInfo{icell}.CoordType |
---|
693 | %case of unstructured coordinates |
---|
694 | case 'scattered' |
---|
695 | % XName= FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)}; |
---|
696 | % YName= FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)}; |
---|
697 | coord_x=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)}); |
---|
698 | coord_y=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)}); |
---|
699 | |
---|
700 | if isequal(ProjMode,'projection') |
---|
701 | if width==0 |
---|
702 | errormsg='range of the projection object is missing'; |
---|
703 | return |
---|
704 | end |
---|
705 | ProjData.ListVarName=[ProjData.ListVarName FieldData.ListVarName(CellInfo{icell}.CoordIndex(end))]; |
---|
706 | ProjData.VarDimName=[ProjData.VarDimName FieldData.ListVarName(CellInfo{icell}.CoordIndex(end))]; |
---|
707 | nbvar=numel(ProjData.ListVarName); |
---|
708 | ProjData.VarAttribute{nbvar}.long_name='abscissa along line'; |
---|
709 | % select the (non false) input data located in the band of projection |
---|
710 | flagsel=(errorflag==0) & ((coord_y -yinf(1))*(xinf(2)-xinf(1))>(coord_x-xinf(1))*(yinf(2)-yinf(1))) ... |
---|
711 | & ((coord_y -ysup(1))*(xsup(2)-xsup(1))<(coord_x-xsup(1))*(ysup(2)-ysup(1))) ... |
---|
712 | & ((coord_y -yinf(2))*(xsup(2)-xinf(2))>(coord_x-xinf(2))*(ysup(2)-yinf(2))) ... |
---|
713 | & ((coord_y -yinf(1))*(xsup(1)-xinf(1))<(coord_x-xinf(1))*(ysup(1)-yinf(1))); |
---|
714 | coord_x=coord_x(flagsel); |
---|
715 | coord_y=coord_y(flagsel); |
---|
716 | costheta=cos(theta); |
---|
717 | sintheta=sin(theta); |
---|
718 | Xproj=(coord_x-ObjectData.Coord(1,1))*costheta + (coord_y-ObjectData.Coord(1,2))*sintheta; %projection on the line |
---|
719 | [Xproj,indsort]=sort(Xproj);% sort points by increasing absissa along the projection line |
---|
720 | ProjData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)})=Xproj; |
---|
721 | for ivar=1:numel(VarIndex) |
---|
722 | ProjData.(VarName{ivar})=FieldData.(VarName{ivar})(flagsel);% restrict variables to the projection band |
---|
723 | ProjData.(VarName{ivar})=ProjData.(VarName{ivar})(indsort);% sort by absissa |
---|
724 | ProjData.ListVarName=[ProjData.ListVarName VarName{ivar}]; |
---|
725 | ProjData.VarDimName=[ProjData.VarDimName FieldData.ListVarName(CellInfo{icell}.CoordIndex(end))]; |
---|
726 | ProjData.VarAttribute{nbvar+ivar}=FieldData.VarAttribute{VarIndex(ivar)};%reproduce var attribute |
---|
727 | if isfield(ProjData.VarAttribute{nbvar+ivar},'Role') |
---|
728 | if strcmp(ProjData.VarAttribute{nbvar+ivar}.Role,'vector_x'); |
---|
729 | ivar_U=nbvar+ivar; |
---|
730 | elseif strcmp(ProjData.VarAttribute{nbvar+ivar}.Role,'vector_y'); |
---|
731 | ivar_V=nbvar+ivar; |
---|
732 | end |
---|
733 | end |
---|
734 | ProjData.VarAttribute{ivar+nbvar}.Role='discrete';% will promote plots of the profiles with continuous lines |
---|
735 | end |
---|
736 | elseif isequal(ProjMode,'interp_lin') %filtering %linear interpolation: |
---|
737 | if ~check_abscissa |
---|
738 | XName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)}; |
---|
739 | ProjData.ListVarName=[ProjData.ListVarName {XName}]; |
---|
740 | ProjData.VarDimName=[ProjData.VarDimName {XName}]; |
---|
741 | nbvar=numel(ProjData.ListVarName); |
---|
742 | ProjData.VarAttribute{nbvar}.long_name='abscissa along line'; |
---|
743 | check_abscissa=1; % define abcissa only once |
---|
744 | end |
---|
745 | if ~isequal(errorflag,0) |
---|
746 | VarName_FF=FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag}; |
---|
747 | indsel=find(FieldData.(VarName_FF)==0); |
---|
748 | coord_x=coord_x(indsel); |
---|
749 | coord_y=coord_y(indsel); |
---|
750 | for ivar=1:numel(CellInfo{icell}.VarIndex) |
---|
751 | VarName=FieldData.ListVarName{CellInfo{icell}.VarIndex(ivar)}; |
---|
752 | FieldData.(VarName)=FieldData.(VarName)(indsel); |
---|
753 | end |
---|
754 | end |
---|
755 | [ProjVar,ListFieldProj,VarAttribute,errormsg]=calc_field_interp([coord_x coord_y],FieldData,CellInfo{icell}.FieldName,XI,YI); |
---|
756 | ProjData.X=Xproj; |
---|
757 | nbvar=numel(ProjData.ListVarName); |
---|
758 | ProjData.ListVarName=[ProjData.ListVarName ListFieldProj]; |
---|
759 | ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute]; |
---|
760 | for ivar=1:numel(VarAttribute) |
---|
761 | ProjData.VarDimName=[ProjData.VarDimName {XName}]; |
---|
762 | if isfield(VarAttribute{ivar},'Role') |
---|
763 | if strcmp(VarAttribute{ivar}.Role,'vector_x'); |
---|
764 | ivar_U=ivar+nbvar; |
---|
765 | elseif strcmp(VarAttribute{ivar}.Role,'vector_y'); |
---|
766 | ivar_V=ivar+nbvar; |
---|
767 | end |
---|
768 | end |
---|
769 | ProjData.VarAttribute{ivar+nbvar}.Role='continuous';% will promote plots of the profiles with continuous lines |
---|
770 | ProjData.(ListFieldProj{ivar})=ProjVar{ivar}; |
---|
771 | end |
---|
772 | end |
---|
773 | case 'tps'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
774 | if strcmp(ProjMode,'interp_tps') |
---|
775 | Coord=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex}); |
---|
776 | NbCentres=FieldData.(FieldData.ListVarName{CellInfo{icell}.NbCentres_tps}); |
---|
777 | SubRange=FieldData.(FieldData.ListVarName{CellInfo{icell}.SubRange_tps}); |
---|
778 | if isfield(CellInfo{icell},'VarIndex_vector_x_tps')&&isfield(CellInfo{icell},'VarIndex_vector_y_tps') |
---|
779 | FieldVar=cat(3,FieldData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x_tps}),FieldData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y_tps})); |
---|
780 | end |
---|
781 | [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord,NbCentres,SubRange,FieldVar,CellInfo{icell}.FieldName,cat(3,XI,YI)); |
---|
782 | ProjData.ListVarName=[ProjData.ListVarName {'X'}]; |
---|
783 | ProjData.VarDimName=[ProjData.VarDimName {'X'}]; |
---|
784 | ProjData.X=Xproj; |
---|
785 | nbvar=numel(ProjData.ListVarName); |
---|
786 | ProjData.VarAttribute{nbvar}.long_name='abscissa along line'; |
---|
787 | ProjVarName=(fieldnames(DataOut))'; |
---|
788 | ProjData.ListVarName=[ProjData.ListVarName ProjVarName]; |
---|
789 | ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute]; |
---|
790 | for ivar=1:numel(VarAttribute) |
---|
791 | ProjData.VarDimName=[ProjData.VarDimName {'X'}]; |
---|
792 | if isfield(VarAttribute{ivar},'Role') |
---|
793 | if strcmp(VarAttribute{ivar}.Role,'vector_x'); |
---|
794 | ivar_U=ivar+nbvar; |
---|
795 | elseif strcmp(VarAttribute{ivar}.Role,'vector_y'); |
---|
796 | ivar_V=ivar+nbvar; |
---|
797 | end |
---|
798 | end |
---|
799 | ProjData.VarAttribute{ivar+nbvar}.Role='continuous';% will promote plots of the profiles with continuous lines |
---|
800 | ProjData.(ProjVarName{ivar})=DataOut.(ProjVarName{ivar}); |
---|
801 | end |
---|
802 | end |
---|
803 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
804 | |
---|
805 | case 'grid' %case of structured coordinates |
---|
806 | if ~isequal(ObjectData.Type,'line')% exclude polyline |
---|
807 | errormsg=['no projection available on ' ObjectData.Type 'for structured coordinates']; |
---|
808 | return |
---|
809 | end% |
---|
810 | test_interp2=0;%default |
---|
811 | AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)}; |
---|
812 | AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)}; |
---|
813 | AX=FieldData.(AXName);% set of x positions |
---|
814 | AY=FieldData.(AYName);% set of y positions |
---|
815 | AName=FieldData.ListVarName{VarIndex(1)}; |
---|
816 | npxy=size(FieldData.(AName)); |
---|
817 | if max(NbDim)==3 % 3D case |
---|
818 | |
---|
819 | else |
---|
820 | npx=npxy(2); |
---|
821 | npy=npxy(1); |
---|
822 | if numel(AX)==2 |
---|
823 | DX=(AX(2)-AX(1))/(npx-1); |
---|
824 | else |
---|
825 | DX_vec=diff(AX); |
---|
826 | DX=max(DX_vec); |
---|
827 | DX_min=min(DX_vec); |
---|
828 | if (DX-DX_min)>0.0001*abs(DX) |
---|
829 | test_interp2=1; |
---|
830 | DX=DX_min; |
---|
831 | end |
---|
832 | end |
---|
833 | if numel(AY)==2 |
---|
834 | DY=(AY(2)-AY(1))/(npy-1); |
---|
835 | else |
---|
836 | DY_vec=diff(AY); |
---|
837 | DY=max(DY_vec); |
---|
838 | DY_min=min(DY_vec); |
---|
839 | if (DY-DY_min)>0.0001*abs(DY) |
---|
840 | test_interp2=1; |
---|
841 | DY=DY_min; |
---|
842 | end |
---|
843 | end |
---|
844 | AXI=linspace(AX(1),AX(end), npx);%set of x positions for the interpolated input data |
---|
845 | AYI=linspace(AY(1),AY(end), npy);%set of x positions for the interpolated input data |
---|
846 | if isfield(ObjectData,'DX') |
---|
847 | DXY_line=ObjectData.DX;%mesh on the projection line |
---|
848 | else |
---|
849 | DXY_line=sqrt(abs(DX*DY));% mesh on the projection line |
---|
850 | end |
---|
851 | dlinx=ObjectData.Coord(2,1)-ObjectData.Coord(1,1); |
---|
852 | dliny=ObjectData.Coord(2,2)-ObjectData.Coord(1,2); |
---|
853 | linelength=sqrt(dlinx*dlinx+dliny*dliny); |
---|
854 | theta=angle(dlinx+i*dliny);%angle of the line |
---|
855 | if isfield(FieldData,'RangeX') |
---|
856 | XMin=min(FieldData.RangeX);%shift of the origin on the line |
---|
857 | else |
---|
858 | XMin=0; |
---|
859 | end |
---|
860 | eval(['ProjData.' AXName '=linspace(XMin,XMin+linelength,linelength/DXY_line+1);'])%abscissa of the new pixels along the line |
---|
861 | y=linspace(-width,width,2*width/DXY_line+1);%ordintes of the new pixels (coordinate across the line) |
---|
862 | eval(['npX=length(ProjData.' AXName ');']) |
---|
863 | npY=length(y); %TODO: utiliser proj_grid |
---|
864 | eval(['[X,Y]=meshgrid(ProjData.' AXName ',y);'])%grid in the line coordinates |
---|
865 | XIMA=ObjectData.Coord(1,1)+(X-XMin)*cos(theta)-Y*sin(theta); |
---|
866 | YIMA=ObjectData.Coord(1,2)+(X-XMin)*sin(theta)+Y*cos(theta); |
---|
867 | XIMA=(XIMA-AX(1))/DX+1;% index of the original image along x |
---|
868 | YIMA=(YIMA-AY(1))/DY+1;% index of the original image along y |
---|
869 | XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line' |
---|
870 | YIMA=reshape(round(YIMA),1,npX*npY); |
---|
871 | flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image |
---|
872 | ind_in=find(flagin); |
---|
873 | ind_out=find(~flagin); |
---|
874 | ICOMB=(XIMA-1)*npy+YIMA; |
---|
875 | ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A |
---|
876 | nbcolor=1; %color images |
---|
877 | if numel(npxy)==2 |
---|
878 | nbcolor=1; |
---|
879 | elseif length(npxy)==3 |
---|
880 | nbcolor=npxy(3); |
---|
881 | else |
---|
882 | errormsg='multicomponent field not projected'; |
---|
883 | display(errormsg) |
---|
884 | return |
---|
885 | end |
---|
886 | nbvar=length(ProjData.ListVarName);% number of var from previous cells |
---|
887 | ProjData.ListVarName=[ProjData.ListVarName {AXName}]; |
---|
888 | ProjData.VarDimName=[ProjData.VarDimName {AXName}]; |
---|
889 | for ivar=VarIndex |
---|
890 | %VarName{ivar}=FieldData.ListVarName{ivar}; |
---|
891 | if test_interp2% interpolate on new grid |
---|
892 | FieldData.(FieldData.ListVarName{ivar})=interp2(FieldData.(AXName),FieldData.(AYName),FieldData.(FieldData.ListVarName{ivar}),AXI,AYI);%TO TEST |
---|
893 | end |
---|
894 | vec_A=reshape(squeeze(FieldData.(FieldData.ListVarName{ivar})),npx*npy,nbcolor); %put the original image in colum |
---|
895 | if nbcolor==1 |
---|
896 | vec_B(ind_in)=vec_A(ICOMB); |
---|
897 | vec_B(ind_out)=zeros(size(ind_out)); |
---|
898 | A_out=reshape(vec_B,npY,npX); |
---|
899 | ProjData.(FieldData.ListVarName{ivar}) =sum(A_out,1)/npY; |
---|
900 | elseif nbcolor==3 |
---|
901 | vec_B(ind_in,1:3)=vec_A(ICOMB,:); |
---|
902 | vec_B(ind_out,1)=zeros(size(ind_out)); |
---|
903 | vec_B(ind_out,2)=zeros(size(ind_out)); |
---|
904 | vec_B(ind_out,3)=zeros(size(ind_out)); |
---|
905 | A_out=reshape(vec_B,npY,npX,nbcolor); |
---|
906 | ProjData.(FieldData.ListVarName{ivar})=squeeze(sum(A_out,1)/npY); |
---|
907 | end |
---|
908 | ProjData.ListVarName=[ProjData.ListVarName FieldData.ListVarName{ivar}]; |
---|
909 | ProjData.VarDimName=[ProjData.VarDimName {AXName}];%to generalize with the initial name of the x coordinate |
---|
910 | ProjData.VarAttribute{ivar}.Role='continuous';% for plot with continuous line |
---|
911 | end |
---|
912 | if nbcolor==3 |
---|
913 | ProjData.VarDimName{end}={AXName,'rgb'}; |
---|
914 | end |
---|
915 | end |
---|
916 | |
---|
917 | end |
---|
918 | if ~isempty(ivar_U) && ~isempty(ivar_V) |
---|
919 | vector_x =ProjData.(ProjData.ListVarName{ivar_U}); |
---|
920 | ProjData.(ProjData.ListVarName{ivar_U}) =cos(theta)*vector_x+sin(theta)*ProjData.(ProjData.ListVarName{ivar_V}); |
---|
921 | ProjData.(ProjData.ListVarName{ivar_V}) =-sin(theta)*vector_x+cos(theta)*ProjData.(ProjData.ListVarName{ivar_V}); |
---|
922 | end |
---|
923 | end |
---|
924 | |
---|
925 | % %shotarter case for horizontal or vertical line (A FAIRE |
---|
926 | % % Rangx=[0.5 npx-0.5];%image coordiantes of corners |
---|
927 | % % Rangy=[npy-0.5 0.5]; |
---|
928 | % % if isfield(Calib,'Pxcmx')&isfield(Calib,'Pxcmy')%old calib |
---|
929 | % % Rangx=Rangx/Calib.Pxcmx; |
---|
930 | % % Rangy=Rangy/Calib.Pxcmy; |
---|
931 | % % else |
---|
932 | % % [Rangx]=phys_XYZ(Calib,Rangx,[0.5 0.5],[0 0]);%case of translations without rotation and quadratic deformation |
---|
933 | % % [xx,Rangy]=phys_XYZ(Calib,[0.5 0.5],Rangy,[0 0]); |
---|
934 | % % end |
---|
935 | % |
---|
936 | % % test_scal=0;%default% 3- 'UserData':(get(handles.Tag,'UserData') |
---|
937 | |
---|
938 | |
---|
939 | %----------------------------------------------------------------- |
---|
940 | %project on a plane |
---|
941 | % AJOUTER flux,circul,error |
---|
942 | function [ProjData,errormsg] = proj_plane(FieldData, ObjectData) |
---|
943 | %----------------------------------------------------------------- |
---|
944 | |
---|
945 | %% rotation angles |
---|
946 | PlaneAngle=[0 0]; |
---|
947 | norm_plane=[0 0 1]; |
---|
948 | %cos_om=1; |
---|
949 | %sin_om=0; |
---|
950 | test90x=0;%=1 for 90 degree rotation alround x axis |
---|
951 | test90y=0;%=1 for 90 degree rotation alround y axis |
---|
952 | % if strcmp(ObjectData.Type,'plane_z') |
---|
953 | % Delta_x=ObjectData.Coord(2,1)-ObjectData.Coord(1,1); |
---|
954 | % Delta_y=ObjectData.Coord(2,2)-ObjectData.Coord(1,2); |
---|
955 | % Delta_mod=sqrt(Delta_x*Delta_x+Delta_y*Delta_y); |
---|
956 | % ObjectData.Angle=[0 0 0]; |
---|
957 | % ObjectData.Angle(1)=90*Delta_x/Delta_mod; |
---|
958 | % ObjectData.0(2)=90*Delta_y/Delta_mod; |
---|
959 | % end |
---|
960 | if isfield(ObjectData,'Angle')&& isequal(size(ObjectData.Angle),[1 2])&& ~isequal(ObjectData.Angle,[0 0]) |
---|
961 | test90y=0;%isequal(ObjectData.Angle,[0 90 0]); |
---|
962 | PlaneAngle=(pi/180)*ObjectData.Angle; |
---|
963 | % om=norm(PlaneAngle);%norm of rotation angle in radians |
---|
964 | % OmAxis=PlaneAngle/om; %unit vector marking the rotation axis |
---|
965 | % cos_om=cos(om); |
---|
966 | % sin_om=sin(om); |
---|
967 | % coeff=OmAxis(3)*(1-cos_om); |
---|
968 | % %components of the unity vector norm_plane normal to the projection plane |
---|
969 | % norm_plane(1)=OmAxis(1)*coeff+OmAxis(2)*sin_om; |
---|
970 | % norm_plane(2)=OmAxis(2)*coeff-OmAxis(1)*sin_om; |
---|
971 | % norm_plane(3)=OmAxis(3)*coeff+cos_om; |
---|
972 | |
---|
973 | M1=[cos(PlaneAngle(1)) sin(PlaneAngle(1)) 0;-sin(PlaneAngle(1)) cos(PlaneAngle(1)) 0;0 0 1]; |
---|
974 | M2=[1 0 0;0 cos(PlaneAngle(2)) sin(PlaneAngle(2));0 -sin(PlaneAngle(2)) cos(PlaneAngle(2))]; |
---|
975 | M=M2*M1;% first rotate in the x,y plane with angle PlaneAngle(1), then slant around the new x axis0 with angle PlaneAngle(2) |
---|
976 | norm_plane=M*[0 0 1]'; |
---|
977 | |
---|
978 | end |
---|
979 | testangle=~isequal(PlaneAngle,[0 0])||~isequal(ObjectData.Coord(1:2),[0 0 ]) ;% && ~test90y && ~test90x;%=1 for slanted plane |
---|
980 | |
---|
981 | %% mesh sizes DX and DY |
---|
982 | DX=[]; |
---|
983 | DY=[];%default |
---|
984 | if isfield(ObjectData,'DX') && ~isempty(ObjectData.DX) |
---|
985 | DX=abs(ObjectData.DX);%mesh of interpolation points |
---|
986 | elseif isfield(FieldData,'CoordMesh') |
---|
987 | DX=FieldData.CoordMesh; |
---|
988 | end |
---|
989 | if isfield(ObjectData,'DY') && ~isempty(ObjectData.DY) |
---|
990 | DY=abs(ObjectData.DY);%mesh of interpolation points |
---|
991 | elseif isfield(FieldData,'CoordMesh') |
---|
992 | DY=FieldData.CoordMesh; |
---|
993 | end |
---|
994 | if ~strcmp(ObjectData.ProjMode,'projection') && (isempty(DX)||isempty(DY)) |
---|
995 | errormsg='DX or DY not defined'; |
---|
996 | return |
---|
997 | end |
---|
998 | InterpMesh=min(DX,DY);%mesh used for interpolation in a slanted plane |
---|
999 | % if strcmp(ObjectData.Type,'plane_z') |
---|
1000 | % InterpMesh=10*InterpMesh;%TODO: temporary, to shorten computation |
---|
1001 | % end |
---|
1002 | |
---|
1003 | %% extrema along each axis |
---|
1004 | testXMin=0;% test if min of X coordinates defined on the projection object, =0 by default |
---|
1005 | testXMax=0;% test if max of X coordinates defined on the projection object, =0 by default |
---|
1006 | testYMin=0;% test if min of Y coordinates defined on the projection object, =0 by default |
---|
1007 | testYMax=0;% test if max of Y coordinates defined on the projection object, =0 by default |
---|
1008 | if isfield(ObjectData,'RangeX') % rangeX defined by the projection object |
---|
1009 | XMin=min(ObjectData.RangeX); |
---|
1010 | XMax=max(ObjectData.RangeX); |
---|
1011 | testXMin=XMax>XMin;%=1 if XMin defined (i.e. RangeY has two distinct elements) |
---|
1012 | testXMax=1;% max of X coordinates defined on the projection object |
---|
1013 | end |
---|
1014 | if isfield(ObjectData,'RangeY') % rangeY defined by the projection object |
---|
1015 | YMin=min(ObjectData.RangeY); |
---|
1016 | YMax=max(ObjectData.RangeY); |
---|
1017 | testYMin=YMax>YMin;%=1 if YMin defined (i.e. RangeY has tow distinct elements) |
---|
1018 | testYMax=1;% max of Y coordinates defined on the projection object |
---|
1019 | end |
---|
1020 | width=0;%default width of the projection band |
---|
1021 | if isfield(ObjectData,'RangeZ') |
---|
1022 | width=max(ObjectData.RangeZ); |
---|
1023 | end |
---|
1024 | |
---|
1025 | %% interpolation range |
---|
1026 | thresh2=[]; |
---|
1027 | if isfield(ObjectData,'RangeInterp') |
---|
1028 | thresh2=ObjectData.RangeInterp*ObjectData.RangeInterp;%square of interpolation range (do not interpolate beyond this range) |
---|
1029 | end |
---|
1030 | |
---|
1031 | %% initiate Matlab structure for physical field |
---|
1032 | [ProjData,errormsg]=proj_heading(FieldData,ObjectData); |
---|
1033 | if ~isempty(errormsg) |
---|
1034 | return |
---|
1035 | end |
---|
1036 | |
---|
1037 | %% reproduce initial plane position and angle |
---|
1038 | if isfield(FieldData,'PlaneCoord')&&length(FieldData.PlaneCoord)==3&& isfield(ProjData,'ProjObjectCoord') |
---|
1039 | if length(ProjData.ProjObjectCoord)==3% if the projection plane has a z coordinate |
---|
1040 | if isfield(ProjData,'.PlaneCoord') && ~isequal(ProjData.PlaneCoord(3),ProjData.ProjObjectCoord) %check the consistency with the z coordinate of the field plane (set by calibration) |
---|
1041 | errormsg='inconsistent z position for field and projection plane'; |
---|
1042 | return |
---|
1043 | end |
---|
1044 | else % the z coordinate is set only by the field plane (by calibration) |
---|
1045 | ProjData.ProjObjectCoord(3)=FieldData.PlaneCoord(3); |
---|
1046 | end |
---|
1047 | if isfield(FieldData,'PlaneAngle') |
---|
1048 | if isfield(ProjData,'ProjObjectAngle') |
---|
1049 | if ~isequal(FieldData.PlaneAngle,ProjData.ProjObjectAngle) %check the consistency with the z coordinate of the field plane (set by calibration) |
---|
1050 | errormsg='inconsistent plane angle for field and projection plane'; |
---|
1051 | return |
---|
1052 | end |
---|
1053 | else |
---|
1054 | ProjData.ProjObjectAngle=FieldData.PlaneAngle; |
---|
1055 | end |
---|
1056 | end |
---|
1057 | end |
---|
1058 | ProjData.NbDim=2; |
---|
1059 | ProjData.ListVarName={}; |
---|
1060 | ProjData.VarDimName={}; |
---|
1061 | ProjData.VarAttribute={}; |
---|
1062 | if ~isempty(DX) && ~isempty(DY) |
---|
1063 | ProjData.CoordMesh=sqrt(DX*DY);%define typical data mesh, useful for mouse selection in plots |
---|
1064 | elseif isfield(FieldData,'CoordMesh') |
---|
1065 | ProjData.CoordMesh=FieldData.CoordMesh; |
---|
1066 | end |
---|
1067 | %error=0;%default |
---|
1068 | %flux=0; |
---|
1069 | %testfalse=0; |
---|
1070 | %ListIndex={}; |
---|
1071 | |
---|
1072 | %% group the variables (fields of 'FieldData') in cells of variables with the same dimensions |
---|
1073 | [CellInfo,NbDimArray,errormsg]=find_field_cells(FieldData); |
---|
1074 | |
---|
1075 | if ~isempty(errormsg) |
---|
1076 | errormsg=['error in proj_field/proj_plane:' errormsg]; |
---|
1077 | return |
---|
1078 | end |
---|
1079 | check_grid=zeros(size(CellInfo));% =1 if a grid is needed , =0 otherwise, for each field cell |
---|
1080 | |
---|
1081 | ProjMode=cell(size(CellInfo)); |
---|
1082 | for icell=1:numel(CellInfo) |
---|
1083 | ProjMode{icell}=ObjectData.ProjMode;% projection mode of the plane object |
---|
1084 | end |
---|
1085 | icell_grid=[];% field cell index which defines the grid |
---|
1086 | if ~strcmp(ObjectData.ProjMode,'projection')&& ~strcmp(ObjectData.Type,'plane_z')% TODO:rationalize |
---|
1087 | %% define the new coordinates in case of interpolation on a imposed grid |
---|
1088 | if ~testYMin |
---|
1089 | errormsg='min Y value not defined for the projection grid';return |
---|
1090 | end |
---|
1091 | if ~testYMax |
---|
1092 | errormsg='max Y value not defined for the projection grid';return |
---|
1093 | end |
---|
1094 | if ~testXMin |
---|
1095 | errormsg='min X value not defined for the projection grid';return |
---|
1096 | end |
---|
1097 | if ~testXMax |
---|
1098 | errormsg='max X value not defined for the projection grid';return |
---|
1099 | end |
---|
1100 | else |
---|
1101 | %% case of a grid requested by the input field |
---|
1102 | for icell=1:numel(CellInfo)% TODO: recalculate coordinates here to get the bounds in the rotated coordinates |
---|
1103 | if isfield(CellInfo{icell},'ProjModeRequest') |
---|
1104 | switch CellInfo{icell}.ProjModeRequest |
---|
1105 | case 'interp_lin' |
---|
1106 | ProjMode{icell}='interp_lin'; |
---|
1107 | case 'interp_tps' |
---|
1108 | ProjMode{icell}='interp_tps'; |
---|
1109 | end |
---|
1110 | end |
---|
1111 | if strcmp(ProjMode{icell},'interp_lin')||strcmp(ProjMode{icell},'interp_tps') |
---|
1112 | check_grid(icell)=1; |
---|
1113 | end |
---|
1114 | if strcmp(CellInfo{icell}.CoordType,'grid')&&NbDimArray(icell)>=2 |
---|
1115 | if ~testangle && isempty(icell_grid)% if the input gridded data is not modified, choose the first one in case of multiple gridded field cells |
---|
1116 | icell_grid=icell; |
---|
1117 | ProjMode{icell}='projection'; |
---|
1118 | end |
---|
1119 | check_grid(icell)=1; |
---|
1120 | end |
---|
1121 | end |
---|
1122 | if ~isempty(find(check_grid))% if a grid is requested by the input field |
---|
1123 | if isempty(icell_grid)% if the grid is not given by cell #icell_grid |
---|
1124 | if ~isfield(FieldData,'XMax') |
---|
1125 | FieldData=find_field_bounds(FieldData); |
---|
1126 | end |
---|
1127 | end |
---|
1128 | end |
---|
1129 | end |
---|
1130 | if ~isempty(find(check_grid))||~strcmp(ObjectData.ProjMode,'projection')%no existing gridded data used |
---|
1131 | if isempty(icell_grid)||~strcmp(ObjectData.ProjMode,'projection')%no existing gridded data used |
---|
1132 | AYName='coord_y'; |
---|
1133 | AXName='coord_x'; |
---|
1134 | if strcmp(ObjectData.ProjMode,'projection')||strcmp(ObjectData.Type,'plane_z') |
---|
1135 | ProjData.coord_y=[FieldData.YMin FieldData.YMax];%note that if projection is done on a grid, the Min and Max along each direction must have been defined |
---|
1136 | ProjData.coord_x=[FieldData.XMin FieldData.XMax]; |
---|
1137 | coord_x_proj=FieldData.XMin:FieldData.CoordMesh:FieldData.XMax; |
---|
1138 | coord_y_proj=FieldData.YMin:FieldData.CoordMesh:FieldData.YMax; |
---|
1139 | else |
---|
1140 | ProjData.coord_y=[ObjectData.RangeY(1) ObjectData.RangeY(2)];%note that if projection is done on a grid, the Min and Max along each direction must have been defined |
---|
1141 | ProjData.coord_x=[ObjectData.RangeX(1) ObjectData.RangeX(2)]; |
---|
1142 | coord_x_proj=ObjectData.RangeX(1):ObjectData.DX:ObjectData.RangeX(2); |
---|
1143 | coord_y_proj=ObjectData.RangeY(1):ObjectData.DY:ObjectData.RangeY(2); |
---|
1144 | end |
---|
1145 | [XI,YI]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates |
---|
1146 | ProjData.VarDimName={AYName,AXName}; |
---|
1147 | % XI=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(3))-YI*sin(PlaneAngle(3));%corresponding coordinates in the original system |
---|
1148 | % YI=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(3))+YI*cos(PlaneAngle(3)); |
---|
1149 | else% we use the existing grid from field cell #icell_grid |
---|
1150 | NbDim=NbDimArray(icell_grid); |
---|
1151 | AYName=FieldData.ListVarName{CellInfo{icell_grid}.CoordIndex(NbDim-1)};%name of input x coordinate (name preserved on projection) |
---|
1152 | AXName=FieldData.ListVarName{CellInfo{icell_grid}.CoordIndex(NbDim)};%name of input y coordinate (name preserved on projection) |
---|
1153 | AYDimName=FieldData.VarDimName{CellInfo{icell_grid}.CoordIndex(NbDim-1)};% |
---|
1154 | AXDimName=FieldData.VarDimName{CellInfo{icell_grid}.CoordIndex(NbDim)};% |
---|
1155 | ProjData.VarDimName={AYDimName,AXDimName}; |
---|
1156 | ProjData.(AYName)=FieldData.(AYName); % new (projected ) y coordinates |
---|
1157 | ProjData.(AXName)=FieldData.(AXName); % new (projected ) y coordinates |
---|
1158 | end |
---|
1159 | ProjData.ListVarName={AYName,AXName}; |
---|
1160 | |
---|
1161 | ProjData.VarAttribute={[],[]}; |
---|
1162 | end |
---|
1163 | |
---|
1164 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
1165 | % LOOP ON FIELD CELLS, PROJECT VARIABLES |
---|
1166 | % CellVarIndex=cells of variable index arrays |
---|
1167 | %ivar_new=0; % index of the current variable in the projected field |
---|
1168 | % icoord=0; |
---|
1169 | nbcoord=0;%number of added coordinate variables brought by projection |
---|
1170 | %nbvar=0; |
---|
1171 | vector_x_proj=[]; |
---|
1172 | vector_y_proj=[]; |
---|
1173 | for icell=1:length(CellInfo) |
---|
1174 | NbDim=NbDimArray(icell); |
---|
1175 | if NbDim<2 |
---|
1176 | continue % only cells represnting 2D or 3D fields are involved |
---|
1177 | end |
---|
1178 | VarIndex=CellInfo{icell}.VarIndex;% indices of the selected variables in the list FieldData.ListVarName |
---|
1179 | %dimensions |
---|
1180 | DimCell=FieldData.VarDimName{VarIndex(1)}; |
---|
1181 | if ischar(DimCell) |
---|
1182 | DimCell={DimCell};%name of dimensions |
---|
1183 | end |
---|
1184 | coord_z=0;%default |
---|
1185 | ListVarName={};% initiate list of projected variables for cell # icell |
---|
1186 | VarDimName={};% initiate coresponding list of dimensions for cell # icell |
---|
1187 | VarAttribute={};% initiate coresponding list of var attributes for cell # icell |
---|
1188 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
1189 | switch CellInfo{icell}.CoordType |
---|
1190 | |
---|
1191 | case 'scattered' |
---|
1192 | %% case of input fields with unstructured coordinates (applies for projMode ='projection' or 'interp_lin') |
---|
1193 | if strcmp(ProjMode{icell},'interp_tps') |
---|
1194 | continue %skip for next cell (needs tps field cell) |
---|
1195 | end |
---|
1196 | coord_x=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)});% initial x coordinates |
---|
1197 | coord_y=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)});% initial y coordinates |
---|
1198 | check3D=(numel(CellInfo{icell}.CoordIndex)==3); |
---|
1199 | if check3D |
---|
1200 | coord_z=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(1)}); |
---|
1201 | end |
---|
1202 | |
---|
1203 | % translate initial coordinates to account for the new origin |
---|
1204 | coord_x=coord_x-ObjectData.Coord(1,1); |
---|
1205 | coord_y=coord_y-ObjectData.Coord(1,2); |
---|
1206 | if check3D |
---|
1207 | coord_z=coord_z-ObjectData.Coord(1,3); |
---|
1208 | end |
---|
1209 | |
---|
1210 | % selection of the vectors in the projection range (3D case) |
---|
1211 | if check3D && width > 0 |
---|
1212 | %components of the unitiy vector normal to the projection plane |
---|
1213 | fieldZ=norm_plane(1)*coord_x + norm_plane(2)*coord_y+ norm_plane(3)*coord_z;% distance to the plane |
---|
1214 | indcut=find(abs(fieldZ) <= width); |
---|
1215 | for ivar=VarIndex |
---|
1216 | VarName=FieldData.ListVarName{ivar}; |
---|
1217 | FieldData.(VarName)=FieldData.(VarName)(indcut); |
---|
1218 | end |
---|
1219 | coord_x=coord_x(indcut); |
---|
1220 | coord_y=coord_y(indcut); |
---|
1221 | coord_z=coord_z(indcut); |
---|
1222 | end |
---|
1223 | |
---|
1224 | %rotate coordinates if needed: coord_X,coord_Y= = coordinates in the new plane |
---|
1225 | Psi=PlaneAngle(1); |
---|
1226 | Theta=PlaneAngle(2); |
---|
1227 | % Phi=PlaneAngle(3); |
---|
1228 | if testangle && ~test90y && ~test90x;%=1 for slanted plane |
---|
1229 | coord_X=(coord_x *cos(Phi) + coord_y* sin(Phi)); |
---|
1230 | coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta); |
---|
1231 | coord_Y=coord_Y+coord_z *sin(Theta); |
---|
1232 | coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER |
---|
1233 | coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi)); |
---|
1234 | else |
---|
1235 | coord_X=coord_x; |
---|
1236 | coord_Y=coord_y; |
---|
1237 | end |
---|
1238 | |
---|
1239 | %restriction to the range of X and Y if imposed by the projection object |
---|
1240 | testin=ones(size(coord_X)); %default |
---|
1241 | testbound=0; |
---|
1242 | if testXMin |
---|
1243 | testin=testin & (coord_X >= XMin); |
---|
1244 | testbound=1; |
---|
1245 | end |
---|
1246 | if testXMax |
---|
1247 | testin=testin & (coord_X <= XMax); |
---|
1248 | testbound=1; |
---|
1249 | end |
---|
1250 | if testYMin |
---|
1251 | testin=testin & (coord_Y >= YMin); |
---|
1252 | testbound=1; |
---|
1253 | end |
---|
1254 | if testYMin |
---|
1255 | testin=testin & (coord_Y <= YMax); |
---|
1256 | testbound=1; |
---|
1257 | end |
---|
1258 | if testbound |
---|
1259 | indcut=find(testin); |
---|
1260 | if isempty(indcut) |
---|
1261 | errormsg='data outside the bounds of the projection object'; |
---|
1262 | return |
---|
1263 | end |
---|
1264 | for ivar=VarIndex |
---|
1265 | VarName=FieldData.ListVarName{ivar}; |
---|
1266 | FieldData.(VarName)=FieldData.(VarName)(indcut); |
---|
1267 | end |
---|
1268 | coord_X=coord_X(indcut); |
---|
1269 | coord_Y=coord_Y(indcut); |
---|
1270 | if check3D |
---|
1271 | coord_Z=coord_Z(indcut); |
---|
1272 | end |
---|
1273 | end |
---|
1274 | |
---|
1275 | % two cases of projection for scattered coordinates |
---|
1276 | switch ProjMode{icell} |
---|
1277 | case 'projection' |
---|
1278 | nbvar=0; |
---|
1279 | %nbvar=numel(ProjData.ListVarName); |
---|
1280 | for ivar=VarIndex %transfer variables to the projection plane |
---|
1281 | VarName=FieldData.ListVarName{ivar}; |
---|
1282 | if ivar==CellInfo{icell}.CoordIndex(end) |
---|
1283 | ProjData.(VarName)=coord_X; |
---|
1284 | elseif ivar==CellInfo{icell}.CoordIndex(end-1) % y coordinate |
---|
1285 | ProjData.(VarName)=coord_Y; |
---|
1286 | elseif ~(check3D && ivar==CellInfo{icell}.CoordIndex(1)) % other variables (except Z coordinate wyhich is not reproduced) |
---|
1287 | ProjData.(VarName)=FieldData.(VarName); |
---|
1288 | end |
---|
1289 | if ~(check3D && ivar==CellInfo{icell}.CoordIndex(1)) |
---|
1290 | ListVarName=[ListVarName VarName]; |
---|
1291 | VarDimName=[VarDimName DimCell]; |
---|
1292 | nbvar=nbvar+1; |
---|
1293 | if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar |
---|
1294 | VarAttribute{nbvar}=FieldData.VarAttribute{ivar}; |
---|
1295 | end |
---|
1296 | end |
---|
1297 | end |
---|
1298 | case 'interp_lin'%interpolate data on a regular grid |
---|
1299 | if isfield(CellInfo{icell},'VarIndex_errorflag') |
---|
1300 | VarName_FF=FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag}; |
---|
1301 | indsel=find(FieldData.(VarName_FF)==0); |
---|
1302 | coord_X=coord_X(indsel); |
---|
1303 | coord_Y=coord_Y(indsel); |
---|
1304 | for ivar=1:numel(CellInfo{icell}.VarIndex) |
---|
1305 | VarName=FieldData.ListVarName{CellInfo{icell}.VarIndex(ivar)}; |
---|
1306 | FieldData.(VarName)=FieldData.(VarName)(indsel); |
---|
1307 | end |
---|
1308 | end |
---|
1309 | % interpolate and calculate field on the grid |
---|
1310 | |
---|
1311 | [VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp([coord_X coord_Y],FieldData,CellInfo{icell}.FieldName,XI,YI); |
---|
1312 | |
---|
1313 | % set to NaN interpolation points which are too far from any initial data (more than 2 CoordMesh) |
---|
1314 | if exist('scatteredInterpolant','file')%recent Matlab versions |
---|
1315 | F=scatteredInterpolant(coord_X, coord_Y,coord_X,'nearest'); |
---|
1316 | G=scatteredInterpolant(coord_X, coord_Y,coord_Y,'nearest'); |
---|
1317 | else |
---|
1318 | F=TriScatteredInterp([coord_X coord_Y],coord_X,'nearest'); |
---|
1319 | G=TriScatteredInterp([coord_X coord_Y],coord_Y,'nearest'); |
---|
1320 | end |
---|
1321 | Distx=F(XI,YI)-XI;% diff of x coordinates with the nearest measurement point |
---|
1322 | Disty=G(XI,YI)-YI;% diff of y coordinates with the nearest measurement point |
---|
1323 | Dist=Distx.*Distx+Disty.*Disty; |
---|
1324 | if ~isempty(thresh2) |
---|
1325 | for ivar=1:numel(VarVal) |
---|
1326 | VarVal{ivar}(Dist>thresh2)=NaN;% % put to NaN interpolated positions further than 4 meshes from initial data |
---|
1327 | end |
---|
1328 | end |
---|
1329 | if isfield(CellInfo{icell},'CheckSub') && CellInfo{icell}.CheckSub && ~isempty(vector_x_proj) |
---|
1330 | ProjData.(FieldData.ListVarName{vector_x_proj})=ProjData.(FieldData.ListVarName{vector_x_proj})-VarVal{1}; |
---|
1331 | ProjData.(FieldData.ListVarName{vector_y_proj})=ProjData.(FieldData.ListVarName{vector_y_proj})-VarVal{2}; |
---|
1332 | ListVarName={};% no new variable |
---|
1333 | VarAttribute={}; |
---|
1334 | else |
---|
1335 | VarDimName=cell(size(ListVarName)); |
---|
1336 | for ilist=1:numel(ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise |
---|
1337 | ListVarName{ilist}=regexprep(ListVarName{ilist},'(.+',''); |
---|
1338 | if ~isempty(find(strcmp(ListVarName{ilist},ProjData.ListVarName))) |
---|
1339 | ListVarName{ilist}=[ListVarName{ilist} '_1']; |
---|
1340 | end |
---|
1341 | ProjData.(ListVarName{ilist})=VarVal{ilist}; |
---|
1342 | VarDimName{ilist}={'coord_y','coord_x'}; |
---|
1343 | end |
---|
1344 | end |
---|
1345 | if isfield (CellInfo{icell},'VarIndex_vector_x')&& isfield (CellInfo{icell},'VarIndex_vector_y') |
---|
1346 | vector_x_proj=CellInfo{icell}.VarIndex_vector_x; %preserve for next cell |
---|
1347 | vector_y_proj=CellInfo{icell}.VarIndex_vector_y; %preserve for next cell |
---|
1348 | end |
---|
1349 | end |
---|
1350 | |
---|
1351 | case 'tps' |
---|
1352 | %% case of tps data (applies only in interp_tps mode) |
---|
1353 | if strcmp(ProjMode{icell},'interp_tps') |
---|
1354 | Coord=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex}); |
---|
1355 | |
---|
1356 | |
---|
1357 | NbCentres=FieldData.(FieldData.ListVarName{CellInfo{icell}.NbCentres_tps}); |
---|
1358 | SubRange=FieldData.(FieldData.ListVarName{CellInfo{icell}.SubRange_tps}); |
---|
1359 | checkUV=0; |
---|
1360 | if isfield(CellInfo{icell},'VarIndex_vector_x_tps')&&isfield(CellInfo{icell},'VarIndex_vector_y_tps') |
---|
1361 | FieldVar=cat(3,FieldData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x_tps}),FieldData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y_tps})); |
---|
1362 | checkUV=1; |
---|
1363 | end |
---|
1364 | |
---|
1365 | %rotate coordinates if needed: coord_X,coord_Y= = coordinates in the new plane |
---|
1366 | Psi=PlaneAngle(1); |
---|
1367 | Theta=PlaneAngle(2); |
---|
1368 | % Phi=PlaneAngle(3); |
---|
1369 | if testangle && ~test90y && ~test90x;%=1 for slanted plane |
---|
1370 | new_XI=XI*cos(Phi) - YI*sin(Phi)+ObjectData.Coord(1); |
---|
1371 | YI=XI *sin(Phi) + YI *cos(Phi)+ObjectData.Coord(2); |
---|
1372 | XI=new_XI; |
---|
1373 | % if checkUV |
---|
1374 | % UValue=cos(Phi)*FieldVar(:,:,1)+ sin(Phi)*FieldVar(:,:,2); |
---|
1375 | % FieldVar(:,:,2)=-sin(Phi)*FieldVar(:,:,1)+ cos(Phi)*FieldVar(:,:,2); |
---|
1376 | % FieldVar(:,:,1)=UValue; |
---|
1377 | % end |
---|
1378 | end |
---|
1379 | |
---|
1380 | % interpolate data using thin plate spline |
---|
1381 | [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord,NbCentres,SubRange,FieldVar,CellInfo{icell}.FieldName,cat(3,XI,YI)); |
---|
1382 | |
---|
1383 | % set to NaN interpolation points which are too far from any initial data (more than 2 CoordMesh) |
---|
1384 | Coord=permute(Coord,[1 3 2]); |
---|
1385 | Coord=reshape(Coord,size(Coord,1)*size(Coord,2),2); |
---|
1386 | if exist('scatteredInterpolant','file')%recent Matlab versions |
---|
1387 | F=scatteredInterpolant(Coord,Coord(:,1),'nearest'); |
---|
1388 | G=scatteredInterpolant(Coord,Coord(:,2),'nearest'); |
---|
1389 | else |
---|
1390 | F=TriScatteredInterp(Coord,Coord(:,1),'nearest'); |
---|
1391 | G=TriScatteredInterp(Coord,Coord(:,2),'nearest'); |
---|
1392 | end |
---|
1393 | Distx=F(XI,YI)-XI;% diff of x coordinates with the nearest measurement point |
---|
1394 | Disty=G(XI,YI)-YI;% diff of y coordinates with the nearest measurement point |
---|
1395 | Dist=Distx.*Distx+Disty.*Disty; |
---|
1396 | ListVarName=(fieldnames(DataOut))'; |
---|
1397 | VarDimName=cell(size(ListVarName)); |
---|
1398 | for ilist=1:numel(ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise |
---|
1399 | VarName=ListVarName{ilist}; |
---|
1400 | VarDimName{ilist}={'coord_y','coord_x'}; |
---|
1401 | ProjData.(VarName)=DataOut.(VarName); |
---|
1402 | if ~isempty(thresh2) |
---|
1403 | ProjData.(VarName)(Dist>thresh2)=NaN;% put to NaN interpolated positions further than RangeInterp from initial data |
---|
1404 | end |
---|
1405 | end |
---|
1406 | end |
---|
1407 | |
---|
1408 | case 'grid' |
---|
1409 | %% case of input fields defined on a structured grid |
---|
1410 | VarName=FieldData.ListVarName{VarIndex(1)};%get the first variable of the cell to get the input matrix dimensions |
---|
1411 | DimValue=size(FieldData.(VarName));%input matrix dimensions |
---|
1412 | DimValue(DimValue==1)=[];%remove singleton dimensions |
---|
1413 | NbDim=numel(DimValue);%update number of space dimensions |
---|
1414 | % nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate |
---|
1415 | if NbDim>=3 |
---|
1416 | if NbDim>3 |
---|
1417 | errormsg='matrices with more than 3 dimensions not handled'; |
---|
1418 | return |
---|
1419 | else |
---|
1420 | if numel(CellInfo{icell}.CoordIndex)==2% the third matrix dimension does not correspond to a space coordinate |
---|
1421 | nbcolor=DimValue(3); |
---|
1422 | DimValue(3)=[]; %number of 'color' components updated |
---|
1423 | NbDim=2;% space dimension set to 2 |
---|
1424 | end |
---|
1425 | end |
---|
1426 | end |
---|
1427 | Coord_z=[]; |
---|
1428 | Coord_y=[]; |
---|
1429 | Coord_x=[]; |
---|
1430 | |
---|
1431 | if testangle |
---|
1432 | ProjMode{icell}='interp_lin'; %request linear interpolation for projection on a tilted plane |
---|
1433 | end |
---|
1434 | |
---|
1435 | if isequal(ProjMode{icell},'projection')% && (~testangle || test90y || test90x) |
---|
1436 | if NbDim==2 && ~testXMin && ~testXMax && ~testYMin && ~testYMax% no range restriction |
---|
1437 | ListVarName=[ListVarName FieldData.ListVarName(VarIndex)]; |
---|
1438 | VarDimName=[VarDimName FieldData.VarDimName(VarIndex)]; |
---|
1439 | if isfield(FieldData,'VarAttribute') |
---|
1440 | VarAttribute=[VarAttribute FieldData.VarAttribute(VarIndex)]; |
---|
1441 | end |
---|
1442 | ProjData.(AYName)=FieldData.(AYName); |
---|
1443 | ProjData.(AXName)=FieldData.(AXName); |
---|
1444 | for ivar=VarIndex |
---|
1445 | VarName=FieldData.ListVarName{ivar}; |
---|
1446 | ProjData.(VarName)=FieldData.(VarName);% no change by projection |
---|
1447 | end |
---|
1448 | else |
---|
1449 | Coord{1}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(1)}); |
---|
1450 | Coord{2}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(2)}); |
---|
1451 | if NbDim==3 |
---|
1452 | Coord{3}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(3)}); |
---|
1453 | end |
---|
1454 | if numel(Coord{NbDim-1})==2% case of coordinate defined only by the first and last values |
---|
1455 | DY=(Coord{NbDim-1}(2)-Coord{NbDim-1}(1))/(DimValue(1)-1); |
---|
1456 | end |
---|
1457 | if numel(Coord{NbDim})==2% case of coordinate defined only by the first and last values |
---|
1458 | DX=(Coord{NbDim}(2)-Coord{NbDim}(1))/(DimValue(2)-1); |
---|
1459 | end |
---|
1460 | if testYMax |
---|
1461 | YIndexMax=(YMax-Coord{NbDim-1}(1))/DY+1;% matrix index corresponding to the max y value for the new field |
---|
1462 | if testYMin%test_direct(indY) |
---|
1463 | YIndexMin=(YMin-Coord{NbDim-1}(1))/DY+1;% matrix index corresponding to the min x value for the new field |
---|
1464 | else |
---|
1465 | YIndexMin=1; |
---|
1466 | end |
---|
1467 | else |
---|
1468 | YIndexMax=numel(Coord{NbDim-1}); |
---|
1469 | YIndexMin=1; |
---|
1470 | end |
---|
1471 | if testXMax |
---|
1472 | XIndexMax=(XMax-Coord{NbDim}(1))/DX+1;% matrix index corresponding to the max y value for the new field |
---|
1473 | if testYMin%test_direct(indY) |
---|
1474 | XIndexMin=(XMin-Coord{NbDim}(1))/DX+1;% matrix index corresponding to the min x value for the new field |
---|
1475 | else |
---|
1476 | XIndexMin=1; |
---|
1477 | end |
---|
1478 | else |
---|
1479 | XIndexMax=numel(Coord{NbDim}); |
---|
1480 | XIndexMin=1; |
---|
1481 | end |
---|
1482 | YIndexRange(1)=ceil(min(YIndexMin,YIndexMax));%first y index to select from the previous field |
---|
1483 | YIndexRange(1)=max(YIndexRange(1),1);% avoid bound lower than the first index |
---|
1484 | YIndexRange(2)=floor(max(YIndexMin,YIndexMax));%last y index to select from the previous field |
---|
1485 | YIndexRange(2)=min(YIndexRange(2),DimValue(NbDim-1));% limit to the last available index |
---|
1486 | XIndexRange(1)=ceil(min(XIndexMin,XIndexMax));%first x index to select from the previous field |
---|
1487 | XIndexRange(1)=max(XIndexRange(1),1);% avoid bound lower than the first index |
---|
1488 | XIndexRange(2)=floor(max(XIndexMin,XIndexMax));%last x index to select from the previous field |
---|
1489 | XIndexRange(2)=min(XIndexRange(2),DimValue(NbDim));% limit to the last available index |
---|
1490 | if test90y |
---|
1491 | ind_new=[3 2 1]; |
---|
1492 | DimCell={AYProjName,AXProjName}; |
---|
1493 | iz=ceil((ObjectData.Coord(1,1)-Coord{3}(1))/DX)+1; |
---|
1494 | for ivar=VarIndex |
---|
1495 | VarName=FieldData.ListVarName{ivar}; |
---|
1496 | ListVarName=[ListVarName VarName]; |
---|
1497 | VarDimName=[VarDimName {DimCell}]; |
---|
1498 | VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes |
---|
1499 | ProjData.(VarName)=permute(FieldData.(VarName),ind_new);% permute x and z indices for 90 degree rotation |
---|
1500 | ProjData.(VarName)=squeeze(ProjData.(VarName)(iz,:,:));% select the z index iz |
---|
1501 | end |
---|
1502 | ProjData.(AYName)=[Ybound(1) Ybound(2)]; %record the new (projected ) y coordinates |
---|
1503 | ProjData.(AXName)=[Coord{1}(end),Coord{1}(1)]; %record the new (projected ) x coordinates |
---|
1504 | else |
---|
1505 | if NbDim==3 |
---|
1506 | DZ=(Coord{1}(end)-Coord{1}(1))/(numel(Coord{1})-1); |
---|
1507 | DimCell(1)=[]; %suppress z variable |
---|
1508 | DimValue(1)=[]; |
---|
1509 | test_direct=1;%TOdo; GENERALIZE, SEE CASE OF points |
---|
1510 | if test_direct(1) |
---|
1511 | iz=ceil((ObjectData.Coord(1,3)-Coord{1}(1))/DZ)+1; |
---|
1512 | else |
---|
1513 | iz=ceil((Coord{1}(1)-ObjectData.Coord(1,3))/DZ)+1; |
---|
1514 | end |
---|
1515 | end |
---|
1516 | for ivar=VarIndex% loop on non coordinate variables |
---|
1517 | VarName=FieldData.ListVarName{ivar}; |
---|
1518 | ListVarName=[ListVarName VarName]; |
---|
1519 | VarDimName=[VarDimName {DimCell}]; |
---|
1520 | if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute)>=ivar |
---|
1521 | VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; |
---|
1522 | end |
---|
1523 | if NbDim==3 |
---|
1524 | ProjData.(VarName)=squeeze(FieldData.(VarName)(iz,YIndexRange(1):YIndexRange(end),XIndexRange(1):XIndexRange(end))); |
---|
1525 | else |
---|
1526 | ProjData.(VarName)=FieldData.(VarName)(YIndexRange(1):YIndexRange(end),XIndexRange(1):XIndexRange(end),:); |
---|
1527 | end |
---|
1528 | end |
---|
1529 | if testXMax |
---|
1530 | ProjData.(AXName)=Coord{NbDim}(1)+DX*(XIndexRange-1); %record the new (projected ) x coordinates |
---|
1531 | else |
---|
1532 | ProjData.(AXName)=FieldData.(AXName); |
---|
1533 | end |
---|
1534 | if testYMax |
---|
1535 | ProjData.(AYName)=Coord{NbDim-1}(1)+DY*(YIndexRange-1); %record the new (projected ) x coordinates |
---|
1536 | else |
---|
1537 | ProjData.(AYName)=FieldData.(AYName); |
---|
1538 | end |
---|
1539 | end |
---|
1540 | end |
---|
1541 | else % case with interpolation on a grid |
---|
1542 | if NbDim==2 %2D case |
---|
1543 | if isequal(ProjMode{icell},'interp_tps') |
---|
1544 | npx_interp_tps=ceil(abs(DX/DAX)); |
---|
1545 | npy_interp_tps=ceil(abs(DY/DAY)); |
---|
1546 | Minterp_tps=ones(npy_interp_tps,npx_interp_tps)/(npx_interp_tps*npy_interp_tps); |
---|
1547 | test_interp_tps=1; |
---|
1548 | else |
---|
1549 | test_interp_tps=0; |
---|
1550 | end |
---|
1551 | Coord{1}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(1)}); |
---|
1552 | Coord{2}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(2)}); |
---|
1553 | if ~(testXMin && testYMin)% % if the range of the projected coordinates is not fully defined by the projection object, find the extrema of the projected field |
---|
1554 | xcorner=[min(Coord{NbDim}) max(Coord{NbDim}) max(Coord{NbDim}) min(Coord{NbDim})]-ObjectData.Coord(1,1);% corner absissa of the original grid with respect to the new origin |
---|
1555 | ycorner=[min(Coord{NbDim-1}) min(Coord{NbDim-1}) max(Coord{NbDim-1}) max(Coord{NbDim-1})]-ObjectData.Coord(1,2);% corner ordinates of the original grid |
---|
1556 | xcor_new=xcorner*cos(PlaneAngle(3))+ycorner*sin(PlaneAngle(3));%coordinates of the corners in new frame |
---|
1557 | ycor_new=-xcorner*sin(PlaneAngle(3))+ycorner*cos(PlaneAngle(3)); |
---|
1558 | if ~testXMin |
---|
1559 | XMin=min(xcor_new); |
---|
1560 | end |
---|
1561 | if ~testXMax |
---|
1562 | XMax=max(xcor_new); |
---|
1563 | end |
---|
1564 | if ~testYMin |
---|
1565 | YMin=min(ycor_new); |
---|
1566 | end |
---|
1567 | if ~testYMax |
---|
1568 | YMax=max(ycor_new); |
---|
1569 | end |
---|
1570 | end |
---|
1571 | coord_x_proj=XMin:DX:XMax; |
---|
1572 | coord_y_proj=YMin:DY:YMax; |
---|
1573 | ProjData.(AYName)=[coord_y_proj(1) coord_y_proj(end)]; %record the new (projected ) y coordinates |
---|
1574 | ProjData.(AXName)=[coord_x_proj(1) coord_x_proj(end)]; %record the new (projected ) x coordinates |
---|
1575 | [X,YI]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates |
---|
1576 | XI=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(3))-YI*sin(PlaneAngle(3));%corresponding coordinates in the original system |
---|
1577 | YI=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(3))+YI*cos(PlaneAngle(3)); |
---|
1578 | if numel(Coord{1})==2% x coordiante defiend by its bounds, get the whole set |
---|
1579 | Coord{1}=linspace(Coord{1}(1),Coord{1}(2),CellInfo{icell}.CoordSize(1)); |
---|
1580 | end |
---|
1581 | if numel(Coord{2})==2% y coordiante defiend by its bounds, get the whole set |
---|
1582 | Coord{2}=linspace(Coord{2}(1),Coord{2}(2),CellInfo{icell}.CoordSize(2)); |
---|
1583 | end |
---|
1584 | [X,Y]=meshgrid(Coord{2},Coord{1});%initial coordinates |
---|
1585 | %name of error flag variable |
---|
1586 | FFName='FF';%default name (if not already used) |
---|
1587 | if isfield(ProjData,'FF') |
---|
1588 | ind=1; |
---|
1589 | while isfield(ProjData,['FF_' num2str(ind)]) |
---|
1590 | ind=ind+1; |
---|
1591 | end |
---|
1592 | FFName=['FF_' num2str(ind)];% append an index to the name of error flag, FF_1,FF_2... |
---|
1593 | end |
---|
1594 | % project all variables in the cell |
---|
1595 | for ivar=VarIndex |
---|
1596 | VarName=FieldData.ListVarName{ivar}; |
---|
1597 | if size(FieldData.(VarName),3)==1 |
---|
1598 | ProjData.(VarName)=interp2(X,Y,double(FieldData.(VarName)),XI,YI,'*linear');%interpolation fct |
---|
1599 | else |
---|
1600 | ProjData.(VarName)=interp2(X,Y,double(FieldData.(VarName)(:,:,1)),XI,YI,'*linear'); |
---|
1601 | for icolor=2:size(FieldData.(VarName),3)% project 'color' components |
---|
1602 | ProjData.(VarName)=cat(3,ProjData.(VarName),interp2(X,Y,double(FieldData.(VarName)(:,:,icolor)),XI,YI,'*linear')); |
---|
1603 | end |
---|
1604 | end |
---|
1605 | if isa(FieldData.(VarName),'uint8') |
---|
1606 | ProjData.(VarName)=uint8(ProjData.(VarName));%put result to integer 8 bits if the initial field is integer (image) |
---|
1607 | elseif isa(FieldData.(VarName),'uint16') |
---|
1608 | ProjData.(VarName)=uint16(ProjData.(VarName));%put result to integer 16 bits if the initial field is integer (image) |
---|
1609 | end |
---|
1610 | ListVarName=[ListVarName VarName]; |
---|
1611 | DimCell(1:2)={AYName,AXName}; |
---|
1612 | VarDimName=[VarDimName {DimCell}]; |
---|
1613 | if isfield(FieldData,'VarAttribute')&&length(FieldData.VarAttribute)>=ivar |
---|
1614 | VarAttribute{length(ListVarName)+nbcoord}=FieldData.VarAttribute{ivar}; |
---|
1615 | end; |
---|
1616 | ProjData.(FFName)=isnan(ProjData.(VarName));%detact NaN (points outside the interpolation range) |
---|
1617 | ProjData.(VarName)(ProjData.(FFName))=0; %set to 0 the NaN data |
---|
1618 | end |
---|
1619 | %update list of variables with error flag |
---|
1620 | ListVarName=[ListVarName FFName]; |
---|
1621 | VarDimName=[VarDimName {DimCell}]; |
---|
1622 | VarAttribute{numel(ListVarName)}.Role='errorflag'; |
---|
1623 | elseif ~testangle |
---|
1624 | % unstructured z coordinate |
---|
1625 | test_sup=(Coord{1}>=ObjectData.Coord(1,3)); |
---|
1626 | iz_sup=find(test_sup); |
---|
1627 | iz=iz_sup(1); |
---|
1628 | if iz>=1 & iz<=npz |
---|
1629 | %ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)]; |
---|
1630 | %ProjData.DimValue=[ProjData.DimValue npY npX]; |
---|
1631 | for ivar=VarIndex |
---|
1632 | VarName=FieldData.ListVarName{ivar}; |
---|
1633 | ListVarName=[ListVarName VarName]; |
---|
1634 | VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes |
---|
1635 | ProjData.(VarName)=squeeze(FieldData.(VarName)(iz,:,:));% select the z index iz |
---|
1636 | %TODO : do a vertical average for a thick plane |
---|
1637 | if test_interp(2) || test_interp(3) |
---|
1638 | ProjData.(VarName)=interp2(Coord{3},Coord{2},ProjData.(VarName),Coord_x,Coord_y); |
---|
1639 | end |
---|
1640 | end |
---|
1641 | end |
---|
1642 | else %projection of structured coordinates on oblique plane |
---|
1643 | % determine the boundaries of the projected field, |
---|
1644 | % first find the 8 summits of the initial volume in the |
---|
1645 | PlaneAngle=ObjectData.Angle*pi/180; |
---|
1646 | % new coordinates |
---|
1647 | Coord{1}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(1)});%initial z coordinates |
---|
1648 | Coord{2}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(2)});%initial y coordinates |
---|
1649 | Coord{3}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(3)});%initial x coordinates |
---|
1650 | summit=zeros(3,8);% initialize summit coordinates |
---|
1651 | summit(1,1:4)=[Coord{3}(1) Coord{3}(end) Coord{3}(1) Coord{3}(end)];%square |
---|
1652 | summit(2,1:4)=[Coord{2}(1) Coord{2}(1) Coord{2}(end) Coord{2}(end)];% square at z= Coord{1}(1) |
---|
1653 | summit(1:2,5:8)=summit(1:2,1:4); |
---|
1654 | summit(3,:)=[Coord{1}(1)*ones(1,4) Coord{1}(end)*ones(1,4)]; |
---|
1655 | %Mrot_inv=rodrigues(-PlaneAngle); |
---|
1656 | newsummit=zeros(3,8);% initialize the rotated summit coordinates |
---|
1657 | ObjectData.Coord=ObjectData.Coord';% set ObjectData.Coord as a vertical vector |
---|
1658 | if size(ObjectData.Coord,1)<3 |
---|
1659 | ObjectData.Coord=[ObjectData.Coord; 0];%add z origin at z=0 by default |
---|
1660 | end |
---|
1661 | |
---|
1662 | M1=[cos(PlaneAngle(1)) sin(PlaneAngle(1)) 0;-sin(PlaneAngle(1)) cos(PlaneAngle(1)) 0;0 0 1]; |
---|
1663 | M2=[1 0 0;0 cos(PlaneAngle(2)) sin(PlaneAngle(2));0 -sin(PlaneAngle(2)) cos(PlaneAngle(2))]; |
---|
1664 | M=M2*M1; |
---|
1665 | M_inv=inv(M); |
---|
1666 | |
---|
1667 | for isummit=1:8% TODO: introduce a function for rotation of n points (to use also for scattered data) |
---|
1668 | newsummit(:,isummit)=M*(summit(:,isummit)-(ObjectData.Coord)); |
---|
1669 | end |
---|
1670 | coord_x_proj=min(newsummit(1,:)):InterpMesh: max(newsummit(1,:));% set of coordinqtes in the projection plane |
---|
1671 | coord_y_proj=min(newsummit(2,:)):InterpMesh: max(newsummit(2,:)); |
---|
1672 | coord_z_proj=-width:width; |
---|
1673 | %Mrot=rodrigues(PlaneAngle);% inverse rotation matrix |
---|
1674 | Origin=M_inv*[coord_x_proj(1);coord_y_proj(1);coord_z_proj(1)]+ObjectData.Coord; |
---|
1675 | npx=numel(coord_x_proj); |
---|
1676 | npy=numel(coord_y_proj); |
---|
1677 | npz=numel(coord_z_proj); |
---|
1678 | |
---|
1679 | %modangle=sqrt(PlaneAngle(1)*PlaneAngle(1)+PlaneAngle(2)*PlaneAngle(2)); |
---|
1680 | % cosphi=PlaneAngle(1)/modangle; |
---|
1681 | % sinphi=PlaneAngle(2)/modangle; |
---|
1682 | iX=[coord_x_proj(end)-coord_x_proj(1);0;0]/(npx-1); |
---|
1683 | iY=[0;coord_y_proj(end)-coord_y_proj(1);0]/(npy-1); |
---|
1684 | iZ=[0;0;coord_z_proj(end)-coord_z_proj(1)]/(npz-1); |
---|
1685 | % iX(1:2)=[cosphi -sinphi;sinphi cosphi]*iX(1:2); |
---|
1686 | % iY(1:2)=[-cosphi -sinphi;sinphi cosphi]*iY(1:2); |
---|
1687 | |
---|
1688 | ix=M_inv*iX;% vector along the new x coordinates transformed into old coordinates |
---|
1689 | iy=M_inv*iY;% vector along y coordinates |
---|
1690 | iz=M_inv*iZ;% vector along z coordinates |
---|
1691 | |
---|
1692 | [Grid_x,Grid_y,Grid_z]=meshgrid(0:npx-1,0:npy-1,0:npz-1); |
---|
1693 | if ismatrix(Grid_x)% add a singleton in case of a single z value |
---|
1694 | Grid_x=shiftdim(Grid_x,-1); |
---|
1695 | Grid_y=shiftdim(Grid_y,-1); |
---|
1696 | Grid_z=shiftdim(Grid_z,-1); |
---|
1697 | end |
---|
1698 | XI=Origin(1)+ix(1)*Grid_x+iy(1)*Grid_y+iz(1)*Grid_z; |
---|
1699 | YI=Origin(2)+ix(2)*Grid_x+iy(2)*Grid_y+iz(2)*Grid_z; |
---|
1700 | ZI=Origin(3)+ix(3)*Grid_x+iy(3)*Grid_y+iz(3)*Grid_z; |
---|
1701 | [X,Y,Z]=meshgrid(Coord{3},Coord{2},Coord{1});% mesh in the initial coordinates |
---|
1702 | for ivar=VarIndex |
---|
1703 | VarName=FieldData.ListVarName{ivar}; |
---|
1704 | ListVarName=[ListVarName VarName]; |
---|
1705 | VarDimName=[VarDimName {{'coord_y','coord_x'}}]; |
---|
1706 | VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes |
---|
1707 | FieldData.(VarName)=permute(FieldData.(VarName),[2 3 1]); |
---|
1708 | ProjData.coord_x=coord_x_proj; |
---|
1709 | ProjData.coord_y=coord_y_proj; |
---|
1710 | ProjData.(VarName)=interp3(X,Y,Z,double(FieldData.(VarName)),XI,YI,ZI,'*linear'); |
---|
1711 | ProjData.(VarName)=nanmean(ProjData.(VarName),3); |
---|
1712 | ProjData.(VarName)=squeeze(ProjData.(VarName)); |
---|
1713 | end |
---|
1714 | end |
---|
1715 | end |
---|
1716 | end |
---|
1717 | % update the global list of projected variables: |
---|
1718 | ProjData.ListVarName=[ProjData.ListVarName ListVarName]; |
---|
1719 | ProjData.VarDimName=[ProjData.VarDimName VarDimName]; |
---|
1720 | ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute]; |
---|
1721 | |
---|
1722 | %% projection of velocity components in the rotated coordinates |
---|
1723 | if testangle |
---|
1724 | ivar_U=[];ivar_V=[];ivar_W=[]; |
---|
1725 | for ivar=1:numel(VarAttribute) |
---|
1726 | if isfield(VarAttribute{ivar},'Role') |
---|
1727 | if strcmp(VarAttribute{ivar}.Role,'vector_x') |
---|
1728 | ivar_U=ivar; |
---|
1729 | elseif strcmp(VarAttribute{ivar}.Role,'vector_y') |
---|
1730 | ivar_V=ivar; |
---|
1731 | elseif strcmp(VarAttribute{ivar}.Role,'vector_z') |
---|
1732 | ivar_W=ivar; |
---|
1733 | end |
---|
1734 | end |
---|
1735 | end |
---|
1736 | if ~isempty(ivar_U) |
---|
1737 | if isempty(ivar_V) |
---|
1738 | msgbox_uvmat('ERROR','v velocity component missing in proj_field.m') |
---|
1739 | return |
---|
1740 | else |
---|
1741 | UName=ListVarName{ivar_U}; |
---|
1742 | VName=ListVarName{ivar_V}; |
---|
1743 | UValue=cos(PlaneAngle(3))*ProjData.(UName)+ sin(PlaneAngle(3))*ProjData.(VName); |
---|
1744 | ProjData.(VName)=(-sin(PlaneAngle(3))*ProjData.(UName)+ cos(PlaneAngle(3))*ProjData.(VName)); |
---|
1745 | ProjData.(UName)=UValue; |
---|
1746 | if ~isempty(ivar_W) |
---|
1747 | WName=FieldData.ListVarName{ivar_W}; |
---|
1748 | VValue=ProjData.(VName)+ ProjData.(WName)*sin(Theta);% |
---|
1749 | ProjData.(WName)=NormVec_X*ProjData.(UName)+ NormVec_Y*ProjData.(VName)+ NormVec_Z* ProjData.(WName); |
---|
1750 | ProjData.(VName)=VValue; |
---|
1751 | end |
---|
1752 | end |
---|
1753 | end |
---|
1754 | end |
---|
1755 | end |
---|
1756 | % %prepare substraction in case of two input fields |
---|
1757 | % SubData.ListVarName={}; |
---|
1758 | % SubData.VarDimName={}; |
---|
1759 | % SubData.VarAttribute={}; |
---|
1760 | % check_remove=zeros(size(ProjData.ListVarName)); |
---|
1761 | % for iproj=1:numel(ProjData.VarAttribute) |
---|
1762 | % if isfield(ProjData.VarAttribute{iproj},'CheckSub')&&isequal(ProjData.VarAttribute{iproj}.CheckSub,1) |
---|
1763 | % VarName=ProjData.ListVarName{iproj}; |
---|
1764 | % SubData.ListVarName=[SubData.ListVarName {VarName}]; |
---|
1765 | % SubData.VarDimName=[SubData.VarDimName ProjData.VarDimName{iproj}]; |
---|
1766 | % SubData.VarAttribute=[SubData.VarAttribute ProjData.VarAttribute{iproj}]; |
---|
1767 | % SubData.(VarName)=ProjData.(VarName); |
---|
1768 | % check_remove(iproj)=1; |
---|
1769 | % end |
---|
1770 | % end |
---|
1771 | % if ~isempty(find(check_remove)) |
---|
1772 | % ind_remove=find(check_remove); |
---|
1773 | % ProjData.ListVarName(ind_remove)=[]; |
---|
1774 | % ProjData.VarDimName(ind_remove)=[]; |
---|
1775 | % ProjData.VarAttribute(ind_remove)=[]; |
---|
1776 | % ProjData=sub_field(ProjData,[],SubData); |
---|
1777 | % end |
---|
1778 | |
---|
1779 | %----------------------------------------------------------------- |
---|
1780 | %projection in a volume |
---|
1781 | function [ProjData,errormsg] = proj_volume(FieldData, ObjectData) |
---|
1782 | |
---|
1783 | %----------------------------------------------------------------- |
---|
1784 | ProjData=FieldData;%default output |
---|
1785 | |
---|
1786 | %% axis origin |
---|
1787 | if isempty(ObjectData.Coord) |
---|
1788 | ObjectData.Coord(1,1)=0;%origin of the plane set to [0 0] by default |
---|
1789 | ObjectData.Coord(1,2)=0; |
---|
1790 | ObjectData.Coord(1,3)=0; |
---|
1791 | end |
---|
1792 | |
---|
1793 | %% rotation angles |
---|
1794 | VolumeAngle=[0 0 0]; |
---|
1795 | norm_plane=[0 0 1]; |
---|
1796 | if isfield(ObjectData,'Angle')&& isequal(size(ObjectData.Angle),[1 3])&& ~isequal(ObjectData.Angle,[0 0 0]) |
---|
1797 | PlaneAngle=ObjectData.Angle; |
---|
1798 | VolumeAngle=ObjectData.Angle; |
---|
1799 | om=norm(VolumeAngle);%norm of rotation angle in radians |
---|
1800 | OmAxis=VolumeAngle/om; %unit vector marking the rotation axis |
---|
1801 | cos_om=cos(pi*om/180); |
---|
1802 | sin_om=sin(pi*om/180); |
---|
1803 | coeff=OmAxis(3)*(1-cos_om); |
---|
1804 | %components of the unity vector norm_plane normal to the projection plane |
---|
1805 | norm_plane(1)=OmAxis(1)*coeff+OmAxis(2)*sin_om; |
---|
1806 | norm_plane(2)=OmAxis(2)*coeff-OmAxis(1)*sin_om; |
---|
1807 | norm_plane(3)=OmAxis(3)*coeff+cos_om; |
---|
1808 | end |
---|
1809 | testangle=~isequal(VolumeAngle,[0 0 0]); |
---|
1810 | |
---|
1811 | %% mesh sizes DX, DY, DZ |
---|
1812 | DX=0; |
---|
1813 | DY=0; %default |
---|
1814 | DZ=0; |
---|
1815 | if isfield(ObjectData,'DX')&~isempty(ObjectData.DX) |
---|
1816 | DX=abs(ObjectData.DX);%mesh of interpolation points |
---|
1817 | end |
---|
1818 | if isfield(ObjectData,'DY')&~isempty(ObjectData.DY) |
---|
1819 | DY=abs(ObjectData.DY);%mesh of interpolation points |
---|
1820 | end |
---|
1821 | if isfield(ObjectData,'DZ')&~isempty(ObjectData.DZ) |
---|
1822 | DZ=abs(ObjectData.DZ);%mesh of interpolation points |
---|
1823 | end |
---|
1824 | if ~strcmp(ProjMode,'projection') && (DX==0||DY==0||DZ==0) |
---|
1825 | errormsg='grid mesh DX , DY or DZ is missing'; |
---|
1826 | return |
---|
1827 | end |
---|
1828 | |
---|
1829 | %% extrema along each axis |
---|
1830 | testXMin=0; |
---|
1831 | testXMax=0; |
---|
1832 | testYMin=0; |
---|
1833 | testYMax=0; |
---|
1834 | if isfield(ObjectData,'RangeX') |
---|
1835 | XMin=min(ObjectData.RangeX); |
---|
1836 | XMax=max(ObjectData.RangeX); |
---|
1837 | testXMin=XMax>XMin; |
---|
1838 | testXMax=1; |
---|
1839 | end |
---|
1840 | if isfield(ObjectData,'RangeY') |
---|
1841 | YMin=min(ObjectData.RangeY); |
---|
1842 | YMax=max(ObjectData.RangeY); |
---|
1843 | testYMin=YMax>YMin; |
---|
1844 | testYMax=1; |
---|
1845 | end |
---|
1846 | width=0;%default width of the projection band |
---|
1847 | if isfield(ObjectData,'RangeZ') |
---|
1848 | ZMin=min(ObjectData.RangeZ); |
---|
1849 | ZMax=max(ObjectData.RangeZ); |
---|
1850 | testZMin=ZMax>ZMin; |
---|
1851 | testZMax=1; |
---|
1852 | end |
---|
1853 | |
---|
1854 | %% initiate Matlab structure for physical field |
---|
1855 | [ProjData,errormsg]=proj_heading(FieldData,ObjectData); |
---|
1856 | if ~isempty(errormsg) |
---|
1857 | return |
---|
1858 | end |
---|
1859 | |
---|
1860 | ProjData.NbDim=3; |
---|
1861 | ProjData.ListVarName={}; |
---|
1862 | ProjData.VarDimName={}; |
---|
1863 | if ~isequal(DX,0)&& ~isequal(DY,0) |
---|
1864 | ProjData.CoordMesh=sqrt(DX*DY);%define typical data mesh, useful for mouse selection in plots |
---|
1865 | elseif isfield(FieldData,'CoordMesh') |
---|
1866 | ProjData.CoordMesh=FieldData.CoordMesh; |
---|
1867 | end |
---|
1868 | |
---|
1869 | error=0;%default |
---|
1870 | flux=0; |
---|
1871 | testfalse=0; |
---|
1872 | ListIndex={}; |
---|
1873 | |
---|
1874 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
1875 | %% group the variables (fields of 'FieldData') in cells of variables with the same dimensions |
---|
1876 | %----------------------------------------------------------------- |
---|
1877 | idimvar=0; |
---|
1878 | % LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS |
---|
1879 | % CellVarIndex=cells of variable index arrays |
---|
1880 | ivar_new=0; % index of the current variable in the projected field |
---|
1881 | icoord=0; |
---|
1882 | nbcoord=0;%number of added coordinate variables brought by projection |
---|
1883 | nbvar=0; |
---|
1884 | for icell=1:length(CellVarIndex) |
---|
1885 | NbDim=NbDimVec(icell); |
---|
1886 | if NbDim<3 |
---|
1887 | continue |
---|
1888 | end |
---|
1889 | VarIndex=CellVarIndex{icell};% indices of the selected variables in the list FieldData.ListVarName |
---|
1890 | VarType=VarTypeCell{icell}; |
---|
1891 | ivar_X=VarType.coord_x; |
---|
1892 | ivar_Y=VarType.coord_y; |
---|
1893 | ivar_Z=VarType.coord_z; |
---|
1894 | ivar_U=VarType.vector_x; |
---|
1895 | ivar_V=VarType.vector_y; |
---|
1896 | ivar_W=VarType.vector_z; |
---|
1897 | ivar_C=VarType.scalar ; |
---|
1898 | ivar_Anc=VarType.ancillary; |
---|
1899 | test_anc=zeros(size(VarIndex)); |
---|
1900 | test_anc(ivar_Anc)=ones(size(ivar_Anc)); |
---|
1901 | ivar_F=VarType.warnflag; |
---|
1902 | ivar_FF=VarType.errorflag; |
---|
1903 | check_unstructured_coord=~isempty(ivar_X) && ~isempty(ivar_Y); |
---|
1904 | DimCell=FieldData.VarDimName{VarIndex(1)}; |
---|
1905 | if ischar(DimCell) |
---|
1906 | DimCell={DimCell};%name of dimensions |
---|
1907 | end |
---|
1908 | |
---|
1909 | %% case of input fields with unstructured coordinates |
---|
1910 | if check_unstructured_coord |
---|
1911 | XName=FieldData.ListVarName{ivar_X}; |
---|
1912 | YName=FieldData.ListVarName{ivar_Y}; |
---|
1913 | eval(['coord_x=FieldData.' XName ';']) |
---|
1914 | eval(['coord_y=FieldData.' YName ';']) |
---|
1915 | if length(ivar_Z)==1 |
---|
1916 | ZName=FieldData.ListVarName{ivar_Z}; |
---|
1917 | eval(['coord_z=FieldData.' ZName ';']) |
---|
1918 | end |
---|
1919 | |
---|
1920 | % translate initial coordinates |
---|
1921 | coord_x=coord_x-ObjectData.Coord(1,1); |
---|
1922 | coord_y=coord_y-ObjectData.Coord(1,2); |
---|
1923 | if ~isempty(ivar_Z) |
---|
1924 | coord_z=coord_z-ObjectData.Coord(1,3); |
---|
1925 | end |
---|
1926 | |
---|
1927 | % selection of the vectors in the projection range |
---|
1928 | % if length(ivar_Z)==1 && width > 0 |
---|
1929 | % %components of the unitiy vector normal to the projection plane |
---|
1930 | % fieldZ=NormVec_X*coord_x + NormVec_Y*coord_y+ NormVec_Z*coord_z;% distance to the plane |
---|
1931 | % indcut=find(abs(fieldZ) <= width); |
---|
1932 | % for ivar=VarIndex |
---|
1933 | % VarName=FieldData.ListVarName{ivar}; |
---|
1934 | % eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);']) |
---|
1935 | % % A VOIR : CAS DE VAR STRUCTUREE MAIS PAS GRILLE REGULIERE : INTERPOLER SUR GRILLE REGULIERE |
---|
1936 | % end |
---|
1937 | % coord_x=coord_x(indcut); |
---|
1938 | % coord_y=coord_y(indcut); |
---|
1939 | % coord_z=coord_z(indcut); |
---|
1940 | % end |
---|
1941 | |
---|
1942 | %rotate coordinates if needed: TODO modify |
---|
1943 | if testangle |
---|
1944 | coord_X=(coord_x *cos(Phi) + coord_y* sin(Phi)); |
---|
1945 | coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta); |
---|
1946 | if ~isempty(ivar_Z) |
---|
1947 | coord_Y=coord_Y+coord_z *sin(Theta); |
---|
1948 | end |
---|
1949 | |
---|
1950 | coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER |
---|
1951 | coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi)); |
---|
1952 | |
---|
1953 | else |
---|
1954 | coord_X=coord_x; |
---|
1955 | coord_Y=coord_y; |
---|
1956 | coord_Z=coord_z; |
---|
1957 | end |
---|
1958 | %restriction to the range of x and y if imposed |
---|
1959 | testin=ones(size(coord_X)); %default |
---|
1960 | testbound=0; |
---|
1961 | if testXMin |
---|
1962 | testin=testin & (coord_X >= XMin); |
---|
1963 | testbound=1; |
---|
1964 | end |
---|
1965 | if testXMax |
---|
1966 | testin=testin & (coord_X <= XMax); |
---|
1967 | testbound=1; |
---|
1968 | end |
---|
1969 | if testYMin |
---|
1970 | testin=testin & (coord_Y >= YMin); |
---|
1971 | testbound=1; |
---|
1972 | end |
---|
1973 | if testYMax |
---|
1974 | testin=testin & (coord_Y <= YMax); |
---|
1975 | testbound=1; |
---|
1976 | end |
---|
1977 | if testbound |
---|
1978 | indcut=find(testin); |
---|
1979 | for ivar=VarIndex |
---|
1980 | VarName=FieldData.ListVarName{ivar}; |
---|
1981 | eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);']) |
---|
1982 | end |
---|
1983 | coord_X=coord_X(indcut); |
---|
1984 | coord_Y=coord_Y(indcut); |
---|
1985 | if length(ivar_Z)==1 |
---|
1986 | coord_Z=coord_Z(indcut); |
---|
1987 | end |
---|
1988 | end |
---|
1989 | % different cases of projection |
---|
1990 | if isequal(ObjectData.ProjMode,'projection')%%%%%%% NOT USED %%%%%%%%%% |
---|
1991 | for ivar=VarIndex %transfer variables to the projection plane |
---|
1992 | VarName=FieldData.ListVarName{ivar}; |
---|
1993 | if ivar==ivar_X %x coordinate |
---|
1994 | eval(['ProjData.' VarName '=coord_X;']) |
---|
1995 | elseif ivar==ivar_Y % y coordinate |
---|
1996 | eval(['ProjData.' VarName '=coord_Y;']) |
---|
1997 | elseif isempty(ivar_Z) || ivar~=ivar_Z % other variables (except Z coordinate wyhich is not reproduced) |
---|
1998 | eval(['ProjData.' VarName '=FieldData.' VarName ';']) |
---|
1999 | end |
---|
2000 | if isempty(ivar_Z) || ivar~=ivar_Z |
---|
2001 | ProjData.ListVarName=[ProjData.ListVarName VarName]; |
---|
2002 | ProjData.VarDimName=[ProjData.VarDimName DimCell]; |
---|
2003 | nbvar=nbvar+1; |
---|
2004 | if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar |
---|
2005 | ProjData.VarAttribute{nbvar}=FieldData.VarAttribute{ivar}; |
---|
2006 | end |
---|
2007 | end |
---|
2008 | end |
---|
2009 | elseif isequal(ObjectData.ProjMode,'interp_lin')||isequal(ObjectData.ProjMode,'interp_tps')%interpolate data on a regular grid |
---|
2010 | coord_x_proj=XMin:DX:XMax; |
---|
2011 | coord_y_proj=YMin:DY:YMax; |
---|
2012 | coord_z_proj=ZMin:DZ:ZMax; |
---|
2013 | DimCell={'coord_z','coord_y','coord_x'}; |
---|
2014 | ProjData.ListVarName={'coord_z','coord_y','coord_x'}; |
---|
2015 | ProjData.VarDimName={'coord_z','coord_y','coord_x'}; |
---|
2016 | nbcoord=2; |
---|
2017 | ProjData.coord_z=[ZMin ZMax]; |
---|
2018 | ProjData.coord_y=[YMin YMax]; |
---|
2019 | ProjData.coord_x=[XMin XMax]; |
---|
2020 | if isempty(ivar_X), ivar_X=0; end; |
---|
2021 | if isempty(ivar_Y), ivar_Y=0; end; |
---|
2022 | if isempty(ivar_Z), ivar_Z=0; end; |
---|
2023 | if isempty(ivar_U), ivar_U=0; end; |
---|
2024 | if isempty(ivar_V), ivar_V=0; end; |
---|
2025 | if isempty(ivar_W), ivar_W=0; end; |
---|
2026 | if isempty(ivar_F), ivar_F=0; end; |
---|
2027 | if isempty(ivar_FF), ivar_FF=0; end; |
---|
2028 | if ~isequal(ivar_FF,0) |
---|
2029 | VarName_FF=FieldData.ListVarName{ivar_FF}; |
---|
2030 | eval(['indsel=find(FieldData.' VarName_FF '==0);']) |
---|
2031 | coord_X=coord_X(indsel); |
---|
2032 | coord_Y=coord_Y(indsel); |
---|
2033 | end |
---|
2034 | FF=zeros(1,length(coord_y_proj)*length(coord_x_proj)); |
---|
2035 | testFF=0; |
---|
2036 | [X,Y,Z]=meshgrid(coord_y_proj,coord_z_proj,coord_x_proj);%grid in the new coordinates |
---|
2037 | for ivar=VarIndex |
---|
2038 | VarName=FieldData.ListVarName{ivar}; |
---|
2039 | if ~( ivar==ivar_X || ivar==ivar_Y || ivar==ivar_Z || ivar==ivar_F || ivar==ivar_FF || test_anc(ivar)==1) |
---|
2040 | ivar_new=ivar_new+1; |
---|
2041 | ProjData.ListVarName=[ProjData.ListVarName {VarName}]; |
---|
2042 | ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; |
---|
2043 | if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar |
---|
2044 | ProjData.VarAttribute{ivar_new+nbcoord}=FieldData.VarAttribute{ivar}; |
---|
2045 | end |
---|
2046 | if ~isequal(ivar_FF,0) |
---|
2047 | eval(['FieldData.' VarName '=FieldData.' VarName '(indsel);']) |
---|
2048 | end |
---|
2049 | % linear interpolation |
---|
2050 | InterpFct=TriScatteredInterp(double(coord_X),double(coord_Y),double(coord_Z),double(FieldData.(VarName))); |
---|
2051 | ProjData.(VarName)=InterpFct(X,Y,Z); |
---|
2052 | % eval(['varline=reshape(ProjData.' VarName ',1,length(coord_y_proj)*length(coord_x_proj));']) |
---|
2053 | % FFlag= isnan(varline); %detect undefined values NaN |
---|
2054 | % indnan=find(FFlag); |
---|
2055 | % if ~isempty(indnan) |
---|
2056 | % varline(indnan)=zeros(size(indnan)); |
---|
2057 | % eval(['ProjData.' VarName '=reshape(varline,length(coord_y_proj),length(coord_x_proj));']) |
---|
2058 | % FF(indnan)=ones(size(indnan)); |
---|
2059 | % testFF=1; |
---|
2060 | % end |
---|
2061 | if ivar==ivar_U |
---|
2062 | ivar_U=ivar_new; |
---|
2063 | end |
---|
2064 | if ivar==ivar_V |
---|
2065 | ivar_V=ivar_new; |
---|
2066 | end |
---|
2067 | if ivar==ivar_W |
---|
2068 | ivar_W=ivar_new; |
---|
2069 | end |
---|
2070 | end |
---|
2071 | end |
---|
2072 | if testFF |
---|
2073 | ProjData.FF=reshape(FF,length(coord_y_proj),length(coord_x_proj)); |
---|
2074 | ProjData.ListVarName=[ProjData.ListVarName {'FF'}]; |
---|
2075 | ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; |
---|
2076 | ProjData.VarAttribute{ivar_new+1+nbcoord}.Role='errorflag'; |
---|
2077 | end |
---|
2078 | end |
---|
2079 | |
---|
2080 | %% case of input fields defined on a structured grid |
---|
2081 | else |
---|
2082 | VarName=FieldData.ListVarName{VarIndex(1)};%get the first variable of the cell to get the input matrix dimensions |
---|
2083 | eval(['DimValue=size(FieldData.' VarName ');'])%input matrix dimensions |
---|
2084 | DimValue(DimValue==1)=[];%remove singleton dimensions |
---|
2085 | NbDim=numel(DimValue);%update number of space dimensions |
---|
2086 | nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate |
---|
2087 | if NbDim>=3 |
---|
2088 | if NbDim>3 |
---|
2089 | errormsg='matrices with more than 3 dimensions not handled'; |
---|
2090 | return |
---|
2091 | else |
---|
2092 | if numel(find(VarType.coord))==2% the third matrix dimension does not correspond to a space coordinate |
---|
2093 | nbcolor=DimValue(3); |
---|
2094 | DimValue(3)=[]; %number of 'color' components updated |
---|
2095 | NbDim=2;% space dimension set to 2 |
---|
2096 | end |
---|
2097 | end |
---|
2098 | end |
---|
2099 | AYName=FieldData.ListVarName{VarType.coord(NbDim-1)};%name of input x coordinate (name preserved on projection) |
---|
2100 | AXName=FieldData.ListVarName{VarType.coord(NbDim)};%name of input y coordinate (name preserved on projection) |
---|
2101 | eval(['AX=FieldData.' AXName ';']) |
---|
2102 | eval(['AY=FieldData.' AYName ';']) |
---|
2103 | ListDimName=FieldData.VarDimName{VarIndex(1)}; |
---|
2104 | ProjData.ListVarName=[ProjData.ListVarName {AYName} {AXName}]; %TODO: check if it already exists in Projdata (several cells) |
---|
2105 | ProjData.VarDimName=[ProjData.VarDimName {AYName} {AXName}]; |
---|
2106 | |
---|
2107 | % for idim=1:length(ListDimName) |
---|
2108 | % DimName=ListDimName{idim}; |
---|
2109 | % if strcmp(DimName,'rgb')||strcmp(DimName,'nb_coord')||strcmp(DimName,'nb_coord_i') |
---|
2110 | % nbcolor=DimValue(idim); |
---|
2111 | % DimValue(idim)=[]; |
---|
2112 | % end |
---|
2113 | % if isequal(DimName,'nb_coord_j')% NOTE: CASE OF TENSOR NOT TREATED |
---|
2114 | % DimValue(idim)=[]; |
---|
2115 | % end |
---|
2116 | % end |
---|
2117 | Coord_z=[]; |
---|
2118 | Coord_y=[]; |
---|
2119 | Coord_x=[]; |
---|
2120 | |
---|
2121 | for idim=1:NbDim %loop on space dimensions |
---|
2122 | test_interp(idim)=0;%test for coordiate interpolation (non regular grid), =0 by default |
---|
2123 | ivar=VarType.coord(idim);% index of the variable corresponding to the current dimension |
---|
2124 | if ~isequal(ivar,0)% a variable corresponds to the dimension #idim |
---|
2125 | eval(['Coord{idim}=FieldData.' FieldData.ListVarName{ivar} ';']) ;% coord values for the input field |
---|
2126 | if numel(Coord{idim})==2 %input array defined on a regular grid |
---|
2127 | DCoord_min(idim)=(Coord{idim}(2)-Coord{idim}(1))/DimValue(idim); |
---|
2128 | else |
---|
2129 | DCoord=diff(Coord{idim});%array of coordinate derivatives for the input field |
---|
2130 | DCoord_min(idim)=min(DCoord); |
---|
2131 | DCoord_max=max(DCoord); |
---|
2132 | % test_direct(idim)=DCoord_max>0;% =1 for increasing values, 0 otherwise |
---|
2133 | if abs(DCoord_max-DCoord_min(idim))>abs(DCoord_max/1000) |
---|
2134 | msgbox_uvmat('ERROR',['non monotonic dimension variable # ' num2str(idim) ' in proj_field.m']) |
---|
2135 | return |
---|
2136 | end |
---|
2137 | test_interp(idim)=(DCoord_max-DCoord_min(idim))> 0.0001*abs(DCoord_max);% test grid regularity |
---|
2138 | end |
---|
2139 | test_direct(idim)=(DCoord_min(idim)>0); |
---|
2140 | else % no variable associated with the dimension #idim, the coordinate value is set equal to the matrix index by default |
---|
2141 | Coord_i_str=['Coord_' num2str(idim)]; |
---|
2142 | DCoord_min(idim)=1;%default |
---|
2143 | Coord{idim}=[0.5 DimValue(idim)-0.5]; |
---|
2144 | test_direct(idim)=1; |
---|
2145 | end |
---|
2146 | end |
---|
2147 | if DY==0 |
---|
2148 | DY=abs(DCoord_min(NbDim-1)); |
---|
2149 | end |
---|
2150 | npY=1+round(abs(Coord{NbDim-1}(end)-Coord{NbDim-1}(1))/DY);%nbre of points after interpol |
---|
2151 | if DX==0 |
---|
2152 | DX=abs(DCoord_min(NbDim)); |
---|
2153 | end |
---|
2154 | npX=1+round(abs(Coord{NbDim}(end)-Coord{NbDim}(1))/DX);%nbre of points after interpol |
---|
2155 | for idim=1:NbDim |
---|
2156 | if test_interp(idim) |
---|
2157 | DimValue(idim)=1+round(abs(Coord{idim}(end)-Coord{idim}(1))/abs(DCoord_min(idim)));%nbre of points after possible interpolation on a regular gri |
---|
2158 | end |
---|
2159 | end |
---|
2160 | Coord_y=linspace(Coord{NbDim-1}(1),Coord{NbDim-1}(end),npY); |
---|
2161 | test_direct_y=test_direct(NbDim-1); |
---|
2162 | Coord_x=linspace(Coord{NbDim}(1),Coord{NbDim}(end),npX); |
---|
2163 | test_direct_x=test_direct(NbDim); |
---|
2164 | DAX=DCoord_min(NbDim); |
---|
2165 | DAY=DCoord_min(NbDim-1); |
---|
2166 | minAX=min(Coord_x); |
---|
2167 | maxAX=max(Coord_x); |
---|
2168 | minAY=min(Coord_y); |
---|
2169 | maxAY=max(Coord_y); |
---|
2170 | xcorner=[minAX maxAX minAX maxAX]-ObjectData.Coord(1,1); |
---|
2171 | ycorner=[maxAY maxAY minAY minAY]-ObjectData.Coord(1,2); |
---|
2172 | xcor_new=xcorner*cos(Phi)+ycorner*sin(Phi);%coord new frame |
---|
2173 | ycor_new=-xcorner*sin(Phi)+ycorner*cos(Phi); |
---|
2174 | if ~testXMax |
---|
2175 | XMax=max(xcor_new); |
---|
2176 | end |
---|
2177 | if ~testXMin |
---|
2178 | XMin=min(xcor_new); |
---|
2179 | end |
---|
2180 | if ~testYMax |
---|
2181 | YMax=max(ycor_new); |
---|
2182 | end |
---|
2183 | if ~testYMin |
---|
2184 | YMin=min(ycor_new); |
---|
2185 | end |
---|
2186 | DXinit=(maxAX-minAX)/(DimValue(NbDim)-1); |
---|
2187 | DYinit=(maxAY-minAY)/(DimValue(NbDim-1)-1); |
---|
2188 | if DX==0 |
---|
2189 | DX=DXinit; |
---|
2190 | end |
---|
2191 | if DY==0 |
---|
2192 | DY=DYinit; |
---|
2193 | end |
---|
2194 | if NbDim==3 |
---|
2195 | DZ=(Coord{1}(end)-Coord{1}(1))/(DimValue(1)-1); |
---|
2196 | if ~test_direct(1) |
---|
2197 | DZ=-DZ; |
---|
2198 | end |
---|
2199 | Coord_z=linspace(Coord{1}(1),Coord{1}(end),DimValue(1)); |
---|
2200 | test_direct_z=test_direct(1); |
---|
2201 | end |
---|
2202 | npX=floor((XMax-XMin)/DX+1); |
---|
2203 | npY=floor((YMax-YMin)/DY+1); |
---|
2204 | if test_direct_y |
---|
2205 | coord_y_proj=linspace(YMin,YMax,npY);%abscissa of the new pixels along the line |
---|
2206 | else |
---|
2207 | coord_y_proj=linspace(YMax,YMin,npY);%abscissa of the new pixels along the line |
---|
2208 | end |
---|
2209 | if test_direct_x |
---|
2210 | coord_x_proj=linspace(XMin,XMax,npX);%abscissa of the new pixels along the line |
---|
2211 | else |
---|
2212 | coord_x_proj=linspace(XMax,XMin,npX);%abscissa of the new pixels along the line |
---|
2213 | end |
---|
2214 | |
---|
2215 | % case with no rotation and interpolation |
---|
2216 | if isequal(ProjMode,'projection') && isequal(Phi,0) && isequal(Theta,0) && isequal(Psi,0) |
---|
2217 | if ~testXMin && ~testXMax && ~testYMin && ~testYMax && NbDim==2 |
---|
2218 | ProjData=FieldData; |
---|
2219 | else |
---|
2220 | indY=NbDim-1; |
---|
2221 | if test_direct(indY) |
---|
2222 | min_indy=ceil((YMin-Coord{indY}(1))/DYinit)+1; |
---|
2223 | YIndexFirst=floor((YMax-Coord{indY}(1))/DYinit)+1; |
---|
2224 | Ybound(1)=Coord{indY}(1)+DYinit*(min_indy-1); |
---|
2225 | Ybound(2)=Coord{indY}(1)+DYinit*(YIndexFirst-1); |
---|
2226 | else |
---|
2227 | min_indy=ceil((Coord{indY}(1)-YMax)/DYinit)+1; |
---|
2228 | max_indy=floor((Coord{indY}(1)-YMin)/DYinit)+1; |
---|
2229 | Ybound(2)=Coord{indY}(1)-DYinit*(max_indy-1); |
---|
2230 | Ybound(1)=Coord{indY}(1)-DYinit*(min_indy-1); |
---|
2231 | end |
---|
2232 | if test_direct(NbDim)==1 |
---|
2233 | min_indx=ceil((XMin-Coord{NbDim}(1))/DXinit)+1; |
---|
2234 | max_indx=floor((XMax-Coord{NbDim}(1))/DXinit)+1; |
---|
2235 | Xbound(1)=Coord{NbDim}(1)+DXinit*(min_indx-1); |
---|
2236 | Xbound(2)=Coord{NbDim}(1)+DXinit*(max_indx-1); |
---|
2237 | else |
---|
2238 | min_indx=ceil((Coord{NbDim}(1)-XMax)/DXinit)+1; |
---|
2239 | max_indx=floor((Coord{NbDim}(1)-XMin)/DXinit)+1; |
---|
2240 | Xbound(2)=Coord{NbDim}(1)+DXinit*(max_indx-1); |
---|
2241 | Xbound(1)=Coord{NbDim}(1)+DXinit*(min_indx-1); |
---|
2242 | end |
---|
2243 | if NbDim==3 |
---|
2244 | DimCell(1)=[]; %suppress z variable |
---|
2245 | DimValue(1)=[]; |
---|
2246 | %structured coordinates |
---|
2247 | if test_direct(1) |
---|
2248 | iz=ceil((ObjectData.Coord(1,3)-Coord{1}(1))/DZ)+1; |
---|
2249 | else |
---|
2250 | iz=ceil((Coord{1}(1)-ObjectData.Coord(1,3))/DZ)+1; |
---|
2251 | end |
---|
2252 | end |
---|
2253 | min_indy=max(min_indy,1);% deals with margin (bound lower than the first index) |
---|
2254 | min_indx=max(min_indx,1); |
---|
2255 | max_indy=min(max_indy,DimValue(1)); |
---|
2256 | max_indx=min(max_indx,DimValue(2)); |
---|
2257 | for ivar=VarIndex% loop on non coordinate variables |
---|
2258 | VarName=FieldData.ListVarName{ivar}; |
---|
2259 | ProjData.ListVarName=[ProjData.ListVarName VarName]; |
---|
2260 | ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; |
---|
2261 | if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute)>=ivar |
---|
2262 | ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; |
---|
2263 | end |
---|
2264 | if NbDim==3 |
---|
2265 | eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,min_indy:max_indy,min_indx:max_indx));']); |
---|
2266 | else |
---|
2267 | eval(['ProjData.' VarName '=FieldData.' VarName '(min_indy:max_indy,min_indx:max_indx,:);']); |
---|
2268 | end |
---|
2269 | end |
---|
2270 | eval(['ProjData.' AYName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates |
---|
2271 | eval(['ProjData.' AXName '=[Xbound(1) Xbound(2)];']) %record the new (projected ) x coordinates |
---|
2272 | end |
---|
2273 | elseif isfield(FieldData,'A') %TO GENERALISE % case with rotation and/or interpolation |
---|
2274 | if NbDim==2 %2D case |
---|
2275 | [X,Y]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates |
---|
2276 | XIMA=ObjectData.Coord(1,1)+(X)*cos(Phi)-Y*sin(Phi);%corresponding coordinates in the original image |
---|
2277 | YIMA=ObjectData.Coord(1,2)+(X)*sin(Phi)+Y*cos(Phi); |
---|
2278 | XIMA=(XIMA-minAX)/DXinit+1;% image index along x |
---|
2279 | YIMA=(-YIMA+maxAY)/DYinit+1;% image index along y |
---|
2280 | XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line' |
---|
2281 | YIMA=reshape(round(YIMA),1,npX*npY); |
---|
2282 | flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image |
---|
2283 | if isequal(ObjectData.ProjMode,'interp_tps') |
---|
2284 | npx_interp_tps=ceil(abs(DX/DAX)); |
---|
2285 | npy_interp_tps=ceil(abs(DY/DAY)); |
---|
2286 | Minterp_tps=ones(npy_interp_tps,npx_interp_tps)/(npx_interp_tps*npy_interp_tps); |
---|
2287 | test_interp_tps=1; |
---|
2288 | else |
---|
2289 | test_interp_tps=0; |
---|
2290 | end |
---|
2291 | eval(['ProjData.' AYName '=[coord_y_proj(1) coord_y_proj(end)];']) %record the new (projected ) y coordinates |
---|
2292 | eval(['ProjData.' AXName '=[coord_x_proj(1) coord_x_proj(end)];']) %record the new (projected ) x coordinates |
---|
2293 | for ivar=VarIndex |
---|
2294 | VarName=FieldData.ListVarName{ivar}; |
---|
2295 | if test_interp(1) || test_interp(2)%interpolate on a regular grid |
---|
2296 | eval(['ProjData.' VarName '=interp2(Coord{2},Coord{1},FieldData.' VarName ',Coord_x,Coord_y'');']) %TO TEST |
---|
2297 | end |
---|
2298 | %filter the field (image) if option 'interp_tps' is used |
---|
2299 | if test_interp_tps |
---|
2300 | Aclass=class(FieldData.A); |
---|
2301 | ProjData.(VarName)=interp_tps2(Minterp_tps,FieldData.(VarName),'valid'); |
---|
2302 | if ~isequal(Aclass,'double') |
---|
2303 | ProjData.(VarName)=Aclass(FieldData.(VarName));%revert to integer values |
---|
2304 | end |
---|
2305 | end |
---|
2306 | eval(['vec_A=reshape(FieldData.' VarName ',[],nbcolor);'])%put the original image in line |
---|
2307 | %ind_in=find(flagin); |
---|
2308 | ind_out=find(~flagin); |
---|
2309 | ICOMB=(XIMA-1)*DimValue(1)+YIMA; |
---|
2310 | ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A |
---|
2311 | vec_B(flagin,1:nbcolor)=vec_A(ICOMB,:); |
---|
2312 | for icolor=1:nbcolor |
---|
2313 | vec_B(ind_out,icolor)=zeros(size(ind_out)); |
---|
2314 | end |
---|
2315 | ProjData.ListVarName=[ProjData.ListVarName VarName]; |
---|
2316 | ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; |
---|
2317 | if isfield(FieldData,'VarAttribute')&&length(FieldData.VarAttribute)>=ivar |
---|
2318 | ProjData.VarAttribute{length(ProjData.ListVarName)+nbcoord}=FieldData.VarAttribute{ivar}; |
---|
2319 | end |
---|
2320 | eval(['ProjData.' VarName '=reshape(vec_B,npY,npX,nbcolor);']); |
---|
2321 | end |
---|
2322 | ProjData.FF=reshape(~flagin,npY,npX);%false flag A FAIRE: tenir compte d'un flga ant???rieur |
---|
2323 | ProjData.ListVarName=[ProjData.ListVarName 'FF']; |
---|
2324 | ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; |
---|
2325 | ProjData.VarAttribute{length(ProjData.ListVarName)}.Role='errorflag'; |
---|
2326 | else %3D case |
---|
2327 | if ~testangle |
---|
2328 | % unstructured z coordinate |
---|
2329 | test_sup=(Coord{1}>=ObjectData.Coord(1,3)); |
---|
2330 | iz_sup=find(test_sup); |
---|
2331 | iz=iz_sup(1); |
---|
2332 | if iz>=1 & iz<=npz |
---|
2333 | %ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)]; |
---|
2334 | %ProjData.DimValue=[ProjData.DimValue npY npX]; |
---|
2335 | for ivar=VarIndex |
---|
2336 | VarName=FieldData.ListVarName{ivar}; |
---|
2337 | ProjData.ListVarName=[ProjData.ListVarName VarName]; |
---|
2338 | ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes |
---|
2339 | eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,:,:));'])% select the z index iz |
---|
2340 | %TODO : do a vertical average for a thick plane |
---|
2341 | if test_interp(2) || test_interp(3) |
---|
2342 | eval(['ProjData.' VarName '=interp2(Coord{3},Coord{2},ProjData.' VarName ',Coord_x,Coord_y'');']) |
---|
2343 | end |
---|
2344 | end |
---|
2345 | end |
---|
2346 | else |
---|
2347 | RotMatrix=rodrigues(om); |
---|
2348 | |
---|
2349 | errormsg='projection of structured coordinates on oblique plane not yet implemented'; |
---|
2350 | %TODO: use interp3 |
---|
2351 | return |
---|
2352 | end |
---|
2353 | end |
---|
2354 | end |
---|
2355 | end |
---|
2356 | |
---|
2357 | %% projection of velocity components in the rotated coordinates |
---|
2358 | if testangle |
---|
2359 | if isempty(ivar_V) |
---|
2360 | msgbox_uvmat('ERROR','v velocity component missing in proj_field.m') |
---|
2361 | return |
---|
2362 | end |
---|
2363 | UName=FieldData.ListVarName{ivar_U}; |
---|
2364 | VName=FieldData.ListVarName{ivar_V}; |
---|
2365 | eval(['ProjData.' UName '=cos(Phi)*ProjData.' UName '+ sin(Phi)*ProjData.' VName ';']) |
---|
2366 | eval(['ProjData.' VName '=cos(Theta)*(-sin(Phi)*ProjData.' UName '+ cos(Phi)*ProjData.' VName ');']) |
---|
2367 | if ~isempty(ivar_W) |
---|
2368 | WName=FieldData.ListVarName{ivar_W}; |
---|
2369 | eval(['ProjData.' VName '=ProjData.' VName '+ ProjData.' WName '*sin(Theta);'])% |
---|
2370 | eval(['ProjData.' WName '=NormVec_X*ProjData.' UName '+ NormVec_Y*ProjData.' VName '+ NormVec_Z* ProjData.' WName ';']); |
---|
2371 | end |
---|
2372 | if ~isequal(Psi,0) |
---|
2373 | eval(['ProjData.' UName '=cos(Psi)* ProjData.' UName '- sin(Psi)*ProjData.' VName ';']); |
---|
2374 | eval(['ProjData.' VName '=sin(Psi)* ProjData.' UName '+ cos(Psi)*ProjData.' VName ';']); |
---|
2375 | end |
---|
2376 | end |
---|
2377 | end |
---|
2378 | |
---|
2379 | %------------------------------------------------------------------------ |
---|
2380 | %--- transfer the global attributes |
---|
2381 | function [ProjData,errormsg]=proj_heading(FieldData,ObjectData) |
---|
2382 | %------------------------------------------------------------------------ |
---|
2383 | ProjData=[];%default |
---|
2384 | errormsg='';%default |
---|
2385 | |
---|
2386 | %% transfer error |
---|
2387 | if isfield(FieldData,'Txt') |
---|
2388 | errormsg=FieldData.Txt; %transmit erreur message |
---|
2389 | return; |
---|
2390 | end |
---|
2391 | |
---|
2392 | %% transfer global attributes |
---|
2393 | if ~isfield(FieldData,'ListGlobalAttribute') |
---|
2394 | ProjData.ListGlobalAttribute={}; |
---|
2395 | else |
---|
2396 | ProjData.ListGlobalAttribute=FieldData.ListGlobalAttribute; |
---|
2397 | end |
---|
2398 | for iattr=1:length(ProjData.ListGlobalAttribute) |
---|
2399 | AttrName=ProjData.ListGlobalAttribute{iattr}; |
---|
2400 | if isfield(FieldData,AttrName) |
---|
2401 | ProjData.(AttrName)=FieldData.(AttrName); |
---|
2402 | end |
---|
2403 | end |
---|
2404 | |
---|
2405 | %% transfer coordinate unit |
---|
2406 | if isfield(ProjData,'CoordUnit') |
---|
2407 | ProjData=rmfield(ProjData,'CoordUnit');% do not transfer by default (to avoid x/y=1 for profiles) |
---|
2408 | end |
---|
2409 | if isfield(FieldData,'CoordUnit') |
---|
2410 | if isfield(ObjectData,'CoordUnit') && ~strcmp(FieldData.CoordUnit,ObjectData.CoordUnit) |
---|
2411 | errormsg=[ObjectData.Type ' in ' ObjectData.CoordUnit ' coordinates, while field in ' FieldData.CoordUnit ]; |
---|
2412 | return |
---|
2413 | elseif strcmp(ObjectData.Type,'plane')|| strcmp(ObjectData.Type,'volume') |
---|
2414 | ProjData.CoordUnit=FieldData.CoordUnit; |
---|
2415 | end |
---|
2416 | end |
---|
2417 | |
---|
2418 | %% store the properties of the projection object |
---|
2419 | ListObject={'Name','Type','ProjMode','angle','RangeX','RangeY','RangeZ','DX','DY','DZ','Coord'}; |
---|
2420 | for ilist=1:length(ListObject) |
---|
2421 | if isfield(ObjectData,ListObject{ilist}) |
---|
2422 | val=ObjectData.(ListObject{ilist}); |
---|
2423 | if ~isempty(val) |
---|
2424 | ProjData.(['ProjObject' ListObject{ilist}])=val; |
---|
2425 | ProjData.ListGlobalAttribute=[ProjData.ListGlobalAttribute {['ProjObject' ListObject{ilist}]}]; |
---|
2426 | end |
---|
2427 | end |
---|
2428 | end |
---|
2429 | |
---|