sisi4s
Loading...
Searching...
No Matches
TricubicInterpolation.hpp
Go to the documentation of this file.
1// David Eberly, Geometric Tools, Redmond WA 98052
2// Copyright (c) 1998-2017
3// Distributed under the Boost Software License, Version 1.0.
4// http://www.boost.org/LICENSE_1_0.txt
5// http://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
6// File Version: 3.0.0 (2016/06/19)
7
8#ifndef TRICUBIC_INTERPOLATION_H
9#define TRICUBIC_INTERPOLATION_H
10
11// different assert routines in sisi4s
12#include <util/Exception.hpp>
13
14
15// The interpolator is for uniformly spaced(x,y z)-values. The input samples
16// must be stored in lexicographical order to represent f(x,y,z); that is,
17// F[c + xBound*(r + yBound*s)] corresponds to f(x,y,z), where c is the index
18// corresponding to x, r is the index corresponding to y, and s is the index
19// corresponding to z. Exact interpolation is achieved by setting catmullRom
20// to 'true', giving you the Catmull-Rom blending matrix. If a smooth
21// interpolation is desired, set catmullRom to 'false' to obtain B-spline
22// blending.
23
24namespace gte
25{
26
27template <typename Real>
28class IntpTricubic3
29{
30public:
31 // Construction.
32 IntpTricubic3(int xBound, int yBound, int zBound, Real xMin,
33 Real xSpacing, Real yMin, Real ySpacing, Real zMin, Real zSpacing,
34 Real const* F, bool catmullRom);
35
36 // Member access.
37 inline int GetXBound() const;
38 inline int GetYBound() const;
39 inline int GetZBound() const;
40 inline int GetQuantity() const;
41 inline Real const* GetF() const;
42 inline Real GetXMin() const;
43 inline Real GetXMax() const;
44 inline Real GetXSpacing() const;
45 inline Real GetYMin() const;
46 inline Real GetYMax() const;
47 inline Real GetYSpacing() const;
48 inline Real GetZMin() const;
49 inline Real GetZMax() const;
50 inline Real GetZSpacing() const;
51
52 // Evaluate the function and its derivatives. The functions clamp the
53 // inputs to xmin <= x <= xmax, ymin <= y <= ymax, and zmin <= z <= zmax.
54 // The first operator is for function evaluation. The second operator is
55 // for function or derivative evaluations. The xOrder argument is the
56 // order of the x-derivative, the yOrder argument is the order of the
57 // y-derivative, and the zOrder argument is the order of the z-derivative.
58 // All orders are zero to get the function value itself.
59 Real operator()(Real x, Real y, Real z) const;
60 Real operator()(int xOrder, int yOrder, int zOrder, Real x, Real y,
61 Real z) const;
62
63private:
64 int mXBound, mYBound, mZBound, mQuantity;
65 Real mXMin, mXMax, mXSpacing, mInvXSpacing;
66 Real mYMin, mYMax, mYSpacing, mInvYSpacing;
67 Real mZMin, mZMax, mZSpacing, mInvZSpacing;
68 Real const* mF;
69 Real mBlend[4][4];
70};
71
72
73template <typename Real>
74IntpTricubic3<Real>::IntpTricubic3(int xBound, int yBound, int zBound,
75 Real xMin, Real xSpacing, Real yMin, Real ySpacing, Real zMin,
76 Real zSpacing, Real const* F, bool catmullRom)
77 :
78 mXBound(xBound),
79 mYBound(yBound),
80 mZBound(zBound),
81 mQuantity(xBound * yBound * zBound),
82 mXMin(xMin),
83 mXSpacing(xSpacing),
84 mYMin(yMin),
85 mYSpacing(ySpacing),
86 mZMin(zMin),
87 mZSpacing(zSpacing),
88 mF(F)
89{
90 // At least a 4x4x4 block of data points are needed to construct the
91 // tricubic interpolation.
92 Assert(mXBound >= 4 && mYBound >= 4 && mZBound >= 4 && mF,
93 "Invalid input Bound.");
94 Assert(mXSpacing > (Real)0 && mYSpacing > (Real)0 &&
95 mZSpacing > (Real)0, "Invalid input Spacing.");
96
97 mXMax = mXMin + mXSpacing * static_cast<Real>(mXBound - 1);
98 mInvXSpacing = ((Real)1) / mXSpacing;
99 mYMax = mYMin + mYSpacing * static_cast<Real>(mYBound - 1);
100 mInvYSpacing = ((Real)1) / mYSpacing;
101 mZMax = mZMin + mZSpacing * static_cast<Real>(mZBound - 1);
102 mInvZSpacing = ((Real)1) / mZSpacing;
103
104 if (catmullRom)
105 {
106 mBlend[0][0] = (Real)0;
107 mBlend[0][1] = (Real)-0.5;
108 mBlend[0][2] = (Real)1;
109 mBlend[0][3] = (Real)-0.5;
110 mBlend[1][0] = (Real)1;
111 mBlend[1][1] = (Real)0;
112 mBlend[1][2] = (Real)-2.5;
113 mBlend[1][3] = (Real)1.5;
114 mBlend[2][0] = (Real)0;
115 mBlend[2][1] = (Real)0.5;
116 mBlend[2][2] = (Real)2;
117 mBlend[2][3] = (Real)-1.5;
118 mBlend[3][0] = (Real)0;
119 mBlend[3][1] = (Real)0;
120 mBlend[3][2] = (Real)-0.5;
121 mBlend[3][3] = (Real)0.5;
122 }
123 else
124 {
125 mBlend[0][0] = (Real)1 / (Real)6;
126 mBlend[0][1] = (Real)-3 / (Real)6;
127 mBlend[0][2] = (Real)3 / (Real)6;
128 mBlend[0][3] = (Real)-1 / (Real)6;;
129 mBlend[1][0] = (Real)4 / (Real)6;
130 mBlend[1][1] = (Real)0 / (Real)6;
131 mBlend[1][2] = (Real)-6 / (Real)6;
132 mBlend[1][3] = (Real)3 / (Real)6;
133 mBlend[2][0] = (Real)1 / (Real)6;
134 mBlend[2][1] = (Real)3 / (Real)6;
135 mBlend[2][2] = (Real)3 / (Real)6;
136 mBlend[2][3] = (Real)-3 / (Real)6;
137 mBlend[3][0] = (Real)0 / (Real)6;
138 mBlend[3][1] = (Real)0 / (Real)6;
139 mBlend[3][2] = (Real)0 / (Real)6;
140 mBlend[3][3] = (Real)1 / (Real)6;
141 }
142}
143
144template <typename Real> inline
145int IntpTricubic3<Real>::GetXBound() const
146{
147 return mXBound;
148}
149
150template <typename Real> inline
151int IntpTricubic3<Real>::GetYBound() const
152{
153 return mYBound;
154}
155
156template <typename Real> inline
157int IntpTricubic3<Real>::GetZBound() const
158{
159 return mZBound;
160}
161
162template <typename Real> inline
163int IntpTricubic3<Real>::GetQuantity() const
164{
165 return mQuantity;
166}
167
168template <typename Real> inline
169Real const* IntpTricubic3<Real>::GetF() const
170{
171 return mF;
172}
173
174template <typename Real> inline
175Real IntpTricubic3<Real>::GetXMin() const
176{
177 return mXMin;
178}
179
180template <typename Real> inline
181Real IntpTricubic3<Real>::GetXMax() const
182{
183 return mXMax;
184}
185
186template <typename Real> inline
187Real IntpTricubic3<Real>::GetXSpacing() const
188{
189 return mXSpacing;
190}
191
192template <typename Real> inline
193Real IntpTricubic3<Real>::GetYMin() const
194{
195 return mYMin;
196}
197
198template <typename Real> inline
199Real IntpTricubic3<Real>::GetYMax() const
200{
201 return mYMax;
202}
203
204template <typename Real> inline
205Real IntpTricubic3<Real>::GetYSpacing() const
206{
207 return mYSpacing;
208}
209
210template <typename Real> inline
211Real IntpTricubic3<Real>::GetZMin() const
212{
213 return mZMin;
214}
215
216template <typename Real> inline
217Real IntpTricubic3<Real>::GetZMax() const
218{
219 return mZMax;
220}
221
222template <typename Real> inline
223Real IntpTricubic3<Real>::GetZSpacing() const
224{
225 return mZSpacing;
226}
227
228template <typename Real>
229Real IntpTricubic3<Real>::operator()(Real x, Real y, Real z) const
230{
231 // Compute x-index and clamp to image.
232 Real xIndex = (x - mXMin) * mInvXSpacing;
233 int ix = static_cast<int>(xIndex);
234 if (ix < 0)
235 {
236 ix = 0;
237 }
238 else if (ix >= mXBound)
239 {
240 ix = mXBound - 1;
241 }
242
243 // Compute y-index and clamp to image.
244 Real yIndex = (y - mYMin) * mInvYSpacing;
245 int iy = static_cast<int>(yIndex);
246 if (iy < 0)
247 {
248 iy = 0;
249 }
250 else if (iy >= mYBound)
251 {
252 iy = mYBound - 1;
253 }
254
255 // Compute z-index and clamp to image.
256 Real zIndex = (z - mZMin) * mInvZSpacing;
257 int iz = static_cast<int>(zIndex);
258 if (iz < 0)
259 {
260 iz = 0;
261 }
262 else if (iz >= mZBound)
263 {
264 iz = mZBound - 1;
265 }
266
267 Real U[4];
268 U[0] = (Real)1;
269 U[1] = xIndex - ix;
270 U[2] = U[1] * U[1];
271 U[3] = U[1] * U[2];
272
273 Real V[4];
274 V[0] = (Real)1;
275 V[1] = yIndex - iy;
276 V[2] = V[1] * V[1];
277 V[3] = V[1] * V[2];
278
279 Real W[4];
280 W[0] = (Real)1;
281 W[1] = zIndex - iz;
282 W[2] = W[1] * W[1];
283 W[3] = W[1] * W[2];
284
285 // Compute P = M*U, Q = M*V, R = M*W.
286 Real P[4], Q[4], R[4];
287 for (int row = 0; row < 4; ++row)
288 {
289 P[row] = (Real)0;
290 Q[row] = (Real)0;
291 R[row] = (Real)0;
292 for (int col = 0; col < 4; ++col)
293 {
294 P[row] += mBlend[row][col] * U[col];
295 Q[row] += mBlend[row][col] * V[col];
296 R[row] += mBlend[row][col] * W[col];
297 }
298 }
299
300 // Compute the tensor product (M*U)(M*V)(M*W)*D where D is the 4x4x4
301 // subimage containing (x,y,z).
302 --ix;
303 --iy;
304 --iz;
305 Real result = (Real)0;
306 for (int slice = 0; slice < 4; ++slice)
307 {
308 int zClamp = iz + slice;
309 if (zClamp < 0)
310 {
311 zClamp = 0;
312 }
313 else if (zClamp > mZBound - 1)
314 {
315 zClamp = mZBound - 1;
316 }
317
318 for (int row = 0; row < 4; ++row)
319 {
320 int yClamp = iy + row;
321 if (yClamp < 0)
322 {
323 yClamp = 0;
324 }
325 else if (yClamp > mYBound - 1)
326 {
327 yClamp = mYBound - 1;
328 }
329
330 for (int col = 0; col < 4; ++col)
331 {
332 int xClamp = ix + col;
333 if (xClamp < 0)
334 {
335 xClamp = 0;
336 }
337 else if (xClamp > mXBound - 1)
338 {
339 xClamp = mXBound - 1;
340 }
341
342 result += P[col] * Q[row] * R[slice] *
343 mF[xClamp + mXBound * (yClamp + mYBound * zClamp)];
344 }
345 }
346 }
347
348 return result;
349}
350
351template <typename Real>
352Real IntpTricubic3<Real>::operator()(int xOrder, int yOrder, int zOrder,
353 Real x, Real y, Real z) const
354{
355 // Compute x-index and clamp to image.
356 Real xIndex = (x - mXMin) * mInvXSpacing;
357 int ix = static_cast<int>(xIndex);
358 if (ix < 0)
359 {
360 ix = 0;
361 }
362 else if (ix >= mXBound)
363 {
364 ix = mXBound - 1;
365 }
366
367 // Compute y-index and clamp to image.
368 Real yIndex = (y - mYMin) * mInvYSpacing;
369 int iy = static_cast<int>(yIndex);
370 if (iy < 0)
371 {
372 iy = 0;
373 }
374 else if (iy >= mYBound)
375 {
376 iy = mYBound - 1;
377 }
378
379 // Compute z-index and clamp to image.
380 Real zIndex = (z - mZMin) * mInvZSpacing;
381 int iz = static_cast<int>(zIndex);
382 if (iz < 0)
383 {
384 iz = 0;
385 }
386 else if (iz >= mZBound)
387 {
388 iz = mZBound - 1;
389 }
390
391 Real U[4], dx, xMult;
392 switch (xOrder)
393 {
394 case 0:
395 dx = xIndex - ix;
396 U[0] = (Real)1;
397 U[1] = dx;
398 U[2] = dx * U[1];
399 U[3] = dx * U[2];
400 xMult = (Real)1;
401 break;
402 case 1:
403 dx = xIndex - ix;
404 U[0] = (Real)0;
405 U[1] = (Real)1;
406 U[2] = ((Real)2) * dx;
407 U[3] = ((Real)3) * dx * dx;
408 xMult = mInvXSpacing;
409 break;
410 case 2:
411 dx = xIndex - ix;
412 U[0] = (Real)0;
413 U[1] = (Real)0;
414 U[2] = (Real)2;
415 U[3] = ((Real)6) * dx;
416 xMult = mInvXSpacing * mInvXSpacing;
417 break;
418 case 3:
419 U[0] = (Real)0;
420 U[1] = (Real)0;
421 U[2] = (Real)0;
422 U[3] = (Real)6;
423 xMult = mInvXSpacing * mInvXSpacing * mInvXSpacing;
424 break;
425 default:
426 return (Real)0;
427 }
428
429 Real V[4], dy, yMult;
430 switch (yOrder)
431 {
432 case 0:
433 dy = yIndex - iy;
434 V[0] = (Real)1;
435 V[1] = dy;
436 V[2] = dy * V[1];
437 V[3] = dy * V[2];
438 yMult = (Real)1;
439 break;
440 case 1:
441 dy = yIndex - iy;
442 V[0] = (Real)0;
443 V[1] = (Real)1;
444 V[2] = ((Real)2) * dy;
445 V[3] = ((Real)3) * dy * dy;
446 yMult = mInvYSpacing;
447 break;
448 case 2:
449 dy = yIndex - iy;
450 V[0] = (Real)0;
451 V[1] = (Real)0;
452 V[2] = (Real)2;
453 V[3] = ((Real)6) * dy;
454 yMult = mInvYSpacing * mInvYSpacing;
455 break;
456 case 3:
457 V[0] = (Real)0;
458 V[1] = (Real)0;
459 V[2] = (Real)0;
460 V[3] = (Real)6;
461 yMult = mInvYSpacing * mInvYSpacing * mInvYSpacing;
462 break;
463 default:
464 return (Real)0;
465 }
466
467 Real W[4], dz, zMult;
468 switch (zOrder)
469 {
470 case 0:
471 dz = zIndex - iz;
472 W[0] = (Real)1;
473 W[1] = dz;
474 W[2] = dz * W[1];
475 W[3] = dz * W[2];
476 zMult = (Real)1;
477 break;
478 case 1:
479 dz = zIndex - iz;
480 W[0] = (Real)0;
481 W[1] = (Real)1;
482 W[2] = ((Real)2) * dz;
483 W[3] = ((Real)3) * dz * dz;
484 zMult = mInvZSpacing;
485 break;
486 case 2:
487 dz = zIndex - iz;
488 W[0] = (Real)0;
489 W[1] = (Real)0;
490 W[2] = (Real)2;
491 W[3] = ((Real)6) * dz;
492 zMult = mInvZSpacing * mInvZSpacing;
493 break;
494 case 3:
495 W[0] = (Real)0;
496 W[1] = (Real)0;
497 W[2] = (Real)0;
498 W[3] = (Real)6;
499 zMult = mInvZSpacing * mInvZSpacing * mInvZSpacing;
500 break;
501 default:
502 return (Real)0;
503 }
504
505 // Compute P = M*U, Q = M*V, and R = M*W.
506 Real P[4], Q[4], R[4];
507 for (int row = 0; row < 4; ++row)
508 {
509 P[row] = (Real)0;
510 Q[row] = (Real)0;
511 R[row] = (Real)0;
512 for (int col = 0; col < 4; ++col)
513 {
514 P[row] += mBlend[row][col] * U[col];
515 Q[row] += mBlend[row][col] * V[col];
516 R[row] += mBlend[row][col] * W[col];
517 }
518 }
519
520 // Compute the tensor product (M*U)(M*V)(M*W)*D where D is the 4x4x4
521 // subimage containing (x,y,z).
522 --ix;
523 --iy;
524 --iz;
525 Real result = (Real)0;
526 for (int slice = 0; slice < 4; ++slice)
527 {
528 int zClamp = iz + slice;
529 if (zClamp < 0)
530 {
531 zClamp = 0;
532 }
533 else if (zClamp > mZBound - 1)
534 {
535 zClamp = mZBound - 1;
536 }
537
538 for (int row = 0; row < 4; ++row)
539 {
540 int yClamp = iy + row;
541 if (yClamp < 0)
542 {
543 yClamp = 0;
544 }
545 else if (yClamp > mYBound - 1)
546 {
547 yClamp = mYBound - 1;
548 }
549
550 for (int col = 0; col < 4; ++col)
551 {
552 int xClamp = ix + col;
553 if (xClamp < 0)
554 {
555 xClamp = 0;
556 }
557 else if (xClamp > mXBound - 1)
558 {
559 xClamp = mXBound - 1;
560 }
561
562 result += P[col] * Q[row] * R[slice] *
563 mF[xClamp + mXBound * (yClamp + mYBound * zClamp)];
564 }
565 }
566 }
567 result *= xMult * yMult * zMult;
568
569 return result;
570}
571
572
573}
574
575#endif
576
Assert
#define Assert(condition, message)
Definition Exception.hpp:10
gte::IntpTricubic3
Definition TricubicInterpolation.hpp:29
gte::IntpTricubic3::GetF
Real const * GetF() const
Definition TricubicInterpolation.hpp:169
gte::IntpTricubic3::GetXMin
Real GetXMin() const
Definition TricubicInterpolation.hpp:175
gte::IntpTricubic3::GetYMax
Real GetYMax() const
Definition TricubicInterpolation.hpp:199
gte::IntpTricubic3::GetXBound
int GetXBound() const
Definition TricubicInterpolation.hpp:145
gte::IntpTricubic3::GetYMin
Real GetYMin() const
Definition TricubicInterpolation.hpp:193
gte::IntpTricubic3::GetYSpacing
Real GetYSpacing() const
Definition TricubicInterpolation.hpp:205
gte::IntpTricubic3::GetZMax
Real GetZMax() const
Definition TricubicInterpolation.hpp:217
gte::IntpTricubic3::GetZBound
int GetZBound() const
Definition TricubicInterpolation.hpp:157
gte::IntpTricubic3::operator()
Real operator()(Real x, Real y, Real z) const
Definition TricubicInterpolation.hpp:229
gte::IntpTricubic3::GetXSpacing
Real GetXSpacing() const
Definition TricubicInterpolation.hpp:187
gte::IntpTricubic3::GetQuantity
int GetQuantity() const
Definition TricubicInterpolation.hpp:163
gte::IntpTricubic3::GetXMax
Real GetXMax() const
Definition TricubicInterpolation.hpp:181
gte::IntpTricubic3::GetYBound
int GetYBound() const
Definition TricubicInterpolation.hpp:151
gte::IntpTricubic3::GetZMin
Real GetZMin() const
Definition TricubicInterpolation.hpp:211
gte::IntpTricubic3::GetZSpacing
Real GetZSpacing() const
Definition TricubicInterpolation.hpp:223
gte::IntpTricubic3::IntpTricubic3
IntpTricubic3(int xBound, int yBound, int zBound, Real xMin, Real xSpacing, Real yMin, Real ySpacing, Real zMin, Real zSpacing, Real const *F, bool catmullRom)
Definition TricubicInterpolation.hpp:74
gte
Definition TricubicInterpolation.hpp:25