iMSTK
Interactive Medical Simulation Toolkit
imstkImplicitFunctionFiniteDifferenceFunctor.h
1 /*
2 ** This file is part of the Interactive Medical Simulation Toolkit (iMSTK)
3 ** iMSTK is distributed under the Apache License, Version 2.0.
4 ** See accompanying NOTICE for details.
5 */
6 
7 #pragma once
8 
9 #include "imstkImplicitGeometry.h"
10 #include "imstkSignedDistanceField.h"
11 
12 namespace imstk
13 {
19 struct ImplicitFunctionGradient
20 {
21  public:
22  virtual ~ImplicitFunctionGradient() = default;
23 
24  virtual Vec3d operator()(const Vec3d& pos) const = 0;
25 
26  virtual void setDx(const Vec3d& dx)
27  {
28  this->m_dx = dx;
29  this->m_invDx = Vec3d(1.0 / dx[0], 1.0 / dx[1], 1.0 / dx[2]);
30  }
31 
32  const Vec3d& getDx() const { return m_dx; }
33 
34  void setFunction(std::shared_ptr<ImplicitGeometry> func) { this->m_func = func; }
35  std::shared_ptr<ImplicitGeometry> getFunction() const { return m_func; }
36 
37  protected:
38  std::shared_ptr<ImplicitGeometry> m_func;
39  Vec3d m_dx = Vec3d(1.0, 1.0, 1.0);
40  Vec3d m_invDx = Vec3d(1.0, 1.0, 1.0);
41 };
42 
48 struct ImplicitFunctionCentralGradient : public ImplicitFunctionGradient
49 {
50  public:
51  ~ImplicitFunctionCentralGradient() override = default;
52 
53  Vec3d operator()(const Vec3d& pos) const override
54  {
55  const ImplicitGeometry& funcRef = *m_func;
56  return Vec3d(
57  funcRef.getFunctionValue(Vec3d(pos[0] + m_dx[0], pos[1], pos[2])) - funcRef.getFunctionValue(Vec3d(pos[0] - m_dx[0], pos[1], pos[2])),
58  funcRef.getFunctionValue(Vec3d(pos[0], pos[1] + m_dx[1], pos[2])) - funcRef.getFunctionValue(Vec3d(pos[0], pos[1] - m_dx[1], pos[2])),
59  funcRef.getFunctionValue(Vec3d(pos[0], pos[1], pos[2] + m_dx[2])) - funcRef.getFunctionValue(Vec3d(pos[0], pos[1], pos[2] - m_dx[2]))).cwiseProduct(m_invDx) * 0.5;
60  }
61 };
62 
68 struct StructuredCentralGradient : public ImplicitFunctionGradient
69 {
70  public:
71  Vec3d operator()(const Vec3d& pos) const override
72  {
73  const SignedDistanceField& funcRef = *static_cast<SignedDistanceField*>(m_func.get());
74  return Vec3d(
75  funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0] + m_dxi[0]), static_cast<int>(pos[1]), static_cast<int>(pos[2]))) -
76  funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0] - m_dxi[0]), static_cast<int>(pos[1]), static_cast<int>(pos[2]))),
77  funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0]), static_cast<int>(pos[1] + m_dxi[1]), static_cast<int>(pos[2]))) -
78  funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0]), static_cast<int>(pos[1] - m_dxi[1]), static_cast<int>(pos[2]))),
79  funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0]), static_cast<int>(pos[1]), static_cast<int>(pos[2] + m_dxi[2]))) -
80  funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0]), static_cast<int>(pos[1]), static_cast<int>(pos[2] - m_dxi[2])))).cwiseProduct(m_invDx) * 0.5;
81  }
82 
83  void setDx(const Vec3i& dx, const Vec3d& dxs)
84  {
85  m_dxi = dx;
86  ImplicitFunctionGradient::setDx(dxs);
87  }
88 
89  protected:
90  using ImplicitFunctionGradient::setDx;
91  Vec3i m_dxi;
92 };
93 
99 struct ImplicitFunctionForwardGradient : public ImplicitFunctionGradient
100 {
101  public:
102  Vec3d operator()(const Vec3d& pos) const override
103  {
104  const ImplicitGeometry& funcRef = *m_func;
105  const double centralValue = funcRef.getFunctionValue(Vec3d(pos[0], pos[1], pos[2]));
106  const double maxValueX = funcRef.getFunctionValue(Vec3d(pos[0] + m_dx[0], pos[1], pos[2]));
107  const double maxValueY = funcRef.getFunctionValue(Vec3d(pos[0], pos[1] + m_dx[1], pos[2]));
108  const double maxValueZ = funcRef.getFunctionValue(Vec3d(pos[0], pos[1], pos[2] + m_dx[2]));
109  return Vec3d(
110  maxValueX - centralValue,
111  maxValueY - centralValue,
112  maxValueZ - centralValue).cwiseProduct(m_invDx);
113  }
114 };
115 
121 struct ImplicitFunctionBackwardGradient : public ImplicitFunctionGradient
122 {
123  public:
124  virtual ~ImplicitFunctionBackwardGradient() = default;
125 
126  Vec3d operator()(const Vec3d& pos) const override
127  {
128  const ImplicitGeometry& funcRef = *m_func;
129  const double centralValue = funcRef.getFunctionValue(Vec3d(pos[0], pos[1], pos[2]));
130  const double minValueX = funcRef.getFunctionValue(Vec3d(pos[0] - m_dx[0], pos[1], pos[2]));
131  const double minValueY = funcRef.getFunctionValue(Vec3d(pos[0], pos[1] - m_dx[1], pos[2]));
132  const double minValueZ = funcRef.getFunctionValue(Vec3d(pos[0], pos[1], pos[2] - m_dx[2]));
133  return Vec3d(
134  centralValue - minValueX,
135  centralValue - minValueY,
136  centralValue - minValueZ).cwiseProduct(m_invDx);
137  }
138 };
139 
145 struct StructuredForwardGradient : public ImplicitFunctionGradient
146 {
147  public:
148  virtual ~StructuredForwardGradient() = default;
149 
150  inline Vec3d operator()(const Vec3d& pos) const override
151  {
152  const SignedDistanceField& funcRef = *static_cast<SignedDistanceField*>(m_func.get());
153  const double centralValue = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0]), static_cast<int>(pos[1]), static_cast<int>(pos[2])));
154  const double maxValueX = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0] + m_dxi[0]), static_cast<int>(pos[1]), static_cast<int>(pos[2])));
155  const double maxValueY = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0]), static_cast<int>(pos[1] + m_dxi[1]), static_cast<int>(pos[2])));
156  const double maxValueZ = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0]), static_cast<int>(pos[1]), static_cast<int>(pos[2] + m_dxi[2])));
157  return Vec3d(
158  maxValueX - centralValue,
159  maxValueY - centralValue,
160  maxValueZ - centralValue).cwiseProduct(m_invDx);
161  }
162 
163  void setDx(const Vec3i& dx, const Vec3d& dxs)
164  {
165  m_dxi = dx;
166  ImplicitFunctionGradient::setDx(dxs.cwiseProduct(m_dxi.cast<double>()));
167  }
168 
169  protected:
170  using ImplicitFunctionGradient::setDx;
171  Vec3i m_dxi;
172 };
173 
179 struct StructuredBackwardGradient : public ImplicitFunctionGradient
180 {
181  public:
182  virtual ~StructuredBackwardGradient() = default;
183 
184  inline Vec3d operator()(const Vec3d& pos) const override
185  {
186  const SignedDistanceField& funcRef = *static_cast<SignedDistanceField*>(m_func.get());
187  const double centralValue = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0]), static_cast<int>(pos[1]), static_cast<int>(pos[2])));
188  const double minValueX = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0] - m_dxi[0]), static_cast<int>(pos[1]), static_cast<int>(pos[2])));
189  const double minValueY = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0]), static_cast<int>(pos[1] - m_dxi[1]), static_cast<int>(pos[2])));
190  const double minValueZ = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0]), static_cast<int>(pos[1]), static_cast<int>(pos[2] - m_dxi[2])));
191  return Vec3d(
192  centralValue - minValueX,
193  centralValue - minValueY,
194  centralValue - minValueZ).cwiseProduct(m_invDx);
195  }
196 
197  void setDx(const Vec3i& dx, const Vec3d& dxs)
198  {
199  m_dxi = dx;
200  ImplicitFunctionGradient::setDx(dxs);
201  }
202 
203  protected:
204  using ImplicitFunctionGradient::setDx;
205  Vec3i m_dxi;
206 };
207 
213 struct ImplicitStructuredCurvature
214 {
215  public:
216  double operator()(const Vec3d& pos) const
217  {
218  const SignedDistanceField& funcRef = *static_cast<SignedDistanceField*>(m_func.get());
219 
220  const double central = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0]), static_cast<int>(pos[1]), static_cast<int>(pos[2])));
221  const double minX = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0] - m_dxi[0]), static_cast<int>(pos[1]), static_cast<int>(pos[2])));
222  const double minY = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0]), static_cast<int>(pos[1] - m_dxi[1]), static_cast<int>(pos[2])));
223  const double minZ = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0]), static_cast<int>(pos[1]), static_cast<int>(pos[2] - m_dxi[2])));
224  const double maxX = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0] + m_dxi[0]), static_cast<int>(pos[1]), static_cast<int>(pos[2])));
225  const double maxY = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0]), static_cast<int>(pos[1] + m_dxi[1]), static_cast<int>(pos[2])));
226  const double maxZ = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0]), static_cast<int>(pos[1]), static_cast<int>(pos[2] + m_dxi[2])));
227 
228  const double minXY = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0] - m_dxi[0]), static_cast<int>(pos[1] - m_dxi[1]), static_cast<int>(pos[2])));
229  const double maxXY = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0] + m_dxi[0]), static_cast<int>(pos[1] + m_dxi[1]), static_cast<int>(pos[2])));
230  const double maxXminY = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0] + m_dxi[0]), static_cast<int>(pos[1] - m_dxi[1]), static_cast<int>(pos[2])));
231  const double minXmaxY = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0] - m_dxi[0]), static_cast<int>(pos[1] + m_dxi[1]), static_cast<int>(pos[2])));
232 
233  const double minXmaxZ = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0] - m_dxi[0]), static_cast<int>(pos[1]), static_cast<int>(pos[2] + m_dxi[2])));
234  const double maxXminZ = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0] + m_dxi[0]), static_cast<int>(pos[1]), static_cast<int>(pos[2] - m_dxi[2])));
235  const double maxXZ = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0] + m_dxi[0]), static_cast<int>(pos[1]), static_cast<int>(pos[2] + m_dxi[2])));
236  const double minXZ = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0] - m_dxi[0]), static_cast<int>(pos[1]), static_cast<int>(pos[2] - m_dxi[2])));
237 
238  const double minYmaxZ = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0]), static_cast<int>(pos[1] - m_dxi[1]), static_cast<int>(pos[2] + m_dxi[2])));
239  const double maxYminZ = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0]), static_cast<int>(pos[1] + m_dxi[1]), static_cast<int>(pos[2] - m_dxi[2])));
240  const double maxYZ = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0]), static_cast<int>(pos[1] + m_dxi[1]), static_cast<int>(pos[2] + m_dxi[2])));
241  const double minYZ = funcRef.getFunctionValueCoord(Vec3i(static_cast<int>(pos[0]), static_cast<int>(pos[1] - m_dxi[1]), static_cast<int>(pos[2] - m_dxi[2])));
242 
243  const double dx = (maxX - minX) * 0.5;
244  const double dxx = maxX - 2.0 * central + minX;
245  const double dx2 = dx * dx;
246 
247  const double dy = (maxY - minY) * 0.5;
248  const double dyy = maxY - 2.0 * central + minY;
249  const double dy2 = dy * dy;
250 
251  const double dz = (maxZ - minZ) * 0.5;
252  const double dzz = maxZ - 2.0 * central + minZ;
253  const double dz2 = dz * dz;
254 
255  const double dxy = (maxXY + minXY - maxXminY - minXmaxY) * 0.25;
256  const double dxz = (minXmaxZ + maxXminZ - maxXZ - minXZ) * 0.25;
257  const double dyz = (minYmaxZ + maxYminZ - maxYZ - minYZ) * 0.25;
258 
259  return ((dxx * (dy2 + dz2) + dyy * (dx2 + dz2) + dzz * (dx2 + dy2) -
260  2.0 * dx * dy * dxy - 2.0 * dx * dz * dxz - 2.0 * dy * dz * dyz) /
261  (dx2 + dy2 + dz2 + std::numeric_limits<double>::epsilon()));
262  }
263 
264  void setDx(const Vec3i& dx, const Vec3d& dxs)
265  {
266  m_dxi = dx;
267  this->m_dx = dxs;
268  this->m_invDx = Vec3d(1.0 / dxs[0], 1.0 / dxs[1], 1.0 / dxs[2]);
269  }
270 
271  public:
272  const Vec3d& getDx() const { return m_dx; }
273 
274  void setFunction(std::shared_ptr<ImplicitGeometry> func) { this->m_func = func; }
275 
276  protected:
277  std::shared_ptr<ImplicitGeometry> m_func;
278  Vec3d m_dx = Vec3d(1.0, 1.0, 1.0);
279  Vec3d m_invDx = Vec3d(1.0, 1.0, 1.0);
280 
281  protected:
282  Vec3i m_dxi;
283 };
284 } // namespace imstk
imstk
Compound Geometry.
Definition: OctreeDebugModel.cpp:10
imstk::StructuredCentralGradient
Gradient given by central finite differences.
Definition: imstkImplicitFunctionFiniteDifferenceFunctor.h:68
imstk::ImplicitGeometry::getFunctionValue
virtual double getFunctionValue(const Vec3d &pos) const =0
Returns function value given position.
imstk::SignedDistanceField::getFunctionValueCoord
double getFunctionValueCoord(const Vec3i &coord) const
Returns signed distance to surface at coordinate inlined for performance.
Definition: imstkSignedDistanceField.h:48
imstk::ImplicitFunctionCentralGradient
Gradient given by central finite differences.
Definition: imstkImplicitFunctionFiniteDifferenceFunctor.h:48
imstk::SignedDistanceField
Structured field of signed distances implemented with ImageData The SDF differs in that when you scal...
Definition: imstkSignedDistanceField.h:26
imstk::StructuredBackwardGradient
Gradient given by backward finite differences.
Definition: imstkImplicitFunctionFiniteDifferenceFunctor.h:179
imstk::ImplicitFunctionForwardGradient
Gradient given by forward finite differences.
Definition: imstkImplicitFunctionFiniteDifferenceFunctor.h:99
imstk::ImplicitGeometry
Class that can represent the geometry of multiple implicit geometries as boolean functions One may su...
Definition: imstkImplicitGeometry.h:19
imstk::StructuredForwardGradient
Gradient given by forward finite differences.
Definition: imstkImplicitFunctionFiniteDifferenceFunctor.h:145
Generated by   doxygen 1.8.13