Mercurial > hg > octave-nkf

comparison scripts/general/bicubic.m @ 11470:eb9e0b597d61

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Use common names for variables in documentation and code for a few more m-script files.
author Rik <octave@nomad.inbox5.com>
date Sun, 09 Jan 2011 13:44:15 -0800
parents fe3c3dfc07eb
children fd0a3ac60b0e
comparison
equal deleted inserted replaced
11469:c776f063fefe 11470:eb9e0b597d61
30 ## @end deftypefn 30 ## @end deftypefn
31 31
32 ## Bicubic interpolation method. 32 ## Bicubic interpolation method.
33 ## Author: Hoxide Ma <hoxide_dirac@yahoo.com.cn> 33 ## Author: Hoxide Ma <hoxide_dirac@yahoo.com.cn>
34 34
35 function F = bicubic (X, Y, Z, XI, YI, extrapval, spline_alpha) 35 function zi = bicubic (x, y, z, xi, yi, extrapval, spline_alpha)
36 36
37 if (nargin < 1 || nargin > 7) 37 if (nargin < 1 || nargin > 7)
38 print_usage (); 38 print_usage ();
39 endif 39 endif
40 40
46 46
47 if (nargin < 6) 47 if (nargin < 6)
48 extrapval = NaN; 48 extrapval = NaN;
49 endif 49 endif
50 50
51 if (isa (X, "single") || isa (Y, "single") || isa (Z, "single") 51 if (isa (x, "single") || isa (y, "single") || isa (z, "single")
52 || isa (XI, "single") || isa (YI, "single")) 52 || isa (xi, "single") || isa (yi, "single"))
53 myeps = eps("single"); 53 myeps = eps("single");
54 else 54 else
55 myeps = eps; 55 myeps = eps;
56 endif 56 endif
57 57
58 if (nargin <= 2) 58 if (nargin <= 2)
59 ## bicubic (Z) or bicubic (Z, 2) 59 ## bicubic (z) or bicubic (z, 2)
60 if (nargin == 1) 60 if (nargin == 1)
61 n = 1; 61 n = 1;
62 else 62 else
63 n = Y; 63 n = y;
64 endif 64 endif
65 Z = X; 65 z = x;
66 X = []; 66 x = [];
67 [rz, cz] = size (Z); 67 [rz, cz] = size (z);
68 s = linspace (1, cz, (cz-1)*pow2(n)+1); 68 s = linspace (1, cz, (cz-1)*pow2(n)+1);
69 t = linspace (1, rz, (rz-1)*pow2(n)+1); 69 t = linspace (1, rz, (rz-1)*pow2(n)+1);
70 elseif (nargin == 3) 70 elseif (nargin == 3)
71 if (! isvector (X) || ! isvector (Y)) 71 if (! isvector (x) || ! isvector (y))
72 error ("bicubic: XI and YI must be vector"); 72 error ("bicubic: XI and YI must be vector");
73 endif 73 endif
74 s = Y; 74 s = y;
75 t = Z; 75 t = z;
76 Z = X; 76 z = x;
77 [rz, cz] = size (Z); 77 [rz, cz] = size (z);
78 elseif (nargin == 5 || nargin == 6) 78 elseif (nargin == 5 || nargin == 6)
79 [rz, cz] = size (Z) ; 79 [rz, cz] = size (z) ;
80 if (isvector (X) && isvector (Y)) 80 if (isvector (x) && isvector (y))
81 if (rz != length (Y) || cz != length (X)) 81 if (rz != length (y) || cz != length (x))
82 error ("bicubic: length of X and Y must match the size of Z"); 82 error ("bicubic: length of X and Y must match the size of Z");
83 endif 83 endif
84 elseif (size_equal (X, Y) && size_equal (X, Z)) 84 elseif (size_equal (x, y) && size_equal (x, z))
85 X = X(1,:); 85 x = x(1,:);
86 Y = Y(:,1); 86 y = y(:,1);
87 else 87 else
88 error ("bicubic: X, Y and Z must be equal size matrices of same size"); 88 error ("bicubic: X, Y and Z must be equal size matrices of same size");
89 endif 89 endif
90 90
91 ## Mark values outside the lookup table. 91 ## Mark values outside the lookup table.
92 xfirst_ind = find (XI < X(1)); 92 xfirst_ind = find (xi < x(1));
93 xlast_ind = find (XI > X(cz)); 93 xlast_ind = find (xi > x(cz));
94 yfirst_ind = find (YI < Y(1)); 94 yfirst_ind = find (yi < y(1));
95 ylast_ind = find (YI > Y(rz)); 95 ylast_ind = find (yi > y(rz));
96 ## Set value outside the table preliminary to min max index. 96 ## Set value outside the table preliminary to min max index.
97 XI(xfirst_ind) = X(1); 97 xi(xfirst_ind) = x(1);
98 XI(xlast_ind) = X(cz); 98 xi(xlast_ind) = x(cz);
99 YI(yfirst_ind) = Y(1); 99 yi(yfirst_ind) = y(1);
100 YI(ylast_ind) = Y(rz); 100 yi(ylast_ind) = y(rz);
101 101
102 102
103 X = reshape (X, 1, cz); 103 x = reshape (x, 1, cz);
104 X(cz) *= 1 + sign (X(cz))*myeps; 104 x(cz) *= 1 + sign (x(cz))*myeps;
105 if (X(cz) == 0) 105 if (x(cz) == 0)
106 X(cz) = myeps; 106 x(cz) = myeps;
107 endif; 107 endif;
108 XI = reshape (XI, 1, length (XI)); 108 xi = reshape (xi, 1, length (xi));
109 [m, i] = sort ([X, XI]); 109 [m, i] = sort ([x, xi]);
110 o = cumsum (i <= cz); 110 o = cumsum (i <= cz);
111 xidx = o(find (i > cz)); 111 xidx = o(find (i > cz));
112 112
113 Y = reshape (Y, rz, 1); 113 y = reshape (y, rz, 1);
114 Y(rz) *= 1 + sign (Y(rz))*myeps; 114 y(rz) *= 1 + sign (y(rz))*myeps;
115 if (Y(rz) == 0) 115 if (y(rz) == 0)
116 Y(rz) = myeps; 116 y(rz) = myeps;
117 endif; 117 endif;
118 YI = reshape (YI, length (YI), 1); 118 yi = reshape (yi, length (yi), 1);
119 [m, i] = sort ([Y; YI]); 119 [m, i] = sort ([y; yi]);
120 o = cumsum (i <= rz); 120 o = cumsum (i <= rz);
121 yidx = o([find(i > rz)]); 121 yidx = o([find(i > rz)]);
122 122
123 ## Set s and t used follow codes. 123 ## Set s and t used follow codes.
124 s = xidx + ((XI .- X(xidx))./(X(xidx+1) .- X(xidx))); 124 s = xidx + ((xi .- x(xidx))./(x(xidx+1) .- x(xidx)));
125 t = yidx + ((YI - Y(yidx))./(Y(yidx+1) - Y(yidx))); 125 t = yidx + ((yi - y(yidx))./(y(yidx+1) - y(yidx)));
126 else 126 else
127 print_usage (); 127 print_usage ();
128 endif 128 endif
129 129
130 if (rz < 3 || cz < 3) 130 if (rz < 3 || cz < 3)
143 t = t - floor (t); 143 t = t - floor (t);
144 indt(d) = rz-1; 144 indt(d) = rz-1;
145 t(d) = 1.0; 145 t(d) = 1.0;
146 d = []; 146 d = [];
147 147
148 p = zeros (size (Z) + 2); 148 p = zeros (size (z) + 2);
149 p(2:rz+1,2:cz+1) = Z; 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,:); 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,:); 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); 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); 153 p(:,cz+2) = (6*(1-a))*p(:,cz+1) - 3*p(:,cz) + (6*a-2)*p(:,cz-1);
154 154
176 cs2 = cs2([1,1,1,1],:); 176 cs2 = cs2([1,1,1,1],:);
177 cs3 = cs3([1,1,1,1],:); 177 cs3 = cs3([1,1,1,1],:);
178 178
179 lent = length (ct0); 179 lent = length (ct0);
180 lens = columns (cs0); 180 lens = columns (cs0);
181 F = zeros (lent, lens); 181 zi = zeros (lent, lens);
182 182
183 for i = 1:lent 183 for i = 1:lent
184 it = indt(i); 184 it = indt(i);
185 int = [it, it+1, it+2, it+3]; 185 int = [it, it+1, it+2, it+3];
186 F(i,:) = ([ct0(i),ct1(i),ct2(i),ct3(i)] 186 zi(i,:) = ([ct0(i),ct1(i),ct2(i),ct3(i)]
187 * (p(int,inds) .* cs0 + p(int,inds+1) .* cs1 187 * (p(int,inds) .* cs0 + p(int,inds+1) .* cs1
188 + p(int,inds+2) .* cs2 + p(int,inds+3) .* cs3)); 188 + p(int,inds+2) .* cs2 + p(int,inds+3) .* cs3));
189 endfor 189 endfor
190 190
191 ## Set points outside the table to extrapval. 191 ## Set points outside the table to extrapval.
192 if (! (isempty (xfirst_ind) && isempty (xlast_ind))) 192 if (! (isempty (xfirst_ind) && isempty (xlast_ind)))
193 F(:, [xfirst_ind, xlast_ind]) = extrapval; 193 zi(:, [xfirst_ind, xlast_ind]) = extrapval;
194 endif 194 endif
195 if (! (isempty (yfirst_ind) && isempty (ylast_ind))) 195 if (! (isempty (yfirst_ind) && isempty (ylast_ind)))
196 F([yfirst_ind; ylast_ind], :) = extrapval; 196 zi([yfirst_ind; ylast_ind], :) = extrapval;
197 endif 197 endif
198 198
199 endfunction 199 endfunction
200 200
201 %!demo 201 %!demo