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