ArraySamplers3-Impl.h
Go to the documentation of this file.
1 /*************************************************************************
2 > File Name: ArraySamplers3-Impl.h
3 > Project Name: CubbyFlow
4 > Author: Chan-Ho Chris Ohk
5 > Purpose: 3-D nearest array sampler class.
6 > Created Time: 2017/05/19
7 > Copyright (c) 2018, Chan-Ho Chris Ohk
8 *************************************************************************/
9 #ifndef CUBBYFLOW_ARRAY_SAMPLERS3_IMPL_H
10 #define CUBBYFLOW_ARRAY_SAMPLERS3_IMPL_H
11 
12 #include <Core/Math/MathUtils.h>
13 
14 namespace CubbyFlow
15 {
16  template <typename T, typename R>
18  const ConstArrayAccessor3<T>& accessor,
19  const Vector3<R>& gridSpacing,
20  const Vector3<R>& gridOrigin)
21  {
22  m_gridSpacing = gridSpacing;
23  m_origin = gridOrigin;
24  m_accessor = accessor;
25  }
26 
27  template <typename T, typename R>
29  {
30  m_gridSpacing = other.m_gridSpacing;
31  m_origin = other.m_origin;
32  m_accessor = other.m_accessor;
33  }
34 
35  template <typename T, typename R>
37  {
38  ssize_t i, j, k;
39  R fx, fy, fz;
40 
41  assert(m_gridSpacing.x > std::numeric_limits<R>::epsilon());
42  assert(m_gridSpacing.y > std::numeric_limits<R>::epsilon());
43  assert(m_gridSpacing.z > std::numeric_limits<R>::epsilon());
44 
45  const Vector3<R> normalizedX = (pt - m_origin) / m_gridSpacing;
46 
47  const ssize_t iSize = static_cast<ssize_t>(m_accessor.size().x);
48  const ssize_t jSize = static_cast<ssize_t>(m_accessor.size().y);
49  const ssize_t kSize = static_cast<ssize_t>(m_accessor.size().z);
50 
51  GetBarycentric(normalizedX.x, 0, iSize - 1, &i, &fx);
52  GetBarycentric(normalizedX.y, 0, jSize - 1, &j, &fy);
53  GetBarycentric(normalizedX.z, 0, kSize - 1, &k, &fz);
54 
55  i = std::min(static_cast<ssize_t>(i + fx + 0.5), iSize - 1);
56  j = std::min(static_cast<ssize_t>(j + fy + 0.5), jSize - 1);
57  k = std::min(static_cast<ssize_t>(k + fz + 0.5), kSize - 1);
58 
59  return m_accessor(i, j, k);
60  }
61 
62  template <typename T, typename R>
64  {
65  ssize_t i, j, k;
66  R fx, fy, fz;
67 
68  assert(m_gridSpacing.x > std::numeric_limits<R>::epsilon());
69  assert(m_gridSpacing.y > std::numeric_limits<R>::epsilon());
70  assert(m_gridSpacing.z > std::numeric_limits<R>::epsilon());
71 
72  const Vector3<R> normalizedX = (pt - m_origin) / m_gridSpacing;
73 
74  const ssize_t iSize = static_cast<ssize_t>(m_accessor.size().x);
75  const ssize_t jSize = static_cast<ssize_t>(m_accessor.size().y);
76  const ssize_t kSize = static_cast<ssize_t>(m_accessor.size().z);
77 
78  GetBarycentric(normalizedX.x, 0, iSize - 1, &i, &fx);
79  GetBarycentric(normalizedX.y, 0, jSize - 1, &j, &fy);
80  GetBarycentric(normalizedX.z, 0, kSize - 1, &k, &fz);
81 
82  index->x = std::min(static_cast<ssize_t>(i + fx + 0.5), iSize - 1);
83  index->y = std::min(static_cast<ssize_t>(j + fy + 0.5), jSize - 1);
84  index->z = std::min(static_cast<ssize_t>(k + fz + 0.5), kSize - 1);
85  }
86 
87  template <typename T, typename R>
88  std::function<T(const Vector3<R>&)> NearestArraySampler<T, R, 3>::Functor() const
89  {
90  NearestArraySampler sampler(*this);
91  return std::bind(&NearestArraySampler::operator(), sampler, std::placeholders::_1);
92  }
93 
94  template <typename T, typename R>
96  const ConstArrayAccessor3<T>& accessor,
97  const Vector3<R>& gridSpacing,
98  const Vector3<R>& gridOrigin)
99  {
100  m_gridSpacing = gridSpacing;
101  m_invGridSpacing = static_cast<R>(1) / m_gridSpacing;
102  m_origin = gridOrigin;
103  m_accessor = accessor;
104  }
105 
106  template <typename T, typename R>
108  {
109  m_gridSpacing = other.m_gridSpacing;
110  m_invGridSpacing = other.m_invGridSpacing;
111  m_origin = other.m_origin;
112  m_accessor = other.m_accessor;
113  }
114 
115  template <typename T, typename R>
117  {
118  ssize_t i, j, k;
119  R fx, fy, fz;
120 
121  assert(m_gridSpacing.x > std::numeric_limits<R>::epsilon());
122  assert(m_gridSpacing.y > std::numeric_limits<R>::epsilon());
123  assert(m_gridSpacing.z > std::numeric_limits<R>::epsilon());
124 
125  const Vector3<R> normalizedX = (pt - m_origin) / m_gridSpacing;
126 
127  const ssize_t iSize = static_cast<ssize_t>(m_accessor.size().x);
128  const ssize_t jSize = static_cast<ssize_t>(m_accessor.size().y);
129  const ssize_t kSize = static_cast<ssize_t>(m_accessor.size().z);
130 
131  GetBarycentric(normalizedX.x, 0, iSize - 1, &i, &fx);
132  GetBarycentric(normalizedX.y, 0, jSize - 1, &j, &fy);
133  GetBarycentric(normalizedX.z, 0, kSize - 1, &k, &fz);
134 
135  const ssize_t ip1 = std::min(i + 1, iSize - 1);
136  const ssize_t jp1 = std::min(j + 1, jSize - 1);
137  const ssize_t kp1 = std::min(k + 1, kSize - 1);
138 
139  return TriLerp(
140  m_accessor(i, j, k), m_accessor(ip1, j, k),
141  m_accessor(i, jp1, k), m_accessor(ip1, jp1, k),
142  m_accessor(i, j, kp1), m_accessor(ip1, j, kp1),
143  m_accessor(i, jp1, kp1), m_accessor(ip1, jp1, kp1),
144  fx, fy, fz);
145  }
146 
147  template <typename T, typename R>
149  const Vector3<R>& pt,
150  std::array<Point3UI, 8>* indices,
151  std::array<R, 8>* weights) const
152  {
153  ssize_t i, j, k;
154  R fx, fy, fz;
155 
156  assert(m_gridSpacing.x > std::numeric_limits<R>::epsilon());
157  assert(m_gridSpacing.y > std::numeric_limits<R>::epsilon());
158  assert(m_gridSpacing.z > std::numeric_limits<R>::epsilon());
159 
160  const Vector3<R> normalizedX = (pt - m_origin) / m_gridSpacing;
161 
162  const ssize_t iSize = static_cast<ssize_t>(m_accessor.size().x);
163  const ssize_t jSize = static_cast<ssize_t>(m_accessor.size().y);
164  const ssize_t kSize = static_cast<ssize_t>(m_accessor.size().z);
165 
166  GetBarycentric(normalizedX.x, 0, iSize - 1, &i, &fx);
167  GetBarycentric(normalizedX.y, 0, jSize - 1, &j, &fy);
168  GetBarycentric(normalizedX.z, 0, kSize - 1, &k, &fz);
169 
170  const ssize_t ip1 = std::min(i + 1, iSize - 1);
171  const ssize_t jp1 = std::min(j + 1, jSize - 1);
172  const ssize_t kp1 = std::min(k + 1, kSize - 1);
173 
174  (*indices)[0] = Point3UI(i, j, k);
175  (*indices)[1] = Point3UI(ip1, j, k);
176  (*indices)[2] = Point3UI(i, jp1, k);
177  (*indices)[3] = Point3UI(ip1, jp1, k);
178  (*indices)[4] = Point3UI(i, j, kp1);
179  (*indices)[5] = Point3UI(ip1, j, kp1);
180  (*indices)[6] = Point3UI(i, jp1, kp1);
181  (*indices)[7] = Point3UI(ip1, jp1, kp1);
182 
183  (*weights)[0] = (1 - fx) * (1 - fy) * (1 - fz);
184  (*weights)[1] = fx * (1 - fy) * (1 - fz);
185  (*weights)[2] = (1 - fx) * fy * (1 - fz);
186  (*weights)[3] = fx * fy * (1 - fz);
187  (*weights)[4] = (1 - fx) * (1 - fy) * fz;
188  (*weights)[5] = fx * (1 - fy) * fz;
189  (*weights)[6] = (1 - fx) * fy * fz;
190  (*weights)[7] = fx * fy * fz;
191  }
192 
193  template <typename T, typename R>
195  const Vector3<R>& x,
196  std::array<Point3UI, 8>* indices,
197  std::array<Vector3<R>, 8>* weights) const
198  {
199  ssize_t i, j, k;
200  R fx, fy, fz;
201 
202  assert(m_gridSpacing.x > 0.0);
203  assert(m_gridSpacing.y > 0.0);
204  assert(m_gridSpacing.z > 0.0);
205 
206  const Vector3<R> normalizedX = (x - m_origin) / m_gridSpacing;
207 
208  const ssize_t iSize = static_cast<ssize_t>(m_accessor.size().x);
209  const ssize_t jSize = static_cast<ssize_t>(m_accessor.size().y);
210  const ssize_t kSize = static_cast<ssize_t>(m_accessor.size().z);
211 
212  GetBarycentric(normalizedX.x, 0, iSize - 1, &i, &fx);
213  GetBarycentric(normalizedX.y, 0, jSize - 1, &j, &fy);
214  GetBarycentric(normalizedX.z, 0, kSize - 1, &k, &fz);
215 
216  const ssize_t ip1 = std::min(i + 1, iSize - 1);
217  const ssize_t jp1 = std::min(j + 1, jSize - 1);
218  const ssize_t kp1 = std::min(k + 1, kSize - 1);
219 
220  (*indices)[0] = Point3UI(i, j, k);
221  (*indices)[1] = Point3UI(ip1, j, k);
222  (*indices)[2] = Point3UI(i, jp1, k);
223  (*indices)[3] = Point3UI(ip1, jp1, k);
224  (*indices)[4] = Point3UI(i, j, kp1);
225  (*indices)[5] = Point3UI(ip1, j, kp1);
226  (*indices)[6] = Point3UI(i, jp1, kp1);
227  (*indices)[7] = Point3UI(ip1, jp1, kp1);
228 
229  (*weights)[0] = Vector3<R>(-m_invGridSpacing.x * (1 - fy) * (1 - fz), -m_invGridSpacing.y * (1 - fx) * (1 - fz), -m_invGridSpacing.z * (1 - fx) * (1 - fy));
230  (*weights)[1] = Vector3<R>(m_invGridSpacing.x * (1 - fy) * (1 - fz), fx * (-m_invGridSpacing.y) * (1 - fz), fx * (1 - fy) * (-m_invGridSpacing.z));
231  (*weights)[2] = Vector3<R>((-m_invGridSpacing.x) * fy * (1 - fz), (1 - fx) * m_invGridSpacing.y * (1 - fz), (1 - fx) * fy * (-m_invGridSpacing.z));
232  (*weights)[3] = Vector3<R>(m_invGridSpacing.x * fy * (1 - fz), fx * m_invGridSpacing.y * (1 - fz), fx * fy * (-m_invGridSpacing.z));
233  (*weights)[4] = Vector3<R>((-m_invGridSpacing.x) * (1 - fy) * fz, (1 - fx) * (-m_invGridSpacing.y) * fz, (1 - fx) * (1 - fy) * m_invGridSpacing.z);
234  (*weights)[5] = Vector3<R>(m_invGridSpacing.x * (1 - fy) * fz, fx * (-m_invGridSpacing.y) * fz, fx * (1 - fy) * m_invGridSpacing.z);
235  (*weights)[6] = Vector3<R>((-m_invGridSpacing.x) * fy * fz, (1 - fx) * m_invGridSpacing.y * fz, (1 - fx) * fy * m_invGridSpacing.z);
236  (*weights)[7] = Vector3<R>(m_invGridSpacing.x * fy * fz, fx * m_invGridSpacing.y * fz, fx * fy * m_invGridSpacing.z);
237  }
238 
239  template <typename T, typename R>
240  std::function<T(const Vector3<R>&)> LinearArraySampler<T, R, 3>::Functor() const
241  {
242  LinearArraySampler sampler(*this);
243  return std::bind(&LinearArraySampler::operator(), sampler, std::placeholders::_1);
244  }
245 
246  template <typename T, typename R>
248  const ConstArrayAccessor3<T>& accessor,
249  const Vector3<R>& gridSpacing,
250  const Vector3<R>& gridOrigin)
251  {
252  m_gridSpacing = gridSpacing;
253  m_origin = gridOrigin;
254  m_accessor = accessor;
255  }
256 
257  template <typename T, typename R>
259  {
260  m_gridSpacing = other.m_gridSpacing;
261  m_origin = other.m_origin;
262  m_accessor = other.m_accessor;
263  }
264 
265  template <typename T, typename R>
267  {
268  ssize_t i, j, k;
269  R fx, fy, fz;
270 
271  assert(m_gridSpacing.x > std::numeric_limits<R>::epsilon());
272  assert(m_gridSpacing.y > std::numeric_limits<R>::epsilon());
273  assert(m_gridSpacing.z > std::numeric_limits<R>::epsilon());
274 
275  const Vector3<R> normalizedX = (pt - m_origin) / m_gridSpacing;
276 
277  const ssize_t iSize = static_cast<ssize_t>(m_accessor.size().x);
278  const ssize_t jSize = static_cast<ssize_t>(m_accessor.size().y);
279  const ssize_t kSize = static_cast<ssize_t>(m_accessor.size().z);
280 
281  GetBarycentric(normalizedX.x, 0, iSize, &i, &fx);
282  GetBarycentric(normalizedX.y, 0, jSize, &j, &fy);
283  GetBarycentric(normalizedX.z, 0, kSize, &k, &fz);
284 
285  const ssize_t is[4] = { std::max(i - 1, ZERO_SSIZE), i, std::min(i + 1, iSize - 1), std::min(i + 2, iSize - 1) };
286  const ssize_t js[4] = { std::max(j - 1, ZERO_SSIZE), j, std::min(j + 1, jSize - 1), std::min(j + 2, jSize - 1) };
287  const ssize_t ks[4] = { std::max(k - 1, ZERO_SSIZE), k, std::min(k + 1, kSize - 1), std::min(k + 2, kSize - 1) };
288 
289  T kValues[4];
290 
291  for (int kk = 0; kk < 4; ++kk)
292  {
293  T jValues[4];
294 
295  for (int jj = 0; jj < 4; ++jj)
296  {
297  jValues[jj] = MonotonicCatmullRom(
298  m_accessor(is[0], js[jj], ks[kk]),
299  m_accessor(is[1], js[jj], ks[kk]),
300  m_accessor(is[2], js[jj], ks[kk]),
301  m_accessor(is[3], js[jj], ks[kk]),
302  fx);
303  }
304 
305  kValues[kk] = MonotonicCatmullRom(jValues[0], jValues[1], jValues[2], jValues[3], fy);
306  }
307 
308  return MonotonicCatmullRom(kValues[0], kValues[1], kValues[2], kValues[3], fz);
309  }
310 
311  template <typename T, typename R>
312  std::function<T(const Vector3<R>&)> CubicArraySampler<T, R, 3>::Functor() const
313  {
314  CubicArraySampler sampler(*this);
315  return std::bind(&CubicArraySampler::operator(), sampler, std::placeholders::_1);
316  }
317 }
318 
319 #endif
3-D vector class.
Definition: Vector3.h:26
T z
Z (or the third) component of the vector.
Definition: Vector3.h:38
3-D read-only array accessor class.
Definition: ArrayAccessor3.h:269
Generic N-D nearest array sampler class.
Definition: ArraySamplers.h:22
3-D point class.
Definition: Point3.h:26
T x
X (or the first) component of the vector.
Definition: Vector3.h:29
Definition: pybind11Utils.h:24
Generic N-D cubic array sampler class.
Definition: ArraySamplers.h:50
T z
Z (or the third) component of the point.
Definition: Point3.h:38
T y
Y (or the second) component of the point.
Definition: Point3.h:35
T x
X (or the first) component of the point.
Definition: Point3.h:29
constexpr ssize_t ZERO_SSIZE
Zero ssize_t.
Definition: Constants.h:20
S TriLerp(const S &f000, const S &f100, const S &f010, const S &f110, const S &f001, const S &f101, const S &f011, const S &f111, T tx, T ty, T fz)
Computes trilinear interpolation.
Definition: MathUtils-Impl.h:199
T MonotonicCatmullRom(const T &f0, const T &f1, const T &f2, const T &f3, T f)
Computes monotonic Catmull-Rom interpolation.
Definition: MathUtils-Impl.h:228
Generic N-D linear array sampler class.
Definition: ArraySamplers.h:36
T y
Y (or the second) component of the vector.
Definition: Vector3.h:35
Point3< size_t > Point3UI
Unsigned integer-type 3D point.
Definition: Point3.h:328
void GetBarycentric(T x, ssize_t iLow, ssize_t iHigh, ssize_t *i, T *f)
Gets the barycentric coordinate.
Definition: MathUtils-Impl.h:151