1 // Copyright 2019-2021 Intel Corporation
2 // SPDX-License-Identifier: Apache-2.0
3
4 #include "UnstructuredVolume.h"
5 #include <algorithm>
6 #include "../common/Data.h"
7 #include "UnstructuredSampler.h"
8 #include "rkcommon/containers/AlignedVector.h"
9 #include "rkcommon/tasking/parallel_for.h"
10
11 // Map cell type to its vertices count
getVerticesCount(uint8_t cellType)12 inline uint32_t getVerticesCount(uint8_t cellType)
13 {
14 switch (cellType) {
15 case VKL_TETRAHEDRON:
16 return 4;
17 case VKL_HEXAHEDRON:
18 return 8;
19 case VKL_WEDGE:
20 return 6;
21 case VKL_PYRAMID:
22 return 5;
23 }
24
25 // Unknown cell type
26 return 1;
27 }
28
29 namespace openvkl {
30 namespace cpu_device {
31
tabIndent(int indent)32 static void tabIndent(int indent)
33 {
34 for (int i = 0; i < indent; i++)
35 std::cerr << "\t";
36 }
37
dumpBVH(Node * root,int indent=0)38 static void dumpBVH(Node *root, int indent = 0)
39 {
40 if (root->nominalLength.x < 0) {
41 auto leaf = (LeafNode *)root;
42 tabIndent(indent);
43 std::cerr << "id: " << leaf->cellID << " bounds: " << leaf->bounds
44 << " range: " << leaf->valueRange
45 << " nom: " << leaf->nominalLength << std::endl;
46 } else {
47 auto inner = (InnerNode *)root;
48 tabIndent(indent);
49 std::cerr << "range: " << inner->valueRange
50 << " nom: " << inner->nominalLength << std::endl;
51
52 tabIndent(indent);
53 std::cerr << "bounds[0]: " << inner->bounds[0] << std::endl;
54 dumpBVH(inner->children[0], indent + 1);
55
56 tabIndent(indent);
57 std::cerr << "bounds[1]: " << inner->bounds[1] << std::endl;
58 dumpBVH(inner->children[1], indent + 1);
59 }
60 }
61
62 template <int W>
~UnstructuredVolume()63 UnstructuredVolume<W>::~UnstructuredVolume()
64 {
65 if (this->ispcEquivalent) {
66 CALL_ISPC(VKLUnstructuredVolume_Destructor, this->ispcEquivalent);
67 }
68
69 if (rtcBVH)
70 rtcReleaseBVH(rtcBVH);
71 if (rtcDevice)
72 rtcReleaseDevice(rtcDevice);
73 }
74
75 template <int W>
commit()76 void UnstructuredVolume<W>::commit()
77 {
78 Volume<W>::commit();
79
80 // hex method planar/nonplanar
81
82 vertexPosition = this->template getParamDataT<vec3f>("vertex.position");
83 vertexValue = this->template getParamDataT<float>("vertex.data", nullptr);
84 indexPrefixed = this->template getParam<bool>("indexPrefixed", false);
85 cellValue = this->template getParamDataT<float>("cell.data", nullptr);
86 cellType = this->template getParamDataT<uint8_t>("cell.type", nullptr);
87
88 background = this->template getParamDataT<float>(
89 "background", 1, VKL_BACKGROUND_UNDEFINED);
90
91 if (!vertexValue && !cellValue) {
92 throw std::runtime_error(
93 "unstructured volume must have 'vertex.data' or 'cell.data'");
94 }
95 if ((!indexPrefixed && !cellType) || (indexPrefixed && cellType)) {
96 throw std::runtime_error(
97 "unstructured volume must have one of 'cell.type' or "
98 "'indexPrefixed'");
99 }
100
101 // index data can be 32 or 64 bit
102 Ref<const Data> index = this->template getParam<Data *>("index");
103
104 switch (index->dataType) {
105 case VKL_UINT:
106 index32Bit = true;
107 index32 = this->template getParamDataT<uint32_t>("index");
108 break;
109 case VKL_ULONG:
110 index32Bit = false;
111 index64 = this->template getParamDataT<uint64_t>("index");
112 break;
113 default:
114 throw std::runtime_error("unstructured volume unsupported index type");
115 }
116
117 // cell.index data can be 32 or 64 bit
118 Ref<const Data> cellIndex = this->template getParam<Data *>("cell.index");
119
120 switch (cellIndex->dataType) {
121 case VKL_UINT:
122 cell32Bit = true;
123 cellIndex32 = this->template getParamDataT<uint32_t>("cell.index");
124 break;
125 case VKL_ULONG:
126 cell32Bit = false;
127 cellIndex64 = this->template getParamDataT<uint64_t>("cell.index");
128 break;
129 default:
130 throw std::runtime_error(
131 "unstructured volume unsupported cell.index type");
132 }
133 nCells = cellIndex->size();
134
135 if (cellType) {
136 if (nCells != cellType->size())
137 throw std::runtime_error(
138 "unstructured volume #cells does not match #cell.type");
139 } else {
140 generatedCellType.resize(nCells);
141
142 for (int i = 0; i < nCells; i++) {
143 auto index = cell32Bit ? (*cellIndex32)[i] : (*cellIndex64)[i];
144 switch (getVertexId(index)) {
145 case 4:
146 generatedCellType[i] = VKL_TETRAHEDRON;
147 break;
148 case 8:
149 generatedCellType[i] = VKL_HEXAHEDRON;
150 break;
151 case 6:
152 generatedCellType[i] = VKL_WEDGE;
153 break;
154 case 5:
155 generatedCellType[i] = VKL_PYRAMID;
156 break;
157 default:
158 throw std::runtime_error(
159 "unstructured volume unsupported cell vertex count");
160 break;
161 }
162 }
163
164 Data *d = new Data(nCells,
165 VKL_UCHAR,
166 generatedCellType.data(),
167 VKL_DATA_SHARED_BUFFER,
168 0);
169 cellType = &(d->as<uint8_t>());
170 d->refDec();
171 }
172
173 hexIterative = this->template getParam<bool>("hexIterative", false);
174
175 bool needTolerances = false;
176 for (int i = 0; i < nCells; i++) {
177 auto cell = (*cellType)[i];
178 if (cell == VKL_WEDGE || cell == VKL_PYRAMID ||
179 (cell == VKL_HEXAHEDRON && hexIterative)) {
180 needTolerances = true;
181 break;
182 }
183 }
184
185 if (needTolerances)
186 calculateIterativeTolerance();
187
188 auto precompute =
189 this->template getParam<bool>("precomputedNormals", false);
190 if (precompute) {
191 if (faceNormals.empty()) {
192 calculateFaceNormals();
193 }
194 } else {
195 if (!faceNormals.empty()) {
196 faceNormals.clear();
197 faceNormals.shrink_to_fit();
198 }
199 }
200
201 buildBvhAndCalculateBounds();
202
203 computeOverlappingNodeMetadata(rtcRoot);
204
205 if (!this->ispcEquivalent) {
206 this->ispcEquivalent = CALL_ISPC(VKLUnstructuredVolume_Constructor);
207 }
208
209 CALL_ISPC(Volume_setBackground, this->ispcEquivalent, background->data());
210
211 CALL_ISPC(
212 VKLUnstructuredVolume_set,
213 this->ispcEquivalent,
214 (const ispc::box3f &)bounds,
215 ispc(vertexPosition),
216 index32Bit ? ispc(index32) : ispc(index64),
217 index32Bit,
218 ispc(vertexValue),
219 ispc(cellValue),
220 cell32Bit ? ispc(cellIndex32) : ispc(cellIndex64),
221 cell32Bit,
222 indexPrefixed,
223 ispc(cellType),
224 (void *)(rtcRoot),
225 faceNormals.empty() ? nullptr
226 : (const ispc::vec3f *)faceNormals.data(),
227 iterativeTolerance.empty() ? nullptr : iterativeTolerance.data(),
228 hexIterative);
229 }
230
231 template <int W>
newSampler()232 Sampler<W> *UnstructuredVolume<W>::newSampler()
233 {
234 return new UnstructuredSampler<W>(this);
235 }
236
237 template <int W>
getCellBBox(size_t id)238 box4f UnstructuredVolume<W>::getCellBBox(size_t id)
239 {
240 // get cell offset in the vertex indices array
241 uint64_t cOffset = getCellOffset(id);
242
243 // iterate through cell vertices
244 box4f bBox;
245 uint32_t maxIdx = getVerticesCount((*cellType)[id]);
246 for (uint32_t i = 0; i < maxIdx; i++) {
247 // get vertex index
248 uint64_t vId = getVertexId(cOffset + i);
249
250 // build 4 dimensional vertex with its position and value
251 const vec3f &v = (*vertexPosition)[vId];
252 float val = cellValue ? (*cellValue)[id] : (*vertexValue)[vId];
253 vec4f p = vec4f(v.x, v.y, v.z, val);
254
255 // extend bounding box
256 if (i == 0)
257 bBox.upper = bBox.lower = p;
258 else
259 bBox.extend(p);
260 }
261
262 return bBox;
263 }
264
errorFunction(void * userPtr,enum RTCError error,const char * str)265 static inline void errorFunction(void *userPtr,
266 enum RTCError error,
267 const char *str)
268 {
269 Device *device = reinterpret_cast<Device *>(userPtr);
270 LogMessageStream(device, VKL_LOG_WARNING)
271 << "error " << error << ": " << str << std::endl;
272 }
273
274 template <int W>
buildBvhAndCalculateBounds()275 void UnstructuredVolume<W>::buildBvhAndCalculateBounds()
276 {
277 rtcDevice = rtcNewDevice(NULL);
278 if (!rtcDevice) {
279 throw std::runtime_error("cannot create device");
280 }
281 rtcSetDeviceErrorFunction(rtcDevice, errorFunction, this->device.ptr);
282
283 containers::AlignedVector<RTCBuildPrimitive> prims;
284 containers::AlignedVector<range1f> range;
285 prims.resize(nCells);
286 range.resize(nCells);
287
288 tasking::parallel_for(nCells, [&](uint64_t taskIndex) {
289 box4f bound = getCellBBox(taskIndex);
290 prims[taskIndex].lower_x = bound.lower.x;
291 prims[taskIndex].lower_y = bound.lower.y;
292 prims[taskIndex].lower_z = bound.lower.z;
293 prims[taskIndex].geomID = taskIndex >> 32;
294 prims[taskIndex].upper_x = bound.upper.x;
295 prims[taskIndex].upper_y = bound.upper.y;
296 prims[taskIndex].upper_z = bound.upper.z;
297 prims[taskIndex].primID = taskIndex & 0xffffffff;
298 range[taskIndex] = range1f(bound.lower.w, bound.upper.w);
299 });
300
301 rtcBVH = rtcNewBVH(rtcDevice);
302 if (!rtcBVH) {
303 throw std::runtime_error("bvh creation failure");
304 }
305
306 RTCBuildArguments arguments = rtcDefaultBuildArguments();
307 arguments.byteSize = sizeof(arguments);
308 arguments.buildFlags = RTC_BUILD_FLAG_NONE;
309 arguments.buildQuality = RTC_BUILD_QUALITY_LOW;
310 arguments.maxBranchingFactor = 2;
311 arguments.maxDepth = 1024;
312 arguments.sahBlockSize = 1;
313 arguments.minLeafSize = 1;
314 arguments.maxLeafSize = 1;
315 arguments.traversalCost = 1.0f;
316 arguments.intersectionCost = 10.0f;
317 arguments.bvh = rtcBVH;
318 arguments.primitives = prims.data();
319 arguments.primitiveCount = prims.size();
320 arguments.primitiveArrayCapacity = prims.size();
321 arguments.createNode = InnerNode::create;
322 arguments.setNodeChildren = InnerNode::setChildren;
323 arguments.setNodeBounds = InnerNode::setBounds;
324 arguments.createLeaf = LeafNode::create;
325 arguments.splitPrimitive = nullptr;
326 arguments.buildProgress = nullptr;
327 arguments.userPtr = range.data();
328
329 rtcRoot = (Node *)rtcBuildBVH(&arguments);
330 if (!rtcRoot) {
331 throw std::runtime_error("bvh build failure");
332 }
333
334 if (rtcRoot->nominalLength.x < 0) {
335 auto &val = ((LeafNode *)rtcRoot)->bounds;
336 bounds = box3f(val.lower, val.upper);
337 } else {
338 auto &vals = ((InnerNode *)rtcRoot)->bounds;
339 bounds = box3f(vals[0].lower, vals[0].upper);
340 bounds.extend(box3f(vals[1].lower, vals[1].upper));
341 }
342 valueRange = rtcRoot->valueRange;
343
344 addLevelToNodes(rtcRoot, 0);
345
346 bvhDepth = getMaxNodeLevel(rtcRoot);
347 }
348
349 template <int W>
calculateIterativeTolerance()350 void UnstructuredVolume<W>::calculateIterativeTolerance()
351 {
352 iterativeTolerance.resize(nCells);
353 const uint32_t wedgeEdges[9][2] = {{0, 1},
354 {1, 2},
355 {2, 0},
356 {3, 4},
357 {4, 5},
358 {5, 3},
359 {0, 3},
360 {1, 4},
361 {2, 5}};
362 const uint32_t pyramidEdges[8][2] = {
363 {0, 1}, {1, 2}, {2, 3}, {3, 0}, {0, 4}, {1, 4}, {2, 4}, {3, 4}};
364 const uint32_t hexDiagonals[4][2] = {{0, 6}, {1, 7}, {2, 4}, {3, 5}};
365
366 // Build all tolerances
367 tasking::parallel_for(nCells, [&](uint64_t taskIndex) {
368 switch ((*cellType)[taskIndex]) {
369 case VKL_HEXAHEDRON:
370 if (!hexIterative)
371 calculateTolerance(taskIndex, hexDiagonals, 4);
372 break;
373 case VKL_WEDGE:
374 calculateTolerance(taskIndex, wedgeEdges, 9);
375 break;
376 case VKL_PYRAMID:
377 calculateTolerance(taskIndex, pyramidEdges, 8);
378 break;
379 default:
380 break;
381 }
382 });
383 }
384
385 template <int W>
calculateTolerance(const uint64_t cellId,const uint32_t edge[][2],const uint32_t count)386 void UnstructuredVolume<W>::calculateTolerance(const uint64_t cellId,
387 const uint32_t edge[][2],
388 const uint32_t count)
389 {
390 uint64_t cOffset = getCellOffset(cellId);
391
392 float longest = 0;
393 for (int i = 0; i < count; i++) {
394 uint64_t vId0 = getVertexId(cOffset + edge[i][0]);
395 uint64_t vId1 = getVertexId(cOffset + edge[i][1]);
396 const vec3f &v0 = (*vertexPosition)[vId0];
397 const vec3f &v1 = (*vertexPosition)[vId1];
398 const float dist = length(v0 - v1);
399 longest = std::max(longest, dist);
400 }
401
402 const float volumeBound = longest * longest * longest;
403 const float determinantTolerance =
404 1e-20 < .00001 * volumeBound ? 1e-20 : .00001 * volumeBound;
405
406 iterativeTolerance[cellId] = determinantTolerance;
407 }
408
409 template <int W>
calculateFaceNormals()410 void UnstructuredVolume<W>::calculateFaceNormals()
411 {
412 // Allocate memory for normal vectors
413 uint64_t numNormals = nCells * 6;
414 faceNormals.resize(numNormals);
415
416 // Define vertices order for normal calculation
417 const uint32_t tetrahedronFaces[4][3] = {
418 {2, 0, 1}, {3, 1, 0}, {3, 2, 1}, {2, 3, 0}};
419 const uint32_t hexahedronFaces[6][3] = {
420 {3, 0, 1}, {5, 1, 0}, {6, 2, 1}, {7, 3, 2}, {7, 4, 0}, {6, 5, 4}};
421 const uint32_t wedgeFaces[5][3] = {
422 {2, 0, 1}, {4, 1, 0}, {5, 2, 1}, {5, 3, 0}, {5, 4, 3}};
423 const uint32_t pyramidFaces[5][3] = {
424 {3, 0, 1}, {4, 1, 0}, {4, 2, 1}, {4, 3, 2}, {3, 4, 0}};
425
426 // Build all normals
427 tasking::parallel_for(nCells, [&](uint64_t taskIndex) {
428 switch ((*cellType)[taskIndex]) {
429 case VKL_TETRAHEDRON:
430 calculateCellNormals(taskIndex, tetrahedronFaces, 4);
431 break;
432 case VKL_HEXAHEDRON:
433 calculateCellNormals(taskIndex, hexahedronFaces, 6);
434 break;
435 case VKL_WEDGE:
436 calculateCellNormals(taskIndex, wedgeFaces, 5);
437 break;
438 case VKL_PYRAMID:
439 calculateCellNormals(taskIndex, pyramidFaces, 5);
440 break;
441 }
442 });
443 }
444
445 // Calculate all normals for arbitrary polyhedron
446 // based on given vertices order
447 template <int W>
calculateCellNormals(const uint64_t cellId,const uint32_t faces[6][3],const uint32_t facesCount)448 void UnstructuredVolume<W>::calculateCellNormals(const uint64_t cellId,
449 const uint32_t faces[6][3],
450 const uint32_t facesCount)
451 {
452 // Get cell offset in the vertex indices array
453 uint64_t cOffset = getCellOffset(cellId);
454
455 // Iterate through all faces
456 for (uint32_t i = 0; i < facesCount; i++) {
457 // Retrieve vertex positions
458 uint64_t vId0 = getVertexId(cOffset + faces[i][0]);
459 uint64_t vId1 = getVertexId(cOffset + faces[i][1]);
460 uint64_t vId2 = getVertexId(cOffset + faces[i][2]);
461 const vec3f &v0 = (*vertexPosition)[vId0];
462 const vec3f &v1 = (*vertexPosition)[vId1];
463 const vec3f &v2 = (*vertexPosition)[vId2];
464
465 // Calculate normal
466 faceNormals[cellId * 6 + i] = normalize(cross(v0 - v1, v2 - v1));
467 }
468 }
469
470 VKL_REGISTER_VOLUME(UnstructuredVolume<VKL_TARGET_WIDTH>,
471 CONCAT1(internal_unstructured_, VKL_TARGET_WIDTH))
472
473 } // namespace cpu_device
474 } // namespace openvkl
475