[651] | 1 | %'tps_coeff_field': calculate the thin plate spline (tps) coefficients within subdomains for a field structure |
---|
[586] | 2 | %--------------------------------------------------------------------- |
---|
| 3 | % DataOut=tps_coeff_field(DataIn,checkall) |
---|
| 4 | % |
---|
| 5 | % OUTPUT: |
---|
[651] | 6 | % DataOut: output field structure, reproducing the input field structure DataIn and adding the fields: |
---|
| 7 | % .Coord_tps |
---|
[822] | 8 | % .[VarName '_tps'] for each eligible input variable VarName (scalar or vector components) |
---|
[651] | 9 | % errormsg: error message, = '' by default |
---|
[586] | 10 | % |
---|
| 11 | % INPUT: |
---|
| 12 | % DataIn: intput field structure |
---|
[822] | 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) |
---|
[586] | 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 | |
---|
[809] | 21 | %======================================================================= |
---|
[1061] | 22 | % Copyright 2008-2019, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France |
---|
[809] | 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 | |
---|
[586] | 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); |
---|
[822] | 56 | if NbDimArray(icell)>=2 && strcmp(CellInfo{icell}.CoordType,'scattered') %if the coordinates are scattered |
---|
[586] | 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 |
---|
[822] | 60 | Y=DataIn.(DataIn.ListVarName{CellInfo{icell}.CoordIndex(end-1)});% value of y coordinate |
---|
[674] | 61 | check_interp_tps=false(numel(DataIn.ListVarName),1); |
---|
[822] | 62 | Index_interp=[];% indices of variables to interpolate |
---|
| 63 | if isfield(CellInfo{icell},'VarIndex_scalar')%interpolate scalar |
---|
[586] | 64 | Index_interp=[Index_interp CellInfo{icell}.VarIndex_scalar]; |
---|
| 65 | end |
---|
[822] | 66 | if isfield(CellInfo{icell},'VarIndex_vector_x')%interpolate vector x component |
---|
[586] | 67 | Index_interp=[Index_interp CellInfo{icell}.VarIndex_vector_x]; |
---|
| 68 | end |
---|
[822] | 69 | if isfield(CellInfo{icell},'VarIndex_vector_y')%interpolate vector y component |
---|
[586] | 70 | Index_interp=[Index_interp CellInfo{icell}.VarIndex_vector_y]; |
---|
| 71 | end |
---|
[674] | 72 | for iselect=1:numel(Index_interp) |
---|
| 73 | Attr=DataIn.VarAttribute{Index_interp(iselect)}; |
---|
[586] | 74 | if ~isfield(Attr,'VarIndex_tps')&& (checkall || (isfield(Attr,'ProjModeRequest')&&strcmp(Attr.ProjModeRequest,'interp_tps'))) |
---|
[674] | 75 | check_interp_tps(Index_interp(iselect))=1; |
---|
[586] | 76 | end |
---|
| 77 | end |
---|
[674] | 78 | ListVarInterp=DataIn.ListVarName(check_interp_tps); |
---|
| 79 | VarIndexInterp=find(check_interp_tps); |
---|
| 80 | if ~isempty(ListVarInterp) |
---|
[586] | 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); |
---|
[674] | 86 | for ilist=1:numel(ListVarInterp) |
---|
[586] | 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 |
---|
[674] | 94 | ListNewVar=cell(1,numel(ListVarInterp)+3); |
---|
[651] | 95 | ListNewVar(1:3)={['SubRange' term],['NbCentre' term],['Coord_tps' term]}; |
---|
[674] | 96 | for ilist=1:numel(ListVarInterp) |
---|
[586] | 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'; |
---|
[651] | 104 | [SubRange,NbCentre,IndSelSubDomain] =set_subdomains([X Y],SubDomainNbPoint);% create subdomains for tps |
---|
[586] | 105 | for isub=1:size(SubRange,3) |
---|
[651] | 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) |
---|
[674] | 109 | end |
---|
| 110 | for ivar=1:numel(ListVarInterp) |
---|
[586] | 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; |
---|
[651] | 127 | DataOut.(['NbCentre' term])=NbCentre; |
---|
[586] | 128 | for ilist=1:numel(VarIndexInterp) |
---|
| 129 | for isub=1:size(SubRange,3) |
---|
[651] | 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 |
---|
[586] | 132 | end |
---|
| 133 | DataOut.(ListNewVar{ilist+3})=Var_tps; |
---|
| 134 | end |
---|
| 135 | end |
---|
| 136 | end |
---|
[674] | 137 | end |
---|