source: trunk/src/tps_coeff_field.m @ 1137

Last change on this file since 1137 was 1137, checked in by sommeria, 8 months ago

tps improved by a cosine window in each subdomain. Default size of subdomains reduced for faster computation in the civ interface

File size: 7.5 KB
RevLine 
[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]40function [DataOut,errormsg]=tps_coeff_field(DataIn,checkall,Smoothing)     
[586]41DataOut=DataIn;%default
[1122]42if ~exist('Smoothing','var')
43Smoothing=0;
44end
[1137]45SubDomainNbPoint=300; %default, estimated nbre of data source points in a subdomain used for tps
[586]46if isfield(DataIn,'SubDomain')
47    SubDomainNbPoint=DataIn.SubDomain;%old convention
48end
49if isfield(DataIn,'SubDomainNbPoint')
50    SubDomainNbPoint=DataIn.SubDomainNbPoint;%
51end
52[CellInfo,NbDimArray,errormsg]=find_field_cells(DataIn);
53if ~isempty(errormsg)
54    errormsg=['tps_coeff_field/find_field_cells/' errormsg];
55    return
56end
57nbtps=0;% indicate the number of tps coordinate sets in the field structure (in general =1)
58
[1122]59for 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]141end
Note: See TracBrowser for help on using the repository browser.