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