Context Navigation

source: trunk/src/proj_field.m @ 656

Last change on this file since 656 was 655, checked in by sommeria, 11 years ago
various bugs corrected
File size: 109.6 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: