GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
__lin_interpn__.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2007-2013 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 3 of the License, or (at your
10 option) any 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, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26 
27 #include "lo-ieee.h"
28 #include "dNDArray.h"
29 #include "oct-locbuf.h"
30 
31 #include "defun.h"
32 #include "error.h"
33 #include "oct-obj.h"
34 
35 // equivalent to isvector.m
36 
37 template <class T>
38 bool
39 isvector (const T& array)
40 {
41  const dim_vector dv = array.dims ();
42  return dv.length () == 2 && (dv(0) == 1 || dv(1) == 1);
43 }
44 
45 // lookup a value in a sorted table (lookup.m)
46 template <class T>
48 lookup (const T *x, octave_idx_type n, T y)
49 {
51 
52  if (x[0] < x[n-1])
53  {
54  // increasing x
55 
56  if (y > x[n-1] || y < x[0])
57  return -1;
58 
59 #ifdef EXHAUSTIF
60  for (j = 0; j < n - 1; j++)
61  {
62  if (x[j] <= y && y <= x[j+1])
63  return j;
64  }
65 #else
66  octave_idx_type j0 = 0;
67  octave_idx_type j1 = n - 1;
68 
69  while (true)
70  {
71  j = (j0+j1)/2;
72 
73  if (y <= x[j+1])
74  {
75  if (x[j] <= y)
76  return j;
77 
78  j1 = j;
79  }
80 
81  if (x[j] <= y)
82  j0 = j;
83  }
84 #endif
85  }
86  else
87  {
88  // decreasing x
89  // previous code with x -> -x and y -> -y
90 
91  if (y > x[0] || y < x[n-1])
92  return -1;
93 
94 #ifdef EXHAUSTIF
95  for (j = 0; j < n - 1; j++)
96  {
97  if (x[j+1] <= y && y <= x[j])
98  return j;
99  }
100 #else
101  octave_idx_type j0 = 0;
102  octave_idx_type j1 = n - 1;
103 
104  while (true)
105  {
106  j = (j0+j1)/2;
107 
108  if (y >= x[j+1])
109  {
110  if (x[j] >= y)
111  return j;
112 
113  j1 = j;
114  }
115 
116  if (x[j] >= y)
117  j0 = j;
118  }
119 #endif
120  }
121 }
122 
123 // n-dimensional linear interpolation
124 
125 template <class T>
126 void
128  octave_idx_type Ni, T extrapval, const T **x,
129  const T *v, const T **y, T *vi)
130 {
131  bool out = false;
132  int bit;
133 
134  OCTAVE_LOCAL_BUFFER (T, coef, 2*n);
136 
137  // loop over all points
138  for (octave_idx_type m = 0; m < Ni; m++)
139  {
140  // loop over all dimensions
141  for (int i = 0; i < n; i++)
142  {
143  index[i] = lookup (x[i], size[i], y[i][m]);
144  out = index[i] == -1;
145 
146  if (out)
147  break;
148  else
149  {
150  octave_idx_type j = index[i];
151  coef[2*i+1] = (y[i][m] - x[i][j])/(x[i][j+1] - x[i][j]);
152  coef[2*i] = 1 - coef[2*i+1];
153  }
154  }
155 
156 
157  if (out)
158  vi[m] = extrapval;
159  else
160  {
161  vi[m] = 0;
162 
163  // loop over all corners of hypercube (1<<n = 2^n)
164  for (int i = 0; i < (1 << n); i++)
165  {
166  T c = 1;
167  octave_idx_type l = 0;
168 
169  // loop over all dimensions
170  for (int j = 0; j < n; j++)
171  {
172  // test if the jth bit in i is set
173  bit = i >> j & 1;
174  l += scale[j] * (index[j] + bit);
175  c *= coef[2*j+bit];
176  }
177 
178  vi[m] += c * v[l];
179  }
180  }
181  }
182 }
183 
184 template <class T, class M>
186 lin_interpn (int n, M *X, const M V, M *Y)
187 {
188  octave_value retval;
189 
190  M Vi = M (Y[0].dims ());
191 
192  OCTAVE_LOCAL_BUFFER (const T *, y, n);
194 
195  for (int i = 0; i < n; i++)
196  {
197  y[i] = Y[i].data ();
198  size[i] = V.dims ()(i);
199  }
200 
201  OCTAVE_LOCAL_BUFFER (const T *, x, n);
203 
204  const T *v = V.data ();
205  T *vi = Vi.fortran_vec ();
206  octave_idx_type Ni = Vi.numel ();
207 
208  T extrapval = octave_NA;
209 
210  // offset in memory of each dimension
211 
212  scale[0] = 1;
213 
214  for (int i = 1; i < n; i++)
215  scale[i] = scale[i-1] * size[i-1];
216 
217  // tests if X[0] is a vector, if yes, assume that all elements of X are
218  // in the ndgrid format.
219 
220  if (! isvector (X[0]))
221  {
222  for (int i = 0; i < n; i++)
223  {
224  if (X[i].dims () != V.dims ())
225  {
226  error ("interpn: incompatible size of argument number %d", i+1);
227  return retval;
228  }
229  else
230  {
231  M tmp = M (dim_vector (size[i], 1));
232 
233  for (octave_idx_type j = 0; j < size[i]; j++)
234  tmp(j) = X[i](scale[i]*j);
235 
236  X[i] = tmp;
237  }
238  }
239  }
240 
241  for (int i = 0; i < n; i++)
242  {
243  if (! isvector (X[i]) && X[i].numel () != size[i])
244  {
245  error ("interpn: incompatible size of argument number %d", i+1);
246  return retval;
247  }
248  else
249  x[i] = X[i].data ();
250  }
251 
252  lin_interpn (n, size, scale, Ni, extrapval, x, v, y, vi);
253 
254  retval = Vi;
255 
256  return retval;
257 }
258 
259 // Perform @var{n}-dimensional interpolation. Each element of then
260 // @var{n}-dimensional array @var{v} represents a value at a location
261 // given by the parameters @var{x1}, @var{x2},...,@var{xn}. The parameters
262 // @var{x1}, @var{x2}, @dots{}, @var{xn} are either @var{n}-dimensional
263 // arrays of the same size as the array @var{v} in the \"ndgrid\" format
264 // or vectors. The parameters @var{y1}, @var{y2}, @dots{}, @var{yn} are
265 // all @var{n}-dimensional arrays of the same size and represent the
266 // points at which the array @var{vi} is interpolated.
267 //
268 //This function only performs linear interpolation.
269 
270 DEFUN (__lin_interpn__, args, ,
271  "-*- texinfo -*-\n\
272 @deftypefn {Built-in Function} {@var{vi} =} __lin_interpn__ (@var{x1}, @var{x2}, @dots{}, @var{xn}, @var{v}, @var{y1}, @var{y2}, @dots{}, @var{yn})\n\
273 Undocumented internal function.\n\
274 @end deftypefn")
275 {
276  octave_value retval;
277 
278  int nargin = args.length ();
279 
280  if (nargin < 2 || nargin % 2 == 0)
281  {
282  print_usage ();
283  return retval;
284  }
285 
286  // dimension of the problem
287  int n = (nargin-1)/2;
288 
289  if (args(n).is_single_type ())
290  {
293 
294  const FloatNDArray V = args(n).float_array_value ();
295 
296  if (error_state)
297  {
298  print_usage ();
299  return retval;
300  }
301 
302  for (int i = 0; i < n; i++)
303  {
304  X[i] = args(i).float_array_value ();
305  Y[i] = args(n+i+1).float_array_value ();
306 
307  if (error_state)
308  {
309  print_usage ();
310  return retval;
311  }
312 
313  if (Y[0].dims () != Y[i].dims ())
314  {
315  error ("interpn: incompatible size of argument number %d", n+i+2);
316  return retval;
317  }
318  }
319 
320  retval = lin_interpn<float, FloatNDArray> (n, X, V, Y);
321  }
322  else
323  {
326 
327  const NDArray V = args(n).array_value ();
328 
329  if (error_state)
330  {
331  print_usage ();
332  return retval;
333  }
334 
335  for (int i = 0; i < n; i++)
336  {
337  X[i] = args(i).array_value ();
338  Y[i] = args(n+i+1).array_value ();
339 
340  if (error_state)
341  {
342  print_usage ();
343  return retval;
344  }
345 
346  if (Y[0].dims () != Y[i].dims ())
347  {
348  error ("interpn: incompatible size of argument number %d", n+i+2);
349  return retval;
350  }
351  }
352 
353  retval = lin_interpn<double, NDArray> (n, X, V, Y);
354  }
355 
356  return retval;
357 }
358 
359 /*
360 ## No test needed for internal helper function.
361 %!assert (1)
362 */