|
5837
|
1 ## Copyright (C) 2005 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 2, or (at your option) |
|
|
8 ## 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, write to the Free |
|
|
17 ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA |
|
|
18 ## 02110-1301, USA. |
|
|
19 |
|
|
20 ## -*- texinfo -*- |
|
|
21 ## @deftypefn {Function File} {@var{zi}=} bicubic (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi}) |
|
|
22 ## |
|
5838
|
23 ## Return a matrix @var{zi} corresponding to the the bicubic |
|
|
24 ## interpolations at @var{xi} and @var{yi} of the data supplied |
|
5837
|
25 ## as @var{x}, @var{y} and @var{z}. |
|
|
26 ## |
|
|
27 ## For further information please see bicubic.pdf available at |
|
|
28 ## @url{http://wiki.woodpecker.org.cn/moin/Octave/Bicubic} |
|
5838
|
29 ## @seealso{interp2} |
|
5837
|
30 ## @end deftypefn |
|
|
31 |
|
|
32 ## Bicubic interpolation method. |
|
|
33 ## Author: Hoxide Ma <hoxide_dirac@yahoo.com.cn> |
|
|
34 |
|
5838
|
35 function F = bicubic (X, Y, Z, XI, YI, spline_alpha) |
|
5837
|
36 |
|
|
37 if (nargin < 1 || nargin > 6) |
|
|
38 print_usage (); |
|
|
39 endif |
|
|
40 |
|
5838
|
41 if (nargin == 6 && prod (size (spline_alpha)) == 1) |
|
5837
|
42 a = spline_alpha |
|
|
43 else |
|
|
44 a = 0.5; |
|
|
45 endif |
|
|
46 |
|
5838
|
47 if (nargin <= 2) |
|
|
48 ## bicubic (Z) or bicubic (Z, 2) |
|
5837
|
49 if (nargin == 1) |
|
|
50 n = 1; |
|
|
51 else |
|
|
52 n = Y; |
|
|
53 endif |
|
5838
|
54 Z = X; |
|
5837
|
55 X = []; |
|
5838
|
56 [rz, cz] = size (Z); |
|
|
57 s = linspace (1, cz, (cz-1)*pow2(n)+1); |
|
|
58 t = linspace (1, rz, (rz-1)*pow2(n)+1); |
|
5837
|
59 elseif (nargin == 3) |
|
5838
|
60 if (! isvector (X) || ! isvector (Y)) |
|
5837
|
61 error ("XI and YI must be vector"); |
|
|
62 endif |
|
|
63 s = Y; |
|
|
64 t = Z; |
|
|
65 Z = X; |
|
5838
|
66 [rz, cz] = size (Z); |
|
5837
|
67 elseif (nargin == 5 || nargin == 6) |
|
|
68 [rz, cz] = size (Z) ; |
|
|
69 if (isvector (X) && isvector (Y)) |
|
|
70 if (rz != length (Y) || cz != length (X)) |
|
|
71 error ("length of X and Y must match the size of Z"); |
|
|
72 endif |
|
6157
|
73 elseif (size_equal (X, Y) && size_equal (X, Z)) |
|
5837
|
74 X = X(1,:); |
|
|
75 Y = Y(:,1); |
|
|
76 else |
|
|
77 error ("X, Y and Z must be martrices of same size"); |
|
|
78 endif |
|
|
79 |
|
|
80 ## mark values outside the lookup table |
|
|
81 xfirst_ind = find (XI < X(1)); |
|
|
82 xlast_ind = find (XI > X(cz)); |
|
|
83 yfirst_ind = find (YI < Y(1)); |
|
|
84 ylast_ind = find (YI > Y(rz)); |
|
|
85 ## set value outside the table preliminary to min max index |
|
|
86 XI(xfirst_ind) = X(1); |
|
|
87 XI(xlast_ind) = X(cz); |
|
|
88 YI(yfirst_ind) = Y(1); |
|
|
89 YI(ylast_ind) = Y(rz); |
|
|
90 |
|
|
91 |
|
5838
|
92 X = reshape (X, 1, cz); |
|
|
93 X(cz) *= 1 + sign (X(cz))*eps; |
|
5837
|
94 if (X(cz) == 0) |
|
|
95 X(cz) = eps; |
|
|
96 endif; |
|
5838
|
97 XI = reshape (XI, 1, length (XI)); |
|
|
98 [m, i] = sort ([X, XI]); |
|
|
99 o = cumsum (i <= cz); |
|
|
100 xidx = o(find (i > cz)); |
|
5837
|
101 |
|
5838
|
102 Y = reshape (Y, rz, 1); |
|
|
103 Y(rz) *= 1 + sign (Y(rz))*eps; |
|
5837
|
104 if (Y(rz) == 0) |
|
|
105 Y(rz) = eps; |
|
|
106 endif; |
|
5838
|
107 YI = reshape (YI, length (YI), 1); |
|
|
108 [m, i] = sort ([Y; YI]); |
|
|
109 o = cumsum (i <= rz); |
|
5837
|
110 yidx = o([find( i> rz)]); |
|
|
111 |
|
|
112 ## set s and t used follow codes |
|
|
113 s = xidx + ((XI .- X(xidx))./(X(xidx+1) .- X(xidx))); |
|
|
114 t = yidx + ((YI - Y(yidx))./(Y(yidx+1) - Y(yidx))); |
|
|
115 else |
|
|
116 print_usage (); |
|
|
117 endif |
|
|
118 |
|
|
119 if (rz < 3 || cz < 3) |
|
|
120 error ("Z at least a 3 by 3 matrices"); |
|
|
121 endif |
|
|
122 |
|
5838
|
123 inds = floor (s); |
|
|
124 d = find (s == cz); |
|
|
125 s = s - floor (s); |
|
5837
|
126 inds(d) = cz-1; |
|
|
127 s(d) = 1.0; |
|
|
128 |
|
|
129 d = []; |
|
5838
|
130 indt = floor (t); |
|
|
131 d = find (t == rz); |
|
|
132 t = t - floor (t); |
|
5837
|
133 indt(d) = rz-1; |
|
|
134 t(d) = 1.0; |
|
|
135 d = []; |
|
|
136 |
|
5838
|
137 p = zeros (size (Z) + 2); |
|
5837
|
138 p(2:rz+1,2:cz+1) = Z; |
|
|
139 p(1,:) = (6*(1-a))*p(2,:) -3*p(3,:) + (6*a-2)*p(4,:); |
|
|
140 p(rz+2,:) = (6*(1-a))*p(rz+1,:) -3*p(rz,:) + (6*a-2)*p(rz-1,:); |
|
|
141 p(:,1) = (6*(1-a))*p(:,2) -3*p(:,3) + (6*a-2)*p(:,4); |
|
|
142 p(:,cz+2) = (6*(1-a))*p(:,cz+1) -3*p(:,cz) + (6*a-2)*p(:,cz-1); |
|
|
143 |
|
|
144 ## calculte the C1(t) C2(t) C3(t) C4(t) and C1(s) C2(s) C3(s) C4(s) |
|
|
145 t2= t.*t; |
|
|
146 t3= t2.*t; |
|
|
147 |
|
|
148 ct0= -a .* t3 + (2 * a) .* t2 - a .* t ; # -a G0 |
|
|
149 ct1 = (2-a) .* t3 + (-3+a) .* t2 + 1 ; # F0 - a G1 |
|
|
150 ct2 = (a-2) .* t3 + (-2 *a + 3) .* t2 + a .* t ; # F1 + a G0 |
|
|
151 ct3 = a .* t3 - a .* t2; # a G1 |
|
|
152 t = [];t2=[]; t3=[]; |
|
|
153 |
|
|
154 s2= s.*s; |
|
|
155 s3= s2.*s; |
|
|
156 |
|
|
157 cs0= -a .* s3 + (2 * a) .* s2 - a .*s ; # -a G0 |
|
|
158 cs1 = (2-a) .* s3 + (-3 + a) .* s2 + 1 ; # F0 - a G1 |
|
|
159 cs2 = (a-2) .* s3 + (-2 *a + 3) .* s2 + a .*s ; # F1 + a G0 |
|
|
160 cs3 = a .* s3 - a .* s2; # a G1 |
|
|
161 s=[] ; s2 = []; s3 = []; |
|
|
162 |
|
|
163 cs0 = cs0([1,1,1,1],:); |
|
|
164 cs1 = cs1([1,1,1,1],:); |
|
|
165 cs2 = cs2([1,1,1,1],:); |
|
|
166 cs3 = cs3([1,1,1,1],:); |
|
|
167 |
|
5838
|
168 lent = length (ct0); |
|
|
169 lens = length (cs0); |
|
|
170 F = zeros (lent, lens); |
|
5837
|
171 |
|
|
172 for i = 1:lent |
|
|
173 it = indt(i); |
|
|
174 int = [it, it+1, it+2, it+3]; |
|
|
175 F(i,:) = [ct0(i),ct1(i),ct2(i),ct3(i)] * ... |
|
5838
|
176 (p(int,inds) .* cs0 + p(int,inds+1) .* cs1 + ... |
|
|
177 p(int,inds+2) .* cs2 + p(int,inds+3) .* cs3); |
|
5837
|
178 endfor |
|
|
179 |
|
|
180 ## set points outside the table to NaN |
|
|
181 if (! (isempty (xfirst_ind) && isempty (xlast_ind))) |
|
|
182 F(:, [xfirst_ind, xlast_ind]) = NaN; |
|
|
183 endif |
|
|
184 if (! (isempty (yfirst_ind) && isempty (ylast_ind))) |
|
|
185 F([yfirst_ind; ylast_ind], :) = NaN; |
|
|
186 endif |
|
|
187 |
|
|
188 endfunction |
|
|
189 |
|
|
190 %!demo |
|
|
191 %! A=[13,-1,12;5,4,3;1,6,2]; |
|
|
192 %! x=[0,1,4]+10; y=[-10,-9,-8]; |
|
|
193 %! xi=linspace(min(x),max(x),17); |
|
|
194 %! yi=linspace(min(y),max(y),26); |
|
|
195 %! mesh(xi,yi,bicubic(x,y,A,xi,yi)); |
|
|
196 %! [x,y] = meshgrid(x,y); |
|
|
197 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; |
