source: trunk/src/proj_field.m @ 402

Last change on this file since 402 was 402, checked in by sommeria, 12 years ago

bugs corrected and improved procedure for object projection

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