1 //=============================================================================
2 //
3 // Copyright (c) Kitware, Inc.
4 // All rights reserved.
5 // See LICENSE.txt for details.
6 //
7 // This software is distributed WITHOUT ANY WARRANTY; without even
8 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
9 // PURPOSE. See the above copyright notice for more information.
10 //
11 // Copyright 2016 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
12 // Copyright 2016 UT-Battelle, LLC.
13 // Copyright 2016 Los Alamos National Security.
14 //
15 // Under the terms of Contract DE-NA0003525 with NTESS,
16 // the U.S. Government retains certain rights in this software.
17 // Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
18 // Laboratory (LANL), the U.S. Government retains certain rights in
19 // this software.
20 //
21 //=============================================================================
22 #include <vtkm/rendering/raytracing/VolumeRendererStructured.h>
23
24 #include <iostream>
25 #include <math.h>
26 #include <stdio.h>
27 #include <vtkm/cont/ArrayHandleCartesianProduct.h>
28 #include <vtkm/cont/ArrayHandleCounting.h>
29 #include <vtkm/cont/ArrayHandleUniformPointCoordinates.h>
30 #include <vtkm/cont/CellSetStructured.h>
31 #include <vtkm/cont/ColorTable.h>
32 #include <vtkm/cont/ErrorBadValue.h>
33 #include <vtkm/cont/Timer.h>
34 #include <vtkm/cont/TryExecute.h>
35 #include <vtkm/rendering/raytracing/Camera.h>
36 #include <vtkm/rendering/raytracing/Logger.h>
37 #include <vtkm/rendering/raytracing/Ray.h>
38 #include <vtkm/rendering/raytracing/RayTracingTypeDefs.h>
39 #include <vtkm/worklet/DispatcherMapField.h>
40 #include <vtkm/worklet/WorkletMapField.h>
41
42 namespace vtkm
43 {
44 namespace rendering
45 {
46 namespace raytracing
47 {
48
49 namespace
50 {
51
52 template <typename Device>
53 class RectilinearLocator
54 {
55 protected:
56 using DefaultHandle = vtkm::cont::ArrayHandle<vtkm::FloatDefault>;
57 using CartesianArrayHandle =
58 vtkm::cont::ArrayHandleCartesianProduct<DefaultHandle, DefaultHandle, DefaultHandle>;
59 using DefaultConstHandle = typename DefaultHandle::ExecutionTypes<Device>::PortalConst;
60 using CartesianConstPortal = typename CartesianArrayHandle::ExecutionTypes<Device>::PortalConst;
61
62 vtkm::Float32 InverseDeltaScalar;
63 DefaultConstHandle CoordPortals[3];
64 CartesianConstPortal Coordinates;
65 vtkm::exec::ConnectivityStructured<vtkm::TopologyElementTagPoint, vtkm::TopologyElementTagCell, 3>
66 Conn;
67 vtkm::Id3 PointDimensions;
68 vtkm::Vec<vtkm::Float32, 3> MinPoint;
69 vtkm::Vec<vtkm::Float32, 3> MaxPoint;
70
71 public:
RectilinearLocator(const CartesianArrayHandle & coordinates,vtkm::cont::CellSetStructured<3> & cellset)72 RectilinearLocator(const CartesianArrayHandle& coordinates,
73 vtkm::cont::CellSetStructured<3>& cellset)
74 : Coordinates(coordinates.PrepareForInput(Device()))
75 , Conn(cellset.PrepareForInput(Device(),
76 vtkm::TopologyElementTagPoint(),
77 vtkm::TopologyElementTagCell()))
78 {
79 CoordPortals[0] = Coordinates.GetFirstPortal();
80 CoordPortals[1] = Coordinates.GetSecondPortal();
81 CoordPortals[2] = Coordinates.GetThirdPortal();
82 PointDimensions = Conn.GetPointDimensions();
83 MinPoint[0] =
84 static_cast<vtkm::Float32>(coordinates.GetPortalConstControl().GetFirstPortal().Get(0));
85 MinPoint[1] =
86 static_cast<vtkm::Float32>(coordinates.GetPortalConstControl().GetSecondPortal().Get(0));
87 MinPoint[2] =
88 static_cast<vtkm::Float32>(coordinates.GetPortalConstControl().GetThirdPortal().Get(0));
89
90 MaxPoint[0] = static_cast<vtkm::Float32>(
91 coordinates.GetPortalConstControl().GetFirstPortal().Get(PointDimensions[0] - 1));
92 MaxPoint[1] = static_cast<vtkm::Float32>(
93 coordinates.GetPortalConstControl().GetSecondPortal().Get(PointDimensions[1] - 1));
94 MaxPoint[2] = static_cast<vtkm::Float32>(
95 coordinates.GetPortalConstControl().GetThirdPortal().Get(PointDimensions[2] - 1));
96 }
97
98 VTKM_EXEC
IsInside(const vtkm::Vec<vtkm::Float32,3> & point) const99 inline bool IsInside(const vtkm::Vec<vtkm::Float32, 3>& point) const
100 {
101 bool inside = true;
102 if (point[0] < MinPoint[0] || point[0] > MaxPoint[0])
103 inside = false;
104 if (point[1] < MinPoint[1] || point[1] > MaxPoint[1])
105 inside = false;
106 if (point[2] < MinPoint[2] || point[2] > MaxPoint[2])
107 inside = false;
108 return inside;
109 }
110
111 VTKM_EXEC
GetCellIndices(const vtkm::Vec<vtkm::Id,3> & cell,vtkm::Vec<vtkm::Id,8> & cellIndices) const112 inline void GetCellIndices(const vtkm::Vec<vtkm::Id, 3>& cell,
113 vtkm::Vec<vtkm::Id, 8>& cellIndices) const
114 {
115 cellIndices[0] = (cell[2] * PointDimensions[1] + cell[1]) * PointDimensions[0] + cell[0];
116 cellIndices[1] = cellIndices[0] + 1;
117 cellIndices[2] = cellIndices[1] + PointDimensions[0];
118 cellIndices[3] = cellIndices[2] - 1;
119 cellIndices[4] = cellIndices[0] + PointDimensions[0] * PointDimensions[1];
120 cellIndices[5] = cellIndices[4] + 1;
121 cellIndices[6] = cellIndices[5] + PointDimensions[0];
122 cellIndices[7] = cellIndices[6] - 1;
123 } // GetCellIndices
124
125 //
126 // Assumes point inside the data set
127 //
128 VTKM_EXEC
LocateCell(vtkm::Vec<vtkm::Id,3> & cell,const vtkm::Vec<vtkm::Float32,3> & point,vtkm::Vec<vtkm::Float32,3> & invSpacing) const129 inline void LocateCell(vtkm::Vec<vtkm::Id, 3>& cell,
130 const vtkm::Vec<vtkm::Float32, 3>& point,
131 vtkm::Vec<vtkm::Float32, 3>& invSpacing) const
132 {
133 for (vtkm::Int32 dim = 0; dim < 3; ++dim)
134 {
135 //
136 // When searching for points, we consider the max value of the cell
137 // to be apart of the next cell. If the point falls on the boundary of the
138 // data set, then it is technically inside a cell. This checks for that case
139 //
140 if (point[dim] == MaxPoint[dim])
141 {
142 cell[dim] = PointDimensions[dim] - 2;
143 continue;
144 }
145
146 bool found = false;
147 vtkm::Float32 minVal = static_cast<vtkm::Float32>(CoordPortals[dim].Get(cell[dim]));
148 const vtkm::Id searchDir = (point[dim] - minVal >= 0.f) ? 1 : -1;
149 vtkm::Float32 maxVal = static_cast<vtkm::Float32>(CoordPortals[dim].Get(cell[dim] + 1));
150
151 while (!found)
152 {
153 if (point[dim] >= minVal && point[dim] < maxVal)
154 {
155 found = true;
156 continue;
157 }
158
159 cell[dim] += searchDir;
160 vtkm::Id nextCellId = searchDir == 1 ? cell[dim] + 1 : cell[dim];
161 BOUNDS_CHECK(CoordPortals[dim], nextCellId);
162 vtkm::Float32 next = static_cast<vtkm::Float32>(CoordPortals[dim].Get(nextCellId));
163 if (searchDir == 1)
164 {
165 minVal = maxVal;
166 maxVal = next;
167 }
168 else
169 {
170 maxVal = minVal;
171 minVal = next;
172 }
173 }
174 invSpacing[dim] = 1.f / (maxVal - minVal);
175 }
176 } // LocateCell
177
178 VTKM_EXEC
GetCellIndex(const vtkm::Vec<vtkm::Id,3> & cell) const179 inline vtkm::Id GetCellIndex(const vtkm::Vec<vtkm::Id, 3>& cell) const
180 {
181 return (cell[2] * (PointDimensions[1] - 1) + cell[1]) * (PointDimensions[0] - 1) + cell[0];
182 }
183
184 VTKM_EXEC
GetPoint(const vtkm::Id & index,vtkm::Vec<vtkm::Float32,3> & point) const185 inline void GetPoint(const vtkm::Id& index, vtkm::Vec<vtkm::Float32, 3>& point) const
186 {
187 BOUNDS_CHECK(Coordinates, index);
188 point = Coordinates.Get(index);
189 }
190
191 VTKM_EXEC
GetMinPoint(const vtkm::Vec<vtkm::Id,3> & cell,vtkm::Vec<vtkm::Float32,3> & point) const192 inline void GetMinPoint(const vtkm::Vec<vtkm::Id, 3>& cell,
193 vtkm::Vec<vtkm::Float32, 3>& point) const
194 {
195 const vtkm::Id pointIndex =
196 (cell[2] * PointDimensions[1] + cell[1]) * PointDimensions[0] + cell[0];
197 point = Coordinates.Get(pointIndex);
198 }
199 }; // class RectilinearLocator
200
201 template <typename Device>
202 class UniformLocator
203 {
204 protected:
205 using UniformArrayHandle = vtkm::cont::ArrayHandleUniformPointCoordinates;
206 using UniformConstPortal = typename UniformArrayHandle::ExecutionTypes<Device>::PortalConst;
207
208 vtkm::Id3 PointDimensions;
209 vtkm::Vec<vtkm::Float32, 3> Origin;
210 vtkm::Vec<vtkm::Float32, 3> InvSpacing;
211 vtkm::Vec<vtkm::Float32, 3> MaxPoint;
212 UniformConstPortal Coordinates;
213 vtkm::exec::ConnectivityStructured<vtkm::TopologyElementTagPoint, vtkm::TopologyElementTagCell, 3>
214 Conn;
215
216 public:
UniformLocator(const UniformArrayHandle & coordinates,vtkm::cont::CellSetStructured<3> & cellset)217 UniformLocator(const UniformArrayHandle& coordinates, vtkm::cont::CellSetStructured<3>& cellset)
218 : Coordinates(coordinates.PrepareForInput(Device()))
219 , Conn(cellset.PrepareForInput(Device(),
220 vtkm::TopologyElementTagPoint(),
221 vtkm::TopologyElementTagCell()))
222 {
223 Origin = Coordinates.GetOrigin();
224 PointDimensions = Conn.GetPointDimensions();
225 vtkm::Vec<vtkm::Float32, 3> spacing = Coordinates.GetSpacing();
226
227 vtkm::Vec<vtkm::Float32, 3> unitLength;
228 unitLength[0] = static_cast<vtkm::Float32>(PointDimensions[0] - 1);
229 unitLength[1] = static_cast<vtkm::Float32>(PointDimensions[1] - 1);
230 unitLength[2] = static_cast<vtkm::Float32>(PointDimensions[2] - 1);
231 MaxPoint = Origin + spacing * unitLength;
232 InvSpacing[0] = 1.f / spacing[0];
233 InvSpacing[1] = 1.f / spacing[1];
234 InvSpacing[2] = 1.f / spacing[2];
235 }
236
237 VTKM_EXEC
IsInside(const vtkm::Vec<vtkm::Float32,3> & point) const238 inline bool IsInside(const vtkm::Vec<vtkm::Float32, 3>& point) const
239 {
240 bool inside = true;
241 if (point[0] < Origin[0] || point[0] > MaxPoint[0])
242 inside = false;
243 if (point[1] < Origin[1] || point[1] > MaxPoint[1])
244 inside = false;
245 if (point[2] < Origin[2] || point[2] > MaxPoint[2])
246 inside = false;
247 return inside;
248 }
249
250 VTKM_EXEC
GetCellIndices(const vtkm::Vec<vtkm::Id,3> & cell,vtkm::Vec<vtkm::Id,8> & cellIndices) const251 inline void GetCellIndices(const vtkm::Vec<vtkm::Id, 3>& cell,
252 vtkm::Vec<vtkm::Id, 8>& cellIndices) const
253 {
254 cellIndices[0] = (cell[2] * PointDimensions[1] + cell[1]) * PointDimensions[0] + cell[0];
255 cellIndices[1] = cellIndices[0] + 1;
256 cellIndices[2] = cellIndices[1] + PointDimensions[0];
257 cellIndices[3] = cellIndices[2] - 1;
258 cellIndices[4] = cellIndices[0] + PointDimensions[0] * PointDimensions[1];
259 cellIndices[5] = cellIndices[4] + 1;
260 cellIndices[6] = cellIndices[5] + PointDimensions[0];
261 cellIndices[7] = cellIndices[6] - 1;
262 } // GetCellIndices
263
264 VTKM_EXEC
GetCellIndex(const vtkm::Vec<vtkm::Id,3> & cell) const265 inline vtkm::Id GetCellIndex(const vtkm::Vec<vtkm::Id, 3>& cell) const
266 {
267 return (cell[2] * (PointDimensions[1] - 1) + cell[1]) * (PointDimensions[0] - 1) + cell[0];
268 }
269
270 VTKM_EXEC
LocateCell(vtkm::Vec<vtkm::Id,3> & cell,const vtkm::Vec<vtkm::Float32,3> & point,vtkm::Vec<vtkm::Float32,3> & invSpacing) const271 inline void LocateCell(vtkm::Vec<vtkm::Id, 3>& cell,
272 const vtkm::Vec<vtkm::Float32, 3>& point,
273 vtkm::Vec<vtkm::Float32, 3>& invSpacing) const
274 {
275 vtkm::Vec<vtkm::Float32, 3> temp = point;
276 temp = temp - Origin;
277 temp = temp * InvSpacing;
278 //make sure that if we border the upper edge, we sample the correct cell
279 if (temp[0] == vtkm::Float32(PointDimensions[0] - 1))
280 temp[0] = vtkm::Float32(PointDimensions[0] - 2);
281 if (temp[1] == vtkm::Float32(PointDimensions[1] - 1))
282 temp[1] = vtkm::Float32(PointDimensions[1] - 2);
283 if (temp[2] == vtkm::Float32(PointDimensions[2] - 1))
284 temp[2] = vtkm::Float32(PointDimensions[2] - 2);
285 cell = temp;
286 invSpacing = InvSpacing;
287 }
288
289 VTKM_EXEC
GetPoint(const vtkm::Id & index,vtkm::Vec<vtkm::Float32,3> & point) const290 inline void GetPoint(const vtkm::Id& index, vtkm::Vec<vtkm::Float32, 3>& point) const
291 {
292 BOUNDS_CHECK(Coordinates, index);
293 point = Coordinates.Get(index);
294 }
295
296 VTKM_EXEC
GetMinPoint(const vtkm::Vec<vtkm::Id,3> & cell,vtkm::Vec<vtkm::Float32,3> & point) const297 inline void GetMinPoint(const vtkm::Vec<vtkm::Id, 3>& cell,
298 vtkm::Vec<vtkm::Float32, 3>& point) const
299 {
300 const vtkm::Id pointIndex =
301 (cell[2] * PointDimensions[1] + cell[1]) * PointDimensions[0] + cell[0];
302 point = Coordinates.Get(pointIndex);
303 }
304
305 }; // class UniformLocator
306
307
308 } //namespace
309
310
311 template <typename Device, typename LocatorType>
312 class Sampler : public vtkm::worklet::WorkletMapField
313 {
314 private:
315 using ColorArrayHandle = typename vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 4>>;
316 using ColorArrayPortal = typename ColorArrayHandle::ExecutionTypes<Device>::PortalConst;
317 ColorArrayPortal ColorMap;
318 vtkm::Id ColorMapSize;
319 vtkm::Float32 MinScalar;
320 vtkm::Float32 SampleDistance;
321 vtkm::Float32 InverseDeltaScalar;
322 LocatorType Locator;
323
324 public:
325 VTKM_CONT
Sampler(const ColorArrayHandle & colorMap,const vtkm::Float32 & minScalar,const vtkm::Float32 & maxScalar,const vtkm::Float32 & sampleDistance,const LocatorType & locator)326 Sampler(const ColorArrayHandle& colorMap,
327 const vtkm::Float32& minScalar,
328 const vtkm::Float32& maxScalar,
329 const vtkm::Float32& sampleDistance,
330 const LocatorType& locator)
331 : ColorMap(colorMap.PrepareForInput(Device()))
332 , MinScalar(minScalar)
333 , SampleDistance(sampleDistance)
334 , Locator(locator)
335 {
336 ColorMapSize = colorMap.GetNumberOfValues() - 1;
337 if ((maxScalar - minScalar) != 0.f)
338 InverseDeltaScalar = 1.f / (maxScalar - minScalar);
339 else
340 InverseDeltaScalar = minScalar;
341 }
342 using ControlSignature = void(FieldIn<>,
343 FieldIn<>,
344 FieldIn<>,
345 FieldIn<>,
346 WholeArrayInOut<>,
347 WholeArrayIn<ScalarRenderingTypes>);
348 using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, WorkIndex);
349
350 template <typename ScalarPortalType, typename ColorBufferType>
operator ()(const vtkm::Vec<vtkm::Float32,3> & rayDir,const vtkm::Vec<vtkm::Float32,3> & rayOrigin,const vtkm::Float32 & minDistance,const vtkm::Float32 & maxDistance,ColorBufferType & colorBuffer,ScalarPortalType & scalars,const vtkm::Id & pixelIndex) const351 VTKM_EXEC void operator()(const vtkm::Vec<vtkm::Float32, 3>& rayDir,
352 const vtkm::Vec<vtkm::Float32, 3>& rayOrigin,
353 const vtkm::Float32& minDistance,
354 const vtkm::Float32& maxDistance,
355 ColorBufferType& colorBuffer,
356 ScalarPortalType& scalars,
357 const vtkm::Id& pixelIndex) const
358 {
359 vtkm::Vec<vtkm::Float32, 4> color;
360 BOUNDS_CHECK(colorBuffer, pixelIndex * 4 + 0);
361 color[0] = colorBuffer.Get(pixelIndex * 4 + 0);
362 BOUNDS_CHECK(colorBuffer, pixelIndex * 4 + 1);
363 color[1] = colorBuffer.Get(pixelIndex * 4 + 1);
364 BOUNDS_CHECK(colorBuffer, pixelIndex * 4 + 2);
365 color[2] = colorBuffer.Get(pixelIndex * 4 + 2);
366 BOUNDS_CHECK(colorBuffer, pixelIndex * 4 + 3);
367 color[3] = colorBuffer.Get(pixelIndex * 4 + 3);
368
369 if (minDistance == -1.f)
370 {
371 return; //TODO: Compact? or just image subset...
372 }
373 //get the initial sample position;
374 vtkm::Vec<vtkm::Float32, 3> sampleLocation;
375 // find the distance to the first sample
376 vtkm::Float32 distance = minDistance + 0.0001f;
377 sampleLocation = rayOrigin + distance * rayDir;
378 // since the calculations are slightly different, we could hit an
379 // edge case where the first sample location may not be in the data set.
380 // Thus, advance to the next sample location
381 while (!Locator.IsInside(sampleLocation) && distance < maxDistance)
382 {
383 distance += SampleDistance;
384 sampleLocation = rayOrigin + distance * rayDir;
385 }
386 /*
387 7----------6
388 /| /|
389 4----------5 |
390 | | | |
391 | 3--------|-2 z y
392 |/ |/ |/
393 0----------1 |__ x
394 */
395 vtkm::Vec<vtkm::Float32, 3> bottomLeft(0, 0, 0);
396 bool newCell = true;
397 //check to see if we left the cell
398 vtkm::Float32 tx = 0.f;
399 vtkm::Float32 ty = 0.f;
400 vtkm::Float32 tz = 0.f;
401 vtkm::Float32 scalar0 = 0.f;
402 vtkm::Float32 scalar1minus0 = 0.f;
403 vtkm::Float32 scalar2minus3 = 0.f;
404 vtkm::Float32 scalar3 = 0.f;
405 vtkm::Float32 scalar4 = 0.f;
406 vtkm::Float32 scalar5minus4 = 0.f;
407 vtkm::Float32 scalar6minus7 = 0.f;
408 vtkm::Float32 scalar7 = 0.f;
409
410 vtkm::Vec<vtkm::Id, 3> cell(0, 0, 0);
411 vtkm::Vec<vtkm::Float32, 3> invSpacing(0.f, 0.f, 0.f);
412
413
414 while (Locator.IsInside(sampleLocation) && distance < maxDistance)
415 {
416 vtkm::Float32 mint = vtkm::Min(tx, vtkm::Min(ty, tz));
417 vtkm::Float32 maxt = vtkm::Max(tx, vtkm::Max(ty, tz));
418 if (maxt > 1.f || mint < 0.f)
419 newCell = true;
420
421 if (newCell)
422 {
423
424 vtkm::Vec<vtkm::Id, 8> cellIndices;
425 Locator.LocateCell(cell, sampleLocation, invSpacing);
426 Locator.GetCellIndices(cell, cellIndices);
427 Locator.GetPoint(cellIndices[0], bottomLeft);
428
429 scalar0 = vtkm::Float32(scalars.Get(cellIndices[0]));
430 vtkm::Float32 scalar1 = vtkm::Float32(scalars.Get(cellIndices[1]));
431 vtkm::Float32 scalar2 = vtkm::Float32(scalars.Get(cellIndices[2]));
432 scalar3 = vtkm::Float32(scalars.Get(cellIndices[3]));
433 scalar4 = vtkm::Float32(scalars.Get(cellIndices[4]));
434 vtkm::Float32 scalar5 = vtkm::Float32(scalars.Get(cellIndices[5]));
435 vtkm::Float32 scalar6 = vtkm::Float32(scalars.Get(cellIndices[6]));
436 scalar7 = vtkm::Float32(scalars.Get(cellIndices[7]));
437
438 // save ourselves a couple extra instructions
439 scalar6minus7 = scalar6 - scalar7;
440 scalar5minus4 = scalar5 - scalar4;
441 scalar1minus0 = scalar1 - scalar0;
442 scalar2minus3 = scalar2 - scalar3;
443
444 tx = (sampleLocation[0] - bottomLeft[0]) * invSpacing[0];
445 ty = (sampleLocation[1] - bottomLeft[1]) * invSpacing[1];
446 tz = (sampleLocation[2] - bottomLeft[2]) * invSpacing[2];
447
448 newCell = false;
449 }
450
451 vtkm::Float32 lerped76 = scalar7 + tx * scalar6minus7;
452 vtkm::Float32 lerped45 = scalar4 + tx * scalar5minus4;
453 vtkm::Float32 lerpedTop = lerped45 + ty * (lerped76 - lerped45);
454
455 vtkm::Float32 lerped01 = scalar0 + tx * scalar1minus0;
456 vtkm::Float32 lerped32 = scalar3 + tx * scalar2minus3;
457 vtkm::Float32 lerpedBottom = lerped01 + ty * (lerped32 - lerped01);
458
459 vtkm::Float32 finalScalar = lerpedBottom + tz * (lerpedTop - lerpedBottom);
460 //normalize scalar
461 finalScalar = (finalScalar - MinScalar) * InverseDeltaScalar;
462
463 vtkm::Id colorIndex =
464 static_cast<vtkm::Id>(finalScalar * static_cast<vtkm::Float32>(ColorMapSize));
465 if (colorIndex < 0)
466 colorIndex = 0;
467 if (colorIndex > ColorMapSize)
468 colorIndex = ColorMapSize;
469
470 vtkm::Vec<vtkm::Float32, 4> sampleColor = ColorMap.Get(colorIndex);
471
472 //composite
473 sampleColor[3] *= (1.f - color[3]);
474 color[0] = color[0] + sampleColor[0] * sampleColor[3];
475 color[1] = color[1] + sampleColor[1] * sampleColor[3];
476 color[2] = color[2] + sampleColor[2] * sampleColor[3];
477 color[3] = sampleColor[3] + color[3];
478 //advance
479 distance += SampleDistance;
480 sampleLocation = sampleLocation + SampleDistance * rayDir;
481
482 //this is linear could just do an addition
483 tx = (sampleLocation[0] - bottomLeft[0]) * invSpacing[0];
484 ty = (sampleLocation[1] - bottomLeft[1]) * invSpacing[1];
485 tz = (sampleLocation[2] - bottomLeft[2]) * invSpacing[2];
486
487 if (color[3] >= 1.f)
488 break;
489 }
490
491 color[0] = vtkm::Min(color[0], 1.f);
492 color[1] = vtkm::Min(color[1], 1.f);
493 color[2] = vtkm::Min(color[2], 1.f);
494 color[3] = vtkm::Min(color[3], 1.f);
495
496 BOUNDS_CHECK(colorBuffer, pixelIndex * 4 + 0);
497 colorBuffer.Set(pixelIndex * 4 + 0, color[0]);
498 BOUNDS_CHECK(colorBuffer, pixelIndex * 4 + 1);
499 colorBuffer.Set(pixelIndex * 4 + 1, color[1]);
500 BOUNDS_CHECK(colorBuffer, pixelIndex * 4 + 2);
501 colorBuffer.Set(pixelIndex * 4 + 2, color[2]);
502 BOUNDS_CHECK(colorBuffer, pixelIndex * 4 + 3);
503 colorBuffer.Set(pixelIndex * 4 + 3, color[3]);
504 }
505 }; //Sampler
506
507 template <typename Device, typename LocatorType>
508 class SamplerCellAssoc : public vtkm::worklet::WorkletMapField
509 {
510 private:
511 using ColorArrayHandle = typename vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 4>>;
512 using ColorArrayPortal = typename ColorArrayHandle::ExecutionTypes<Device>::PortalConst;
513 ColorArrayPortal ColorMap;
514 vtkm::Id ColorMapSize;
515 vtkm::Float32 MinScalar;
516 vtkm::Float32 SampleDistance;
517 vtkm::Float32 InverseDeltaScalar;
518 LocatorType Locator;
519
520 public:
521 VTKM_CONT
SamplerCellAssoc(const ColorArrayHandle & colorMap,const vtkm::Float32 & minScalar,const vtkm::Float32 & maxScalar,const vtkm::Float32 & sampleDistance,const LocatorType & locator)522 SamplerCellAssoc(const ColorArrayHandle& colorMap,
523 const vtkm::Float32& minScalar,
524 const vtkm::Float32& maxScalar,
525 const vtkm::Float32& sampleDistance,
526 const LocatorType& locator)
527 : ColorMap(colorMap.PrepareForInput(Device()))
528 , MinScalar(minScalar)
529 , SampleDistance(sampleDistance)
530 , Locator(locator)
531 {
532 ColorMapSize = colorMap.GetNumberOfValues() - 1;
533 if ((maxScalar - minScalar) != 0.f)
534 InverseDeltaScalar = 1.f / (maxScalar - minScalar);
535 else
536 InverseDeltaScalar = minScalar;
537 }
538 using ControlSignature = void(FieldIn<>,
539 FieldIn<>,
540 FieldIn<>,
541 FieldIn<>,
542 WholeArrayInOut<>,
543 WholeArrayIn<ScalarRenderingTypes>);
544 using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, WorkIndex);
545
546 template <typename ScalarPortalType, typename ColorBufferType>
operator ()(const vtkm::Vec<vtkm::Float32,3> & rayDir,const vtkm::Vec<vtkm::Float32,3> & rayOrigin,const vtkm::Float32 & minDistance,const vtkm::Float32 & maxDistance,ColorBufferType & colorBuffer,const ScalarPortalType & scalars,const vtkm::Id & pixelIndex) const547 VTKM_EXEC void operator()(const vtkm::Vec<vtkm::Float32, 3>& rayDir,
548 const vtkm::Vec<vtkm::Float32, 3>& rayOrigin,
549 const vtkm::Float32& minDistance,
550 const vtkm::Float32& maxDistance,
551 ColorBufferType& colorBuffer,
552 const ScalarPortalType& scalars,
553 const vtkm::Id& pixelIndex) const
554 {
555 vtkm::Vec<vtkm::Float32, 4> color;
556 BOUNDS_CHECK(colorBuffer, pixelIndex * 4 + 0);
557 color[0] = colorBuffer.Get(pixelIndex * 4 + 0);
558 BOUNDS_CHECK(colorBuffer, pixelIndex * 4 + 1);
559 color[1] = colorBuffer.Get(pixelIndex * 4 + 1);
560 BOUNDS_CHECK(colorBuffer, pixelIndex * 4 + 2);
561 color[2] = colorBuffer.Get(pixelIndex * 4 + 2);
562 BOUNDS_CHECK(colorBuffer, pixelIndex * 4 + 3);
563 color[3] = colorBuffer.Get(pixelIndex * 4 + 3);
564
565 if (minDistance == -1.f)
566 return; //TODO: Compact? or just image subset...
567 //get the initial sample position;
568 vtkm::Vec<vtkm::Float32, 3> sampleLocation;
569 // find the distance to the first sample
570 vtkm::Float32 distance = minDistance + 0.0001f;
571 sampleLocation = rayOrigin + distance * rayDir;
572 // since the calculations are slightly different, we could hit an
573 // edge case where the first sample location may not be in the data set.
574 // Thus, advance to the next sample location
575 while (!Locator.IsInside(sampleLocation) && distance < maxDistance)
576 {
577 distance += SampleDistance;
578 sampleLocation = rayOrigin + distance * rayDir;
579 }
580
581 /*
582 7----------6
583 /| /|
584 4----------5 |
585 | | | |
586 | 3--------|-2 z y
587 |/ |/ |/
588 0----------1 |__ x
589 */
590 bool newCell = true;
591 vtkm::Float32 tx = 2.f;
592 vtkm::Float32 ty = 2.f;
593 vtkm::Float32 tz = 2.f;
594 vtkm::Float32 scalar0 = 0.f;
595 vtkm::Vec<vtkm::Float32, 4> sampleColor(0.f, 0.f, 0.f, 0.f);
596 vtkm::Vec<vtkm::Float32, 3> bottomLeft(0.f, 0.f, 0.f);
597 vtkm::Vec<vtkm::Float32, 3> invSpacing(0.f, 0.f, 0.f);
598 vtkm::Vec<vtkm::Id, 3> cell(0, 0, 0);
599 while (Locator.IsInside(sampleLocation) && distance < maxDistance)
600 {
601 vtkm::Float32 mint = vtkm::Min(tx, vtkm::Min(ty, tz));
602 vtkm::Float32 maxt = vtkm::Max(tx, vtkm::Max(ty, tz));
603 if (maxt > 1.f || mint < 0.f)
604 newCell = true;
605 if (newCell)
606 {
607 Locator.LocateCell(cell, sampleLocation, invSpacing);
608 vtkm::Id cellId = Locator.GetCellIndex(cell);
609
610 scalar0 = vtkm::Float32(scalars.Get(cellId));
611 vtkm::Float32 normalizedScalar = (scalar0 - MinScalar) * InverseDeltaScalar;
612 vtkm::Id colorIndex =
613 static_cast<vtkm::Id>(normalizedScalar * static_cast<vtkm::Float32>(ColorMapSize));
614 if (colorIndex < 0)
615 colorIndex = 0;
616 if (colorIndex > ColorMapSize)
617 colorIndex = ColorMapSize;
618 sampleColor = ColorMap.Get(colorIndex);
619 Locator.GetMinPoint(cell, bottomLeft);
620 tx = (sampleLocation[0] - bottomLeft[0]) * invSpacing[0];
621 ty = (sampleLocation[1] - bottomLeft[1]) * invSpacing[1];
622 tz = (sampleLocation[2] - bottomLeft[2]) * invSpacing[2];
623 newCell = false;
624 }
625
626 // just repeatably composite
627 vtkm::Float32 alpha = sampleColor[3] * (1.f - color[3]);
628 color[0] = color[0] + sampleColor[0] * alpha;
629 color[1] = color[1] + sampleColor[1] * alpha;
630 color[2] = color[2] + sampleColor[2] * alpha;
631 color[3] = alpha + color[3];
632 //advance
633 distance += SampleDistance;
634 sampleLocation = sampleLocation + SampleDistance * rayDir;
635
636 if (color[3] >= 1.f)
637 break;
638 tx = (sampleLocation[0] - bottomLeft[0]) * invSpacing[0];
639 ty = (sampleLocation[1] - bottomLeft[1]) * invSpacing[1];
640 tz = (sampleLocation[2] - bottomLeft[2]) * invSpacing[2];
641 }
642 color[0] = vtkm::Min(color[0], 1.f);
643 color[1] = vtkm::Min(color[1], 1.f);
644 color[2] = vtkm::Min(color[2], 1.f);
645 color[3] = vtkm::Min(color[3], 1.f);
646
647 BOUNDS_CHECK(colorBuffer, pixelIndex * 4 + 0);
648 colorBuffer.Set(pixelIndex * 4 + 0, color[0]);
649 BOUNDS_CHECK(colorBuffer, pixelIndex * 4 + 1);
650 colorBuffer.Set(pixelIndex * 4 + 1, color[1]);
651 BOUNDS_CHECK(colorBuffer, pixelIndex * 4 + 2);
652 colorBuffer.Set(pixelIndex * 4 + 2, color[2]);
653 BOUNDS_CHECK(colorBuffer, pixelIndex * 4 + 3);
654 colorBuffer.Set(pixelIndex * 4 + 3, color[3]);
655 }
656 }; //SamplerCell
657
658 class CalcRayStart : public vtkm::worklet::WorkletMapField
659 {
660 vtkm::Float32 Xmin;
661 vtkm::Float32 Ymin;
662 vtkm::Float32 Zmin;
663 vtkm::Float32 Xmax;
664 vtkm::Float32 Ymax;
665 vtkm::Float32 Zmax;
666
667 public:
668 VTKM_CONT
CalcRayStart(const vtkm::Bounds boundingBox)669 CalcRayStart(const vtkm::Bounds boundingBox)
670 {
671 Xmin = static_cast<vtkm::Float32>(boundingBox.X.Min);
672 Xmax = static_cast<vtkm::Float32>(boundingBox.X.Max);
673 Ymin = static_cast<vtkm::Float32>(boundingBox.Y.Min);
674 Ymax = static_cast<vtkm::Float32>(boundingBox.Y.Max);
675 Zmin = static_cast<vtkm::Float32>(boundingBox.Z.Min);
676 Zmax = static_cast<vtkm::Float32>(boundingBox.Z.Max);
677 }
678
679 VTKM_EXEC
rcp(vtkm::Float32 f) const680 vtkm::Float32 rcp(vtkm::Float32 f) const { return 1.0f / f; }
681
682 VTKM_EXEC
rcp_safe(vtkm::Float32 f) const683 vtkm::Float32 rcp_safe(vtkm::Float32 f) const { return rcp((fabs(f) < 1e-8f) ? 1e-8f : f); }
684
685 using ControlSignature = void(FieldIn<>, FieldOut<>, FieldInOut<>, FieldInOut<>, FieldIn<>);
686 using ExecutionSignature = void(_1, _2, _3, _4, _5);
687 template <typename Precision>
operator ()(const vtkm::Vec<Precision,3> & rayDir,vtkm::Float32 & minDistance,vtkm::Float32 & distance,vtkm::Float32 & maxDistance,const vtkm::Vec<Precision,3> & rayOrigin) const688 VTKM_EXEC void operator()(const vtkm::Vec<Precision, 3>& rayDir,
689 vtkm::Float32& minDistance,
690 vtkm::Float32& distance,
691 vtkm::Float32& maxDistance,
692 const vtkm::Vec<Precision, 3>& rayOrigin) const
693 {
694 vtkm::Float32 dirx = static_cast<vtkm::Float32>(rayDir[0]);
695 vtkm::Float32 diry = static_cast<vtkm::Float32>(rayDir[1]);
696 vtkm::Float32 dirz = static_cast<vtkm::Float32>(rayDir[2]);
697 vtkm::Float32 origx = static_cast<vtkm::Float32>(rayOrigin[0]);
698 vtkm::Float32 origy = static_cast<vtkm::Float32>(rayOrigin[1]);
699 vtkm::Float32 origz = static_cast<vtkm::Float32>(rayOrigin[2]);
700
701 vtkm::Float32 invDirx = rcp_safe(dirx);
702 vtkm::Float32 invDiry = rcp_safe(diry);
703 vtkm::Float32 invDirz = rcp_safe(dirz);
704
705 vtkm::Float32 odirx = origx * invDirx;
706 vtkm::Float32 odiry = origy * invDiry;
707 vtkm::Float32 odirz = origz * invDirz;
708
709 vtkm::Float32 xmin = Xmin * invDirx - odirx;
710 vtkm::Float32 ymin = Ymin * invDiry - odiry;
711 vtkm::Float32 zmin = Zmin * invDirz - odirz;
712 vtkm::Float32 xmax = Xmax * invDirx - odirx;
713 vtkm::Float32 ymax = Ymax * invDiry - odiry;
714 vtkm::Float32 zmax = Zmax * invDirz - odirz;
715
716
717 minDistance = vtkm::Max(
718 vtkm::Max(vtkm::Max(vtkm::Min(ymin, ymax), vtkm::Min(xmin, xmax)), vtkm::Min(zmin, zmax)),
719 minDistance);
720 vtkm::Float32 exitDistance =
721 vtkm::Min(vtkm::Min(vtkm::Max(ymin, ymax), vtkm::Max(xmin, xmax)), vtkm::Max(zmin, zmax));
722 maxDistance = vtkm::Min(maxDistance, exitDistance);
723 if (maxDistance < minDistance)
724 {
725 minDistance = -1.f; //flag for miss
726 }
727 else
728 {
729 distance = minDistance;
730 }
731 }
732 }; //class CalcRayStart
733
VolumeRendererStructured()734 VolumeRendererStructured::VolumeRendererStructured()
735 {
736 IsSceneDirty = false;
737 IsUniformDataSet = true;
738 SampleDistance = -1.f;
739 }
740
SetColorMap(const vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32,4>> & colorMap)741 void VolumeRendererStructured::SetColorMap(
742 const vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 4>>& colorMap)
743 {
744 ColorMap = colorMap;
745 }
746
SetData(const vtkm::cont::CoordinateSystem & coords,const vtkm::cont::Field & scalarField,const vtkm::cont::CellSetStructured<3> & cellset,const vtkm::Range & scalarRange)747 void VolumeRendererStructured::SetData(const vtkm::cont::CoordinateSystem& coords,
748 const vtkm::cont::Field& scalarField,
749 const vtkm::cont::CellSetStructured<3>& cellset,
750 const vtkm::Range& scalarRange)
751 {
752 if (coords.GetData().IsSameType(CartesianArrayHandle()))
753 IsUniformDataSet = false;
754 IsSceneDirty = true;
755 SpatialExtent = coords.GetBounds();
756 Coordinates = coords.GetData();
757 ScalarField = &scalarField;
758 Cellset = cellset;
759 ScalarRange = scalarRange;
760 }
761
762 template <typename Precision>
763 struct VolumeRendererStructured::RenderFunctor
764 {
765 protected:
766 vtkm::rendering::raytracing::VolumeRendererStructured* Self;
767 vtkm::rendering::raytracing::Ray<Precision>& Rays;
768
769 public:
770 VTKM_CONT
RenderFunctorvtkm::rendering::raytracing::VolumeRendererStructured::RenderFunctor771 RenderFunctor(vtkm::rendering::raytracing::VolumeRendererStructured* self,
772 vtkm::rendering::raytracing::Ray<Precision>& rays)
773 : Self(self)
774 , Rays(rays)
775 {
776 }
777
778 template <typename Device>
operator ()vtkm::rendering::raytracing::VolumeRendererStructured::RenderFunctor779 VTKM_CONT bool operator()(Device)
780 {
781 VTKM_IS_DEVICE_ADAPTER_TAG(Device);
782
783 this->Self->RenderOnDevice(this->Rays, Device());
784 return true;
785 }
786 };
787
Render(vtkm::rendering::raytracing::Ray<vtkm::Float32> & rays)788 void VolumeRendererStructured::Render(vtkm::rendering::raytracing::Ray<vtkm::Float32>& rays)
789 {
790 RenderFunctor<vtkm::Float32> functor(this, rays);
791 vtkm::cont::TryExecute(functor);
792 }
793
794 //void
795 //VolumeRendererStructured::Render(vtkm::rendering::raytracing::Ray<vtkm::Float64>& rays)
796 //{
797 // RenderFunctor<vtkm::Float64> functor(this, rays);
798 // vtkm::cont::TryExecute(functor);
799 //}
800
801 template <typename Precision, typename Device>
RenderOnDevice(vtkm::rendering::raytracing::Ray<Precision> & rays,Device)802 void VolumeRendererStructured::RenderOnDevice(vtkm::rendering::raytracing::Ray<Precision>& rays,
803 Device)
804 {
805 vtkm::cont::Timer<Device> renderTimer;
806 Logger* logger = Logger::GetInstance();
807 logger->OpenLogEntry("volume_render_structured");
808 logger->AddLogData("device", GetDeviceString(Device()));
809
810 if (SampleDistance <= 0.f)
811 {
812 vtkm::Vec<vtkm::Float32, 3> extent;
813 extent[0] = static_cast<vtkm::Float32>(this->SpatialExtent.X.Length());
814 extent[1] = static_cast<vtkm::Float32>(this->SpatialExtent.Y.Length());
815 extent[2] = static_cast<vtkm::Float32>(this->SpatialExtent.Z.Length());
816 const vtkm::Float32 defaultNumberOfSamples = 200.f;
817 SampleDistance = vtkm::Magnitude(extent) / defaultNumberOfSamples;
818 }
819
820 vtkm::cont::Timer<Device> timer;
821 vtkm::worklet::DispatcherMapField<CalcRayStart> calcRayStartDispatcher(
822 CalcRayStart(this->SpatialExtent));
823 calcRayStartDispatcher.SetDevice(Device());
824 calcRayStartDispatcher.Invoke(
825 rays.Dir, rays.MinDistance, rays.Distance, rays.MaxDistance, rays.Origin);
826
827 vtkm::Float64 time = timer.GetElapsedTime();
828 logger->AddLogData("calc_ray_start", time);
829 timer.Reset();
830
831 bool isSupportedField =
832 (ScalarField->GetAssociation() == vtkm::cont::Field::Association::POINTS ||
833 ScalarField->GetAssociation() == vtkm::cont::Field::Association::CELL_SET);
834 if (!isSupportedField)
835 throw vtkm::cont::ErrorBadValue("Field not accociated with cell set or points");
836 bool isAssocPoints = ScalarField->GetAssociation() == vtkm::cont::Field::Association::POINTS;
837
838 if (IsUniformDataSet)
839 {
840 vtkm::cont::ArrayHandleUniformPointCoordinates vertices;
841 vertices = Coordinates.Cast<vtkm::cont::ArrayHandleUniformPointCoordinates>();
842 UniformLocator<Device> locator(vertices, Cellset);
843
844 if (isAssocPoints)
845 {
846 vtkm::worklet::DispatcherMapField<Sampler<Device, UniformLocator<Device>>> samplerDispatcher(
847 Sampler<Device, UniformLocator<Device>>(ColorMap,
848 vtkm::Float32(ScalarRange.Min),
849 vtkm::Float32(ScalarRange.Max),
850 SampleDistance,
851 locator));
852 samplerDispatcher.SetDevice(Device());
853 samplerDispatcher.Invoke(rays.Dir,
854 rays.Origin,
855 rays.MinDistance,
856 rays.MaxDistance,
857 rays.Buffers.at(0).Buffer,
858 *ScalarField);
859 }
860 else
861 {
862 vtkm::worklet::DispatcherMapField<SamplerCellAssoc<Device, UniformLocator<Device>>>(
863 SamplerCellAssoc<Device, UniformLocator<Device>>(ColorMap,
864 vtkm::Float32(ScalarRange.Min),
865 vtkm::Float32(ScalarRange.Max),
866 SampleDistance,
867 locator))
868 .Invoke(rays.Dir,
869 rays.Origin,
870 rays.MinDistance,
871 rays.MaxDistance,
872 rays.Buffers.at(0).Buffer,
873 *ScalarField);
874 }
875 }
876 else
877 {
878 CartesianArrayHandle vertices;
879 vertices = Coordinates.Cast<CartesianArrayHandle>();
880 RectilinearLocator<Device> locator(vertices, Cellset);
881 if (isAssocPoints)
882 {
883 vtkm::worklet::DispatcherMapField<Sampler<Device, RectilinearLocator<Device>>>
884 samplerDispatcher(
885 Sampler<Device, RectilinearLocator<Device>>(ColorMap,
886 vtkm::Float32(ScalarRange.Min),
887 vtkm::Float32(ScalarRange.Max),
888 SampleDistance,
889 locator));
890 samplerDispatcher.SetDevice(Device());
891 samplerDispatcher.Invoke(rays.Dir,
892 rays.Origin,
893 rays.MinDistance,
894 rays.MaxDistance,
895 rays.Buffers.at(0).Buffer,
896 *ScalarField);
897 }
898 else
899 {
900 vtkm::worklet::DispatcherMapField<SamplerCellAssoc<Device, RectilinearLocator<Device>>>
901 rectilinearLocatorDispatcher(
902 SamplerCellAssoc<Device, RectilinearLocator<Device>>(ColorMap,
903 vtkm::Float32(ScalarRange.Min),
904 vtkm::Float32(ScalarRange.Max),
905 SampleDistance,
906 locator));
907 rectilinearLocatorDispatcher.SetDevice(Device());
908 rectilinearLocatorDispatcher.Invoke(rays.Dir,
909 rays.Origin,
910 rays.MinDistance,
911 rays.MaxDistance,
912 rays.Buffers.at(0).Buffer,
913 *ScalarField);
914 }
915 }
916
917 time = timer.GetElapsedTime();
918 logger->AddLogData("sample", time);
919 timer.Reset();
920
921 time = renderTimer.GetElapsedTime();
922 logger->CloseLogEntry(time);
923 } //Render
924
SetSampleDistance(const vtkm::Float32 & distance)925 void VolumeRendererStructured::SetSampleDistance(const vtkm::Float32& distance)
926 {
927 if (distance <= 0.f)
928 throw vtkm::cont::ErrorBadValue("Sample distance must be positive.");
929 SampleDistance = distance;
930 }
931 }
932 }
933 } //namespace vtkm::rendering::raytracing
934