1## Copyright (C) 2013 Nir Krakauer 2## 3## This program is free software; you can redistribute it and/or modify 4## it under the terms of the GNU General Public License as published by 5## the Free Software Foundation; either version 2 of the License, or 6## (at your option) any later version. 7## 8## This program is distributed in the hope that it will be useful, 9## but WITHOUT ANY WARRANTY; without even the implied warranty of 10## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11## GNU General Public License for more details. 12## 13## You should have received a copy of the GNU General Public License 14## along with this program; If not, see <http://www.gnu.org/licenses/>. 15 16## -*- texinfo -*- 17## @deftypefn{Function File}{[@var{yi}] =} tps_val(@var{x}, @var{coefs}, @var{xi}, @var{vectorize}=true) 18## 19## Evaluates a thin plate spline at given points @* 20## @var{xi} 21## 22## @var{coefs} should be the vector of fitted coefficients returned from @code{tpaps(x, y, [p])} 23## 24## @var{x} should be @var{n} by @var{d} in size, where @var{n} is the number of points and @var{d} the number of dimensions; @var{coefs} should be @var{n} + @var{d} + 1 by 1; @var{xi} should be @var{k} by @var{d} 25## 26## The logical argument @var{vectorize} controls whether an @var{k} by @var{n} by @var{d} intermediate array is formed to speed up computation (the default) or whether looping is used to economize on memory 27## 28## The returned @var{yi} will be @var{k} by 1 29## 30## See the documentation to @code{tpaps} for more information 31## 32## @end deftypefn 33## @seealso{tpaps, tps_val_der} 34 35## Author: Nir Krakauer <nkrakauer@ccny.cuny.edu> 36 37function [yi]=tps_val(x,coefs,xi,vectorize=true) 38 39 40 [n d] = size(x); #d: number of dimensions; n: number of points 41 k = size(xi, 1); #number of points for which to find the spline function value 42 43 #form of the Green's function for solutions 44 G = @(r) merge(r == 0, 0, r .^ 2 .* log(r)); 45 46 a = coefs(1:n); 47 b = coefs((n+1):end); 48 49 50 yi = [ones(k, 1) xi] * b; 51 if vectorize 52 if d == 1 53 yi = yi + G(abs(x' - xi)) * a; 54 else 55 yi = yi + G(sqrt(sumsq((reshape(x, 1, n, d) - reshape(xi, k, 1, d)), 3))) * a; 56 endif 57 else 58 dist = @(x1, x2) norm(x2 - x1, 2, "rows"); #Euclidian distance between points in d-dimensional space 59 warn_state = warning ("query", "Octave:broadcast").state; 60 warning ("off", "Octave:broadcast"); #turn off warning message for automatic broadcasting when dist is called 61 unwind_protect 62 if k > n ##choose from either of two ways of computing the values of the thin plate spline at xi 63 for i = 1:n 64 yi = yi + a(i)*G(dist(x(i, :), xi)); 65 endfor 66 else 67 for i = 1:k 68 yi(i) = yi(i) + dot(a, G(dist(x, xi(i, :)))); 69 endfor 70 endif 71 unwind_protect_cleanup 72 warning (warn_state, "Octave:broadcast"); 73 end_unwind_protect 74 endif 75 76 77endfunction 78 79%!shared x,y,c,xi 80%! x = ([1:10 10.5 11.3])'; y = sin(x); 81%! c = tpaps(x,y,1); 82%!assert (tpaps(x,y,1,x), tps_val(x,c,x), 100*eps); 83%! x = ((1 ./ (1:100))' - 0.5) * ([0.2 0.6]); 84%! y = x(:, 1) .^ 2 + x(:, 2) .^ 2; 85%! c = tpaps(x,y,1); 86%!assert (tpaps(x,y,1,x), tps_val(x,c,x), 100*eps); 87%!assert (tps_val(x,c,x,true), tps_val(x,c,x,false), 100*eps); 88