1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 2007-2021 The Octave Project Developers
4 //
5 // See the file COPYRIGHT.md in the top-level directory of this
6 // distribution or <https://octave.org/copyright/>.
7 //
8 // This file is part of Octave.
9 //
10 // Octave is free software: you can redistribute it and/or modify it
11 // under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // Octave is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with Octave; see the file COPYING. If not, see
22 // <https://www.gnu.org/licenses/>.
23 //
24 ////////////////////////////////////////////////////////////////////////
25
26 #if defined (HAVE_CONFIG_H)
27 # include "config.h"
28 #endif
29
30 #include "lo-ieee.h"
31 #include "dNDArray.h"
32 #include "oct-locbuf.h"
33
34 #include "defun.h"
35 #include "error.h"
36 #include "ovl.h"
37
38 // equivalent to isvector.m
39
40 template <typename T>
41 bool
isvector(const T & array)42 isvector (const T& array)
43 {
44 const dim_vector dv = array.dims ();
45 return dv.ndims () == 2 && (dv(0) == 1 || dv(1) == 1);
46 }
47
48 // lookup a value in a sorted table (lookup.m)
49 template <typename T>
50 octave_idx_type
lookup(const T * x,octave_idx_type n,T y)51 lookup (const T *x, octave_idx_type n, T y)
52 {
53 octave_idx_type j;
54
55 if (x[0] < x[n-1])
56 {
57 // increasing x
58
59 if (y > x[n-1] || y < x[0])
60 return -1;
61
62 #if defined (EXHAUSTIF)
63 for (j = 0; j < n - 1; j++)
64 {
65 if (x[j] <= y && y <= x[j+1])
66 return j;
67 }
68 #else
69 octave_idx_type j0 = 0;
70 octave_idx_type j1 = n - 1;
71
72 while (true)
73 {
74 j = (j0+j1)/2;
75
76 if (y <= x[j+1])
77 {
78 if (x[j] <= y)
79 return j;
80
81 j1 = j;
82 }
83
84 if (x[j] <= y)
85 j0 = j;
86 }
87 #endif
88 }
89 else
90 {
91 // decreasing x
92 // previous code with x -> -x and y -> -y
93
94 if (y > x[0] || y < x[n-1])
95 return -1;
96
97 #if defined (EXHAUSTIF)
98 for (j = 0; j < n - 1; j++)
99 {
100 if (x[j+1] <= y && y <= x[j])
101 return j;
102 }
103 #else
104 octave_idx_type j0 = 0;
105 octave_idx_type j1 = n - 1;
106
107 while (true)
108 {
109 j = (j0+j1)/2;
110
111 if (y >= x[j+1])
112 {
113 if (x[j] >= y)
114 return j;
115
116 j1 = j;
117 }
118
119 if (x[j] >= y)
120 j0 = j;
121 }
122 #endif
123 }
124 }
125
126 // n-dimensional linear interpolation
127
128 template <typename T>
129 void
lin_interpn(int n,const octave_idx_type * size,const octave_idx_type * scale,octave_idx_type Ni,T extrapval,const T ** x,const T * v,const T ** y,T * vi)130 lin_interpn (int n, const octave_idx_type *size, const octave_idx_type *scale,
131 octave_idx_type Ni, T extrapval, const T **x,
132 const T *v, const T **y, T *vi)
133 {
134 bool out = false;
135 int bit;
136
137 OCTAVE_LOCAL_BUFFER (T, coef, 2*n);
138 OCTAVE_LOCAL_BUFFER (octave_idx_type, index, n);
139
140 // loop over all points
141 for (octave_idx_type m = 0; m < Ni; m++)
142 {
143 // loop over all dimensions
144 for (int i = 0; i < n; i++)
145 {
146 index[i] = lookup (x[i], size[i], y[i][m]);
147 out = index[i] == -1;
148
149 if (out)
150 break;
151 else
152 {
153 octave_idx_type j = index[i];
154 coef[2*i+1] = (y[i][m] - x[i][j])/(x[i][j+1] - x[i][j]);
155 coef[2*i] = 1 - coef[2*i+1];
156 }
157 }
158
159 if (out)
160 vi[m] = extrapval;
161 else
162 {
163 vi[m] = 0;
164
165 // loop over all corners of hypercube (1<<n = 2^n)
166 for (int i = 0; i < (1 << n); i++)
167 {
168 T c = 1;
169 octave_idx_type l = 0;
170
171 // loop over all dimensions
172 for (int j = 0; j < n; j++)
173 {
174 // test if the jth bit in i is set
175 bit = i >> j & 1;
176 l += scale[j] * (index[j] + bit);
177 c *= coef[2*j+bit];
178 }
179
180 vi[m] += c * v[l];
181 }
182 }
183 }
184 }
185
186 template <typename T, typename M>
187 octave_value
lin_interpn(int n,M * X,const M V,M * Y)188 lin_interpn (int n, M *X, const M V, M *Y)
189 {
190 octave_value retval;
191
192 M Vi = M (Y[0].dims ());
193
194 OCTAVE_LOCAL_BUFFER (const T *, y, n);
195 OCTAVE_LOCAL_BUFFER (octave_idx_type, size, n);
196
197 for (int i = 0; i < n; i++)
198 {
199 y[i] = Y[i].data ();
200 size[i] = V.dims ()(i);
201 }
202
203 OCTAVE_LOCAL_BUFFER (const T *, x, n);
204 OCTAVE_LOCAL_BUFFER (octave_idx_type, scale, n);
205
206 const T *v = V.data ();
207 T *vi = Vi.fortran_vec ();
208 octave_idx_type Ni = Vi.numel ();
209
210 T extrapval = octave_NA;
211
212 // offset in memory of each dimension
213
214 scale[0] = 1;
215
216 for (int i = 1; i < n; i++)
217 scale[i] = scale[i-1] * size[i-1];
218
219 // tests if X[0] is a vector, if yes, assume that all elements of X are
220 // in the ndgrid format.
221
222 if (! isvector (X[0]))
223 {
224 for (int i = 0; i < n; i++)
225 {
226 if (X[i].dims () != V.dims ())
227 error ("interpn: incompatible size of argument number %d", i+1);
228
229 M tmp = M (dim_vector (size[i], 1));
230
231 for (octave_idx_type j = 0; j < size[i]; j++)
232 tmp(j) = X[i](scale[i]*j);
233
234 X[i] = tmp;
235 }
236 }
237
238 for (int i = 0; i < n; i++)
239 {
240 if (! isvector (X[i]) && X[i].numel () != size[i])
241 error ("interpn: incompatible size of argument number %d", i+1);
242
243 x[i] = X[i].data ();
244 }
245
246 lin_interpn (n, size, scale, Ni, extrapval, x, v, y, vi);
247
248 retval = Vi;
249
250 return retval;
251 }
252
253 // Perform @var{n}-dimensional interpolation. Each element of then
254 // @var{n}-dimensional array @var{v} represents a value at a location
255 // given by the parameters @var{x1}, @var{x2},...,@var{xn}. The parameters
256 // @var{x1}, @var{x2}, @dots{}, @var{xn} are either @var{n}-dimensional
257 // arrays of the same size as the array @var{v} in the "ndgrid" format
258 // or vectors. The parameters @var{y1}, @var{y2}, @dots{}, @var{yn} are
259 // all @var{n}-dimensional arrays of the same size and represent the
260 // points at which the array @var{vi} is interpolated.
261 //
262 //This function only performs linear interpolation.
263
264 DEFUN (__lin_interpn__, args, ,
265 doc: /* -*- texinfo -*-
266 @deftypefn {} {@var{vi} =} __lin_interpn__ (@var{x1}, @var{x2}, @dots{}, @var{xn}, @var{v}, @var{y1}, @var{y2}, @dots{}, @var{yn})
267 Undocumented internal function.
268 @end deftypefn */)
269 {
270 int nargin = args.length ();
271
272 if (nargin < 2 || nargin % 2 == 0)
273 print_usage ();
274
275 octave_value retval;
276
277 // dimension of the problem
278 int n = (nargin-1)/2;
279
280 if (args(n).is_single_type ())
281 {
282 OCTAVE_LOCAL_BUFFER (FloatNDArray, X, n);
283 OCTAVE_LOCAL_BUFFER (FloatNDArray, Y, n);
284
285 const FloatNDArray V = args(n).float_array_value ();
286
287 for (int i = 0; i < n; i++)
288 {
289 X[i] = args(i).float_array_value ();
290 Y[i] = args(n+i+1).float_array_value ();
291
292 if (Y[0].dims () != Y[i].dims ())
293 error ("interpn: incompatible size of argument number %d", n+i+2);
294 }
295
296 retval = lin_interpn<float, FloatNDArray> (n, X, V, Y);
297 }
298 else
299 {
300 OCTAVE_LOCAL_BUFFER (NDArray, X, n);
301 OCTAVE_LOCAL_BUFFER (NDArray, Y, n);
302
303 const NDArray V = args(n).array_value ();
304
305 for (int i = 0; i < n; i++)
306 {
307 X[i] = args(i).array_value ();
308 Y[i] = args(n+i+1).array_value ();
309
310 if (Y[0].dims () != Y[i].dims ())
311 error ("interpn: incompatible size of argument number %d", n+i+2);
312 }
313
314 retval = lin_interpn<double, NDArray> (n, X, V, Y);
315 }
316
317 return retval;
318 }
319
320 /*
321 ## No test needed for internal helper function.
322 %!assert (1)
323 */
324