|
6702
|
1 /* |
|
|
2 |
|
|
3 Copyright (C) 2007 Alexander Barth |
|
|
4 |
|
|
5 This file is part of Octave. |
|
|
6 |
|
|
7 Octave is free software; you can redistribute it and/or modify it |
|
|
8 under the terms of the GNU General Public License as published by the |
|
|
9 Free Software Foundation; either version 2, or (at your option) any |
|
|
10 later version. |
|
|
11 |
|
|
12 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
15 for more details. |
|
|
16 |
|
|
17 You should have received a copy of the GNU General Public License |
|
|
18 along with Octave; see the file COPYING. If not, write to the Free |
|
|
19 Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA |
|
|
20 02110-1301, USA. |
|
|
21 |
|
|
22 */ |
|
|
23 |
|
|
24 #ifdef HAVE_CONFIG_H |
|
|
25 #include <config.h> |
|
|
26 #endif |
|
|
27 |
|
|
28 #include "dNDArray.h" |
|
|
29 |
|
|
30 #include "defun-dld.h" |
|
|
31 #include "error.h" |
|
|
32 #include "oct-obj.h" |
|
|
33 |
|
|
34 // equivalent to isvector.m |
|
|
35 |
|
|
36 bool |
|
|
37 isvector (const NDArray& array) |
|
|
38 { |
|
|
39 const dim_vector dv = array.dims (); |
|
|
40 return dv.length () == 2 && (dv(0) == 1 || dv(1) == 1); |
|
|
41 } |
|
|
42 |
|
|
43 // lookup a value in a sorted table (lookup.m) |
|
|
44 octave_idx_type |
|
|
45 lookup (const double *x, octave_idx_type n, double y) |
|
|
46 { |
|
|
47 octave_idx_type j, j0, j1; |
|
|
48 |
|
|
49 if (y > x[n-1] || y < x[0]) |
|
|
50 return -1; |
|
|
51 |
|
|
52 #ifdef EXHAUSTIF |
|
|
53 for (j = 0; j < n - 1; j++) |
|
|
54 { |
|
|
55 if (x[j] <= y && y <= x[j+1]) |
|
|
56 return j; |
|
|
57 } |
|
|
58 #else |
|
|
59 j0 = 0; |
|
|
60 j1 = n - 1; |
|
|
61 |
|
|
62 while (true) |
|
|
63 { |
|
|
64 j = (j0+j1)/2; |
|
|
65 |
|
|
66 if (y <= x[j+1]) |
|
|
67 { |
|
|
68 if (x[j] <= y) |
|
|
69 return j; |
|
|
70 |
|
|
71 j1 = j; |
|
|
72 } |
|
|
73 |
|
|
74 if (x[j] <= y) |
|
|
75 j0 = j; |
|
|
76 } |
|
|
77 |
|
|
78 #endif |
|
|
79 } |
|
|
80 |
|
|
81 // n-dimensional linear interpolation |
|
|
82 |
|
|
83 void |
|
|
84 lin_interpn (int n, const octave_idx_type *size, const octave_idx_type *scale, |
|
|
85 octave_idx_type Ni, double extrapval, const double **x, |
|
|
86 const double *v, const double **y, double *vi) |
|
|
87 { |
|
|
88 bool out = false; |
|
|
89 int bit; |
|
|
90 |
|
|
91 OCTAVE_LOCAL_BUFFER (double, coef, 2*n); |
|
|
92 OCTAVE_LOCAL_BUFFER (octave_idx_type, index, n); |
|
|
93 |
|
|
94 // loop over all points |
|
|
95 for (octave_idx_type m = 0; m < Ni; m++) |
|
|
96 { |
|
|
97 // loop over all dimensions |
|
|
98 for (int i = 0; i < n; i++) |
|
|
99 { |
|
|
100 index[i] = lookup (x[i], size[i], y[i][m]); |
|
|
101 out = index[i] == -1; |
|
|
102 |
|
|
103 if (out) |
|
|
104 break; |
|
|
105 else |
|
|
106 { |
|
|
107 octave_idx_type j = index[i]; |
|
|
108 coef[2*i+1] = (y[i][m] - x[i][j])/(x[i][j+1] - x[i][j]); |
|
|
109 coef[2*i] = 1 - coef[2*i+1]; |
|
|
110 } |
|
|
111 } |
|
|
112 |
|
|
113 |
|
|
114 if (out) |
|
|
115 vi[m] = extrapval; |
|
|
116 else |
|
|
117 { |
|
|
118 vi[m] = 0; |
|
|
119 |
|
|
120 // loop over all corners of hypercube (1<<n = 2^n) |
|
|
121 for (int i = 0; i < (1 << n); i++) |
|
|
122 { |
|
|
123 double c = 1; |
|
|
124 octave_idx_type l = 0; |
|
|
125 |
|
|
126 // loop over all dimensions |
|
|
127 for (int j = 0; j < n; j++) |
|
|
128 { |
|
|
129 // test if the jth bit in i is set |
|
|
130 bit = i >> j & 1; |
|
|
131 l += scale[j] * (index[j] + bit); |
|
|
132 c *= coef[2*j+bit]; |
|
|
133 } |
|
|
134 |
|
|
135 vi[m] += c * v[l]; |
|
|
136 } |
|
|
137 } |
|
|
138 } |
|
|
139 } |
|
|
140 |
|
|
141 DEFUN_DLD (__lin_interpn__, args, , |
|
|
142 "-*- texinfo -*-\n\ |
|
|
143 @deftypefn {Loadable Function} {@var{vi} =} __lin_interpn__ (@var{x1}, @var{x2}, @dots{}, @var{xn}, @var{v}, @var{y1}, @var{y2}, @dots{}, @var{yn})\n\ |
|
|
144 Perform @var{n}-dimensional interpolation. Each element of then\n\ |
|
|
145 @var{n}-dimensional array @var{v} represents a value at a location\n\ |
|
|
146 given by the parameters @var{x1}, @var{x2},...,@var{xn}. The parameters\n\ |
|
|
147 @var{x1}, @var{x2}, @dots{}, @var{xn} are either @var{n}-dimensional\n\ |
|
|
148 arrays of the same size as the array @var{v} in the \"ndgrid\" format\n\ |
|
|
149 or vectors. The parameters @var{y1}, @var{y2}, @dots{}, @var{yn} are\n\ |
|
|
150 all @var{n}-dimensional arrays of the same size and represent the\n\ |
|
|
151 points at which the array @var{vi} is interpolated.\n\ |
|
|
152 \n\ |
|
|
153 This function only performs linear interpolation.\n\ |
|
|
154 @seealso{interp1, interp2, ndgrid}\n\ |
|
|
155 @end deftypefn") |
|
|
156 { |
|
|
157 octave_value retval; |
|
|
158 |
|
|
159 int nargin = args.length (); |
|
|
160 |
|
|
161 if (nargin < 2 || nargin % 2 == 0) |
|
|
162 { |
|
|
163 print_usage (); |
|
|
164 return retval; |
|
|
165 } |
|
|
166 |
|
|
167 // dimension of the problem |
|
|
168 int n = (nargin-1)/2; |
|
|
169 |
|
|
170 OCTAVE_LOCAL_BUFFER (NDArray, X, n); |
|
|
171 OCTAVE_LOCAL_BUFFER (NDArray, Y, n); |
|
|
172 |
|
|
173 OCTAVE_LOCAL_BUFFER (const double *, x, n); |
|
|
174 OCTAVE_LOCAL_BUFFER (const double *, y, n); |
|
|
175 OCTAVE_LOCAL_BUFFER (octave_idx_type, scale, n); |
|
|
176 OCTAVE_LOCAL_BUFFER (octave_idx_type, size, n); |
|
|
177 |
|
|
178 const NDArray V = args(n).array_value (); |
|
|
179 NDArray Vi = NDArray (args(n+1).dims ()); |
|
|
180 |
|
|
181 if (error_state) |
|
|
182 { |
|
|
183 print_usage (); |
|
|
184 return retval; |
|
|
185 } |
|
|
186 |
|
|
187 const double *v = V.data (); |
|
|
188 double *vi = Vi.fortran_vec (); |
|
|
189 octave_idx_type Ni = Vi.numel (); |
|
|
190 |
|
6742
|
191 double extrapval = octave_NA; |
|
6702
|
192 |
|
|
193 for (int i = 0; i < n; i++) |
|
|
194 { |
|
|
195 X[i] = args(i).array_value (); |
|
|
196 Y[i] = args(n+i+1).array_value (); |
|
|
197 |
|
|
198 if (error_state) |
|
|
199 { |
|
|
200 print_usage (); |
|
|
201 return retval; |
|
|
202 } |
|
|
203 |
|
|
204 y[i] = Y[i].data (); |
|
|
205 size[i] = V.dims()(i); |
|
|
206 |
|
|
207 if (Y[0].dims () != Y[i].dims ()) |
|
|
208 { |
|
|
209 error ("interpn: incompatible size of argument number %d", n+i+2); |
|
|
210 return retval; |
|
|
211 } |
|
|
212 } |
|
|
213 |
|
|
214 // offset in memory of each dimension |
|
|
215 |
|
|
216 scale[0] = 1; |
|
|
217 |
|
|
218 for (int i = 1; i < n; i++) |
|
|
219 scale[i] = scale[i-1] * size[i-1]; |
|
|
220 |
|
|
221 // tests if X[0] is a vector, if yes, assume that all elements of X are |
|
|
222 // in the ndgrid format. |
|
|
223 |
|
|
224 if (! isvector (X[0])) |
|
|
225 { |
|
|
226 for (int i = 0; i < n; i++) |
|
|
227 { |
|
|
228 if (X[i].dims () != V.dims ()) |
|
|
229 { |
|
|
230 error ("interpn: incompatible size of argument number %d", i+1); |
|
|
231 return retval; |
|
|
232 } |
|
|
233 else |
|
|
234 { |
|
|
235 NDArray tmp = NDArray (dim_vector (size[i], 1)); |
|
|
236 |
|
|
237 for (octave_idx_type j = 0; j < size[i]; j++) |
|
|
238 tmp(j) = X[i](scale[i]*j); |
|
|
239 |
|
|
240 X[i] = tmp; |
|
|
241 } |
|
|
242 } |
|
|
243 } |
|
|
244 |
|
|
245 for (int i = 0; i < n; i++) |
|
|
246 { |
|
|
247 if (! isvector (X[i]) && X[i].numel () != size[i]) |
|
|
248 { |
|
|
249 error ("interpn: incompatible size of argument number %d", i+1); |
|
|
250 return retval; |
|
|
251 } |
|
|
252 else |
|
|
253 x[i] = X[i].data (); |
|
|
254 } |
|
|
255 |
|
|
256 lin_interpn (n, size, scale, Ni, extrapval, x, v, y, vi); |
|
|
257 |
|
|
258 retval = Vi; |
|
|
259 |
|
|
260 return retval; |
|
|
261 } |
