Line data Source code
1 : #ifndef CRPROPA_GRID_H
2 : #define CRPROPA_GRID_H
3 :
4 : #include "crpropa/Referenced.h"
5 : #include "crpropa/Vector3.h"
6 :
7 : #include "kiss/string.h"
8 : #include "kiss/logger.h"
9 :
10 : #include <vector>
11 : #include <type_traits>
12 : #include <string>
13 : #include <sstream>
14 : #if HAVE_SIMD
15 : #include <immintrin.h>
16 : #include <smmintrin.h>
17 : #endif // HAVE_SIMD
18 :
19 : namespace crpropa {
20 :
21 : /** If set to TRILINEAR, use trilinear interpolation (standard)
22 : If set to TRICUBIC, use tricubic interpolation instead of trilinear interpolation
23 : If set to NEAREST_NEIGHBOUR , use nearest neighbour interpolation instead of trilinear interpolation */
24 : enum interpolationType {
25 : TRILINEAR = 0,
26 : TRICUBIC,
27 : NEAREST_NEIGHBOUR
28 : };
29 :
30 : /** Lower and upper neighbour in a periodically continued unit grid */
31 : inline void periodicClamp(double x, int n, int &lo, int &hi) {
32 100051 : lo = ((int(floor(x)) % (n)) + (n)) % (n);
33 10 : hi = (lo + 1) % (n);
34 : }
35 :
36 : /** grid index in a reflective continued unit grid */
37 : inline int reflectiveBoundary(int index, int n) {
38 1963 : while ((index < -0.5) or (index > (n-0.5)))
39 806 : index = 2 * n * (index > (n-0.5)) - index-1;
40 : return index;
41 : }
42 :
43 : /** grid index in a periodically continued unit grid */
44 : inline int periodicBoundary(int index, int n) {
45 768 : return ((index % (n)) + (n)) % (n);
46 : }
47 :
48 : /** Lower and upper neighbour in a reflectively repeated unit grid */
49 20 : inline void reflectiveClamp(double x, int n, int &lo, int &hi, double &res) {
50 36 : while ((x < -0.5) or (x > (n-0.5)))
51 16 : x = 2 * n * (x > (n-0.5)) -x-1;
52 20 : res = x;
53 20 : lo = floor(x);
54 20 : hi = lo + (lo < n-1);
55 20 : if (x<0) {
56 10 : lo=0;
57 10 : hi=0;
58 : }
59 20 : }
60 :
61 : /** Symmetrical round */
62 51 : inline double round(double r) {
63 51 : return (r > 0.0) ? floor(r + 0.5) : ceil(r - 0.5);
64 : }
65 :
66 : /**
67 : * \addtogroup Core
68 : * @{
69 : */
70 :
71 : /**
72 : @class GridProperties
73 : @brief Combines parameters that uniquely define Grid class
74 : */
75 : class GridProperties: public Referenced {
76 : public:
77 : size_t Nx, Ny, Nz; // Number of grid points
78 : Vector3d origin; // Position of the lower left front corner of the volume
79 : Vector3d spacing; // Spacing vector between gridpoints
80 : bool reflective; // using reflective repetition of the grid instead of periodic
81 : interpolationType ipol; // Interpolation type used between grid points
82 : bool clipVolume; // Set grid values to 0 outside the volume if true
83 :
84 : /** Constructor for cubic grid
85 : @param origin Position of the lower left front corner of the volume
86 : @param N Number of grid points in one direction
87 : @param spacing Spacing between grid points
88 : */
89 6 : GridProperties(Vector3d origin, size_t N, double spacing) :
90 6 : origin(origin), Nx(N), Ny(N), Nz(N), spacing(Vector3d(spacing)), reflective(false), ipol(TRILINEAR), clipVolume(false) {
91 : }
92 :
93 : /** Constructor for non-cubic grid
94 : @param origin Position of the lower left front corner of the volume
95 : @param Nx Number of grid points in x-direction
96 : @param Ny Number of grid points in y-direction
97 : @param Nz Number of grid points in z-direction
98 : @param spacing Spacing between grid points
99 : */
100 2 : GridProperties(Vector3d origin, size_t Nx, size_t Ny, size_t Nz, double spacing) :
101 2 : origin(origin), Nx(Nx), Ny(Ny), Nz(Nz), spacing(Vector3d(spacing)), reflective(false), ipol(TRILINEAR), clipVolume(false) {
102 : }
103 :
104 : /** Constructor for non-cubic grid with spacing vector
105 : @param origin Position of the lower left front corner of the volume
106 : @param Nx Number of grid points in x-direction
107 : @param Ny Number of grid points in y-direction
108 : @param Nz Number of grid points in z-direction
109 : @param spacing Spacing vector between grid points
110 : */
111 1 : GridProperties(Vector3d origin, size_t Nx, size_t Ny, size_t Nz, Vector3d spacing) :
112 1 : origin(origin), Nx(Nx), Ny(Ny), Nz(Nz), spacing(spacing), reflective(false), ipol(TRILINEAR), clipVolume(false) {
113 : }
114 :
115 2 : virtual ~GridProperties() {
116 4 : }
117 :
118 : /** If True, the repetition of the grid is refletive instead of periodic. */
119 : void setReflective(bool b) {
120 0 : reflective = b;
121 : }
122 :
123 : /** set the type of interpolation between grid points.
124 : * @param i: interpolationType (TRILINEAR, TRICUBIC, NEAREST_NEIGHBOUR) */
125 : void setInterpolationType(interpolationType i) {
126 2 : ipol = i;
127 2 : }
128 :
129 : /** If True, the grid is set to zero outside of the volume. */
130 : void setClipVolume(bool b) {
131 2 : clipVolume = b;
132 : }
133 :
134 : /** show all GridProperty parameters
135 : * @param unit unit for the lengthscale (origin, spacing). Default is 1 = SI units
136 : */
137 0 : std::string getDescription(double unit = 1) const {
138 0 : std::stringstream ss;
139 0 : ss << "GridProperties:\torigin: " << origin / unit
140 0 : << "\t" << "Nx: " << Nx << " Ny: " << Ny << " Nz: " << Nz
141 0 : << "\t" << "spacing: " << spacing / unit
142 0 : << "\t" << "refletive: " << reflective
143 0 : << "\t" << "interpolation: " << ipol
144 0 : << "\n";
145 0 : return ss.str();
146 0 : }
147 : };
148 :
149 : /**
150 : @class Grid
151 : @brief Template class for fields on a periodic grid with trilinear interpolation
152 :
153 : The grid spacing is constant with diffrent resolution along all three axes.
154 : Values are calculated by trilinear interpolation of the surrounding 8 grid points.
155 : The grid is periodically (default) or reflectively extended.
156 : The grid sample positions are at 1/2 * size/N, 3/2 * size/N ... (2N-1)/2 * size/N.
157 : */
158 : template<typename T>
159 9 : class Grid: public Referenced {
160 : std::vector<T> grid;
161 : size_t Nx, Ny, Nz; /**< Number of grid points */
162 : Vector3d origin; /**< Origin of the volume that is represented by the grid. */
163 : Vector3d gridOrigin; /**< Grid origin */
164 : Vector3d spacing; /**< Distance between grid points, determines the extension of the grid */
165 : bool clipVolume; /**< If set to true, all values outside of the grid will be 0*/
166 : bool reflective; /**< If set to true, the grid is repeated reflectively instead of periodically */
167 : interpolationType ipolType; /**< Type of interpolation between the grid points */
168 :
169 : public:
170 : /** Constructor for cubic grid
171 : @param origin Position of the lower left front corner of the volume
172 : @param N Number of grid points in one direction
173 : @param spacing Spacing between grid points
174 : */
175 11 : Grid(Vector3d origin, size_t N, double spacing) {
176 : setOrigin(origin);
177 11 : setGridSize(N, N, N);
178 : setSpacing(Vector3d(spacing));
179 : setReflective(false);
180 : setClipVolume(false);
181 : setInterpolationType(TRILINEAR);
182 0 : }
183 :
184 : /** Constructor for non-cubic grid
185 : @param origin Position of the lower left front corner of the volume
186 : @param Nx Number of grid points in x-direction
187 : @param Ny Number of grid points in y-direction
188 : @param Nz Number of grid points in z-direction
189 : @param spacing Spacing between grid points
190 : */
191 7 : Grid(Vector3d origin, size_t Nx, size_t Ny, size_t Nz, double spacing) {
192 : setOrigin(origin);
193 7 : setGridSize(Nx, Ny, Nz);
194 : setSpacing(Vector3d(spacing));
195 : setReflective(false);
196 : setClipVolume(false);
197 : setInterpolationType(TRILINEAR);
198 0 : }
199 :
200 : /** Constructor for non-cubic grid with spacing vector
201 : @param origin Position of the lower left front corner of the volume
202 : @param Nx Number of grid points in x-direction
203 : @param Ny Number of grid points in y-direction
204 : @param Nz Number of grid points in z-direction
205 : @param spacing Spacing vector between grid points
206 : */
207 3 : Grid(Vector3d origin, size_t Nx, size_t Ny, size_t Nz, Vector3d spacing) {
208 : setOrigin(origin);
209 3 : setGridSize(Nx, Ny, Nz);
210 : setSpacing(spacing);
211 : setReflective(false);
212 : setClipVolume(false);
213 : setInterpolationType(TRILINEAR);
214 0 : }
215 :
216 : /** Constructor for GridProperties
217 : @param p GridProperties instance
218 : */
219 11 : Grid(const GridProperties &p) :
220 11 : origin(p.origin), spacing(p.spacing), reflective(p.reflective), ipolType(p.ipol) {
221 11 : setGridSize(p.Nx, p.Ny, p.Nz);
222 11 : setClipVolume(p.clipVolume);
223 0 : }
224 :
225 : void setOrigin(Vector3d origin) {
226 : this->origin = origin;
227 : this->gridOrigin = origin + spacing/2;
228 : }
229 :
230 : /** Resize grid, also enlarges the volume as the spacing stays constant */
231 32 : void setGridSize(size_t Nx, size_t Ny, size_t Nz) {
232 32 : this->Nx = Nx;
233 32 : this->Ny = Ny;
234 32 : this->Nz = Nz;
235 32 : grid.resize(Nx * Ny * Nz);
236 : setOrigin(origin);
237 32 : }
238 :
239 : void setSpacing(Vector3d spacing) {
240 : this->spacing = spacing;
241 : setOrigin(origin);
242 : }
243 :
244 : void setReflective(bool b) {
245 8 : reflective = b;
246 : }
247 :
248 : // If set to true, all values outside of the grid will be 0.
249 : void setClipVolume(bool b) {
250 23 : clipVolume = b;
251 : }
252 :
253 : /** Change the interpolation type to the routine specified by the user. Check if this routine is
254 : contained in the enum interpolationType and thus supported by CRPropa.*/
255 0 : void setInterpolationType(interpolationType ipolType) {
256 0 : if (ipolType == TRILINEAR || ipolType == TRICUBIC || ipolType == NEAREST_NEIGHBOUR) {
257 34 : this->ipolType = ipolType;
258 0 : if ((ipolType == TRICUBIC) && (std::is_same<T, Vector3d>::value)) {
259 0 : KISS_LOG_WARNING << "Tricubic interpolation on Grid3d works only with float-precision, doubles will be downcasted";
260 : }
261 : } else {
262 0 : throw std::runtime_error("InterpolationType: unknown interpolation type");
263 : }
264 0 : }
265 :
266 : interpolationType getInterpolationType() {
267 2 : return ipolType;
268 : }
269 :
270 2 : std::string getInterpolationTypeName() {
271 2 : if (ipolType == TRILINEAR)
272 0 : return "TRILINEAR";
273 2 : if (ipolType == TRICUBIC)
274 1 : return "TRICUBIC";
275 1 : if (ipolType == NEAREST_NEIGHBOUR)
276 1 : return "NEAREST_NEIGHBOUR";
277 :
278 0 : return "NOT_UNDERSTOOD";
279 : }
280 :
281 : /** returns the position of the lower left front corner of the volume */
282 : Vector3d getOrigin() const {
283 : return origin;
284 : }
285 :
286 : bool getClipVolume() const {
287 4 : return clipVolume;
288 : }
289 :
290 : size_t getNx() const {
291 764 : return Nx;
292 : }
293 :
294 : size_t getNy() const {
295 41868 : return Ny;
296 : }
297 :
298 : size_t getNz() const {
299 2663965 : return Nz;
300 : }
301 :
302 : /** Calculates the total size of the grid in bytes */
303 : size_t getSizeOf() const {
304 0 : return sizeof(grid) + (sizeof(grid[0]) * grid.size());
305 : }
306 :
307 : Vector3d getSpacing() const {
308 : return spacing;
309 : }
310 :
311 : bool isReflective() const {
312 4 : return reflective;
313 : }
314 :
315 : /** Choose the interpolation algorithm based on the set interpolation type.
316 : By default this it the trilinear interpolation. The user can change the
317 : routine with the setInterpolationType function.*/
318 100091 : T interpolate(const Vector3d &position) {
319 : // check for volume
320 100091 : if (clipVolume) {
321 5 : Vector3d edge = origin + Vector3d(Nx, Ny, Nz) * spacing;
322 5 : bool isInVolume = (position.x >= origin.x) && (position.x <= edge.x);
323 5 : isInVolume &= (position.y >= origin.y) && (position.y <= edge.y);
324 5 : isInVolume &= (position.z >= origin.z) && (position.z <= edge.z);
325 5 : if (!isInVolume)
326 : return T(0.);
327 : }
328 :
329 100086 : if (ipolType == TRICUBIC)
330 18 : return tricubicInterpolate(T(), position);
331 100068 : else if (ipolType == NEAREST_NEIGHBOUR)
332 13 : return closestValue(position);
333 : else
334 100055 : return trilinearInterpolate(position);
335 : }
336 :
337 : /** Inspector & Mutator */
338 : T &get(size_t ix, size_t iy, size_t iz) {
339 8391214 : return grid[ix * Ny * Nz + iy * Nz + iz];
340 : }
341 :
342 : /** Inspector */
343 : const T &get(size_t ix, size_t iy, size_t iz) const {
344 100072 : return grid[ix * Ny * Nz + iy * Nz + iz];
345 : }
346 :
347 : const T &periodicGet(size_t ix, size_t iy, size_t iz) const {
348 192 : ix = periodicBoundary(ix, Nx);
349 192 : iy = periodicBoundary(iy, Ny);
350 192 : iz = periodicBoundary(iz, Nz);
351 192 : return grid[ix * Ny * Nz + iy * Nz + iz];
352 : }
353 :
354 0 : const T &reflectiveGet(size_t ix, size_t iy, size_t iz) const {
355 0 : ix = reflectiveBoundary(ix, Nx);
356 0 : iy = reflectiveBoundary(iy, Ny);
357 0 : iz = reflectiveBoundary(iz, Nz);
358 0 : return grid[ix * Ny * Nz + iy * Nz + iz];
359 : }
360 :
361 : T getValue(size_t ix, size_t iy, size_t iz) {
362 0 : return grid[ix * Ny * Nz + iy * Nz + iz];
363 : }
364 :
365 : void setValue(size_t ix, size_t iy, size_t iz, T value) {
366 0 : grid[ix * Ny * Nz + iy * Nz + iz] = value;
367 : }
368 :
369 : /** Return a reference to the grid values */
370 : std::vector<T> &getGrid() {
371 10102 : return grid;
372 : }
373 :
374 : /** Position of the grid point of a given index */
375 10103 : Vector3d positionFromIndex(int index) const {
376 10103 : int ix = index / (Ny * Nz);
377 10103 : int iy = (index / Nz) % Ny;
378 10103 : int iz = index % Nz;
379 10103 : return Vector3d(ix, iy, iz) * spacing + gridOrigin;
380 : }
381 :
382 : /** Value of a grid point that is closest to a given position / nearest neighbour interpolation */
383 17 : T closestValue(const Vector3d &position) const {
384 : Vector3d r = (position - gridOrigin) / spacing;
385 : int ix, iy, iz;
386 17 : if (reflective) {
387 6 : ix = round(r.x);
388 6 : iy = round(r.y);
389 6 : iz = round(r.z);
390 9 : while ((ix < -0.5) or (ix > (Nx-0.5)))
391 3 : ix = 2 * Nx * (ix > (Nx-0.5)) - ix-1;
392 11 : while ((iy < -0.5) or (iy > (Ny-0.5)))
393 5 : iy = 2 * Ny * (iy > (Ny-0.5)) - iy-1;
394 9 : while ((iz < -0.5) or (iz > (Nz-0.5)))
395 3 : iz = 2 * Nz * (iz > (Nz-0.5)) - iz-1;
396 : } else {
397 11 : ix = round(fmod(r.x, Nx));
398 11 : iy = round(fmod(r.y, Ny));
399 11 : iz = round(fmod(r.z, Nz));
400 11 : ix = (ix + Nx * (ix < 0)) % Nx;
401 11 : iy = (iy + Ny * (iy < 0)) % Ny;
402 11 : iz = (iz + Nz * (iz < 0)) % Nz;
403 : }
404 17 : return get(ix, iy, iz);
405 : }
406 :
407 : private:
408 : #ifdef HAVE_SIMD
409 : __m128 simdperiodicGet(size_t ix, size_t iy, size_t iz) const {
410 576 : ix = periodicBoundary(ix, Nx);
411 576 : iy = periodicBoundary(iy, Ny);
412 576 : iz = periodicBoundary(iz, Nz);
413 576 : return convertVector3fToSimd(grid[ix * Ny * Nz + iy * Nz + iz]);
414 : }
415 :
416 384 : __m128 simdreflectiveGet(size_t ix, size_t iy, size_t iz) const {
417 384 : ix = reflectiveBoundary(ix, Nx);
418 384 : iy = reflectiveBoundary(iy, Ny);
419 384 : iz = reflectiveBoundary(iz, Nz);
420 384 : return convertVector3fToSimd(grid[ix * Ny * Nz + iy * Nz + iz]);
421 : }
422 :
423 : __m128 convertVector3fToSimd(const Vector3f v) const {
424 : __m128 simdVar = _mm_set_ps(0,v.z,v.y,v.x);
425 : return simdVar;
426 : }
427 :
428 : Vector3f convertSimdToVector3f(__m128 res) const {
429 : float vec[4];
430 : _mm_store_ps(&vec[0], res);
431 : Vector3f result = Vector3f(vec[0], vec[1], vec[2]);
432 : return result;
433 : }
434 :
435 : /** Vectorized cubic Interpolator in 1D */
436 : __m128 CubicInterpolate(__m128 p0,__m128 p1,__m128 p2,__m128 p3,double position) const {
437 : __m128 c1 = _mm_set1_ps (1/2.);
438 : __m128 c2 = _mm_set1_ps (3/2.);
439 : __m128 c3 = _mm_set1_ps (2.);
440 : __m128 c4 = _mm_set1_ps (5/2.);
441 :
442 315 : __m128 pos = _mm_set1_ps (position);
443 315 : __m128 pos2 = _mm_set1_ps (position*position);
444 300 : __m128 pos3 = _mm_set1_ps (position*position*position);
445 :
446 : /** SIMD optimized routine to calculate 'res = ((-0.5*p0+3/2.*p1-3/2.*p2+0.5*p3)*pos*pos*pos+(p0-5/2.*p1+p2*2-0.5*p3)*pos*pos+(-0.5*p0+0.5*p2)*pos+p1);'
447 : where terms are used as:
448 : term = (-0.5*p0+0.5*p2)*pos
449 : term2 = (p0-5/2.*p1+p2*2-0.5*p3)*pos*pos;
450 : term3 = (-0.5*p0+3/2.*p1-3/2.*p2+0.5*p3)*pos*pos*pos; */
451 : __m128 term = _mm_mul_ps(_mm_sub_ps(_mm_mul_ps(c1,p2),_mm_mul_ps(c1,p0)),pos);
452 : __m128 term2 = _mm_mul_ps(_mm_sub_ps(_mm_add_ps(p0,_mm_mul_ps(c3,p2)),_mm_add_ps(_mm_mul_ps(c4,p1),_mm_mul_ps(c1,p3))),pos2);
453 : __m128 term3 = _mm_mul_ps(_mm_sub_ps(_mm_add_ps(_mm_mul_ps(c2,p1),_mm_mul_ps(c1,p3)),_mm_add_ps(_mm_mul_ps(c1,p0),_mm_mul_ps(c2,p2))),pos3);
454 : __m128 res = _mm_add_ps(_mm_add_ps(_mm_add_ps(term3,term2),term),p1);
455 : return res;
456 : }
457 : #endif // HAVE_SIMD
458 : /** Interpolate the grid tricubic at a given position (see https://www.paulinternet.nl/?page=bicubic, http://graphics.cs.cmu.edu/nsp/course/15-462/Fall04/assts/catmullRom.pdf) */
459 15 : Vector3f tricubicInterpolate(Vector3f, const Vector3d &position) const {
460 : #ifdef HAVE_SIMD
461 : // position on a unit grid
462 : Vector3d r = (position - gridOrigin) / spacing;
463 :
464 : int iX0, iY0, iZ0;
465 15 : iX0 = floor(r.x);
466 15 : iY0 = floor(r.y);
467 15 : iZ0 = floor(r.z);
468 :
469 : double fX, fY, fZ;
470 15 : fX = r.x - iX0;
471 15 : fY = r.y - iY0;
472 15 : fZ = r.z - iZ0;
473 :
474 : int nrCubicInterpolations = 4;
475 15 : __m128 interpolateVaryX[nrCubicInterpolations];
476 15 : __m128 interpolateVaryY[nrCubicInterpolations];
477 15 : __m128 interpolateVaryZ[nrCubicInterpolations];
478 : /** Perform 1D interpolations while iterating in each for loop over the index of another direction */
479 75 : for (int iLoopX = -1; iLoopX < nrCubicInterpolations-1; iLoopX++) {
480 300 : for (int iLoopY = -1; iLoopY < nrCubicInterpolations-1; iLoopY++) {
481 1200 : for (int iLoopZ = -1; iLoopZ < nrCubicInterpolations-1; iLoopZ++) {
482 960 : if (reflective)
483 384 : interpolateVaryZ[iLoopZ+1] = simdreflectiveGet(iX0+iLoopX, iY0+iLoopY, iZ0+iLoopZ);
484 : else
485 576 : interpolateVaryZ[iLoopZ+1] = simdperiodicGet(iX0+iLoopX, iY0+iLoopY, iZ0+iLoopZ);
486 : }
487 240 : interpolateVaryY[iLoopY+1] = CubicInterpolate(interpolateVaryZ[0], interpolateVaryZ[1], interpolateVaryZ[2], interpolateVaryZ[3], fZ);
488 : }
489 60 : interpolateVaryX[iLoopX+1] = CubicInterpolate(interpolateVaryY[0], interpolateVaryY[1], interpolateVaryY[2], interpolateVaryY[3], fY);
490 : }
491 15 : __m128 result = CubicInterpolate(interpolateVaryX[0], interpolateVaryX[1], interpolateVaryX[2], interpolateVaryX[3], fX);
492 15 : return convertSimdToVector3f(result);
493 : #else // HAVE_SIMD
494 : throw std::runtime_error( "Tried to use tricubic Interpolation without SIMD_EXTENSION. SIMD Optimization is necessary for tricubic interpolation of vector grids.\n");
495 : #endif // HAVE_SIMD
496 15 : }
497 :
498 : /** Vectorized cubic Interpolator in 1D that returns a scalar (see https://www.paulinternet.nl/?page=bicubic, http://graphics.cs.cmu.edu/nsp/course/15-462/Fall04/assts/catmullRom.pdf) */
499 63 : double CubicInterpolateScalar(double p0,double p1,double p2,double p3,double pos) const {
500 63 : return((-0.5*p0+3/2.*p1-3/2.*p2+0.5*p3)*pos*pos*pos+(p0-5/2.*p1+p2*2-0.5*p3)*pos*pos+(-0.5*p0+0.5*p2)*pos+p1);
501 : }
502 :
503 : /** Interpolate the grid tricubic at a given position (see https://www.paulinternet.nl/?page=bicubic, http://graphics.cs.cmu.edu/nsp/course/15-462/Fall04/assts/catmullRom.pdf) */
504 3 : double tricubicInterpolate(double, const Vector3d &position) const {
505 : /** position on a unit grid */
506 : Vector3d r = (position - gridOrigin) / spacing;
507 :
508 : int iX0, iY0, iZ0;
509 3 : iX0 = floor(r.x);
510 3 : iY0 = floor(r.y);
511 3 : iZ0 = floor(r.z);
512 :
513 : double fX, fY, fZ;
514 3 : fX = r.x - iX0;
515 3 : fY = r.y - iY0;
516 3 : fZ = r.z - iZ0;
517 :
518 : int nrCubicInterpolations = 4;
519 3 : double interpolateVaryX[nrCubicInterpolations];
520 3 : double interpolateVaryY[nrCubicInterpolations];
521 3 : double interpolateVaryZ[nrCubicInterpolations];
522 : /** Perform 1D interpolations while iterating in each for loop over the index of another direction */
523 15 : for (int iLoopX = -1; iLoopX < nrCubicInterpolations-1; iLoopX++) {
524 60 : for (int iLoopY = -1; iLoopY < nrCubicInterpolations-1; iLoopY++) {
525 240 : for (int iLoopZ = -1; iLoopZ < nrCubicInterpolations-1; iLoopZ++) {
526 192 : if (reflective)
527 0 : interpolateVaryZ[iLoopZ+1] = reflectiveGet(iX0+iLoopX, iY0+iLoopY, iZ0+iLoopZ);
528 : else
529 192 : interpolateVaryZ[iLoopZ+1] = periodicGet(iX0+iLoopX, iY0+iLoopY, iZ0+iLoopZ);
530 : }
531 48 : interpolateVaryY[iLoopY+1] = CubicInterpolateScalar(interpolateVaryZ[0], interpolateVaryZ[1], interpolateVaryZ[2], interpolateVaryZ[3], fZ);
532 : }
533 12 : interpolateVaryX[iLoopX+1] = CubicInterpolateScalar(interpolateVaryY[0], interpolateVaryY[1], interpolateVaryY[2], interpolateVaryY[3], fY);
534 : }
535 3 : double result = CubicInterpolateScalar(interpolateVaryX[0], interpolateVaryX[1], interpolateVaryX[2], interpolateVaryX[3], fX);
536 3 : return result;
537 3 : }
538 :
539 : /** Interpolate the grid trilinear at a given position */
540 100055 : T trilinearInterpolate(const Vector3d &position) const {
541 : /** position on a unit grid */
542 : Vector3d r = (position - gridOrigin) / spacing;
543 :
544 : /** indices of lower (0) and upper (1) neighbours. The neighbours span a grid
545 : with the origin at [iX0, iY0, iZ0] and the most distant corner [iX1, iY1, iZ1]. */
546 : int iX0, iX1, iY0, iY1, iZ0, iZ1;
547 : double resX, resY, resZ, fX0, fY0, fZ0;
548 :
549 100055 : if (reflective) {
550 6 : reflectiveClamp(r.x, Nx, iX0, iX1, resX);
551 6 : reflectiveClamp(r.y, Ny, iY0, iY1, resY);
552 6 : reflectiveClamp(r.z, Nz, iZ0, iZ1, resZ);
553 6 : fX0 = resX - floor(resX);
554 6 : fY0 = resY - floor(resY);
555 6 : fZ0 = resZ - floor(resZ);
556 : } else {
557 100049 : periodicClamp(r.x, Nx, iX0, iX1);
558 100049 : periodicClamp(r.y, Ny, iY0, iY1);
559 100049 : periodicClamp(r.z, Nz, iZ0, iZ1);
560 100049 : fX0 = r.x - floor(r.x);
561 100049 : fY0 = r.y - floor(r.y);
562 100049 : fZ0 = r.z - floor(r.z);
563 : }
564 :
565 : /** linear fraction to upper neighbours based on lower neighbours calculated above */
566 100055 : double fX1 = 1 - fX0;
567 100055 : double fY1 = 1 - fY0;
568 100055 : double fZ1 = 1 - fZ0;
569 :
570 : /** trilinear interpolation (see http://paulbourke.net/miscellaneous/interpolation) */
571 : T b(0.);
572 100055 : b += get(iX0, iY0, iZ0) * fX1 * fY1 * fZ1;
573 100055 : b += get(iX1, iY0, iZ0) * fX0 * fY1 * fZ1;
574 100055 : b += get(iX0, iY1, iZ0) * fX1 * fY0 * fZ1;
575 100055 : b += get(iX0, iY0, iZ1) * fX1 * fY1 * fZ0;
576 9 : b += get(iX1, iY0, iZ1) * fX0 * fY1 * fZ0;
577 9 : b += get(iX0, iY1, iZ1) * fX1 * fY0 * fZ0;
578 9 : b += get(iX1, iY1, iZ0) * fX0 * fY0 * fZ1;
579 9 : b += get(iX1, iY1, iZ1) * fX0 * fY0 * fZ0;
580 :
581 100055 : return b;
582 : }
583 :
584 : }; // class Grid
585 :
586 : typedef Grid<double> Grid1d;
587 : typedef Grid<float> Grid1f;
588 : typedef Grid<Vector3f> Grid3f;
589 : typedef Grid<Vector3d> Grid3d;
590 :
591 : /** @}*/
592 :
593 : } // namespace crpropa
594 :
595 : #endif // CRPROPA_GRID_H
|
