octave_packages/m/general/bicubic.m

   1 ## Copyright (C) 2005-2012 Hoxide Ma
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn {Function File} {@var{zi} =} bicubic (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi}, @var{extrapval})
  21 ##
  22 ## Return a matrix @var{zi} corresponding to the bicubic
  23 ## interpolations at @var{xi} and @var{yi} of the data supplied
  24 ## as @var{x}, @var{y} and @var{z}.  Points outside the grid are set
  25 ## to @var{extrapval}.
  26 ##
  27 ## See @url{http://wiki.woodpecker.org.cn/moin/Octave/Bicubic}
  28 ## for further information.
  29 ## @seealso{interp2}
  30 ## @end deftypefn
  31
  32 ## Bicubic interpolation method.
  33 ## Author: Hoxide Ma <hoxide_dirac@yahoo.com.cn>
  34
  35 function zi = bicubic (x, y, z, xi, yi, extrapval, spline_alpha)
  36
  37   if (nargin < 1 || nargin > 7)
  38     print_usage ();
  39   endif
  40
  41   if (nargin == 7 && isscalar(spline_alpha))
  42     a = spline_alpha;
  43   else
  44     a = 0.5;
  45   endif
  46
  47   if (nargin < 6)
  48     extrapval = NaN;
  49   endif
  50
  51   if (isa (x, "single") || isa (y, "single") || isa (z, "single")
  52       || isa (xi, "single") || isa (yi, "single"))
  53     myeps = eps("single");
  54   else
  55     myeps = eps;
  56   endif
  57
  58   if (nargin <= 2)
  59     ## bicubic (z) or bicubic (z, 2)
  60     if (nargin == 1)
  61       n = 1;
  62     else
  63       n = y;
  64     endif
  65     z = x;
  66     x = [];
  67     [rz, cz] = size (z);
  68     s = linspace (1, cz, (cz-1)*pow2(n)+1);
  69     t = linspace (1, rz, (rz-1)*pow2(n)+1);
  70   elseif (nargin == 3)
  71     if (! isvector (x) || ! isvector (y))
  72       error ("bicubic: XI and YI must be vector");
  73     endif
  74     s = y;
  75     t = z;
  76     z = x;
  77     [rz, cz] = size (z);
  78   elseif (nargin == 5 || nargin == 6)
  79     [rz, cz] = size (z) ;
  80     if (isvector (x) && isvector (y))
  81       if (rz != length (y) || cz != length (x))
  82         error ("bicubic: length of X and Y must match the size of Z");
  83       endif
  84     elseif (size_equal (x, y) && size_equal (x, z))
  85       x = x(1,:);
  86       y = y(:,1);
  87     else
  88       error ("bicubic: X, Y and Z must be equal size matrices of same size");
  89     endif
  90
  91     ## Mark values outside the lookup table.
  92     xfirst_ind = find (xi < x(1));
  93     xlast_ind  = find (xi > x(cz));
  94     yfirst_ind = find (yi < y(1));
  95     ylast_ind  = find (yi > y(rz));
  96     ## Set value outside the table preliminary to min max index.
  97     xi(xfirst_ind) = x(1);
  98     xi(xlast_ind) = x(cz);
  99     yi(yfirst_ind) = y(1);
 100     yi(ylast_ind) = y(rz);
 101
 102
 103     x = reshape (x, 1, cz);
 104     x(cz) *= 1 + sign (x(cz))*myeps;
 105     if (x(cz) == 0)
 106       x(cz) = myeps;
 107     endif;
 108     xi = reshape (xi, 1, length (xi));
 109     [m, i] = sort ([x, xi]);
 110     o = cumsum (i <= cz);
 111     xidx = o(find (i > cz));
 112
 113     y = reshape (y, rz, 1);
 114     y(rz) *= 1 + sign (y(rz))*myeps;
 115     if (y(rz) == 0)
 116       y(rz) = myeps;
 117     endif;
 118     yi = reshape (yi, length (yi), 1);
 119     [m, i] = sort ([y; yi]);
 120     o = cumsum (i <= rz);
 121     yidx = o([find(i > rz)]);
 122
 123     ## Set s and t used follow codes.
 124     s = xidx + ((xi .- x(xidx))./(x(xidx+1) .- x(xidx)));
 125     t = yidx + ((yi - y(yidx))./(y(yidx+1) - y(yidx)));
 126   else
 127     print_usage ();
 128   endif
 129
 130   if (rz < 3 || cz < 3)
 131     error ("bicubic: Z at least a 3 by 3 matrices");
 132   endif
 133
 134   inds = floor (s);
 135   d = find (s == cz);
 136   s = s - floor (s);
 137   inds(d) = cz-1;
 138   s(d) = 1.0;
 139
 140   d = [];
 141   indt = floor (t);
 142   d = find (t == rz);
 143   t = t - floor (t);
 144   indt(d) = rz-1;
 145   t(d) = 1.0;
 146   d = [];
 147
 148   p = zeros (size (z) + 2);
 149   p(2:rz+1,2:cz+1) = z;
 150   p(1,:) =    (6*(1-a))*p(2,:)    - 3*p(3,:)  + (6*a-2)*p(4,:);
 151   p(rz+2,:) = (6*(1-a))*p(rz+1,:) - 3*p(rz,:) + (6*a-2)*p(rz-1,:);
 152   p(:,1) =    (6*(1-a))*p(:,2)    - 3*p(:,3)  + (6*a-2)*p(:,4);
 153   p(:,cz+2) = (6*(1-a))*p(:,cz+1) - 3*p(:,cz) + (6*a-2)*p(:,cz-1);
 154
 155   ## Calculte the C1(t) C2(t) C3(t) C4(t) and C1(s) C2(s) C3(s) C4(s).
 156   t2 = t.*t;
 157   t3 = t2.*t;
 158
 159   ct0 =    -a .* t3 +     (2 * a) .* t2 - a .* t ;      # -a G0
 160   ct1 = (2-a) .* t3 +      (-3+a) .* t2          + 1 ;  # F0 - a G1
 161   ct2 = (a-2) .* t3 + (-2 *a + 3) .* t2 + a .* t ;      # F1 + a G0
 162   ct3 =     a .* t3 -           a .* t2;                # a G1
 163   t = []; t2 = []; t3 = [];
 164
 165   s2 = s.*s;
 166   s3 = s2.*s;
 167
 168   cs0 =    -a .* s3 +     (2 * a) .* s2 - a .*s ;      # -a G0
 169   cs1 = (2-a) .* s3 +    (-3 + a) .* s2         + 1 ;  # F0 - a G1
 170   cs2 = (a-2) .* s3 + (-2 *a + 3) .* s2 + a .*s ;      # F1 + a G0
 171   cs3 =     a .* s3 -           a .* s2;               # a G1
 172   s = []; s2 = []; s3 = [];
 173
 174   cs0 = cs0([1,1,1,1],:);
 175   cs1 = cs1([1,1,1,1],:);
 176   cs2 = cs2([1,1,1,1],:);
 177   cs3 = cs3([1,1,1,1],:);
 178
 179   lent = length (ct0);
 180   lens = columns (cs0);
 181   zi = zeros (lent, lens);
 182
 183   for i = 1:lent
 184     it = indt(i);
 185     int = [it, it+1, it+2, it+3];
 186     zi(i,:) = ([ct0(i),ct1(i),ct2(i),ct3(i)]
 187               * (p(int,inds) .* cs0 + p(int,inds+1) .* cs1
 188                  + p(int,inds+2) .* cs2 + p(int,inds+3) .* cs3));
 189   endfor
 190
 191   ## Set points outside the table to extrapval.
 192   if (! (isempty (xfirst_ind) && isempty (xlast_ind)))
 193     zi(:, [xfirst_ind, xlast_ind]) = extrapval;
 194   endif
 195   if (! (isempty (yfirst_ind) && isempty (ylast_ind)))
 196     zi([yfirst_ind; ylast_ind], :) = extrapval;
 197   endif
 198
 199 endfunction
 200
 201 %!demo
 202 %! A=[13,-1,12;5,4,3;1,6,2];
 203 %! x=[0,1,4]+10; y=[-10,-9,-8];
 204 %! xi=linspace(min(x),max(x),17);
 205 %! yi=linspace(min(y),max(y),26)';
 206 %! mesh(xi,yi,bicubic(x,y,A,xi,yi));
 207 %! [x,y] = meshgrid(x,y);
 208 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;