source: trunk/src/tps_coeff_field.m @ 643

Last change on this file since 643 was 586, checked in by sommeria, 11 years ago

tps_coeff_field introduced and several bugs corrected

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