Context Navigation

source: trunk/src/proj_field.m @ 1054

Last change on this file since 1054 was 1054, checked in by sommeria, 6 years ago
bed_scan updated
File size: 120.2 KB

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

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

Download in other formats: