Mercurial > forge

--- a/main/splines/INDEX	Fri Apr 05 11:39:35 2013 +0000
+++ b/main/splines/INDEX	Fri Apr 05 19:24:09 2013 +0000
@@ -8,4 +8,5 @@
  fnder
  fnval
  catmullrom
-
+ tpaps
+ tps_val
--- a/main/splines/NEWS	Fri Apr 05 11:39:35 2013 +0000
+++ b/main/splines/NEWS	Fri Apr 05 19:24:09 2013 +0000
@@ -1,3 +1,10 @@
+Summary of important user-visible changes for splines 1.2.0:
+-------------------------------------------------------------------
+
+ ** The following functions are new:
+
+      tpaps      tps_val
+
 Summary of important user-visible changes for splines 1.1.2:
 -------------------------------------------------------------------
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/splines/inst/tpaps.m	Fri Apr 05 19:24:09 2013 +0000
@@ -0,0 +1,158 @@
+## Copyright (C) 2013 Nir Krakauer
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn{Function File}{[@var{yi} @var{p}] =} tpaps(@var{x}, @var{y}, @var{p}, @var{xi})
+## @deftypefnx{Function File}{[@var{coefs} @var{p}] =} tpaps(@var{x}, @var{y}, @var{p}, [])
+##
+## Thin plate smoothing of scattered values in multi-D @*
+## approximately interpolate [@var{x},@var{y}] at @var{xi}
+##
+## The chosen thin plate spline minimizes the sum of squared deviations from the given points plus a penalty term proportional to the curvature of the spline function
+##
+## @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{y} and @var{w} should be @var{n} by 1; @var{xi} should be @var{k} by @var{d}; the points in @var{x} should be distinct
+##
+## @table @asis
+## @item @var{p}=0
+##       maximum smoothing: flat surface
+## @item @var{p}=1
+##       no smoothing: interpolation
+## @item @var{p}<0 or not given
+##       an intermediate amount of smoothing is chosen (such that the smoothing term and the interpolation term are of the same magnitude)
+## @end table
+##
+## If @var{xi} is not specified, returns a vector @var{coefs} of the @var{n} + @var{d} + 1 fitted thin plate spline coefficients.
+## Given @var{coefs}, the value of the thin-plate spline at any @var{xi} can be determined with @code{tps_val}
+##
+## Note: Computes the pseudoinverse of an @var{n} by @var{n} matrix, so not recommended for very large @var{n}
+##
+## Example usages:
+## @example
+## x = ([1:10 10.5 11.3])'; y = sin(x); xi = (0:0.1:12)';
+## yi = tpaps(x, y, 0.5, xi);
+## plot(x, y, xi, yi)
+## @end example
+##
+## @example
+## x = rand(100, 2)*2 - 1;
+## y = x(:, 1) .^ 2 + x(:, 2) .^ 2;
+## scatter(x(:, 1), x(:, 2), 10, y, "filled")
+## [x1 y1] = meshgrid((-1:0.2:1)', (-1:0.2:1)');
+## xi = [x1(:) y1(:)];
+## yi = tpaps(x, y, 1, xi);
+## contourf(x1, y1, reshape(yi, 11, 11))
+## @end example
+##
+## Reference:
+## David Eberly (2011), Thin-Plate Splines, www.geometrictools.com/Documentation/ThinPlateSplines.pdf
+## Bouhamidi, A. (2005) Weighted thin plate splines, Analysis and Applications, 3: 297-324
+##
+## @end deftypefn
+## @seealso{csaps, tps_val}
+
+## Author: Nir Krakauer <nkrakauer@ccny.cuny.edu>
+
+function [ret, p]=tpaps(x,y,p,xi)
+
+
+  if(nargin < 4)
+    xi = [];
+    if(nargin < 3) || p < 0
+      p = [];
+    endif
+  endif
+
+
+  [n d] = size(x); #d: number of dimensions; n: number of points [y should be n*1]
+
+  dist = @(x1, x2) norm(x2 - x1, 2, "rows"); #Euclidian distance between points in d-dimensional space
+
+  #form of the Green's function for solutions
+  G = @(r) merge(r == 0, 0, r .^ 2 .* log(r));
+
+  N = [ones(n, 1) x];
+
+  if p == 0 #infinite regularization (no curvature allowed), solution is a regression plane
+
+    b = N \ y;
+
+    a = zeros(n, 1);
+
+  else
+
+    #coefficient matrices
+    #need pairwise distances between points
+    M = zeros(n);
+     warn_state = warning ("query", "Octave:broadcast").state;
+     warning ("off", "Octave:broadcast"); #turn off warning message for automatic broadcasting when dist is called
+     unwind_protect
+    for i = 1:(n-1) #M is symmetric, so only need to compute half
+      M(i, (i+1):n) = G(dist(x(i, :), x((i+1):n, :)));
+    endfor
+     unwind_protect_cleanup
+     warning (warn_state, "Octave:broadcast");
+     end_unwind_protect
+    M = M + M';
+    if isempty(p) #choose an intermediate value for the regularization parameter
+      lambda = mean(spdiags(M, 1:min(n, 3))(:));
+      p = 1 / (lambda + 1);
+    else #use the given value
+      lambda = (1 - p) / p;
+    endif
+
+    M = M + lambda*eye(n); #add regularization term
+
+    M_inv = pinv(M);
+
+    b = pinv(N' * M_inv * N) * N' * M_inv * y;
+
+    a = M_inv * (y - N*b);
+
+  endif
+
+
+  if isempty(xi) #return the coefficients
+      ret = [a' b']';
+  else #return the thin plate spline values at xi
+      k = size(xi, 1);
+      ret = [ones(k, 1) xi] * b;
+      if p ~= 0
+         warn_state = warning ("query", "Octave:broadcast").state;
+         warning ("off", "Octave:broadcast"); #turn off warning message for automatic broadcasting when dist is called
+         unwind_protect
+        if k > n ##choose from either of two ways of computing the values of the thin plate spline at xi
+          for i = 1:n
+            ret = ret + a(i)*G(dist(x(i, :), xi));
+          endfor
+        else
+          for i = 1:k
+            ret(i) = ret(i) + dot(a, G(dist(x, xi(i, :))));
+          endfor
+        endif
+         unwind_protect_cleanup
+         warning (warn_state, "Octave:broadcast");
+         end_unwind_protect
+      endif
+  endif
+
+
+endfunction
+
+%!shared x,y
+%! x = ([1:10 10.5 11.3])'; y = sin(x);
+%!assert (tpaps(x,y,1,x), y, 1E-13);
+%! x = rand(100, 2)*2 - 1;
+%! y = x(:, 1) .^ 2 + x(:, 2) .^ 2;
+%!assert (tpaps(x,y,1,x), y, 1E-11);
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/splines/inst/tps_val.m	Fri Apr 05 19:24:09 2013 +0000
@@ -0,0 +1,77 @@
+## Copyright (C) 2013 Nir Krakauer
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn{Function File}{[@var{yi}] =} tps_val(@var{x}, @var{coefs}, @var{xi})
+##
+## Evaluates a thin plate spline at given points @*
+## @var{xi}
+##
+## @var{coefs} should be the vector of fitted coefficients returned from @code{tpaps(x, y, [p])}
+##
+## @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}
+##
+## See the documentation to @code{tpaps} for more information
+##
+## @end deftypefn
+## @seealso{tpaps}
+
+## Author: Nir Krakauer <nkrakauer@ccny.cuny.edu>
+
+function [yi]=tps_val(x,coefs,xi)
+
+
+  [n d] = size(x); #d: number of dimensions; n: number of points [y should be n*1]
+  k = size(xi, 1); #number of points for which to find the spline function value
+
+  dist = @(x1, x2) norm(x2 - x1, 2, "rows"); #Euclidian distance between points in d-dimensional space
+
+  #form of the Green's function for solutions
+  G = @(r) merge(r == 0, 0, r .^ 2 .* log(r));
+
+ a = coefs(1:n);
+ b = coefs((n+1):end);
+
+
+ yi = [ones(k, 1) xi] * b;
+
+   warn_state = warning ("query", "Octave:broadcast").state;
+   warning ("off", "Octave:broadcast"); #turn off warning message for automatic broadcasting when dist is called
+   unwind_protect
+  if k > n ##choose from either of two ways of computing the values of the thin plate spline at xi
+    for i = 1:n
+      yi = yi + a(i)*G(dist(x(i, :), xi));
+    endfor
+  else
+    for i = 1:k
+      yi(i) = yi(i) + dot(a, G(dist(x, xi(i, :))));
+    endfor
+  endif
+   unwind_protect_cleanup
+   warning (warn_state, "Octave:broadcast");
+   end_unwind_protect
+
+
+
+endfunction
+
+%!shared x,y,c
+%! x = ([1:10 10.5 11.3])'; y = sin(x);
+%! c = tpaps(x,y,1);
+%!assert (tpaps(x,y,1,x), tps_val(x,c,x));
+%! x = rand(100, 2)*2 - 1;
+%! y = x(:, 1) .^ 2 + x(:, 2) .^ 2;
+%! c = tpaps(x,y,1);
+%!assert (tpaps(x,y,1,x), tps_val(x,c,x));