source: trunk/src/tps_coeff_field.m @ 690

Last change on this file since 690 was 674, checked in by sommeria, 11 years ago

various bugs repaired, in particula timing

File size: 7.0 KB
Line 
1%'tps_coeff_field': calculate the thin plate spline (tps) coefficients within subdomains for a field structure
2%---------------------------------------------------------------------
3% DataOut=tps_coeff_field(DataIn,checkall)
4%
5% OUTPUT:
6% DataOut: output field structure, reproducing the input field structure DataIn and adding the fields:
7%         .Coord_tps
8%         .[VarName '_tps'] for each eligible input variable VarName (scalar ofr vector components)
9% errormsg: error message, = '' by default
10%
11% INPUT:
12% DataIn: intput field structure
13% checkall:=1 if tps is needed for all fields (a projection mode interp_tps is needed), =0 otherwise
14%
15% called functions:
16% 'find_field_cells': analyse the input field structure, grouping the variables  into 'fields' with common coordinates
17% 'set_subdomains': sort a set of points defined by scattered coordinates in subdomains, as needed for tps interpolation
18% 'tps_coeff': calculate the thin plate spline (tps) coefficients for a single domain.
19
20function [DataOut,errormsg]=tps_coeff_field(DataIn,checkall)     
21DataOut=DataIn;%default
22SubDomainNbPoint=1000; %default, estimated nbre of data source points in a subdomain used for tps
23if isfield(DataIn,'SubDomain')
24    SubDomainNbPoint=DataIn.SubDomain;%old convention
25end
26if isfield(DataIn,'SubDomainNbPoint')
27    SubDomainNbPoint=DataIn.SubDomainNbPoint;%
28end
29[CellInfo,NbDimArray,errormsg]=find_field_cells(DataIn);
30if ~isempty(errormsg)
31    errormsg=['tps_coeff_field/find_field_cells/' errormsg];
32    return
33end
34nbtps=0;% indicate the number of tps coordinate sets in the field structure (in general =1)
35
36for icell=1:numel(CellInfo);
37    if NbDimArray(icell)>=2 && strcmp(CellInfo{icell}.CoordType,'scattered')% %if the coordinates are scattered
38        NbCoord=NbDimArray(icell);% dimension of space
39        nbtps=nbtps+1;% indicate the number of tps coordinate sets in the field structure (in general =1)
40        X=DataIn.(DataIn.ListVarName{CellInfo{icell}.CoordIndex(end)});% value of x coordinate
41        Y=DataIn.(DataIn.ListVarName{CellInfo{icell}.CoordIndex(end-1)});% value of x coordinate
42        check_interp_tps=false(numel(DataIn.ListVarName),1);
43        %for ivar=1:numel(CellInfo{icell}.VarIndex)
44        Index_interp=[];
45        if isfield(CellInfo{icell},'VarIndex_scalar')
46            Index_interp=[Index_interp CellInfo{icell}.VarIndex_scalar];
47        end
48        if isfield(CellInfo{icell},'VarIndex_vector_x')
49            Index_interp=[Index_interp CellInfo{icell}.VarIndex_vector_x];
50        end
51        if isfield(CellInfo{icell},'VarIndex_vector_y')
52            Index_interp=[Index_interp CellInfo{icell}.VarIndex_vector_y];
53        end
54        for iselect=1:numel(Index_interp)
55            Attr=DataIn.VarAttribute{Index_interp(iselect)};
56            if ~isfield(Attr,'VarIndex_tps')&& (checkall || (isfield(Attr,'ProjModeRequest')&&strcmp(Attr.ProjModeRequest,'interp_tps')))
57                check_interp_tps(Index_interp(iselect))=1;
58            end
59        end
60       
61        %VarIndexInterp=CellInfo{icell}.VarIndex(check_interp_tps);% indices of variables to interpolate through tps
62        %         if ~isempty(VarIndexInterp)
63        ListVarInterp=DataIn.ListVarName(check_interp_tps);
64        VarIndexInterp=find(check_interp_tps);
65        if ~isempty(ListVarInterp)
66            % exclude data points marked 'false' for interpolation
67            if isfield(CellInfo{icell},'VarIndex_errorflag')
68                FF=DataIn.(DataIn.ListVarName{CellInfo{icell}.VarIndex_errorflag});% error flag
69                X=X(FF==0);
70                Y=Y(FF==0);
71                for ilist=1:numel(ListVarInterp)
72                    DataIn.(ListVarInterp{ilist})=DataIn.(ListVarInterp{ilist})(FF==0);
73                end
74            end
75            term='';
76            if nbtps>1
77                term=['_' num2str(nbtps-1)];
78            end
79            ListNewVar=cell(1,numel(ListVarInterp)+3);
80            ListNewVar(1:3)={['SubRange' term],['NbCentre' term],['Coord_tps' term]};
81            for ilist=1:numel(ListVarInterp)
82                ListNewVar{ilist+3}=[ListVarInterp{ilist} '_tps' term];
83            end
84            nbvar=numel(DataIn.ListVarName);
85            DataOut.ListVarName=[DataIn.ListVarName ListNewVar];
86            DataOut.VarDimName=[DataIn.VarDimName {{'nb_coord','nb_bounds',['nb_subdomain' term]}} {['nb_subdomain' term]} ...
87                {{['nb_tps' term],'nb_coord',['nb_subdomain' term]}}];
88            DataOut.VarAttribute{nbvar+3}.Role='coord_tps';
89            [SubRange,NbCentre,IndSelSubDomain] =set_subdomains([X Y],SubDomainNbPoint);% create subdomains for tps
90            for isub=1:size(SubRange,3)
91                ind_sel=IndSelSubDomain(1:NbCentre(isub),isub);% array indices selected for the subdomain
92                DataOut.(['Coord_tps' term])(1:NbCentre(isub),1:2,isub)=[X(ind_sel) Y(ind_sel)];
93                DataOut.(['Coord_tps' term])(NbCentre(isub)+1:NbCentre(isub)+3,1:2,isub)=0;%matrix of zeros to complement the matrix Coord_tps (conveninent for file storage)
94                %fill=zeros(NbCoord+1,NbCoord,size(SubRange,3)); %matrix of zeros to complement the matrix Data.Civ1_Coord_tps (conveninent for file storage)
95                %                 fill=zeros(NbCoord+1,NbCoord+1,NbCoord); %matrix of zeros to complement the matrix Data.Civ1_Coord_tps (conveninent for file storage)
96                %                 Coord_tps=cat(1,Coord_tps,fill);
97            end
98            for ivar=1:numel(ListVarInterp)
99                DataOut.VarDimName{nbvar+3+ivar}={['nb_tps' term],['nb_subdomain' term]};
100                DataOut.VarAttribute{nbvar+3+ivar}=DataIn.VarAttribute{CellInfo{icell}.VarIndex_vector_x};%reproduce attributes of velocity
101                if ~isfield(DataIn.VarAttribute{VarIndexInterp(ivar)},'Role')
102                    DataOut.VarAttribute{nbvar+3+ivar}.Role='scalar_tps';
103                else
104                    DataOut.VarAttribute{nbvar+3+ivar}.Role=[DataIn.VarAttribute{VarIndexInterp(ivar)}.Role '_tps'];
105                end
106                DataOut.VarAttribute{VarIndexInterp(ivar)}.VarIndex_tps=nbvar+3+ivar;% indicate the tps correspondance in the source data
107            end
108            if isfield(DataOut,'ListDimName')%cleaning'
109                DataOut=rmfield(DataOut,'ListDimName');
110            end
111            if isfield(DataOut,'DimValue')%cleaning
112                DataOut=rmfield(DataOut,'DimValue');
113            end
114            DataOut.(['SubRange' term])=SubRange;
115            DataOut.(['NbCentre' term])=NbCentre;
116            %             DataOut.(['Coord_tps' term])=Coord_tps;
117            for ilist=1:numel(VarIndexInterp)
118                %                 Var_tps=zeros(size(IndSelSubDomain)+[NbCoord+1 0]);%default spline
119                for isub=1:size(SubRange,3)
120                    ind_sel=IndSelSubDomain(1:NbCentre(isub),isub);% array indices selected for the subdomain
121                    [tild,Var_tps(1:NbCentre(isub)+NbCoord+1,isub)]=tps_coeff([X(ind_sel) Y(ind_sel)],DataIn.(ListVarInterp{ilist})(ind_sel),0);%calculate the tps coeff in the subdomain
122                end
123                DataOut.(ListNewVar{ilist+3})=Var_tps;
124            end
125        end
126    end
127end
Note: See TracBrowser for help on using the repository browser.