Mercurial > octave
annotate src/DLD-FUNCTIONS/__lin_interpn__.cc @ 7561:a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
author | Alexander Barth |
---|---|
date | Thu, 06 Mar 2008 01:56:55 -0500 |
parents | 93c65f2a5668 |
children | 82be108cc558 72830070a17b |
rev | line source |
---|---|
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 | |
7016 | 9 Free Software Foundation; either version 3 of the License, or (at your |
10 option) any later version. | |
6702 | 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 | |
7016 | 18 along with Octave; see the file COPYING. If not, see |
19 <http://www.gnu.org/licenses/>. | |
6702 | 20 |
21 */ | |
22 | |
23 #ifdef HAVE_CONFIG_H | |
24 #include <config.h> | |
25 #endif | |
26 | |
27 #include "dNDArray.h" | |
28 | |
29 #include "defun-dld.h" | |
30 #include "error.h" | |
31 #include "oct-obj.h" | |
32 | |
33 // equivalent to isvector.m | |
34 | |
35 bool | |
36 isvector (const NDArray& array) | |
37 { | |
38 const dim_vector dv = array.dims (); | |
39 return dv.length () == 2 && (dv(0) == 1 || dv(1) == 1); | |
40 } | |
41 | |
42 // lookup a value in a sorted table (lookup.m) | |
43 octave_idx_type | |
44 lookup (const double *x, octave_idx_type n, double y) | |
45 { | |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
46 octave_idx_type j; |
6702 | 47 |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
48 if (x[0] < x[n-1]) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
49 { |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
50 // increasing x |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
51 |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
52 if (y > x[n-1] || y < x[0]) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
53 return -1; |
6702 | 54 |
55 #ifdef EXHAUSTIF | |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
56 for (j = 0; j < n - 1; j++) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
57 { |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
58 if (x[j] <= y && y <= x[j+1]) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
59 return j; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
60 } |
6702 | 61 #else |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
62 octave_idx_type j0 = 0; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
63 octave_idx_type j1 = n - 1; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
64 |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
65 while (true) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
66 { |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
67 j = (j0+j1)/2; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
68 |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
69 if (y <= x[j+1]) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
70 { |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
71 if (x[j] <= y) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
72 return j; |
6702 | 73 |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
74 j1 = j; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
75 } |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
76 |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
77 if (x[j] <= y) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
78 j0 = j; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
79 } |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
80 #endif |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
81 } |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
82 else |
6702 | 83 { |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
84 // decreasing x |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
85 // previous code with x -> -x and y -> -y |
6702 | 86 |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
87 if (y > x[0] || y < x[n-1]) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
88 return -1; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
89 |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
90 #ifdef EXHAUSTIF |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
91 for (j = 0; j < n - 1; j++) |
6702 | 92 { |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
93 if (x[j+1] <= y && y <= x[j]) |
6702 | 94 return j; |
95 } | |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
96 #else |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
97 octave_idx_type j0 = 0; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
98 octave_idx_type j1 = n - 1; |
6702 | 99 |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
100 while (true) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
101 { |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
102 j = (j0+j1)/2; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
103 |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
104 if (y >= x[j+1]) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
105 { |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
106 if (x[j] >= y) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
107 return j; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
108 |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
109 j1 = j; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
110 } |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
111 |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
112 if (x[j] >= y) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
113 j0 = j; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
114 } |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
115 #endif |
6702 | 116 } |
117 } | |
118 | |
119 // n-dimensional linear interpolation | |
120 | |
121 void | |
122 lin_interpn (int n, const octave_idx_type *size, const octave_idx_type *scale, | |
123 octave_idx_type Ni, double extrapval, const double **x, | |
124 const double *v, const double **y, double *vi) | |
125 { | |
126 bool out = false; | |
127 int bit; | |
128 | |
129 OCTAVE_LOCAL_BUFFER (double, coef, 2*n); | |
130 OCTAVE_LOCAL_BUFFER (octave_idx_type, index, n); | |
131 | |
132 // loop over all points | |
133 for (octave_idx_type m = 0; m < Ni; m++) | |
134 { | |
135 // loop over all dimensions | |
136 for (int i = 0; i < n; i++) | |
137 { | |
138 index[i] = lookup (x[i], size[i], y[i][m]); | |
139 out = index[i] == -1; | |
140 | |
141 if (out) | |
142 break; | |
143 else | |
144 { | |
145 octave_idx_type j = index[i]; | |
146 coef[2*i+1] = (y[i][m] - x[i][j])/(x[i][j+1] - x[i][j]); | |
147 coef[2*i] = 1 - coef[2*i+1]; | |
148 } | |
149 } | |
150 | |
151 | |
152 if (out) | |
153 vi[m] = extrapval; | |
154 else | |
155 { | |
156 vi[m] = 0; | |
157 | |
158 // loop over all corners of hypercube (1<<n = 2^n) | |
159 for (int i = 0; i < (1 << n); i++) | |
160 { | |
161 double c = 1; | |
162 octave_idx_type l = 0; | |
163 | |
164 // loop over all dimensions | |
165 for (int j = 0; j < n; j++) | |
166 { | |
167 // test if the jth bit in i is set | |
168 bit = i >> j & 1; | |
169 l += scale[j] * (index[j] + bit); | |
170 c *= coef[2*j+bit]; | |
171 } | |
172 | |
173 vi[m] += c * v[l]; | |
174 } | |
175 } | |
176 } | |
177 } | |
178 | |
6945 | 179 // Perform @var{n}-dimensional interpolation. Each element of then |
180 // @var{n}-dimensional array @var{v} represents a value at a location | |
181 // given by the parameters @var{x1}, @var{x2},...,@var{xn}. The parameters | |
182 // @var{x1}, @var{x2}, @dots{}, @var{xn} are either @var{n}-dimensional | |
183 // arrays of the same size as the array @var{v} in the \"ndgrid\" format | |
184 // or vectors. The parameters @var{y1}, @var{y2}, @dots{}, @var{yn} are | |
185 // all @var{n}-dimensional arrays of the same size and represent the | |
186 // points at which the array @var{vi} is interpolated. | |
187 // | |
188 //This function only performs linear interpolation. | |
189 | |
6702 | 190 DEFUN_DLD (__lin_interpn__, args, , |
191 "-*- texinfo -*-\n\ | |
192 @deftypefn {Loadable Function} {@var{vi} =} __lin_interpn__ (@var{x1}, @var{x2}, @dots{}, @var{xn}, @var{v}, @var{y1}, @var{y2}, @dots{}, @var{yn})\n\ | |
6945 | 193 Undocumented internal function.\n\ |
6702 | 194 @end deftypefn") |
195 { | |
196 octave_value retval; | |
197 | |
198 int nargin = args.length (); | |
199 | |
200 if (nargin < 2 || nargin % 2 == 0) | |
201 { | |
202 print_usage (); | |
203 return retval; | |
204 } | |
205 | |
206 // dimension of the problem | |
207 int n = (nargin-1)/2; | |
208 | |
209 OCTAVE_LOCAL_BUFFER (NDArray, X, n); | |
210 OCTAVE_LOCAL_BUFFER (NDArray, Y, n); | |
211 | |
212 OCTAVE_LOCAL_BUFFER (const double *, x, n); | |
213 OCTAVE_LOCAL_BUFFER (const double *, y, n); | |
214 OCTAVE_LOCAL_BUFFER (octave_idx_type, scale, n); | |
215 OCTAVE_LOCAL_BUFFER (octave_idx_type, size, n); | |
216 | |
217 const NDArray V = args(n).array_value (); | |
218 NDArray Vi = NDArray (args(n+1).dims ()); | |
219 | |
220 if (error_state) | |
221 { | |
222 print_usage (); | |
223 return retval; | |
224 } | |
225 | |
226 const double *v = V.data (); | |
227 double *vi = Vi.fortran_vec (); | |
228 octave_idx_type Ni = Vi.numel (); | |
229 | |
6742 | 230 double extrapval = octave_NA; |
6702 | 231 |
232 for (int i = 0; i < n; i++) | |
233 { | |
234 X[i] = args(i).array_value (); | |
235 Y[i] = args(n+i+1).array_value (); | |
236 | |
237 if (error_state) | |
238 { | |
239 print_usage (); | |
240 return retval; | |
241 } | |
242 | |
243 y[i] = Y[i].data (); | |
244 size[i] = V.dims()(i); | |
245 | |
246 if (Y[0].dims () != Y[i].dims ()) | |
247 { | |
248 error ("interpn: incompatible size of argument number %d", n+i+2); | |
249 return retval; | |
250 } | |
251 } | |
252 | |
253 // offset in memory of each dimension | |
254 | |
255 scale[0] = 1; | |
256 | |
257 for (int i = 1; i < n; i++) | |
258 scale[i] = scale[i-1] * size[i-1]; | |
259 | |
260 // tests if X[0] is a vector, if yes, assume that all elements of X are | |
261 // in the ndgrid format. | |
262 | |
263 if (! isvector (X[0])) | |
264 { | |
265 for (int i = 0; i < n; i++) | |
266 { | |
267 if (X[i].dims () != V.dims ()) | |
268 { | |
269 error ("interpn: incompatible size of argument number %d", i+1); | |
270 return retval; | |
271 } | |
272 else | |
273 { | |
274 NDArray tmp = NDArray (dim_vector (size[i], 1)); | |
275 | |
276 for (octave_idx_type j = 0; j < size[i]; j++) | |
277 tmp(j) = X[i](scale[i]*j); | |
278 | |
279 X[i] = tmp; | |
280 } | |
281 } | |
282 } | |
283 | |
284 for (int i = 0; i < n; i++) | |
285 { | |
286 if (! isvector (X[i]) && X[i].numel () != size[i]) | |
287 { | |
288 error ("interpn: incompatible size of argument number %d", i+1); | |
289 return retval; | |
290 } | |
291 else | |
292 x[i] = X[i].data (); | |
293 } | |
294 | |
295 lin_interpn (n, size, scale, Ni, extrapval, x, v, y, vi); | |
296 | |
297 retval = Vi; | |
298 | |
299 return retval; | |
300 } |