1 //============================================================================
2 // Copyright (c) Kitware, Inc.
3 // All rights reserved.
4 // See LICENSE.txt for details.
5 // This software is distributed WITHOUT ANY WARRANTY; without even
6 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
7 // PURPOSE. See the above copyright notice for more information.
8 //
9 // Copyright 2015 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
10 // Copyright 2015 UT-Battelle, LLC.
11 // Copyright 2015 Los Alamos National Security.
12 //
13 // Under the terms of Contract DE-NA0003525 with NTESS,
14 // the U.S. Government retains certain rights in this software.
15 //
16 // Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
17 // Laboratory (LANL), the U.S. Government retains certain rights in
18 // this software.
19 //============================================================================
20
21 #include <vtkm/VectorAnalysis.h>
22 #include <vtkm/cont/Algorithm.h>
23 #include <vtkm/rendering/raytracing/BVHTraverser.h>
24 #include <vtkm/rendering/raytracing/QuadIntersector.h>
25 #include <vtkm/rendering/raytracing/RayOperations.h>
26 #include <vtkm/worklet/DispatcherMapTopology.h>
27
28 namespace vtkm
29 {
30 namespace rendering
31 {
32 namespace raytracing
33 {
34 namespace detail
35 {
36
37 template <typename Device>
38 class QuadLeafIntersector : public vtkm::cont::ExecutionObjectBase
39 {
40 public:
41 using IdType = vtkm::Vec<vtkm::Id, 5>;
42 using IdHandle = vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 5>>;
43 using IdArrayPortal = typename IdHandle::ExecutionTypes<Device>::PortalConst;
44 using FloatHandle = vtkm::cont::ArrayHandle<vtkm::Float32>;
45 using FloatPortal = typename FloatHandle::ExecutionTypes<Device>::PortalConst;
46 IdArrayPortal QuadIds;
47
QuadLeafIntersector(const IdHandle & quadIds)48 QuadLeafIntersector(const IdHandle& quadIds)
49 : QuadIds(quadIds.PrepareForInput(Device()))
50 {
51 }
52
53 template <typename vec3, typename Precision>
quad(const vec3 & ray_origin,const vec3 & ray_direction,const vec3 & v00,const vec3 & v10,const vec3 & v11,const vec3 & v01,Precision & u,Precision & v,Precision & t) const54 VTKM_EXEC bool quad(const vec3& ray_origin,
55 const vec3& ray_direction,
56 const vec3& v00,
57 const vec3& v10,
58 const vec3& v11,
59 const vec3& v01,
60 Precision& u,
61 Precision& v,
62 Precision& t) const
63 {
64
65 /* An Efficient Ray-Quadrilateral Intersection Test
66 Ares Lagae Philip Dutr´e
67 http://graphics.cs.kuleuven.be/publications/LD05ERQIT/index.html
68
69 v01 *------------ * v11
70 |\ |
71 | \ |
72 | \ |
73 | \ |
74 | \ |
75 | \ |
76 v00 *------------* v10
77 */
78 // Rejects rays that are parallel to Q, and rays that intersect the plane of
79 // Q either on the left of the line V00V01 or on the right of the line V00V10.
80
81 vec3 E03 = v01 - v00;
82 vec3 P = vtkm::Cross(ray_direction, E03);
83 vec3 E01 = v10 - v00;
84 Precision det = vtkm::dot(E01, P);
85
86 if (vtkm::Abs(det) < vtkm::Epsilon<Precision>())
87 return false;
88 Precision inv_det = 1.0f / det;
89 vec3 T = ray_origin - v00;
90 Precision alpha = vtkm::dot(T, P) * inv_det;
91 if (alpha < 0.0)
92 return false;
93 vec3 Q = vtkm::Cross(T, E01);
94 Precision beta = vtkm::dot(ray_direction, Q) * inv_det;
95 if (beta < 0.0)
96 return false;
97
98 if ((alpha + beta) > 1.0f)
99 {
100
101 // Rejects rays that intersect the plane of Q either on the
102 // left of the line V11V10 or on the right of the line V11V01.
103
104 vec3 E23 = v01 - v11;
105 vec3 E21 = v10 - v11;
106 vec3 P_prime = vtkm::Cross(ray_direction, E21);
107 Precision det_prime = vtkm::dot(E23, P_prime);
108 if (vtkm::Abs(det_prime) < vtkm::Epsilon<Precision>())
109 return false;
110 Precision inv_det_prime = 1.0f / det_prime;
111 vec3 T_prime = ray_origin - v11;
112 Precision alpha_prime = vtkm::dot(T_prime, P_prime) * inv_det_prime;
113 if (alpha_prime < 0.0f)
114 return false;
115 vec3 Q_prime = vtkm::Cross(T_prime, E23);
116 Precision beta_prime = vtkm::dot(ray_direction, Q_prime) * inv_det_prime;
117 if (beta_prime < 0.0f)
118 return false;
119 }
120
121 // Compute the ray parameter of the intersection point, and
122 // reject the ray if it does not hit Q.
123
124 t = vtkm::dot(E03, Q) * inv_det;
125 if (t < 0.0)
126 return false;
127
128
129 // Compute the barycentric coordinates of V11
130 Precision alpha_11, beta_11;
131 vec3 E02 = v11 - v00;
132 vec3 n = vtkm::Cross(E01, E02);
133
134 if ((vtkm::Abs(n[0]) >= vtkm::Abs(n[1])) && (vtkm::Abs(n[0]) >= vtkm::Abs(n[2])))
135 {
136
137 alpha_11 = ((E02[1] * E03[2]) - (E02[2] * E03[1])) / n[0];
138 beta_11 = ((E01[1] * E02[2]) - (E01[2] * E02[1])) / n[0];
139 }
140 else if ((vtkm::Abs(n[1]) >= vtkm::Abs(n[0])) && (vtkm::Abs(n[1]) >= vtkm::Abs(n[2])))
141 {
142
143 alpha_11 = ((E02[2] * E03[0]) - (E02[0] * E03[2])) / n[1];
144 beta_11 = ((E01[2] * E02[0]) - (E01[0] * E02[2])) / n[1];
145 }
146 else
147 {
148
149 alpha_11 = ((E02[0] * E03[1]) - (E02[1] * E03[0])) / n[2];
150 beta_11 = ((E01[0] * E02[1]) - (E01[1] * E02[0])) / n[2];
151 }
152
153 // Compute the bilinear coordinates of the intersection point.
154 if (vtkm::Abs(alpha_11 - 1.0f) < vtkm::Epsilon<Precision>())
155 {
156
157 u = alpha;
158 if (vtkm::Abs(beta_11 - 1.0f) < vtkm::Epsilon<Precision>())
159 v = beta;
160 else
161 v = beta / ((u * (beta_11 - 1.0f)) + 1.0f);
162 }
163 else if (vtkm::Abs(beta_11 - 1.0) < vtkm::Epsilon<Precision>())
164 {
165
166 v = beta;
167 u = alpha / ((v * (alpha_11 - 1.0f)) + 1.0f);
168 }
169 else
170 {
171
172 Precision A = 1.0f - beta_11;
173 Precision B = (alpha * (beta_11 - 1.0f)) - (beta * (alpha_11 - 1.0f)) - 1.0f;
174 Precision C = alpha;
175 Precision D = (B * B) - (4.0f * A * C);
176 Precision QQ = -0.5f * (B + ((B < 0.0f ? -1.0f : 1.0f) * vtkm::Sqrt(D)));
177 u = QQ / A;
178 if ((u < 0.0f) || (u > 1.0f))
179 u = C / QQ;
180 v = beta / ((u * (beta_11 - 1.0f)) + 1.0f);
181 }
182
183 return true;
184 }
185
186 template <typename PointPortalType, typename LeafPortalType, typename Precision>
IntersectLeaf(const vtkm::Int32 & currentNode,const vtkm::Vec<Precision,3> & origin,const vtkm::Vec<Precision,3> & dir,const PointPortalType & points,vtkm::Id & hitIndex,Precision & closestDistance,Precision & minU,Precision & minV,LeafPortalType leafs,const Precision & minDistance) const187 VTKM_EXEC inline void IntersectLeaf(
188 const vtkm::Int32& currentNode,
189 const vtkm::Vec<Precision, 3>& origin,
190 const vtkm::Vec<Precision, 3>& dir,
191 const PointPortalType& points,
192 vtkm::Id& hitIndex,
193 Precision& closestDistance, // closest distance in this set of primitives
194 Precision& minU,
195 Precision& minV,
196 LeafPortalType leafs,
197 const Precision& minDistance) const // report intesections past this distance
198 {
199 const vtkm::Id quadCount = leafs.Get(currentNode);
200 for (vtkm::Id i = 1; i <= quadCount; ++i)
201 {
202 const vtkm::Id quadIndex = leafs.Get(currentNode + i);
203 if (quadIndex < QuadIds.GetNumberOfValues())
204 {
205 IdType pointIndex = QuadIds.Get(quadIndex);
206 Precision dist;
207 vtkm::Vec<Precision, 3> q, r, s, t;
208 q = vtkm::Vec<Precision, 3>(points.Get(pointIndex[1]));
209 r = vtkm::Vec<Precision, 3>(points.Get(pointIndex[2]));
210 s = vtkm::Vec<Precision, 3>(points.Get(pointIndex[3]));
211 t = vtkm::Vec<Precision, 3>(points.Get(pointIndex[4]));
212 Precision u, v;
213
214 bool ret = quad(origin, dir, q, r, s, t, u, v, dist);
215 if (ret)
216 {
217 if (dist < closestDistance && dist > minDistance)
218 {
219 //matid = vtkm::Vec<, 3>(points.Get(cur_offset + 2))[0];
220 closestDistance = dist;
221 hitIndex = quadIndex;
222 minU = u;
223 minV = v;
224 }
225 }
226 }
227 } // for
228 }
229 };
230
231 struct IntersectFunctor
232 {
233 template <typename Device, typename Precision>
operator ()vtkm::rendering::raytracing::detail::IntersectFunctor234 VTKM_CONT bool operator()(Device,
235 QuadIntersector* self,
236 Ray<Precision>& rays,
237 bool returnCellIndex)
238 {
239 VTKM_IS_DEVICE_ADAPTER_TAG(Device);
240 self->IntersectRaysImp(Device(), rays, returnCellIndex);
241 return true;
242 }
243 };
244
245 class CalculateNormals : public vtkm::worklet::WorkletMapField
246 {
247 public:
248 VTKM_CONT
CalculateNormals()249 CalculateNormals() {}
250 typedef void ControlSignature(FieldIn<>,
251 FieldIn<>,
252 FieldOut<>,
253 FieldOut<>,
254 FieldOut<>,
255 WholeArrayIn<Vec3RenderingTypes>,
256 WholeArrayIn<>);
257 typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6, _7);
258 template <typename Precision, typename PointPortalType, typename IndicesPortalType>
operator ()(const vtkm::Id & hitIndex,const vtkm::Vec<Precision,3> & rayDir,Precision & normalX,Precision & normalY,Precision & normalZ,const PointPortalType & points,const IndicesPortalType & indicesPortal) const259 VTKM_EXEC inline void operator()(const vtkm::Id& hitIndex,
260 const vtkm::Vec<Precision, 3>& rayDir,
261 Precision& normalX,
262 Precision& normalY,
263 Precision& normalZ,
264 const PointPortalType& points,
265 const IndicesPortalType& indicesPortal) const
266 {
267 if (hitIndex < 0)
268 return;
269
270 vtkm::Vec<vtkm::Id, 5> quadId = indicesPortal.Get(hitIndex);
271
272 vtkm::Vec<Precision, 3> a, b, c;
273 a = points.Get(quadId[1]);
274 b = points.Get(quadId[2]);
275 c = points.Get(quadId[3]);
276
277 vtkm::Vec<Precision, 3> normal = vtkm::TriangleNormal(a, b, c);
278 vtkm::Normalize(normal);
279
280 //flip the normal if its pointing the wrong way
281 if (vtkm::dot(normal, rayDir) > 0.f)
282 normal = -normal;
283
284 normalX = normal[0];
285 normalY = normal[1];
286 normalZ = normal[2];
287 }
288 }; //class CalculateNormals
289
290 template <typename Precision>
291 class GetScalar : public vtkm::worklet::WorkletMapField
292 {
293 private:
294 Precision MinScalar;
295 Precision invDeltaScalar;
296
297 public:
298 VTKM_CONT
GetScalar(const vtkm::Float32 & minScalar,const vtkm::Float32 & maxScalar)299 GetScalar(const vtkm::Float32& minScalar, const vtkm::Float32& maxScalar)
300 : MinScalar(minScalar)
301 {
302 //Make sure the we don't divide by zero on
303 //something like an iso-surface
304 if (maxScalar - MinScalar != 0.f)
305 invDeltaScalar = 1.f / (maxScalar - MinScalar);
306 else
307 invDeltaScalar = 1.f / minScalar;
308 }
309 typedef void ControlSignature(FieldIn<>,
310 FieldInOut<>,
311 WholeArrayIn<ScalarRenderingTypes>,
312 WholeArrayIn<>);
313 typedef void ExecutionSignature(_1, _2, _3, _4);
314 template <typename ScalarPortalType, typename IndicesPortalType>
operator ()(const vtkm::Id & hitIndex,Precision & scalar,const ScalarPortalType & scalars,const IndicesPortalType & indicesPortal) const315 VTKM_EXEC void operator()(const vtkm::Id& hitIndex,
316 Precision& scalar,
317 const ScalarPortalType& scalars,
318 const IndicesPortalType& indicesPortal) const
319 {
320 if (hitIndex < 0)
321 return;
322
323 //TODO: this should be interpolated?
324 vtkm::Vec<vtkm::Id, 5> pointId = indicesPortal.Get(hitIndex);
325
326 scalar = Precision(scalars.Get(pointId[0]));
327 //normalize
328 scalar = (scalar - MinScalar) * invDeltaScalar;
329 }
330 }; //class GetScalar
331
332 } // namespace detail
333
QuadIntersector()334 QuadIntersector::QuadIntersector()
335 : ShapeIntersector()
336 {
337 }
338
~QuadIntersector()339 QuadIntersector::~QuadIntersector()
340 {
341 }
342
343
IntersectRays(Ray<vtkm::Float32> & rays,bool returnCellIndex)344 void QuadIntersector::IntersectRays(Ray<vtkm::Float32>& rays, bool returnCellIndex)
345 {
346 vtkm::cont::TryExecute(detail::IntersectFunctor(), this, rays, returnCellIndex);
347 }
348
IntersectRays(Ray<vtkm::Float64> & rays,bool returnCellIndex)349 void QuadIntersector::IntersectRays(Ray<vtkm::Float64>& rays, bool returnCellIndex)
350 {
351 vtkm::cont::TryExecute(detail::IntersectFunctor(), this, rays, returnCellIndex);
352 }
353
354 template <typename Device, typename Precision>
IntersectRaysImp(Device,Ray<Precision> & rays,bool vtkmNotUsed (returnCellIndex))355 void QuadIntersector::IntersectRaysImp(Device,
356 Ray<Precision>& rays,
357 bool vtkmNotUsed(returnCellIndex))
358 {
359
360 detail::QuadLeafIntersector<Device> leafIntersector(this->QuadIds);
361
362 BVHTraverser<detail::QuadLeafIntersector> traverser;
363 traverser.IntersectRays(rays, this->BVH, leafIntersector, this->CoordsHandle, Device());
364
365 RayOperations::UpdateRayStatus(rays);
366 }
367
368 template <typename Precision>
IntersectionDataImp(Ray<Precision> & rays,const vtkm::cont::Field * scalarField,const vtkm::Range & scalarRange)369 void QuadIntersector::IntersectionDataImp(Ray<Precision>& rays,
370 const vtkm::cont::Field* scalarField,
371 const vtkm::Range& scalarRange)
372 {
373 ShapeIntersector::IntersectionPoint(rays);
374
375 // TODO: if this is nodes of a mesh, support points
376 bool isSupportedField =
377 (scalarField->GetAssociation() == vtkm::cont::Field::Association::POINTS ||
378 scalarField->GetAssociation() == vtkm::cont::Field::Association::CELL_SET);
379 if (!isSupportedField)
380 throw vtkm::cont::ErrorBadValue("Field not accociated with a cell set");
381
382 vtkm::worklet::DispatcherMapField<detail::CalculateNormals>(detail::CalculateNormals())
383 .Invoke(rays.HitIdx, rays.Dir, rays.NormalX, rays.NormalY, rays.NormalZ, CoordsHandle, QuadIds);
384
385 vtkm::worklet::DispatcherMapField<detail::GetScalar<Precision>>(
386 detail::GetScalar<Precision>(vtkm::Float32(scalarRange.Min), vtkm::Float32(scalarRange.Max)))
387 .Invoke(rays.HitIdx, rays.Scalar, *scalarField, QuadIds);
388 }
389
IntersectionData(Ray<vtkm::Float32> & rays,const vtkm::cont::Field * scalarField,const vtkm::Range & scalarRange)390 void QuadIntersector::IntersectionData(Ray<vtkm::Float32>& rays,
391 const vtkm::cont::Field* scalarField,
392 const vtkm::Range& scalarRange)
393 {
394 IntersectionDataImp(rays, scalarField, scalarRange);
395 }
396
IntersectionData(Ray<vtkm::Float64> & rays,const vtkm::cont::Field * scalarField,const vtkm::Range & scalarRange)397 void QuadIntersector::IntersectionData(Ray<vtkm::Float64>& rays,
398 const vtkm::cont::Field* scalarField,
399 const vtkm::Range& scalarRange)
400 {
401 IntersectionDataImp(rays, scalarField, scalarRange);
402 }
403
GetNumberOfShapes() const404 vtkm::Id QuadIntersector::GetNumberOfShapes() const
405 {
406 return QuadIds.GetNumberOfValues();
407 }
408 }
409 }
410 } //namespace vtkm::rendering::raytracing
411