1 ## Copyright (C) 2005-2012 Hoxide Ma
3 ## This file is part of Octave.
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.
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.
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/>.
20 ## @deftypefn {Function File} {@var{zi} =} bicubic (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi}, @var{extrapval})
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}.
27 ## See @url{http://wiki.woodpecker.org.cn/moin/Octave/Bicubic}
28 ## for further information.
32 ## Bicubic interpolation method.
33 ## Author: Hoxide Ma <hoxide_dirac@yahoo.com.cn>
35 function zi = bicubic (x, y, z, xi, yi, extrapval, spline_alpha)
37 if (nargin < 1 || nargin > 7)
41 if (nargin == 7 && isscalar(spline_alpha))
51 if (isa (x, "single") || isa (y, "single") || isa (z, "single")
52 || isa (xi, "single") || isa (yi, "single"))
53 myeps = eps("single");
59 ## bicubic (z) or bicubic (z, 2)
68 s = linspace (1, cz, (cz-1)*pow2(n)+1);
69 t = linspace (1, rz, (rz-1)*pow2(n)+1);
71 if (! isvector (x) || ! isvector (y))
72 error ("bicubic: XI and YI must be vector");
78 elseif (nargin == 5 || nargin == 6)
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");
84 elseif (size_equal (x, y) && size_equal (x, z))
88 error ("bicubic: X, Y and Z must be equal size matrices of same size");
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);
103 x = reshape (x, 1, cz);
104 x(cz) *= 1 + sign (x(cz))*myeps;
108 xi = reshape (xi, 1, length (xi));
109 [m, i] = sort ([x, xi]);
110 o = cumsum (i <= cz);
111 xidx = o(find (i > cz));
113 y = reshape (y, rz, 1);
114 y(rz) *= 1 + sign (y(rz))*myeps;
118 yi = reshape (yi, length (yi), 1);
119 [m, i] = sort ([y; yi]);
120 o = cumsum (i <= rz);
121 yidx = o([find(i > rz)]);
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)));
130 if (rz < 3 || cz < 3)
131 error ("bicubic: Z at least a 3 by 3 matrices");
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);
155 ## Calculte the C1(t) C2(t) C3(t) C4(t) and C1(s) C2(s) C3(s) C4(s).
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 = [];
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 = [];
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],:);
180 lens = columns (cs0);
181 zi = zeros (lent, lens);
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));
191 ## Set points outside the table to extrapval.
192 if (! (isempty (xfirst_ind) && isempty (xlast_ind)))
193 zi(:, [xfirst_ind, xlast_ind]) = extrapval;
195 if (! (isempty (yfirst_ind) && isempty (ylast_ind)))
196 zi([yfirst_ind; ylast_ind], :) = extrapval;
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;