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 or 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 has been chosen), |
---|
14 | % =0 otherwise (tps only needed to get spatial derivatives of scattered data) |
---|
15 | % |
---|
16 | % called functions: |
---|
17 | % 'find_field_cells': analyse the input field structure, grouping the variables into 'fields' with common coordinates |
---|
18 | % 'set_subdomains': sort a set of points defined by scattered coordinates in subdomains, as needed for tps interpolation |
---|
19 | % 'tps_coeff': calculate the thin plate spline (tps) coefficients for a single domain. |
---|
20 | |
---|
21 | %======================================================================= |
---|
22 | % Copyright 2008-2019, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France |
---|
23 | % http://www.legi.grenoble-inp.fr |
---|
24 | % Joel.Sommeria - Joel.Sommeria (A) legi.cnrs.fr |
---|
25 | % |
---|
26 | % This file is part of the toolbox UVMAT. |
---|
27 | % |
---|
28 | % UVMAT is free software; you can redistribute it and/or modify |
---|
29 | % it under the terms of the GNU General Public License as published |
---|
30 | % by the Free Software Foundation; either version 2 of the license, |
---|
31 | % or (at your option) any later version. |
---|
32 | % |
---|
33 | % UVMAT is distributed in the hope that it will be useful, |
---|
34 | % but WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
35 | % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
---|
36 | % GNU General Public License (see LICENSE.txt) for more details. |
---|
37 | %======================================================================= |
---|
38 | |
---|
39 | function [DataOut,errormsg]=tps_coeff_field(DataIn,checkall) |
---|
40 | DataOut=DataIn;%default |
---|
41 | SubDomainNbPoint=1000; %default, estimated nbre of data source points in a subdomain used for tps |
---|
42 | if isfield(DataIn,'SubDomain') |
---|
43 | SubDomainNbPoint=DataIn.SubDomain;%old convention |
---|
44 | end |
---|
45 | if isfield(DataIn,'SubDomainNbPoint') |
---|
46 | SubDomainNbPoint=DataIn.SubDomainNbPoint;% |
---|
47 | end |
---|
48 | [CellInfo,NbDimArray,errormsg]=find_field_cells(DataIn); |
---|
49 | if ~isempty(errormsg) |
---|
50 | errormsg=['tps_coeff_field/find_field_cells/' errormsg]; |
---|
51 | return |
---|
52 | end |
---|
53 | nbtps=0;% indicate the number of tps coordinate sets in the field structure (in general =1) |
---|
54 | |
---|
55 | for icell=1:numel(CellInfo); |
---|
56 | if NbDimArray(icell)>=2 && strcmp(CellInfo{icell}.CoordType,'scattered') %if the coordinates are scattered |
---|
57 | NbCoord=NbDimArray(icell);% dimension of space |
---|
58 | nbtps=nbtps+1;% indicate the number of tps coordinate sets in the field structure (in general =1) |
---|
59 | X=DataIn.(DataIn.ListVarName{CellInfo{icell}.CoordIndex(end)});% value of x coordinate |
---|
60 | Y=DataIn.(DataIn.ListVarName{CellInfo{icell}.CoordIndex(end-1)});% value of y coordinate |
---|
61 | check_interp_tps=false(numel(DataIn.ListVarName),1); |
---|
62 | Index_interp=[];% indices of variables to interpolate |
---|
63 | if isfield(CellInfo{icell},'VarIndex_scalar')%interpolate scalar |
---|
64 | Index_interp=[Index_interp CellInfo{icell}.VarIndex_scalar]; |
---|
65 | end |
---|
66 | if isfield(CellInfo{icell},'VarIndex_vector_x')%interpolate vector x component |
---|
67 | Index_interp=[Index_interp CellInfo{icell}.VarIndex_vector_x]; |
---|
68 | end |
---|
69 | if isfield(CellInfo{icell},'VarIndex_vector_y')%interpolate vector y component |
---|
70 | Index_interp=[Index_interp CellInfo{icell}.VarIndex_vector_y]; |
---|
71 | end |
---|
72 | for iselect=1:numel(Index_interp) |
---|
73 | Attr=DataIn.VarAttribute{Index_interp(iselect)}; |
---|
74 | if ~isfield(Attr,'VarIndex_tps')&& (checkall || (isfield(Attr,'ProjModeRequest')&&strcmp(Attr.ProjModeRequest,'interp_tps'))) |
---|
75 | check_interp_tps(Index_interp(iselect))=1; |
---|
76 | end |
---|
77 | end |
---|
78 | ListVarInterp=DataIn.ListVarName(check_interp_tps); |
---|
79 | VarIndexInterp=find(check_interp_tps); |
---|
80 | if ~isempty(ListVarInterp) |
---|
81 | % exclude data points marked 'false' for interpolation |
---|
82 | if isfield(CellInfo{icell},'VarIndex_errorflag') |
---|
83 | FF=DataIn.(DataIn.ListVarName{CellInfo{icell}.VarIndex_errorflag});% error flag |
---|
84 | X=X(FF==0); |
---|
85 | Y=Y(FF==0); |
---|
86 | for ilist=1:numel(ListVarInterp) |
---|
87 | DataIn.(ListVarInterp{ilist})=DataIn.(ListVarInterp{ilist})(FF==0); |
---|
88 | end |
---|
89 | end |
---|
90 | term=''; |
---|
91 | if nbtps>1 |
---|
92 | term=['_' num2str(nbtps-1)]; |
---|
93 | end |
---|
94 | ListNewVar=cell(1,numel(ListVarInterp)+3); |
---|
95 | ListNewVar(1:3)={['SubRange' term],['NbCentre' term],['Coord_tps' term]}; |
---|
96 | for ilist=1:numel(ListVarInterp) |
---|
97 | ListNewVar{ilist+3}=[ListVarInterp{ilist} '_tps' term]; |
---|
98 | end |
---|
99 | nbvar=numel(DataIn.ListVarName); |
---|
100 | DataOut.ListVarName=[DataIn.ListVarName ListNewVar]; |
---|
101 | DataOut.VarDimName=[DataIn.VarDimName {{'nb_coord','nb_bounds',['nb_subdomain' term]}} {['nb_subdomain' term]} ... |
---|
102 | {{['nb_tps' term],'nb_coord',['nb_subdomain' term]}}]; |
---|
103 | DataOut.VarAttribute{nbvar+3}.Role='coord_tps'; |
---|
104 | [SubRange,NbCentre,IndSelSubDomain] =set_subdomains([X Y],SubDomainNbPoint);% create subdomains for tps |
---|
105 | for isub=1:size(SubRange,3) |
---|
106 | ind_sel=IndSelSubDomain(1:NbCentre(isub),isub);% array indices selected for the subdomain |
---|
107 | DataOut.(['Coord_tps' term])(1:NbCentre(isub),1:2,isub)=[X(ind_sel) Y(ind_sel)]; |
---|
108 | 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) |
---|
109 | end |
---|
110 | for ivar=1:numel(ListVarInterp) |
---|
111 | DataOut.VarDimName{nbvar+3+ivar}={['nb_tps' term],['nb_subdomain' term]}; |
---|
112 | DataOut.VarAttribute{nbvar+3+ivar}=DataIn.VarAttribute{CellInfo{icell}.VarIndex_vector_x};%reproduce attributes of velocity |
---|
113 | if ~isfield(DataIn.VarAttribute{VarIndexInterp(ivar)},'Role') |
---|
114 | DataOut.VarAttribute{nbvar+3+ivar}.Role='scalar_tps'; |
---|
115 | else |
---|
116 | DataOut.VarAttribute{nbvar+3+ivar}.Role=[DataIn.VarAttribute{VarIndexInterp(ivar)}.Role '_tps']; |
---|
117 | end |
---|
118 | DataOut.VarAttribute{VarIndexInterp(ivar)}.VarIndex_tps=nbvar+3+ivar;% indicate the tps correspondance in the source data |
---|
119 | end |
---|
120 | if isfield(DataOut,'ListDimName')%cleaning' |
---|
121 | DataOut=rmfield(DataOut,'ListDimName'); |
---|
122 | end |
---|
123 | if isfield(DataOut,'DimValue')%cleaning |
---|
124 | DataOut=rmfield(DataOut,'DimValue'); |
---|
125 | end |
---|
126 | DataOut.(['SubRange' term])=SubRange; |
---|
127 | DataOut.(['NbCentre' term])=NbCentre; |
---|
128 | for ilist=1:numel(VarIndexInterp) |
---|
129 | for isub=1:size(SubRange,3) |
---|
130 | ind_sel=IndSelSubDomain(1:NbCentre(isub),isub);% array indices selected for the subdomain |
---|
131 | [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 |
---|
132 | end |
---|
133 | DataOut.(ListNewVar{ilist+3})=Var_tps; |
---|
134 | end |
---|
135 | end |
---|
136 | end |
---|
137 | end |
---|