Context Navigation

source: trunk/src/proj_field.m @ 889

Last change on this file since 889 was 888, checked in by sommeria, 9 years ago
'proj_field
File size: 111.7 KB

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

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

Download in other formats: