Context Navigation

source: trunk/src/proj_field.m @ 672

Last change on this file since 672 was 672, checked in by sommeria, 11 years ago
various bugs corrected
File size: 110.8 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	for icell=1:numel(CellInfo)% TODO: recalculate coordinates here to get the bounds in the rotated coordinates
1015	if isfield(CellInfo{icell},'ProjModeRequest')
1016	switch CellInfo{icell}.ProjModeRequest
1017	case 'interp_lin'
1018	ProjMode{icell}='interp_lin';
1019	case 'interp_tps'
1020	ProjMode{icell}='interp_tps';
1021	end
1022	end
1023	if strcmp(ProjMode{icell},'interp_lin')\|\|strcmp(ProjMode{icell},'interp_tps')
1024	check_grid(icell)=1;
1025	end
1026	if strcmp(CellInfo{icell}.CoordType,'grid')&&NbDimArray(icell)>=2
1027	if ~testangle && isempty(icell_grid)% if the input gridded data is not modified, choose the first one in case of multiple gridded field cells
1028	icell_grid=icell;
1029	ProjMode{icell}='projection';
1030	end
1031	check_grid(icell)=1;
1032	end
1033	end
1034	if ~isempty(find(check_grid))% if a grid is requested by the input field
1035	if isempty(icell_grid)% if the grid is not given by cell #icell_grid
1036	if ~isfield(FieldData,'XMax')
1037	FieldData=find_field_bounds(FieldData);
1038	% CellIndex=find(check_grid);
1039	% for igrid=1:numel(CellIndex)
1040	% icell=CellIndex(igrid);% TODO: recalculate coordinates here to get the bounds in the rotated coordinates
1041	% NbDim=NbDimArray(icell);
1042	% for idim=1:NbDim %loop on space dimensions
1043	% ivar=CellInfo{icell}.CoordIndex(idim);% index of the variable corresponding to the current dimension
1044	% Coord{idim}=FieldData.(FieldData.ListVarName{ivar});% coord values for the input field
1045	% if ~isequal(ivar,0)% a variable corresponds to the dimension #idim
1046	%
1047	%
1048	%
1049	% else
1050	% Coord_i_str=['Coord_' num2str(idim)];
1051	% DCoord_min(idim)=1;%default
1052	% Coord{idim}=[0.5 DimValue(idim)-0.5];
1053	% test_direct(idim)=1;
1054	% end
1055	% % default bounds, case of gridded data
1056	% if strcmp(CellInfo{icell}.CoordType,'grid')
1057	% % find default mesh
1058	% for idim=1:NbDim %loop on space dimensions
1059	% test_interp(idim)=0;%test for coordiate interpolation (non regular grid), =0 by default
1060	% ivar=CellInfo{icell}.CoordIndex(idim);% index of the variable corresponding to the current dimension
1061	% DimValue=size(FieldData.(FieldData.ListVarName{ivar}));
1062	% if ~isequal(ivar,0)% a variable corresponds to the dimension #idim
1063	% Coord{idim}=FieldData.(FieldData.ListVarName{ivar});% coord values for the input field
1064	% if numel(Coord{idim})==2 %input array defined on a regular grid
1065	% DCoord_min(idim)=(Coord{idim}(2)-Coord{idim}(1))/DimValue(idim);
1066	% Coord{idim}=linspace(Coord{idim}(1),Coord{idim}(2),DimValue(idim));
1067	% else
1068	% DCoord=diff(Coord{idim});%array of coordinate derivatives for the input field
1069	% DCoord_min(idim)=min(DCoord);
1070	% DCoord_max=max(DCoord);
1071	% if abs(DCoord_max-DCoord_min(idim))>abs(DCoord_max/1000)
1072	% msgbox_uvmat('ERROR',['non monotonic dimension variable # ' num2str(idim) ' in proj_field.m'])
1073	% return
1074	% end
1075	% test_interp(idim)=(DCoord_max-DCoord_min(idim))> 0.0001*abs(DCoord_max);% test grid regularity
1076	% end
1077	% test_direct(idim)=(DCoord_min(idim)>0);
1078	% else % no variable associated with the dimension #idim, the coordinate value is set equal to the matrix index by default
1079	% % Coord_i_str=['Coord_' num2str(idim)];
1080	% % DCoord_min(idim)=1;%default
1081	% % Coord{idim}=[0.5 DimValue(idim)-0.5];
1082	% % test_direct(idim)=1;
1083	% end
1084	% end
1085	% if isempty(DY)
1086	% DY=abs(DCoord_min(NbDim-1));
1087	% end
1088	% npY=1+round(abs(Coord{NbDim-1}(end)-Coord{NbDim-1}(1))/DY);%nbre of points after interpol
1089	% if isempty(DX)
1090	% DX=abs(DCoord_min(NbDim));
1091	% end
1092	% npX=1+round(abs(Coord{NbDim}(end)-Coord{NbDim}(1))/DX);%nbre of points after interpol
1093	% for idim=1:NbDim
1094	% if test_interp(idim)
1095	% 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
1096	% end
1097	% end
1098	% Coord_y=linspace(Coord{NbDim-1}(1),Coord{NbDim-1}(end),npY);
1099	% test_direct_y=test_direct(NbDim-1);
1100	% Coord_x=linspace(Coord{NbDim}(1),Coord{NbDim}(end),npX);
1101	% test_direct_x=test_direct(NbDim);
1102	% DAX=DCoord_min(NbDim);
1103	% DAY=DCoord_min(NbDim-1);
1104	% minAX=min(Coord_x);
1105	% maxAX=max(Coord_x);
1106	% minAY=min(Coord_y);
1107	% maxAY=max(Coord_y);
1108	% xcorner=[minAX maxAX minAX maxAX]-ObjectData.Coord(1,1);% image corners in the new coordiantes
1109	% ycorner=[maxAY maxAY minAY minAY]-ObjectData.Coord(1,2);
1110	% xcor_new=xcornercos_om+ycornersin_om;%coord new frame
1111	% ycor_new=-xcornersin_om+ycornercos_om;
1112	% if ~testXMax
1113	% XMax(igrid)=max(xcor_new);
1114	% end
1115	% if ~testXMin
1116	% XMin(igrid)=min(xcor_new);
1117	% end
1118	% if ~testYMax
1119	% YMax(igrid)=max(ycor_new);
1120	% end
1121	% if ~testYMin
1122	% YMin(igrid)=min(ycor_new);
1123	% end
1124	% DX(igrid)=(maxAX-minAX)/(DimValue(NbDim)-1);
1125	% DY(igrid)=(maxAY-minAY)/(DimValue(NbDim-1)-1);
1126	% if NbDim==3
1127	% DZ(igrid)=(Coord{1}(end)-Coord{1}(1))/(DimValue(1)-1);
1128	% if ~test_direct(1)
1129	% DZ=-DZ;
1130	% end
1131	% Coord_z=linspace(Coord{1}(1),Coord{1}(end),DimValue(1));
1132	% test_direct_z=test_direct(1);
1133	% end
1134	% else
1135	%
1136	% end
1137	% YMax=max(YMax);
1138	% YMin=min(YMin);
1139	% XMax=max(XMax);
1140	% XMin=min(XMin);
1141	end
1142	end
1143	end
1144	end
1145	if ~isempty(find(check_grid))\|\|~strcmp(ObjectData.ProjMode,'projection')%no existing gridded data used
1146	if isempty(icell_grid)\|\|~strcmp(ObjectData.ProjMode,'projection')%no existing gridded data used
1147	AYName='coord_y';
1148	AXName='coord_x';
1149	if ~strcmp(ObjectData.ProjMode,'projection')
1150	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
1151	ProjData.coord_x=[ObjectData.RangeX(1) ObjectData.RangeX(2)];
1152	coord_x_proj=ObjectData.RangeX(1):ObjectData.DX:ObjectData.RangeX(2);
1153	coord_y_proj=ObjectData.RangeY(1):ObjectData.DY:ObjectData.RangeY(2);
1154	else
1155	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
1156	ProjData.coord_x=[FieldData.XMin FieldData.XMax];
1157	coord_x_proj=FieldData.XMin:FieldData.CoordMesh:FieldData.XMax;
1158	coord_y_proj=FieldData.YMin:FieldData.CoordMesh:FieldData.YMax;
1159	end
1160	[X,YI]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates
1161	XI=ObjectData.Coord(1,1)+(X)cos(PlaneAngle(3))-YIsin(PlaneAngle(3));%corresponding coordinates in the original system
1162	YI=ObjectData.Coord(1,2)+(X)sin(PlaneAngle(3))+YIcos(PlaneAngle(3));
1163	else% we use the existing grid from field cell #icell_grid
1164	NbDim=NbDimArray(icell_grid);
1165	AYName=FieldData.ListVarName{CellInfo{icell_grid}.CoordIndex(NbDim-1)};%name of input x coordinate (name preserved on projection)
1166	AXName=FieldData.ListVarName{CellInfo{icell_grid}.CoordIndex(NbDim)};%name of input y coordinate (name preserved on projection)
1167	ProjData.(AYName)=FieldData.(AYName); % new (projected ) y coordinates
1168	ProjData.(AXName)=FieldData.(AXName); % new (projected ) y coordinates
1169	end
1170	ProjData.ListVarName={AYName,AXName};
1171	ProjData.VarDimName={AYName,AXName};
1172	ProjData.VarAttribute={[],[]};
1173	end
1174
1175	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1176	% LOOP ON FIELD CELLS, PROJECT VARIABLES
1177	% CellVarIndex=cells of variable index arrays
1178	%ivar_new=0; % index of the current variable in the projected field
1179	% icoord=0;
1180	nbcoord=0;%number of added coordinate variables brought by projection
1181	%nbvar=0;
1182	vector_x_proj=[];
1183	vector_y_proj=[];
1184	for icell=1:length(CellInfo)
1185	NbDim=NbDimArray(icell);
1186	if NbDim<2
1187	continue % only cells represnting 2D or 3D fields are involved
1188	end
1189	VarIndex=CellInfo{icell}.VarIndex;% indices of the selected variables in the list FieldData.ListVarName
1190	% ivar_U=[];ivar_V=[];ivar_W=[];
1191	% if isfield(CellInfo{icell},'VarIndex_vector_x_tps')&&isfield(CellInfo{icell},'VarIndex_vector_y_tps')
1192	% ivar_U=CellInfo{icell}.VarIndex_vector_x_tps;
1193	% ivar_V=CellInfo{icell}.VarIndex_vector_y_tps;
1194	% elseif isfield(CellInfo{icell},'VarIndex_vector_x')&&isfield(CellInfo{icell},'VarIndex_vector_y')
1195	% ivar_U=CellInfo{icell}.VarIndex_vector_x;
1196	% ivar_V=CellInfo{icell}.VarIndex_vector_y;
1197	% end
1198	% if isfield(CellInfo{icell},'VarIndex_vector_z')
1199	% ivar_W=CellInfo{icell}.VarIndex_vector_z;
1200	% end
1201	%dimensions
1202	DimCell=FieldData.VarDimName{VarIndex(1)};
1203	if ischar(DimCell)
1204	DimCell={DimCell};%name of dimensions
1205	end
1206	coord_z=0;%default
1207	ListVarName={};% initiate list of projected variables for cell # icell
1208	VarDimName={};% initiate coresponding list of dimensions for cell # icell
1209	VarAttribute={};% initiate coresponding list of var attributes for cell # icell
1210	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1211	switch CellInfo{icell}.CoordType
1212
1213	case 'scattered'
1214	%% case of input fields with unstructured coordinates (applies for projMode ='projection' or 'interp_lin')
1215	if strcmp(ProjMode{icell},'interp_tps')
1216	continue %skip for next cell (needs tps field cell)
1217	end
1218	coord_x=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)});% initial x coordinates
1219	coord_y=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)});% initial y coordinates
1220	check3D=(numel(CellInfo{icell}.CoordIndex)==3);
1221	if check3D
1222	coord_z=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(1)});
1223	end
1224
1225	% translate initial coordinates to account for the new origin
1226	coord_x=coord_x-ObjectData.Coord(1,1);
1227	coord_y=coord_y-ObjectData.Coord(1,2);
1228	if check3D
1229	coord_z=coord_z-ObjectData.Coord(1,3);
1230	end
1231
1232	% selection of the vectors in the projection range (3D case)
1233	if check3D && width > 0
1234	%components of the unitiy vector normal to the projection plane
1235	fieldZ=norm_plane(1)coord_x + norm_plane(2)coord_y+ norm_plane(3)*coord_z;% distance to the plane
1236	indcut=find(abs(fieldZ) <= width);
1237	for ivar=VarIndex
1238	VarName=FieldData.ListVarName{ivar};
1239	FieldData.(VarName)=FieldData.(VarName)(indcut);
1240	end
1241	coord_x=coord_x(indcut);
1242	coord_y=coord_y(indcut);
1243	coord_z=coord_z(indcut);
1244	end
1245
1246	%rotate coordinates if needed: coord_X,coord_Y= = coordinates in the new plane
1247	Psi=PlaneAngle(1);
1248	Theta=PlaneAngle(2);
1249	Phi=PlaneAngle(3);
1250	if testangle && ~test90y && ~test90x;%=1 for slanted plane
1251	coord_X=(coord_x cos(Phi) + coord_y sin(Phi));
1252	coord_Y=(-coord_x sin(Phi) + coord_y cos(Phi))*cos(Theta);
1253	coord_Y=coord_Y+coord_z *sin(Theta);
1254	coord_X=(coord_X cos(Psi) - coord_Y sin(Psi));%A VERIFIER
1255	coord_Y=(coord_X sin(Psi) + coord_Y cos(Psi));
1256	else
1257	coord_X=coord_x;
1258	coord_Y=coord_y;
1259	end
1260
1261	%restriction to the range of X and Y if imposed by the projection object
1262	testin=ones(size(coord_X)); %default
1263	testbound=0;
1264	if testXMin
1265	testin=testin & (coord_X >= XMin);
1266	testbound=1;
1267	end
1268	if testXMax
1269	testin=testin & (coord_X <= XMax);
1270	testbound=1;
1271	end
1272	if testYMin
1273	testin=testin & (coord_Y >= YMin);
1274	testbound=1;
1275	end
1276	if testYMin
1277	testin=testin & (coord_Y <= YMax);
1278	testbound=1;
1279	end
1280	if testbound
1281	indcut=find(testin);
1282	if isempty(indcut)
1283	errormsg='data outside the bounds of the projection object';
1284	return
1285	end
1286	for ivar=VarIndex
1287	VarName=FieldData.ListVarName{ivar};
1288	FieldData.(VarName)=FieldData.(VarName)(indcut);
1289	end
1290	coord_X=coord_X(indcut);
1291	coord_Y=coord_Y(indcut);
1292	if check3D
1293	coord_Z=coord_Z(indcut);
1294	end
1295	end
1296
1297	% two cases of projection for scattered coordinates
1298	switch ProjMode{icell}
1299	case 'projection'
1300	nbvar=0;
1301	%nbvar=numel(ProjData.ListVarName);
1302	for ivar=VarIndex %transfer variables to the projection plane
1303	VarName=FieldData.ListVarName{ivar};
1304	if ivar==CellInfo{icell}.CoordIndex(end)
1305	ProjData.(VarName)=coord_X;
1306	elseif ivar==CellInfo{icell}.CoordIndex(end-1) % y coordinate
1307	ProjData.(VarName)=coord_Y;
1308	elseif ~(check3D && ivar==CellInfo{icell}.CoordIndex(1)) % other variables (except Z coordinate wyhich is not reproduced)
1309	ProjData.(VarName)=FieldData.(VarName);
1310	end
1311	if ~(check3D && ivar==CellInfo{icell}.CoordIndex(1))
1312	ListVarName=[ListVarName VarName];
1313	VarDimName=[VarDimName DimCell];
1314	nbvar=nbvar+1;
1315	if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar
1316	VarAttribute{nbvar}=FieldData.VarAttribute{ivar};
1317	end
1318	end
1319	end
1320	case 'interp_lin'%interpolate data on a regular grid
1321	if isfield(CellInfo{icell},'VarIndex_errorflag')
1322	VarName_FF=FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag};
1323	indsel=find(FieldData.(VarName_FF)==0);
1324	coord_X=coord_X(indsel);
1325	coord_Y=coord_Y(indsel);
1326	for ivar=1:numel(CellInfo{icell}.VarIndex)
1327	VarName=FieldData.ListVarName{CellInfo{icell}.VarIndex(ivar)};
1328	FieldData.(VarName)=FieldData.(VarName)(indsel);
1329	end
1330	end
1331	% interpolate and calculate field on the grid
1332	[VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp([coord_X coord_Y],FieldData,CellInfo{icell}.FieldName,XI,YI);
1333
1334	if isfield(CellInfo{icell},'CheckSub') && CellInfo{icell}.CheckSub && ~isempty(vector_x_proj)
1335	ProjData.(ListVarName{vector_x_proj})=ProjData.(ListVarName{vector_x_proj})-VarVal{1};
1336	ProjData.(ListVarName{vector_y_proj})=ProjData.(ListVarName{vector_y_proj})-VarVal{2};
1337	else
1338	VarDimName=cell(size(ListVarName));
1339	for ilist=1:numel(ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise
1340	ListVarName{ilist}=regexprep(ListVarName{ilist},'(.+','');
1341	if ~isempty(find(strcmp(ListVarName{ilist},ProjData.ListVarName)))
1342	ListVarName{ilist}=[ListVarName{ilist} '_1'];
1343	end
1344	ProjData.(ListVarName{ilist})=VarVal{ilist};
1345	VarDimName{ilist}={'coord_y','coord_x'};
1346	end
1347	end
1348	end
1349
1350	case 'tps'
1351	%% case of tps data (applies only in interp_tps mode)
1352	if strcmp(ProjMode{icell},'interp_tps')
1353	Coord=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex});
1354	NbCentres=FieldData.(FieldData.ListVarName{CellInfo{icell}.NbCentres_tps});
1355	SubRange=FieldData.(FieldData.ListVarName{CellInfo{icell}.SubRange_tps});
1356	if isfield(CellInfo{icell},'VarIndex_vector_x_tps')&&isfield(CellInfo{icell},'VarIndex_vector_y_tps')
1357	FieldVar=cat(3,FieldData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x_tps}),FieldData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y_tps}));
1358	end
1359	[DataOut,VarAttribute,errormsg]=calc_field_tps(Coord,NbCentres,SubRange,FieldVar,CellInfo{icell}.FieldName,cat(3,XI,YI));
1360	ListVarName=(fieldnames(DataOut))';
1361	VarDimName=cell(size(ListVarName));
1362	for ilist=1:numel(ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise
1363	VarName=ListVarName{ilist};
1364	ProjData.(VarName)=DataOut.(VarName);
1365	VarDimName{ilist}={'coord_y','coord_x'};
1366	end
1367	end
1368
1369	case 'grid'
1370	%% case of input fields defined on a structured grid
1371	VarName=FieldData.ListVarName{VarIndex(1)};%get the first variable of the cell to get the input matrix dimensions
1372	DimValue=size(FieldData.(VarName));%input matrix dimensions
1373	DimValue(DimValue==1)=[];%remove singleton dimensions
1374	NbDim=numel(DimValue);%update number of space dimensions
1375	nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate
1376	if NbDim>=3
1377	if NbDim>3
1378	errormsg='matrices with more than 3 dimensions not handled';
1379	return
1380	else
1381	if numel(CellInfo{icell}.CoordIndex)==2% the third matrix dimension does not correspond to a space coordinate
1382	nbcolor=DimValue(3);
1383	DimValue(3)=[]; %number of 'color' components updated
1384	NbDim=2;% space dimension set to 2
1385	end
1386	end
1387	end
1388	Coord_z=[];
1389	Coord_y=[];
1390	Coord_x=[];
1391
1392
1393
1394	% % default bounds
1395	% minAX=min(Coord_x);
1396	% maxAX=max(Coord_x);
1397	% minAY=min(Coord_y);
1398	% maxAY=max(Coord_y);
1399	% xcorner=[minAX maxAX minAX maxAX]-ObjectData.Coord(1,1);% image corners in the new coordiantes
1400	% ycorner=[maxAY maxAY minAY minAY]-ObjectData.Coord(1,2);
1401	% xcor_new=xcornercos_om+ycornersin_om;%coord new frame
1402	% ycor_new=-xcornersin_om+ycornercos_om;
1403	% if ~testXMax
1404	% XMax=max(xcor_new);
1405	% end
1406	% if ~testXMin
1407	% XMin=min(xcor_new);
1408	% end
1409	% if ~testYMax
1410	% YMax=max(ycor_new);
1411	% end
1412	% if ~testYMin
1413	% YMin=min(ycor_new);
1414	% end
1415	% DXinit=(maxAX-minAX)/(DimValue(NbDim)-1);
1416	% DYinit=(maxAY-minAY)/(DimValue(NbDim-1)-1);
1417	% if DX==0
1418	% DX=DXinit;
1419	% end
1420	% if DY==0
1421	% DY=DYinit;
1422	% end
1423	% if NbDim==3
1424	% DZ=(Coord{1}(end)-Coord{1}(1))/(DimValue(1)-1);
1425	% if ~test_direct(1)
1426	% DZ=-DZ;
1427	% end
1428	% Coord_z=linspace(Coord{1}(1),Coord{1}(end),DimValue(1));
1429	% test_direct_z=test_direct(1);
1430	% end
1431	%define new coordiantes if not already done by the definition of a projection grid
1432	if ~check_grid
1433
1434
1435	end
1436
1437	if testangle
1438	ProjMode{icell}='interp_lin'; %request linear interpolation for projection on a tilted plane
1439	end
1440
1441	% npX=floor((XMax-XMin)/DX+1);
1442	% npY=floor((YMax-YMin)/DY+1);
1443	% if test_direct_y
1444	% coord_y_proj=linspace(YMin,YMax,npY);%abscissa of the new pixels along the line
1445	% else
1446	% coord_y_proj=linspace(YMax,YMin,npY);%abscissa of the new pixels along the line
1447	% end
1448	% if test_direct_x
1449	% coord_x_proj=linspace(XMin,XMax,npX);%abscissa of the new pixels along the line
1450	% else
1451	% coord_x_proj=linspace(XMax,XMin,npX);%abscissa of the new pixels along the line
1452	% end
1453	if isequal(ProjMode{icell},'projection')% && (~testangle \|\| test90y \|\| test90x)
1454	if NbDim==2 && ~testXMin && ~testXMax && ~testYMin && ~testYMax% no range restriction
1455	ListVarName=[ListVarName FieldData.ListVarName(VarIndex)];
1456	VarDimName=[VarDimName FieldData.VarDimName(VarIndex)];
1457	if isfield(FieldData,'VarAttribute')
1458	VarAttribute=[VarAttribute FieldData.VarAttribute(VarIndex)];
1459	end
1460
1461	ProjData.(AYName)=FieldData.(AYName);
1462	ProjData.(AXName)=FieldData.(AXName);
1463	for ivar=VarIndex
1464	VarName=FieldData.ListVarName{ivar};
1465	ProjData.(VarName)=FieldData.(VarName);% no change by projection
1466	end
1467	else
1468	indY=NbDim-1;
1469	if test_direct(indY)
1470	min_indy=ceil((YMin-Coord{indY}(1))/DYinit)+1;
1471	max_indy=floor((YMax-Coord{indY}(1))/DYinit)+1;
1472	Ybound(1)=Coord{indY}(1)+DYinit*(min_indy-1);
1473	Ybound(2)=Coord{indY}(1)+DYinit*(max_indy-1);
1474	else
1475	min_indy=ceil((Coord{indY}(1)-YMax)/DYinit)+1;
1476	max_indy=floor((Coord{indY}(1)-YMin)/DYinit)+1;
1477	Ybound(2)=Coord{indY}(1)-DYinit*(max_indy-1);
1478	Ybound(1)=Coord{indY}(1)-DYinit*(min_indy-1);
1479	end
1480	if test_direct(NbDim)==1
1481	min_indx=ceil((XMin-Coord{NbDim}(1))/DXinit)+1;
1482	max_indx=floor((XMax-Coord{NbDim}(1))/DXinit)+1;
1483	Xbound(1)=Coord{NbDim}(1)+DXinit*(min_indx-1);
1484	Xbound(2)=Coord{NbDim}(1)+DXinit*(max_indx-1);
1485	else
1486	min_indx=ceil((Coord{NbDim}(1)-XMax)/DXinit)+1;
1487	max_indx=floor((Coord{NbDim}(1)-XMin)/DXinit)+1;
1488	Xbound(2)=Coord{NbDim}(1)+DXinit*(max_indx-1);
1489	Xbound(1)=Coord{NbDim}(1)+DXinit*(min_indx-1);
1490	end
1491	min_indy=max(min_indy,1);% deals with margin (bound lower than the first index)
1492	min_indx=max(min_indx,1);
1493	if test90y
1494	ind_new=[3 2 1];
1495	DimCell={AYProjName,AXProjName};
1496	iz=ceil((ObjectData.Coord(1,1)-Coord{3}(1))/DX)+1;
1497	for ivar=VarIndex
1498	VarName=FieldData.ListVarName{ivar};
1499	ListVarName=[ListVarName VarName];
1500	VarDimName=[VarDimName {DimCell}];
1501	VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes
1502	eval(['ProjData.' VarName '=permute(FieldData.' VarName ',ind_new);'])% permute x and z indices for 90 degree rotation
1503	eval(['ProjData.' VarName '=squeeze(ProjData.' VarName '(iz,:,:));'])% select the z index iz
1504	end
1505	eval(['ProjData.' AYName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates
1506	eval(['ProjData.' AXName '=[Coord{1}(end),Coord{1}(1)];']) %record the new (projected ) x coordinates
1507	else
1508	if NbDim==3
1509	DimCell(1)=[]; %suppress z variable
1510	DimValue(1)=[];
1511	if test_direct(1)
1512	iz=ceil((ObjectData.Coord(1,3)-Coord{1}(1))/DZ)+1;
1513	else
1514	iz=ceil((Coord{1}(1)-ObjectData.Coord(1,3))/DZ)+1;
1515	end
1516	end
1517	max_indy=min(max_indy,DimValue(1));%introduce bounds in y and x indices
1518	max_indx=min(max_indx,DimValue(2));
1519	for ivar=VarIndex% loop on non coordinate variables
1520	VarName=FieldData.ListVarName{ivar};
1521	ListVarName=[ListVarName VarName];
1522	VarDimName=[VarDimName {DimCell}];
1523	if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute)>=ivar
1524	VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar};
1525	end
1526	if NbDim==3
1527	eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,min_indy:max_indy,min_indx:max_indx));']);
1528	else
1529	eval(['ProjData.' VarName '=FieldData.' VarName '(min_indy:max_indy,min_indx:max_indx,:);']);
1530	end
1531	end
1532	eval(['ProjData.' AYName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates
1533	eval(['ProjData.' AXName '=[Xbound(1) Xbound(2)];']) %record the new (projected ) x coordinates
1534	end
1535	end
1536	else % case with interpolation
1537	if NbDim==2 %2D case
1538	if isequal(ProjMode{icell},'interp_tps')
1539	npx_interp_tps=ceil(abs(DX/DAX));
1540	npy_interp_tps=ceil(abs(DY/DAY));
1541	Minterp_tps=ones(npy_interp_tps,npx_interp_tps)/(npx_interp_tps*npy_interp_tps);
1542	test_interp_tps=1;
1543	else
1544	test_interp_tps=0;
1545	end
1546	coord_x_proj=XMin:DX:XMax;
1547	coord_y_proj=YMin:DY:YMax;
1548	[X,YI]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates
1549	XI=ObjectData.Coord(1,1)+(X)cos(PlaneAngle(3))-YIsin(PlaneAngle(3));%corresponding coordinates in the original system
1550	YI=ObjectData.Coord(1,2)+(X)sin(PlaneAngle(3))+YIcos(PlaneAngle(3));
1551	Coord{1}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(1)});
1552	Coord{2}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(2)});
1553	if numel(Coord{1})==2% x coordiante defiend by its bounds, get the whole set
1554	Coord{1}=linspace(Coord{1}(1),Coord{1}(2),CellInfo{icell}.CoordSize(1));
1555	end
1556	if numel(Coord{2})==2% y coordiante defiend by its bounds, get the whole set
1557	Coord{2}=linspace(Coord{2}(1),Coord{2}(2),CellInfo{icell}.CoordSize(2));
1558	end
1559	[X,Y]=meshgrid(Coord{2},Coord{1});%initial coordinates
1560	%name of error flag variable
1561	FFName='FF';%default name (if not already used)
1562	if isfield(ProjData,'FF')
1563	ind=1;
1564	while isfield(ProjData,['FF_' num2str(ind)])
1565	ind=ind+1;
1566	end
1567	FFName=['FF_' num2str(ind)];% append an index to the name of error flag, FF_1,FF_2...
1568	end
1569	% project all variables in the cell
1570	for ivar=VarIndex
1571	VarName=FieldData.ListVarName{ivar};
1572	if size(FieldData.(VarName),3)==1
1573	ProjData.(VarName)=interp2(X,Y,double(FieldData.(VarName)),XI,YI,'*linear');
1574	else
1575	ProjData.(VarName)=interp2(X,Y,double(FieldData.(VarName)(:,:,1)),XI,YI,'*linear');
1576	for icolor=2:size(FieldData.(VarName),3)% project 'color' components
1577	ProjData.(VarName)=cat(3,ProjData.(VarName),interp2(X,Y,double(FieldData.(VarName)(:,:,icolor)),XI,YI,'*linear')); %TO TEST
1578	end
1579	end
1580	ListVarName=[ListVarName VarName];
1581	DimCell(1:2)={AYName,AXName};
1582	VarDimName=[VarDimName {DimCell}];
1583	if isfield(FieldData,'VarAttribute')&&length(FieldData.VarAttribute)>=ivar
1584	VarAttribute{length(ListVarName)+nbcoord}=FieldData.VarAttribute{ivar};
1585	end;
1586	ProjData.(FFName)=isnan(ProjData.(VarName));%detact NaN (points outside the interpolation range)
1587	ProjData.(VarName)(ProjData.(FFName))=0; %set to 0 the NaN data
1588	end
1589	%update list of variables with error flag
1590	ListVarName=[ListVarName FFName];
1591	VarDimName=[VarDimName {DimCell}];
1592	VarAttribute{numel(ListVarName)}.Role='errorflag';
1593	elseif ~testangle
1594	% unstructured z coordinate
1595	test_sup=(Coord{1}>=ObjectData.Coord(1,3));
1596	iz_sup=find(test_sup);
1597	iz=iz_sup(1);
1598	if iz>=1 & iz<=npz
1599	%ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)];
1600	%ProjData.DimValue=[ProjData.DimValue npY npX];
1601	for ivar=VarIndex
1602	VarName=FieldData.ListVarName{ivar};
1603	ListVarName=[ListVarName VarName];
1604	VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes
1605	eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,:,:));'])% select the z index iz
1606	%TODO : do a vertical average for a thick plane
1607	if test_interp(2) \|\| test_interp(3)
1608	eval(['ProjData.' VarName '=interp2(Coord{3},Coord{2},ProjData.' VarName ',Coord_x,Coord_y'');'])
1609	end
1610	end
1611	end
1612	else
1613	errormsg='projection of structured coordinates on oblique plane not yet implemented';
1614	%TODO: use interp_lin3
1615	return
1616	end
1617	end
1618	end
1619	% update the global list of projected variables:
1620	ProjData.ListVarName=[ProjData.ListVarName ListVarName];
1621	ProjData.VarDimName=[ProjData.VarDimName VarDimName];
1622	ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute];
1623
1624	%% projection of velocity components in the rotated coordinates
1625	if testangle
1626	ivar_U=[];ivar_V=[];ivar_W=[];
1627	for ivar=1:numel(VarAttribute)
1628	if isfield(VarAttribute{ivar},'Role')
1629	if strcmp(VarAttribute{ivar}.Role,'vector_x')
1630	ivar_U=ivar;
1631	elseif strcmp(VarAttribute{ivar}.Role,'vector_y')
1632	ivar_V=ivar;
1633	elseif strcmp(VarAttribute{ivar}.Role,'vector_z')
1634	ivar_W=ivar;
1635	end
1636	end
1637	end
1638	if ~isempty(ivar_U)
1639	if isempty(ivar_V)
1640	msgbox_uvmat('ERROR','v velocity component missing in proj_field.m')
1641	return
1642	else
1643	UName=ListVarName{ivar_U};
1644	VName=ListVarName{ivar_V};
1645	ProjData.(UName)=cos(PlaneAngle(3))ProjData.(UName)+ sin(PlaneAngle(3))ProjData.(VName);
1646	ProjData.(VName)=(-sin(PlaneAngle(3))ProjData.(UName)+ cos(PlaneAngle(3))ProjData.(VName));
1647	if ~isempty(ivar_W)
1648	WName=FieldData.ListVarName{ivar_W};
1649	eval(['ProjData.' VName '=ProjData.' VName '+ ProjData.' WName '*sin(Theta);'])%
1650	eval(['ProjData.' WName '=NormVec_XProjData.' UName '+ NormVec_YProjData.' VName '+ NormVec_Z* ProjData.' WName ';']);
1651	end
1652	% if ~isequal(Psi,0)
1653	% eval(['ProjData.' UName '=cos(Psi)* ProjData.' UName '- sin(Psi)*ProjData.' VName ';']);
1654	% eval(['ProjData.' VName '=sin(Psi)* ProjData.' UName '+ cos(Psi)*ProjData.' VName ';']);
1655	% end
1656	end
1657	end
1658	end
1659	end
1660	% %prepare substraction in case of two input fields
1661	% SubData.ListVarName={};
1662	% SubData.VarDimName={};
1663	% SubData.VarAttribute={};
1664	% check_remove=zeros(size(ProjData.ListVarName));
1665	% for iproj=1:numel(ProjData.VarAttribute)
1666	% if isfield(ProjData.VarAttribute{iproj},'CheckSub')&&isequal(ProjData.VarAttribute{iproj}.CheckSub,1)
1667	% VarName=ProjData.ListVarName{iproj};
1668	% SubData.ListVarName=[SubData.ListVarName {VarName}];
1669	% SubData.VarDimName=[SubData.VarDimName ProjData.VarDimName{iproj}];
1670	% SubData.VarAttribute=[SubData.VarAttribute ProjData.VarAttribute{iproj}];
1671	% SubData.(VarName)=ProjData.(VarName);
1672	% check_remove(iproj)=1;
1673	% end
1674	% end
1675	% if ~isempty(find(check_remove))
1676	% ind_remove=find(check_remove);
1677	% ProjData.ListVarName(ind_remove)=[];
1678	% ProjData.VarDimName(ind_remove)=[];
1679	% ProjData.VarAttribute(ind_remove)=[];
1680	% ProjData=sub_field(ProjData,[],SubData);
1681	% end
1682
1683	%-----------------------------------------------------------------
1684	%projection in a volume
1685	function [ProjData,errormsg] = proj_volume(FieldData, ObjectData)
1686
1687	%-----------------------------------------------------------------
1688	ProjData=FieldData;%default output
1689
1690	%% axis origin
1691	if isempty(ObjectData.Coord)
1692	ObjectData.Coord(1,1)=0;%origin of the plane set to [0 0] by default
1693	ObjectData.Coord(1,2)=0;
1694	ObjectData.Coord(1,3)=0;
1695	end
1696
1697	%% rotation angles
1698	VolumeAngle=[0 0 0];
1699	norm_plane=[0 0 1];
1700	if isfield(ObjectData,'Angle')&& isequal(size(ObjectData.Angle),[1 3])&& ~isequal(ObjectData.Angle,[0 0 0])
1701	PlaneAngle=ObjectData.Angle;
1702	VolumeAngle=ObjectData.Angle;
1703	om=norm(VolumeAngle);%norm of rotation angle in radians
1704	OmAxis=VolumeAngle/om; %unit vector marking the rotation axis
1705	cos_om=cos(pi*om/180);
1706	sin_om=sin(pi*om/180);
1707	coeff=OmAxis(3)*(1-cos_om);
1708	%components of the unity vector norm_plane normal to the projection plane
1709	norm_plane(1)=OmAxis(1)coeff+OmAxis(2)sin_om;
1710	norm_plane(2)=OmAxis(2)coeff-OmAxis(1)sin_om;
1711	norm_plane(3)=OmAxis(3)*coeff+cos_om;
1712	end
1713	testangle=~isequal(VolumeAngle,[0 0 0]);
1714
1715	%% mesh sizes DX, DY, DZ
1716	DX=0;
1717	DY=0; %default
1718	DZ=0;
1719	if isfield(ObjectData,'DX')&~isempty(ObjectData.DX)
1720	DX=abs(ObjectData.DX);%mesh of interpolation points
1721	end
1722	if isfield(ObjectData,'DY')&~isempty(ObjectData.DY)
1723	DY=abs(ObjectData.DY);%mesh of interpolation points
1724	end
1725	if isfield(ObjectData,'DZ')&~isempty(ObjectData.DZ)
1726	DZ=abs(ObjectData.DZ);%mesh of interpolation points
1727	end
1728	if ~strcmp(ProjMode,'projection') && (DX==0\|\|DY==0\|\|DZ==0)
1729	errormsg='grid mesh DX , DY or DZ is missing';
1730	return
1731	end
1732
1733	%% extrema along each axis
1734	testXMin=0;
1735	testXMax=0;
1736	testYMin=0;
1737	testYMax=0;
1738	if isfield(ObjectData,'RangeX')
1739	XMin=min(ObjectData.RangeX);
1740	XMax=max(ObjectData.RangeX);
1741	testXMin=XMax>XMin;
1742	testXMax=1;
1743	end
1744	if isfield(ObjectData,'RangeY')
1745	YMin=min(ObjectData.RangeY);
1746	YMax=max(ObjectData.RangeY);
1747	testYMin=YMax>YMin;
1748	testYMax=1;
1749	end
1750	width=0;%default width of the projection band
1751	if isfield(ObjectData,'RangeZ')
1752	ZMin=min(ObjectData.RangeZ);
1753	ZMax=max(ObjectData.RangeZ);
1754	testZMin=ZMax>ZMin;
1755	testZMax=1;
1756	end
1757
1758	%% initiate Matlab structure for physical field
1759	[ProjData,errormsg]=proj_heading(FieldData,ObjectData);
1760	ProjData.NbDim=3;
1761	ProjData.ListVarName={};
1762	ProjData.VarDimName={};
1763	if ~isequal(DX,0)&& ~isequal(DY,0)
1764	ProjData.CoordMesh=sqrt(DX*DY);%define typical data mesh, useful for mouse selection in plots
1765	elseif isfield(FieldData,'CoordMesh')
1766	ProjData.CoordMesh=FieldData.CoordMesh;
1767	end
1768
1769	error=0;%default
1770	flux=0;
1771	testfalse=0;
1772	ListIndex={};
1773
1774	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1775	%% group the variables (fields of 'FieldData') in cells of variables with the same dimensions
1776	%-----------------------------------------------------------------
1777	idimvar=0;
1778	% LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS
1779	% CellVarIndex=cells of variable index arrays
1780	ivar_new=0; % index of the current variable in the projected field
1781	icoord=0;
1782	nbcoord=0;%number of added coordinate variables brought by projection
1783	nbvar=0;
1784	for icell=1:length(CellVarIndex)
1785	NbDim=NbDimVec(icell);
1786	if NbDim<3
1787	continue
1788	end
1789	VarIndex=CellVarIndex{icell};% indices of the selected variables in the list FieldData.ListVarName
1790	VarType=VarTypeCell{icell};
1791	ivar_X=VarType.coord_x;
1792	ivar_Y=VarType.coord_y;
1793	ivar_Z=VarType.coord_z;
1794	ivar_U=VarType.vector_x;
1795	ivar_V=VarType.vector_y;
1796	ivar_W=VarType.vector_z;
1797	ivar_C=VarType.scalar ;
1798	ivar_Anc=VarType.ancillary;
1799	test_anc=zeros(size(VarIndex));
1800	test_anc(ivar_Anc)=ones(size(ivar_Anc));
1801	ivar_F=VarType.warnflag;
1802	ivar_FF=VarType.errorflag;
1803	check_unstructured_coord=~isempty(ivar_X) && ~isempty(ivar_Y);
1804	DimCell=FieldData.VarDimName{VarIndex(1)};
1805	if ischar(DimCell)
1806	DimCell={DimCell};%name of dimensions
1807	end
1808
1809	%% case of input fields with unstructured coordinates
1810	if check_unstructured_coord
1811	XName=FieldData.ListVarName{ivar_X};
1812	YName=FieldData.ListVarName{ivar_Y};
1813	eval(['coord_x=FieldData.' XName ';'])
1814	eval(['coord_y=FieldData.' YName ';'])
1815	if length(ivar_Z)==1
1816	ZName=FieldData.ListVarName{ivar_Z};
1817	eval(['coord_z=FieldData.' ZName ';'])
1818	end
1819
1820	% translate initial coordinates
1821	coord_x=coord_x-ObjectData.Coord(1,1);
1822	coord_y=coord_y-ObjectData.Coord(1,2);
1823	if ~isempty(ivar_Z)
1824	coord_z=coord_z-ObjectData.Coord(1,3);
1825	end
1826
1827	% selection of the vectors in the projection range
1828	% if length(ivar_Z)==1 && width > 0
1829	% %components of the unitiy vector normal to the projection plane
1830	% fieldZ=NormVec_Xcoord_x + NormVec_Ycoord_y+ NormVec_Z*coord_z;% distance to the plane
1831	% indcut=find(abs(fieldZ) <= width);
1832	% for ivar=VarIndex
1833	% VarName=FieldData.ListVarName{ivar};
1834	% eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);'])
1835	% % A VOIR : CAS DE VAR STRUCTUREE MAIS PAS GRILLE REGULIERE : INTERPOLER SUR GRILLE REGULIERE
1836	% end
1837	% coord_x=coord_x(indcut);
1838	% coord_y=coord_y(indcut);
1839	% coord_z=coord_z(indcut);
1840	% end
1841
1842	%rotate coordinates if needed: TODO modify
1843	if testangle
1844	coord_X=(coord_x cos(Phi) + coord_y sin(Phi));
1845	coord_Y=(-coord_x sin(Phi) + coord_y cos(Phi))*cos(Theta);
1846	if ~isempty(ivar_Z)
1847	coord_Y=coord_Y+coord_z *sin(Theta);
1848	end
1849
1850	coord_X=(coord_X cos(Psi) - coord_Y sin(Psi));%A VERIFIER
1851	coord_Y=(coord_X sin(Psi) + coord_Y cos(Psi));
1852
1853	else
1854	coord_X=coord_x;
1855	coord_Y=coord_y;
1856	coord_Z=coord_z;
1857	end
1858	%restriction to the range of x and y if imposed
1859	testin=ones(size(coord_X)); %default
1860	testbound=0;
1861	if testXMin
1862	testin=testin & (coord_X >= XMin);
1863	testbound=1;
1864	end
1865	if testXMax
1866	testin=testin & (coord_X <= XMax);
1867	testbound=1;
1868	end
1869	if testYMin
1870	testin=testin & (coord_Y >= YMin);
1871	testbound=1;
1872	end
1873	if testYMax
1874	testin=testin & (coord_Y <= YMax);
1875	testbound=1;
1876	end
1877	if testbound
1878	indcut=find(testin);
1879	for ivar=VarIndex
1880	VarName=FieldData.ListVarName{ivar};
1881	eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);'])
1882	end
1883	coord_X=coord_X(indcut);
1884	coord_Y=coord_Y(indcut);
1885	if length(ivar_Z)==1
1886	coord_Z=coord_Z(indcut);
1887	end
1888	end
1889	% different cases of projection
1890	if isequal(ObjectData.ProjMode,'projection')%%%%%%% NOT USED %%%%%%%%%%
1891	for ivar=VarIndex %transfer variables to the projection plane
1892	VarName=FieldData.ListVarName{ivar};
1893	if ivar==ivar_X %x coordinate
1894	eval(['ProjData.' VarName '=coord_X;'])
1895	elseif ivar==ivar_Y % y coordinate
1896	eval(['ProjData.' VarName '=coord_Y;'])
1897	elseif isempty(ivar_Z) \|\| ivar~=ivar_Z % other variables (except Z coordinate wyhich is not reproduced)
1898	eval(['ProjData.' VarName '=FieldData.' VarName ';'])
1899	end
1900	if isempty(ivar_Z) \|\| ivar~=ivar_Z
1901	ProjData.ListVarName=[ProjData.ListVarName VarName];
1902	ProjData.VarDimName=[ProjData.VarDimName DimCell];
1903	nbvar=nbvar+1;
1904	if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar
1905	ProjData.VarAttribute{nbvar}=FieldData.VarAttribute{ivar};
1906	end
1907	end
1908	end
1909	elseif isequal(ObjectData.ProjMode,'interp_lin')\|\|isequal(ObjectData.ProjMode,'interp_tps')%interpolate data on a regular grid
1910	coord_x_proj=XMin:DX:XMax;
1911	coord_y_proj=YMin:DY:YMax;
1912	coord_z_proj=ZMin:DZ:ZMax;
1913	DimCell={'coord_z','coord_y','coord_x'};
1914	ProjData.ListVarName={'coord_z','coord_y','coord_x'};
1915	ProjData.VarDimName={'coord_z','coord_y','coord_x'};
1916	nbcoord=2;
1917	ProjData.coord_z=[ZMin ZMax];
1918	ProjData.coord_y=[YMin YMax];
1919	ProjData.coord_x=[XMin XMax];
1920	if isempty(ivar_X), ivar_X=0; end;
1921	if isempty(ivar_Y), ivar_Y=0; end;
1922	if isempty(ivar_Z), ivar_Z=0; end;
1923	if isempty(ivar_U), ivar_U=0; end;
1924	if isempty(ivar_V), ivar_V=0; end;
1925	if isempty(ivar_W), ivar_W=0; end;
1926	if isempty(ivar_F), ivar_F=0; end;
1927	if isempty(ivar_FF), ivar_FF=0; end;
1928	if ~isequal(ivar_FF,0)
1929	VarName_FF=FieldData.ListVarName{ivar_FF};
1930	eval(['indsel=find(FieldData.' VarName_FF '==0);'])
1931	coord_X=coord_X(indsel);
1932	coord_Y=coord_Y(indsel);
1933	end
1934	FF=zeros(1,length(coord_y_proj)*length(coord_x_proj));
1935	testFF=0;
1936	[X,Y,Z]=meshgrid(coord_y_proj,coord_z_proj,coord_x_proj);%grid in the new coordinates
1937	for ivar=VarIndex
1938	VarName=FieldData.ListVarName{ivar};
1939	if ~( ivar==ivar_X \|\| ivar==ivar_Y \|\| ivar==ivar_Z \|\| ivar==ivar_F \|\| ivar==ivar_FF \|\| test_anc(ivar)==1)
1940	ivar_new=ivar_new+1;
1941	ProjData.ListVarName=[ProjData.ListVarName {VarName}];
1942	ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
1943	if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar
1944	ProjData.VarAttribute{ivar_new+nbcoord}=FieldData.VarAttribute{ivar};
1945	end
1946	if ~isequal(ivar_FF,0)
1947	eval(['FieldData.' VarName '=FieldData.' VarName '(indsel);'])
1948	end
1949	% linear interpolation
1950	InterpFct=TriScatteredInterp(double(coord_X),double(coord_Y),double(coord_Z),double(FieldData.(VarName)));
1951	ProjData.(VarName)=InterpFct(X,Y,Z);
1952	% eval(['varline=reshape(ProjData.' VarName ',1,length(coord_y_proj)*length(coord_x_proj));'])
1953	% FFlag= isnan(varline); %detect undefined values NaN
1954	% indnan=find(FFlag);
1955	% if~isempty(indnan)
1956	% varline(indnan)=zeros(size(indnan));
1957	% eval(['ProjData.' VarName '=reshape(varline,length(coord_y_proj),length(coord_x_proj));'])
1958	% FF(indnan)=ones(size(indnan));
1959	% testFF=1;
1960	% end
1961	if ivar==ivar_U
1962	ivar_U=ivar_new;
1963	end
1964	if ivar==ivar_V
1965	ivar_V=ivar_new;
1966	end
1967	if ivar==ivar_W
1968	ivar_W=ivar_new;
1969	end
1970	end
1971	end
1972	if testFF
1973	ProjData.FF=reshape(FF,length(coord_y_proj),length(coord_x_proj));
1974	ProjData.ListVarName=[ProjData.ListVarName {'FF'}];
1975	ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
1976	ProjData.VarAttribute{ivar_new+1+nbcoord}.Role='errorflag';
1977	end
1978	end
1979
1980	%% case of input fields defined on a structured grid
1981	else
1982	VarName=FieldData.ListVarName{VarIndex(1)};%get the first variable of the cell to get the input matrix dimensions
1983	eval(['DimValue=size(FieldData.' VarName ');'])%input matrix dimensions
1984	DimValue(DimValue==1)=[];%remove singleton dimensions
1985	NbDim=numel(DimValue);%update number of space dimensions
1986	nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate
1987	if NbDim>=3
1988	if NbDim>3
1989	errormsg='matrices with more than 3 dimensions not handled';
1990	return
1991	else
1992	if numel(find(VarType.coord))==2% the third matrix dimension does not correspond to a space coordinate
1993	nbcolor=DimValue(3);
1994	DimValue(3)=[]; %number of 'color' components updated
1995	NbDim=2;% space dimension set to 2
1996	end
1997	end
1998	end
1999	AYName=FieldData.ListVarName{VarType.coord(NbDim-1)};%name of input x coordinate (name preserved on projection)
2000	AXName=FieldData.ListVarName{VarType.coord(NbDim)};%name of input y coordinate (name preserved on projection)
2001	eval(['AX=FieldData.' AXName ';'])
2002	eval(['AY=FieldData.' AYName ';'])
2003	ListDimName=FieldData.VarDimName{VarIndex(1)};
2004	ProjData.ListVarName=[ProjData.ListVarName {AYName} {AXName}]; %TODO: check if it already exists in Projdata (several cells)
2005	ProjData.VarDimName=[ProjData.VarDimName {AYName} {AXName}];
2006
2007	% for idim=1:length(ListDimName)
2008	% DimName=ListDimName{idim};
2009	% if strcmp(DimName,'rgb')\|\|strcmp(DimName,'nb_coord')\|\|strcmp(DimName,'nb_coord_i')
2010	% nbcolor=DimValue(idim);
2011	% DimValue(idim)=[];
2012	% end
2013	% if isequal(DimName,'nb_coord_j')% NOTE: CASE OF TENSOR NOT TREATED
2014	% DimValue(idim)=[];
2015	% end
2016	% end
2017	Coord_z=[];
2018	Coord_y=[];
2019	Coord_x=[];
2020
2021	for idim=1:NbDim %loop on space dimensions
2022	test_interp(idim)=0;%test for coordiate interpolation (non regular grid), =0 by default
2023	ivar=VarType.coord(idim);% index of the variable corresponding to the current dimension
2024	if ~isequal(ivar,0)% a variable corresponds to the dimension #idim
2025	eval(['Coord{idim}=FieldData.' FieldData.ListVarName{ivar} ';']) ;% coord values for the input field
2026	if numel(Coord{idim})==2 %input array defined on a regular grid
2027	DCoord_min(idim)=(Coord{idim}(2)-Coord{idim}(1))/DimValue(idim);
2028	else
2029	DCoord=diff(Coord{idim});%array of coordinate derivatives for the input field
2030	DCoord_min(idim)=min(DCoord);
2031	DCoord_max=max(DCoord);
2032	% test_direct(idim)=DCoord_max>0;% =1 for increasing values, 0 otherwise
2033	if abs(DCoord_max-DCoord_min(idim))>abs(DCoord_max/1000)
2034	msgbox_uvmat('ERROR',['non monotonic dimension variable # ' num2str(idim) ' in proj_field.m'])
2035	return
2036	end
2037	test_interp(idim)=(DCoord_max-DCoord_min(idim))> 0.0001*abs(DCoord_max);% test grid regularity
2038	end
2039	test_direct(idim)=(DCoord_min(idim)>0);
2040	else % no variable associated with the dimension #idim, the coordinate value is set equal to the matrix index by default
2041	Coord_i_str=['Coord_' num2str(idim)];
2042	DCoord_min(idim)=1;%default
2043	Coord{idim}=[0.5 DimValue(idim)-0.5];
2044	test_direct(idim)=1;
2045	end
2046	end
2047	if DY==0
2048	DY=abs(DCoord_min(NbDim-1));
2049	end
2050	npY=1+round(abs(Coord{NbDim-1}(end)-Coord{NbDim-1}(1))/DY);%nbre of points after interpol
2051	if DX==0
2052	DX=abs(DCoord_min(NbDim));
2053	end
2054	npX=1+round(abs(Coord{NbDim}(end)-Coord{NbDim}(1))/DX);%nbre of points after interpol
2055	for idim=1:NbDim
2056	if test_interp(idim)
2057	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
2058	end
2059	end
2060	Coord_y=linspace(Coord{NbDim-1}(1),Coord{NbDim-1}(end),npY);
2061	test_direct_y=test_direct(NbDim-1);
2062	Coord_x=linspace(Coord{NbDim}(1),Coord{NbDim}(end),npX);
2063	test_direct_x=test_direct(NbDim);
2064	DAX=DCoord_min(NbDim);
2065	DAY=DCoord_min(NbDim-1);
2066	minAX=min(Coord_x);
2067	maxAX=max(Coord_x);
2068	minAY=min(Coord_y);
2069	maxAY=max(Coord_y);
2070	xcorner=[minAX maxAX minAX maxAX]-ObjectData.Coord(1,1);
2071	ycorner=[maxAY maxAY minAY minAY]-ObjectData.Coord(1,2);
2072	xcor_new=xcornercos(Phi)+ycornersin(Phi);%coord new frame
2073	ycor_new=-xcornersin(Phi)+ycornercos(Phi);
2074	if ~testXMax
2075	XMax=max(xcor_new);
2076	end
2077	if ~testXMin
2078	XMin=min(xcor_new);
2079	end
2080	if ~testYMax
2081	YMax=max(ycor_new);
2082	end
2083	if ~testYMin
2084	YMin=min(ycor_new);
2085	end
2086	DXinit=(maxAX-minAX)/(DimValue(NbDim)-1);
2087	DYinit=(maxAY-minAY)/(DimValue(NbDim-1)-1);
2088	if DX==0
2089	DX=DXinit;
2090	end
2091	if DY==0
2092	DY=DYinit;
2093	end
2094	if NbDim==3
2095	DZ=(Coord{1}(end)-Coord{1}(1))/(DimValue(1)-1);
2096	if ~test_direct(1)
2097	DZ=-DZ;
2098	end
2099	Coord_z=linspace(Coord{1}(1),Coord{1}(end),DimValue(1));
2100	test_direct_z=test_direct(1);
2101	end
2102	npX=floor((XMax-XMin)/DX+1);
2103	npY=floor((YMax-YMin)/DY+1);
2104	if test_direct_y
2105	coord_y_proj=linspace(YMin,YMax,npY);%abscissa of the new pixels along the line
2106	else
2107	coord_y_proj=linspace(YMax,YMin,npY);%abscissa of the new pixels along the line
2108	end
2109	if test_direct_x
2110	coord_x_proj=linspace(XMin,XMax,npX);%abscissa of the new pixels along the line
2111	else
2112	coord_x_proj=linspace(XMax,XMin,npX);%abscissa of the new pixels along the line
2113	end
2114
2115	% case with no rotation and interpolation
2116	if isequal(ProjMode,'projection') && isequal(Phi,0) && isequal(Theta,0) && isequal(Psi,0)
2117	if ~testXMin && ~testXMax && ~testYMin && ~testYMax && NbDim==2
2118	ProjData=FieldData;
2119	else
2120	indY=NbDim-1;
2121	if test_direct(indY)
2122	min_indy=ceil((YMin-Coord{indY}(1))/DYinit)+1;
2123	max_indy=floor((YMax-Coord{indY}(1))/DYinit)+1;
2124	Ybound(1)=Coord{indY}(1)+DYinit*(min_indy-1);
2125	Ybound(2)=Coord{indY}(1)+DYinit*(max_indy-1);
2126	else
2127	min_indy=ceil((Coord{indY}(1)-YMax)/DYinit)+1;
2128	max_indy=floor((Coord{indY}(1)-YMin)/DYinit)+1;
2129	Ybound(2)=Coord{indY}(1)-DYinit*(max_indy-1);
2130	Ybound(1)=Coord{indY}(1)-DYinit*(min_indy-1);
2131	end
2132	if test_direct(NbDim)==1
2133	min_indx=ceil((XMin-Coord{NbDim}(1))/DXinit)+1;
2134	max_indx=floor((XMax-Coord{NbDim}(1))/DXinit)+1;
2135	Xbound(1)=Coord{NbDim}(1)+DXinit*(min_indx-1);
2136	Xbound(2)=Coord{NbDim}(1)+DXinit*(max_indx-1);
2137	else
2138	min_indx=ceil((Coord{NbDim}(1)-XMax)/DXinit)+1;
2139	max_indx=floor((Coord{NbDim}(1)-XMin)/DXinit)+1;
2140	Xbound(2)=Coord{NbDim}(1)+DXinit*(max_indx-1);
2141	Xbound(1)=Coord{NbDim}(1)+DXinit*(min_indx-1);
2142	end
2143	if NbDim==3
2144	DimCell(1)=[]; %suppress z variable
2145	DimValue(1)=[];
2146	%structured coordinates
2147	if test_direct(1)
2148	iz=ceil((ObjectData.Coord(1,3)-Coord{1}(1))/DZ)+1;
2149	else
2150	iz=ceil((Coord{1}(1)-ObjectData.Coord(1,3))/DZ)+1;
2151	end
2152	end
2153	min_indy=max(min_indy,1);% deals with margin (bound lower than the first index)
2154	min_indx=max(min_indx,1);
2155	max_indy=min(max_indy,DimValue(1));
2156	max_indx=min(max_indx,DimValue(2));
2157	for ivar=VarIndex% loop on non coordinate variables
2158	VarName=FieldData.ListVarName{ivar};
2159	ProjData.ListVarName=[ProjData.ListVarName VarName];
2160	ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
2161	if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute)>=ivar
2162	ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar};
2163	end
2164	if NbDim==3
2165	eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,min_indy:max_indy,min_indx:max_indx));']);
2166	else
2167	eval(['ProjData.' VarName '=FieldData.' VarName '(min_indy:max_indy,min_indx:max_indx,:);']);
2168	end
2169	end
2170	eval(['ProjData.' AYName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates
2171	eval(['ProjData.' AXName '=[Xbound(1) Xbound(2)];']) %record the new (projected ) x coordinates
2172	end
2173	elseif isfield(FieldData,'A') %TO GENERALISE % case with rotation and/or interpolation
2174	if NbDim==2 %2D case
2175	[X,Y]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates
2176	XIMA=ObjectData.Coord(1,1)+(X)cos(Phi)-Ysin(Phi);%corresponding coordinates in the original image
2177	YIMA=ObjectData.Coord(1,2)+(X)sin(Phi)+Ycos(Phi);
2178	XIMA=(XIMA-minAX)/DXinit+1;% image index along x
2179	YIMA=(-YIMA+maxAY)/DYinit+1;% image index along y
2180	XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line'
2181	YIMA=reshape(round(YIMA),1,npX*npY);
2182	flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image
2183	if isequal(ObjectData.ProjMode,'interp_tps')
2184	npx_interp_tps=ceil(abs(DX/DAX));
2185	npy_interp_tps=ceil(abs(DY/DAY));
2186	Minterp_tps=ones(npy_interp_tps,npx_interp_tps)/(npx_interp_tps*npy_interp_tps);
2187	test_interp_tps=1;
2188	else
2189	test_interp_tps=0;
2190	end
2191	eval(['ProjData.' AYName '=[coord_y_proj(1) coord_y_proj(end)];']) %record the new (projected ) y coordinates
2192	eval(['ProjData.' AXName '=[coord_x_proj(1) coord_x_proj(end)];']) %record the new (projected ) x coordinates
2193	for ivar=VarIndex
2194	VarName=FieldData.ListVarName{ivar};
2195	if test_interp(1) \|\| test_interp(2)%interpolate on a regular grid
2196	eval(['ProjData.' VarName '=interp2(Coord{2},Coord{1},FieldData.' VarName ',Coord_x,Coord_y'');']) %TO TEST
2197	end
2198	%filter the field (image) if option 'interp_tps' is used
2199	if test_interp_tps
2200	Aclass=class(FieldData.A);
2201	ProjData.(VarName)=interp_tps2(Minterp_tps,FieldData.(VarName),'valid');
2202	if ~isequal(Aclass,'double')
2203	ProjData.(VarName)=Aclass(FieldData.(VarName));%revert to integer values
2204	end
2205	end
2206	eval(['vec_A=reshape(FieldData.' VarName ',[],nbcolor);'])%put the original image in line
2207	%ind_in=find(flagin);
2208	ind_out=find(~flagin);
2209	ICOMB=(XIMA-1)*DimValue(1)+YIMA;
2210	ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
2211	vec_B(flagin,1:nbcolor)=vec_A(ICOMB,:);
2212	for icolor=1:nbcolor
2213	vec_B(ind_out,icolor)=zeros(size(ind_out));
2214	end
2215	ProjData.ListVarName=[ProjData.ListVarName VarName];
2216	ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
2217	if isfield(FieldData,'VarAttribute')&&length(FieldData.VarAttribute)>=ivar
2218	ProjData.VarAttribute{length(ProjData.ListVarName)+nbcoord}=FieldData.VarAttribute{ivar};
2219	end
2220	eval(['ProjData.' VarName '=reshape(vec_B,npY,npX,nbcolor);']);
2221	end
2222	ProjData.FF=reshape(~flagin,npY,npX);%false flag A FAIRE: tenir compte d'un flga antï¿œrieur
2223	ProjData.ListVarName=[ProjData.ListVarName 'FF'];
2224	ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
2225	ProjData.VarAttribute{length(ProjData.ListVarName)}.Role='errorflag';
2226	else %3D case
2227	if ~testangle
2228	% unstructured z coordinate
2229	test_sup=(Coord{1}>=ObjectData.Coord(1,3));
2230	iz_sup=find(test_sup);
2231	iz=iz_sup(1);
2232	if iz>=1 & iz<=npz
2233	%ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)];
2234	%ProjData.DimValue=[ProjData.DimValue npY npX];
2235	for ivar=VarIndex
2236	VarName=FieldData.ListVarName{ivar};
2237	ProjData.ListVarName=[ProjData.ListVarName VarName];
2238	ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes
2239	eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,:,:));'])% select the z index iz
2240	%TODO : do a vertical average for a thick plane
2241	if test_interp(2) \|\| test_interp(3)
2242	eval(['ProjData.' VarName '=interp2(Coord{3},Coord{2},ProjData.' VarName ',Coord_x,Coord_y'');'])
2243	end
2244	end
2245	end
2246	else
2247	errormsg='projection of structured coordinates on oblique plane not yet implemented';
2248	%TODO: use interp3
2249	return
2250	end
2251	end
2252	end
2253	end
2254
2255	%% projection of velocity components in the rotated coordinates
2256	if testangle
2257	if isempty(ivar_V)
2258	msgbox_uvmat('ERROR','v velocity component missing in proj_field.m')
2259	return
2260	end
2261	UName=FieldData.ListVarName{ivar_U};
2262	VName=FieldData.ListVarName{ivar_V};
2263	eval(['ProjData.' UName '=cos(Phi)ProjData.' UName '+ sin(Phi)ProjData.' VName ';'])
2264	eval(['ProjData.' VName '=cos(Theta)(-sin(Phi)ProjData.' UName '+ cos(Phi)*ProjData.' VName ');'])
2265	if ~isempty(ivar_W)
2266	WName=FieldData.ListVarName{ivar_W};
2267	eval(['ProjData.' VName '=ProjData.' VName '+ ProjData.' WName '*sin(Theta);'])%
2268	eval(['ProjData.' WName '=NormVec_XProjData.' UName '+ NormVec_YProjData.' VName '+ NormVec_Z* ProjData.' WName ';']);
2269	end
2270	if ~isequal(Psi,0)
2271	eval(['ProjData.' UName '=cos(Psi)* ProjData.' UName '- sin(Psi)*ProjData.' VName ';']);
2272	eval(['ProjData.' VName '=sin(Psi)* ProjData.' UName '+ cos(Psi)*ProjData.' VName ';']);
2273	end
2274	end
2275	end
2276
2277	%------------------------------------------------------------------------
2278	%--- transfer the global attributes
2279	function [ProjData,errormsg]=proj_heading(FieldData,ObjectData)
2280	%------------------------------------------------------------------------
2281	ProjData=[];%default
2282	errormsg='';%default
2283
2284	%% transfer error
2285	if isfield(FieldData,'Txt')
2286	errormsg=FieldData.Txt; %transmit erreur message
2287	return;
2288	end
2289
2290	%% transfer global attributes
2291	if ~isfield(FieldData,'ListGlobalAttribute')
2292	ProjData.ListGlobalAttribute={};
2293	else
2294	ProjData.ListGlobalAttribute=FieldData.ListGlobalAttribute;
2295	end
2296	for iattr=1:length(ProjData.ListGlobalAttribute)
2297	AttrName=ProjData.ListGlobalAttribute{iattr};
2298	if isfield(FieldData,AttrName)
2299	ProjData.(AttrName)=FieldData.(AttrName);
2300	end
2301	end
2302
2303	%% transfer coordinate unit
2304	if isfield(ProjData,'CoordUnit')
2305	ProjData=rmfield(ProjData,'CoordUnit');% do not transfer by default (to avoid x/y=1 for profiles)
2306	end
2307	if isfield(FieldData,'CoordUnit')
2308	if isfield(ObjectData,'CoordUnit') && ~strcmp(FieldData.CoordUnit,ObjectData.CoordUnit)
2309	errormsg=[ObjectData.Type ' in ' ObjectData.CoordUnit ' coordinates, while field in ' FieldData.CoordUnit ];
2310	return
2311	elseif strcmp(ObjectData.Type,'plane')\|\| strcmp(ObjectData.Type,'volume')
2312	ProjData.CoordUnit=FieldData.CoordUnit;
2313	end
2314	end
2315
2316	%% store the properties of the projection object
2317	ListObject={'Name','Type','ProjMode','angle','RangeX','RangeY','RangeZ','DX','DY','DZ','Coord'};
2318	for ilist=1:length(ListObject)
2319	if isfield(ObjectData,ListObject{ilist})
2320	val=ObjectData.(ListObject{ilist});
2321	if ~isempty(val)
2322	ProjData.(['ProjObject' ListObject{ilist}])=val;
2323	ProjData.ListGlobalAttribute=[ProjData.ListGlobalAttribute {['ProjObject' ListObject{ilist}]}];
2324	end
2325	end
2326	end
2327

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

Download in other formats: