1 // Copyright 2009-2021 Intel Corporation
2 // SPDX-License-Identifier: Apache-2.0
3 
4 #include "scene_subdiv_mesh.h"
5 #include "scene.h"
6 #include "../subdiv/patch_eval.h"
7 #include "../subdiv/patch_eval_simd.h"
8 
9 #include "../../common/algorithms/parallel_sort.h"
10 #include "../../common/algorithms/parallel_prefix_sum.h"
11 #include "../../common/algorithms/parallel_for.h"
12 
13 /*! maximum number of user vertex buffers for subdivision surfaces */
14 #define RTC_MAX_USER_VERTEX_BUFFERS 65536
15 
16 namespace embree
17 {
18 #if defined(EMBREE_LOWEST_ISA)
19 
SubdivMesh(Device * device)20   SubdivMesh::SubdivMesh (Device* device)
21     : Geometry(device,GTY_SUBDIV_MESH,0,1),
22       displFunc(nullptr),
23       tessellationRate(2.0f),
24       numHalfEdges(0),
25       faceStartEdge(device,0),
26       halfEdgeFace(device,0),
27       invalid_face(device,0),
28       commitCounter(0)
29   {
30 
31     vertices.resize(numTimeSteps);
32     vertex_buffer_tags.resize(numTimeSteps);
33     topology.resize(1);
34     topology[0] = Topology(this);
35   }
36 
addElementsToCount(GeometryCounts & counts) const37   void SubdivMesh::addElementsToCount (GeometryCounts & counts) const
38   {
39     if (numTimeSteps == 1) counts.numSubdivPatches += numPrimitives;
40     else                   counts.numMBSubdivPatches += numPrimitives;
41   }
42 
setMask(unsigned mask)43   void SubdivMesh::setMask (unsigned mask)
44   {
45     this->mask = mask;
46     Geometry::update();
47   }
48 
setSubdivisionMode(unsigned topologyID,RTCSubdivisionMode mode)49   void SubdivMesh::setSubdivisionMode (unsigned topologyID, RTCSubdivisionMode mode)
50   {
51     if (topologyID >= topology.size())
52       throw_RTCError(RTC_ERROR_INVALID_OPERATION,"invalid topology ID");
53     topology[topologyID].setSubdivisionMode(mode);
54     Geometry::update();
55   }
56 
setVertexAttributeTopology(unsigned int vertexAttribID,unsigned int topologyID)57   void SubdivMesh::setVertexAttributeTopology(unsigned int vertexAttribID, unsigned int topologyID)
58   {
59     if (vertexAttribID < vertexAttribs.size()){
60       if (topologyID < topology.size()) {
61         if ((unsigned)vertexAttribs[vertexAttribID].userData != topologyID) {
62           vertexAttribs[vertexAttribID].userData = topologyID;
63           commitCounter++; // triggers recalculation of cached interpolation data
64         }
65       } else {
66         throw_RTCError(RTC_ERROR_INVALID_OPERATION, "invalid topology specified");
67       }
68     } else {
69       throw_RTCError(RTC_ERROR_INVALID_OPERATION, "invalid vertex attribute specified");
70     }
71   }
72 
setNumTimeSteps(unsigned int numTimeSteps)73   void SubdivMesh::setNumTimeSteps (unsigned int numTimeSteps)
74   {
75     vertices.resize(numTimeSteps);
76     vertex_buffer_tags.resize(numTimeSteps);
77     Geometry::setNumTimeSteps(numTimeSteps);
78   }
79 
setVertexAttributeCount(unsigned int N)80   void SubdivMesh::setVertexAttributeCount (unsigned int N)
81   {
82     vertexAttribs.resize(N);
83     vertex_attrib_buffer_tags.resize(N);
84     Geometry::update();
85   }
86 
setTopologyCount(unsigned int N)87   void SubdivMesh::setTopologyCount (unsigned int N)
88   {
89     if (N == 0)
90       throw_RTCError(RTC_ERROR_INVALID_ARGUMENT,"at least one topology has to exist")
91 
92     size_t begin = topology.size();
93     topology.resize(N);
94     for (size_t i = begin; i < topology.size(); i++)
95       topology[i] = Topology(this);
96   }
97 
setBuffer(RTCBufferType type,unsigned int slot,RTCFormat format,const Ref<Buffer> & buffer,size_t offset,size_t stride,unsigned int num)98   void SubdivMesh::setBuffer(RTCBufferType type, unsigned int slot, RTCFormat format, const Ref<Buffer>& buffer, size_t offset, size_t stride, unsigned int num)
99   {
100     /* verify that all accesses are 4 bytes aligned */
101     if (((size_t(buffer->getPtr()) + offset) & 0x3) || (stride & 0x3))
102       throw_RTCError(RTC_ERROR_INVALID_OPERATION, "data must be 4 bytes aligned");
103 
104     if (type != RTC_BUFFER_TYPE_LEVEL)
105       commitCounter++;
106 
107     if (type == RTC_BUFFER_TYPE_VERTEX)
108     {
109       if (format != RTC_FORMAT_FLOAT3)
110         throw_RTCError(RTC_ERROR_INVALID_OPERATION, "invalid vertex buffer format");
111 
112       if (slot >= vertices.size())
113         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid vertex buffer slot");
114 
115       vertices[slot].set(buffer, offset, stride, num, format);
116       vertices[slot].checkPadding16();
117     }
118     else if (type == RTC_BUFFER_TYPE_VERTEX_ATTRIBUTE)
119     {
120       if (format < RTC_FORMAT_FLOAT || format > RTC_FORMAT_FLOAT16)
121         throw_RTCError(RTC_ERROR_INVALID_OPERATION, "invalid vertex attribute buffer format");
122 
123       if (slot >= vertexAttribs.size())
124         throw_RTCError(RTC_ERROR_INVALID_OPERATION, "invalid vertex attribute buffer slot");
125 
126       vertexAttribs[slot].set(buffer, offset, stride, num, format);
127       vertexAttribs[slot].checkPadding16();
128     }
129     else if (type == RTC_BUFFER_TYPE_FACE)
130     {
131       if (slot != 0)
132         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
133       if (format != RTC_FORMAT_UINT)
134         throw_RTCError(RTC_ERROR_INVALID_OPERATION, "invalid face buffer format");
135 
136       faceVertices.set(buffer, offset, stride, num, format);
137       setNumPrimitives(num);
138     }
139     else if (type == RTC_BUFFER_TYPE_INDEX)
140     {
141       if (format != RTC_FORMAT_UINT)
142         throw_RTCError(RTC_ERROR_INVALID_OPERATION, "invalid face buffer format");
143 
144       if (slot >= topology.size())
145         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid index buffer slot");
146 
147       topology[slot].vertexIndices.set(buffer, offset, stride, num, format);
148     }
149     else if (type == RTC_BUFFER_TYPE_EDGE_CREASE_INDEX)
150     {
151       if (slot != 0)
152         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
153       if (format != RTC_FORMAT_UINT2)
154         throw_RTCError(RTC_ERROR_INVALID_OPERATION, "invalid edge crease index buffer format");
155 
156       edge_creases.set(buffer, offset, stride, num, format);
157     }
158     else if (type == RTC_BUFFER_TYPE_EDGE_CREASE_WEIGHT)
159     {
160       if (slot != 0)
161         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
162       if (format != RTC_FORMAT_FLOAT)
163         throw_RTCError(RTC_ERROR_INVALID_OPERATION, "invalid edge crease weight buffer format");
164 
165       edge_crease_weights.set(buffer, offset, stride, num, format);
166     }
167     else if (type == RTC_BUFFER_TYPE_VERTEX_CREASE_INDEX)
168     {
169       if (slot != 0)
170         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
171       if (format != RTC_FORMAT_UINT)
172         throw_RTCError(RTC_ERROR_INVALID_OPERATION, "invalid vertex crease index buffer format");
173 
174       vertex_creases.set(buffer, offset, stride, num, format);
175     }
176     else if (type == RTC_BUFFER_TYPE_VERTEX_CREASE_WEIGHT)
177     {
178       if (slot != 0)
179         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
180       if (format != RTC_FORMAT_FLOAT)
181         throw_RTCError(RTC_ERROR_INVALID_OPERATION, "invalid vertex crease weight buffer format");
182 
183       vertex_crease_weights.set(buffer, offset, stride, num, format);
184     }
185     else if (type == RTC_BUFFER_TYPE_HOLE)
186     {
187       if (slot != 0)
188         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
189       if (format != RTC_FORMAT_UINT)
190         throw_RTCError(RTC_ERROR_INVALID_OPERATION, "invalid hole buffer format");
191 
192       holes.set(buffer, offset, stride, num, format);
193     }
194     else if (type == RTC_BUFFER_TYPE_LEVEL)
195     {
196       if (slot != 0)
197         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
198       if (format != RTC_FORMAT_FLOAT)
199         throw_RTCError(RTC_ERROR_INVALID_OPERATION, "invalid level buffer format");
200 
201       levels.set(buffer, offset, stride, num, format);
202     }
203     else
204     {
205       throw_RTCError(RTC_ERROR_INVALID_ARGUMENT,"unknown buffer type");
206     }
207   }
208 
getBuffer(RTCBufferType type,unsigned int slot)209   void* SubdivMesh::getBuffer(RTCBufferType type, unsigned int slot)
210   {
211     if (type == RTC_BUFFER_TYPE_VERTEX)
212     {
213       if (slot >= vertices.size())
214         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
215       return vertices[slot].getPtr();
216     }
217     else if (type == RTC_BUFFER_TYPE_VERTEX_ATTRIBUTE)
218     {
219       if (slot >= vertexAttribs.size())
220         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
221       return vertexAttribs[slot].getPtr();
222     }
223     else if (type == RTC_BUFFER_TYPE_FACE)
224     {
225       if (slot != 0)
226         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
227       return faceVertices.getPtr();
228     }
229     else if (type == RTC_BUFFER_TYPE_INDEX)
230     {
231       if (slot >= topology.size())
232         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
233       return topology[slot].vertexIndices.getPtr();
234     }
235     else if (type == RTC_BUFFER_TYPE_EDGE_CREASE_INDEX)
236     {
237       if (slot != 0)
238         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
239       return edge_creases.getPtr();
240     }
241     else if (type == RTC_BUFFER_TYPE_EDGE_CREASE_WEIGHT)
242     {
243       if (slot != 0)
244         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
245       return edge_crease_weights.getPtr();
246     }
247     else if (type == RTC_BUFFER_TYPE_VERTEX_CREASE_INDEX)
248     {
249       if (slot != 0)
250         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
251       return vertex_creases.getPtr();
252     }
253     else if (type == RTC_BUFFER_TYPE_VERTEX_CREASE_WEIGHT)
254     {
255       if (slot != 0)
256         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
257       return vertex_crease_weights.getPtr();
258     }
259     else if (type == RTC_BUFFER_TYPE_HOLE)
260     {
261       if (slot != 0)
262         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
263       return holes.getPtr();
264     }
265     else if (type == RTC_BUFFER_TYPE_LEVEL)
266     {
267       if (slot != 0)
268         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
269       return levels.getPtr();
270     }
271     else
272     {
273       throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "unknown buffer type");
274       return nullptr;
275     }
276   }
277 
updateBuffer(RTCBufferType type,unsigned int slot)278   void SubdivMesh::updateBuffer(RTCBufferType type, unsigned int slot)
279   {
280     if (type != RTC_BUFFER_TYPE_LEVEL)
281       commitCounter++;
282 
283     if (type == RTC_BUFFER_TYPE_VERTEX)
284     {
285       if (slot >= vertices.size())
286         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
287       vertices[slot].setModified();
288     }
289     else if (type == RTC_BUFFER_TYPE_VERTEX_ATTRIBUTE)
290     {
291       if (slot >= vertexAttribs.size())
292         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
293       vertexAttribs[slot].setModified();
294     }
295     else if (type == RTC_BUFFER_TYPE_FACE)
296     {
297       if (slot != 0)
298         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
299       faceVertices.setModified();
300     }
301     else if (type == RTC_BUFFER_TYPE_INDEX)
302     {
303       if (slot >= topology.size())
304         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
305       topology[slot].vertexIndices.setModified();
306     }
307     else if (type == RTC_BUFFER_TYPE_EDGE_CREASE_INDEX)
308     {
309       if (slot != 0)
310         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
311       edge_creases.setModified();
312     }
313     else if (type == RTC_BUFFER_TYPE_EDGE_CREASE_WEIGHT)
314     {
315       if (slot != 0)
316         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
317       edge_crease_weights.setModified();
318     }
319     else if (type == RTC_BUFFER_TYPE_VERTEX_CREASE_INDEX)
320     {
321       if (slot != 0)
322         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
323       vertex_creases.setModified();
324     }
325     else if (type == RTC_BUFFER_TYPE_VERTEX_CREASE_WEIGHT)
326     {
327       if (slot != 0)
328         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
329       vertex_crease_weights.setModified();
330     }
331     else if (type == RTC_BUFFER_TYPE_HOLE)
332     {
333       if (slot != 0)
334         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
335       holes.setModified();
336     }
337     else if (type == RTC_BUFFER_TYPE_LEVEL)
338     {
339       if (slot != 0)
340         throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid buffer slot");
341       levels.setModified();
342     }
343     else
344     {
345       throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "unknown buffer type");
346     }
347 
348     Geometry::update();
349   }
350 
setDisplacementFunction(RTCDisplacementFunctionN func)351   void SubdivMesh::setDisplacementFunction (RTCDisplacementFunctionN func)
352   {
353     this->displFunc = func;
354   }
355 
setTessellationRate(float N)356   void SubdivMesh::setTessellationRate(float N)
357   {
358     tessellationRate = N;
359     levels.setModified();
360   }
361 
pair64(unsigned int x,unsigned int y)362   __forceinline uint64_t pair64(unsigned int x, unsigned int y)
363   {
364     if (x<y) std::swap(x,y);
365     return (((uint64_t)x) << 32) | (uint64_t)y;
366   }
367 
Topology(SubdivMesh * mesh)368   SubdivMesh::Topology::Topology(SubdivMesh* mesh)
369     : mesh(mesh), subdiv_mode(RTC_SUBDIVISION_MODE_SMOOTH_BOUNDARY), halfEdges(mesh->device,0)
370   {
371   }
372 
setSubdivisionMode(RTCSubdivisionMode mode)373   void SubdivMesh::Topology::setSubdivisionMode (RTCSubdivisionMode mode)
374   {
375     if (subdiv_mode == mode) return;
376     subdiv_mode = mode;
377     mesh->updateBuffer(RTC_BUFFER_TYPE_VERTEX_CREASE_WEIGHT, 0);
378   }
379 
update()380   void SubdivMesh::Topology::update () {
381     vertexIndices.setModified();
382   }
383 
verify(size_t numVertices)384   bool SubdivMesh::Topology::verify (size_t numVertices)
385   {
386     size_t ofs = 0;
387     for (size_t i=0; i<mesh->size(); i++)
388     {
389       int valence = mesh->faceVertices[i];
390       for (size_t j=ofs; j<ofs+valence; j++)
391       {
392         if (j >= vertexIndices.size())
393           return false;
394 
395         if (vertexIndices[j] >= numVertices)
396           return false;
397       }
398       ofs += valence;
399     }
400     return true;
401   }
402 
calculateHalfEdges()403   void SubdivMesh::Topology::calculateHalfEdges()
404   {
405     const size_t blockSize = 4096;
406     const size_t numEdges = mesh->numEdges();
407     const size_t numFaces = mesh->numFaces();
408     const size_t numHalfEdges = mesh->numHalfEdges;
409 
410     /* allocate temporary array */
411     halfEdges0.resize(numEdges);
412     halfEdges1.resize(numEdges);
413 
414     /* create all half edges */
415     parallel_for( size_t(0), numFaces, blockSize, [&](const range<size_t>& r)
416     {
417       for (size_t f=r.begin(); f<r.end(); f++)
418       {
419 	const unsigned N = mesh->faceVertices[f];
420 	const unsigned e = mesh->faceStartEdge[f];
421 
422 	for (unsigned de=0; de<N; de++)
423 	{
424 	  HalfEdge* edge = &halfEdges[e+de];
425           int nextOfs = (de == (N-1)) ? -int(N-1) : +1;
426           int prevOfs = (de ==     0) ? +int(N-1) : -1;
427 
428 	  const unsigned int startVertex = vertexIndices[e+de];
429           const unsigned int endVertex = vertexIndices[e+de+nextOfs];
430 	  const uint64_t key = SubdivMesh::Edge(startVertex,endVertex);
431 
432           /* we always have to use the geometry topology to lookup creases */
433           const unsigned int startVertex0 = mesh->topology[0].vertexIndices[e+de];
434           const unsigned int endVertex0 = mesh->topology[0].vertexIndices[e+de+nextOfs];
435 	  const uint64_t key0 = SubdivMesh::Edge(startVertex0,endVertex0);
436 
437 	  edge->vtx_index              = startVertex;
438 	  edge->next_half_edge_ofs     = nextOfs;
439 	  edge->prev_half_edge_ofs     = prevOfs;
440 	  edge->opposite_half_edge_ofs = 0;
441 	  edge->edge_crease_weight     = mesh->edgeCreaseMap.lookup(key0,0.0f);
442 	  edge->vertex_crease_weight   = mesh->vertexCreaseMap.lookup(startVertex0,0.0f);
443 	  edge->edge_level             = mesh->getEdgeLevel(e+de);
444           edge->patch_type             = HalfEdge::COMPLEX_PATCH; // type gets updated below
445           edge->vertex_type            = HalfEdge::REGULAR_VERTEX;
446 
447           if (unlikely(mesh->holeSet.lookup(unsigned(f))))
448 	    halfEdges1[e+de] = SubdivMesh::KeyHalfEdge(std::numeric_limits<uint64_t>::max(),edge);
449 	  else
450 	    halfEdges1[e+de] = SubdivMesh::KeyHalfEdge(key,edge);
451 	}
452       }
453     });
454 
455     /* sort half edges to find adjacent edges */
456     radix_sort_u64(halfEdges1.data(),halfEdges0.data(),numHalfEdges);
457 
458     /* link all adjacent pairs of edges */
459     parallel_for( size_t(0), numHalfEdges, blockSize, [&](const range<size_t>& r)
460     {
461       /* skip if start of adjacent edges was not in our range */
462       size_t e=r.begin();
463       if (e != 0 && (halfEdges1[e].key == halfEdges1[e-1].key)) {
464 	const uint64_t key = halfEdges1[e].key;
465 	while (e<r.end() && halfEdges1[e].key == key) e++;
466       }
467 
468       /* process all adjacent edges starting in our range */
469       while (e<r.end())
470       {
471 	const uint64_t key = halfEdges1[e].key;
472 	if (key == std::numeric_limits<uint64_t>::max()) break;
473 	size_t N=1; while (e+N<numHalfEdges && halfEdges1[e+N].key == key) N++;
474 
475         /* border edges are identified by not having an opposite edge set */
476 	if (N == 1) {
477           halfEdges1[e].edge->edge_crease_weight = float(inf);
478 	}
479 
480         /* standard edge shared between two faces */
481         else if (N == 2)
482         {
483           /* create edge crease if winding order mismatches between neighboring patches */
484           if (halfEdges1[e+0].edge->next()->vtx_index != halfEdges1[e+1].edge->vtx_index)
485           {
486             halfEdges1[e+0].edge->edge_crease_weight = float(inf);
487             halfEdges1[e+1].edge->edge_crease_weight = float(inf);
488           }
489           /* otherwise mark edges as opposites of each other */
490           else {
491             halfEdges1[e+0].edge->setOpposite(halfEdges1[e+1].edge);
492             halfEdges1[e+1].edge->setOpposite(halfEdges1[e+0].edge);
493           }
494 	}
495 
496         /* non-manifold geometry is handled by keeping vertices fixed during subdivision */
497         else {
498 	  for (size_t i=0; i<N; i++) {
499 	    halfEdges1[e+i].edge->vertex_crease_weight = inf;
500             halfEdges1[e+i].edge->vertex_type = HalfEdge::NON_MANIFOLD_EDGE_VERTEX;
501             halfEdges1[e+i].edge->edge_crease_weight = inf;
502 
503 	    halfEdges1[e+i].edge->next()->vertex_crease_weight = inf;
504             halfEdges1[e+i].edge->next()->vertex_type = HalfEdge::NON_MANIFOLD_EDGE_VERTEX;
505             halfEdges1[e+i].edge->next()->edge_crease_weight = inf;
506 	  }
507 	}
508 	e+=N;
509       }
510     });
511 
512     /* set subdivision mode and calculate patch types */
513     parallel_for( size_t(0), numFaces, blockSize, [&](const range<size_t>& r)
514     {
515       for (size_t f=r.begin(); f<r.end(); f++)
516       {
517         HalfEdge* edge = &halfEdges[mesh->faceStartEdge[f]];
518 
519         /* for vertex topology we also test if vertices are valid */
520         if (this == &mesh->topology[0])
521         {
522           /* calculate if face is valid */
523           for (size_t t=0; t<mesh->numTimeSteps; t++)
524             mesh->invalidFace(f,t) = !edge->valid(mesh->vertices[t]) || mesh->holeSet.lookup(unsigned(f));
525         }
526 
527         /* pin some edges and vertices */
528         for (size_t i=0; i<mesh->faceVertices[f]; i++)
529         {
530           /* pin corner vertices when requested by user */
531           if (subdiv_mode == RTC_SUBDIVISION_MODE_PIN_CORNERS && edge[i].isCorner())
532             edge[i].vertex_crease_weight = float(inf);
533 
534           /* pin all border vertices when requested by user */
535           else if (subdiv_mode == RTC_SUBDIVISION_MODE_PIN_BOUNDARY && edge[i].vertexHasBorder())
536             edge[i].vertex_crease_weight = float(inf);
537 
538           /* pin all edges and vertices when requested by user */
539           else if (subdiv_mode == RTC_SUBDIVISION_MODE_PIN_ALL) {
540             edge[i].edge_crease_weight = float(inf);
541             edge[i].vertex_crease_weight = float(inf);
542           }
543         }
544 
545         /* we have to calculate patch_type last! */
546         HalfEdge::PatchType patch_type = edge->patchType();
547         for (size_t i=0; i<mesh->faceVertices[f]; i++)
548           edge[i].patch_type = patch_type;
549       }
550     });
551   }
552 
updateHalfEdges()553   void SubdivMesh::Topology::updateHalfEdges()
554   {
555     /* we always use the geometry topology to lookup creases */
556     mvector<HalfEdge>& halfEdgesGeom = mesh->topology[0].halfEdges;
557 
558     /* assume we do no longer recalculate in the future and clear these arrays */
559     halfEdges0.clear();
560     halfEdges1.clear();
561 
562     /* calculate which data to update */
563     const bool updateEdgeCreases   = mesh->topology[0].vertexIndices.isLocalModified() || mesh->edge_creases.isLocalModified()   || mesh->edge_crease_weights.isLocalModified();
564     const bool updateVertexCreases = mesh->topology[0].vertexIndices.isLocalModified() || mesh->vertex_creases.isLocalModified() || mesh->vertex_crease_weights.isLocalModified();
565     const bool updateLevels = mesh->levels.isLocalModified();
566 
567     /* parallel loop over all half edges */
568     parallel_for( size_t(0), mesh->numHalfEdges, size_t(4096), [&](const range<size_t>& r)
569     {
570       for (size_t i=r.begin(); i!=r.end(); i++)
571       {
572 	HalfEdge& edge = halfEdges[i];
573 
574 	if (updateLevels)
575 	  edge.edge_level = mesh->getEdgeLevel(i);
576 
577 	if (updateEdgeCreases) {
578 	  if (edge.hasOpposite()) // leave weight at inf for borders
579             edge.edge_crease_weight = mesh->edgeCreaseMap.lookup((uint64_t)halfEdgesGeom[i].getEdge(),0.0f);
580 	}
581 
582         /* we only use user specified vertex_crease_weight if the vertex is manifold */
583         if (updateVertexCreases && edge.vertex_type != HalfEdge::NON_MANIFOLD_EDGE_VERTEX)
584         {
585 	  edge.vertex_crease_weight = mesh->vertexCreaseMap.lookup(halfEdgesGeom[i].vtx_index,0.0f);
586 
587           /* pin corner vertices when requested by user */
588           if (subdiv_mode == RTC_SUBDIVISION_MODE_PIN_CORNERS && edge.isCorner())
589             edge.vertex_crease_weight = float(inf);
590 
591           /* pin all border vertices when requested by user */
592           else if (subdiv_mode == RTC_SUBDIVISION_MODE_PIN_BOUNDARY && edge.vertexHasBorder())
593             edge.vertex_crease_weight = float(inf);
594 
595           /* pin every vertex when requested by user */
596           else if (subdiv_mode == RTC_SUBDIVISION_MODE_PIN_ALL) {
597             edge.edge_crease_weight = float(inf);
598             edge.vertex_crease_weight = float(inf);
599           }
600         }
601 
602         /* update patch type */
603         if (updateEdgeCreases || updateVertexCreases) {
604           edge.patch_type = edge.patchType();
605         }
606       }
607     });
608   }
609 
initializeHalfEdgeStructures()610   void SubdivMesh::Topology::initializeHalfEdgeStructures ()
611   {
612     /* if vertex indices not set we ignore this topology */
613     if (!vertexIndices)
614       return;
615 
616     /* allocate half edge array */
617     halfEdges.resize(mesh->numEdges());
618 
619     /* check if we have to recalculate the half edges */
620     bool recalculate = false;
621     recalculate |= vertexIndices.isLocalModified();
622     recalculate |= mesh->faceVertices.isLocalModified();
623     recalculate |= mesh->holes.isLocalModified();
624 
625     /* check if we can simply update the half edges */
626     bool update = false;
627     update |= mesh->topology[0].vertexIndices.isLocalModified(); // we use this buffer to copy creases to interpolation topologies
628     update |= mesh->edge_creases.isLocalModified();
629     update |= mesh->edge_crease_weights.isLocalModified();
630     update |= mesh->vertex_creases.isLocalModified();
631     update |= mesh->vertex_crease_weights.isLocalModified();
632     update |= mesh->levels.isLocalModified();
633 
634     /* now either recalculate or update the half edges */
635     if (recalculate) calculateHalfEdges();
636     else if (update) updateHalfEdges();
637 
638     /* cleanup some state for static scenes */
639     /* if (mesh->scene_ == nullptr || mesh->scene_->isStaticAccel())
640     {
641       halfEdges0.clear();
642       halfEdges1.clear();
643     } */
644 
645     /* clear modified state of all buffers */
646     vertexIndices.clearLocalModified();
647   }
648 
printStatistics()649   void SubdivMesh::printStatistics()
650   {
651     size_t numBilinearFaces = 0;
652     size_t numRegularQuadFaces = 0;
653     size_t numIrregularQuadFaces = 0;
654     size_t numComplexFaces = 0;
655 
656     for (size_t e=0, f=0; f<numFaces(); e+=faceVertices[f++])
657     {
658       switch (topology[0].halfEdges[e].patch_type) {
659       case HalfEdge::BILINEAR_PATCH      : numBilinearFaces++;   break;
660       case HalfEdge::REGULAR_QUAD_PATCH  : numRegularQuadFaces++;   break;
661       case HalfEdge::IRREGULAR_QUAD_PATCH: numIrregularQuadFaces++; break;
662       case HalfEdge::COMPLEX_PATCH       : numComplexFaces++;   break;
663       }
664     }
665 
666     std::cout << "numFaces = " << numFaces() << ", "
667               << "numBilinearFaces = " << numBilinearFaces << " (" << 100.0f * numBilinearFaces / numFaces() << "%), "
668               << "numRegularQuadFaces = " << numRegularQuadFaces << " (" << 100.0f * numRegularQuadFaces / numFaces() << "%), "
669               << "numIrregularQuadFaces " << numIrregularQuadFaces << " (" << 100.0f * numIrregularQuadFaces / numFaces() << "%) "
670               << "numComplexFaces " << numComplexFaces << " (" << 100.0f * numComplexFaces / numFaces() << "%) "
671               << std::endl;
672   }
673 
initializeHalfEdgeStructures()674   void SubdivMesh::initializeHalfEdgeStructures ()
675   {
676     double t0 = getSeconds();
677 
678     invalid_face.resize(numFaces()*numTimeSteps);
679 
680     /* calculate start edge of each face */
681     faceStartEdge.resize(numFaces());
682 
683     if (faceVertices.isLocalModified())
684     {
685       numHalfEdges = parallel_prefix_sum(faceVertices,faceStartEdge,numFaces(),0,std::plus<unsigned>());
686 
687       /* calculate face of each half edge */
688       halfEdgeFace.resize(numHalfEdges);
689       for (size_t f=0, h=0; f<numFaces(); f++)
690         for (size_t e=0; e<faceVertices[f]; e++)
691           halfEdgeFace[h++] = (unsigned int) f;
692     }
693 
694     /* create set with all vertex creases */
695     if (vertex_creases.isLocalModified() || vertex_crease_weights.isLocalModified())
696       vertexCreaseMap.init(vertex_creases,vertex_crease_weights);
697 
698     /* create map with all edge creases */
699     if (edge_creases.isLocalModified() || edge_crease_weights.isLocalModified())
700       edgeCreaseMap.init(edge_creases,edge_crease_weights);
701 
702     /* create set with all holes */
703     if (holes.isLocalModified())
704       holeSet.init(holes);
705 
706     /* create topology */
707     for (auto& t: topology)
708       t.initializeHalfEdgeStructures();
709 
710     /* create interpolation cache mapping for interpolatable meshes */
711     for (size_t i=0; i<vertex_buffer_tags.size(); i++)
712       vertex_buffer_tags[i].resize(numFaces()*numInterpolationSlots4(vertices[i].getStride()));
713     for (size_t i=0; i<vertexAttribs.size(); i++)
714       if (vertexAttribs[i]) vertex_attrib_buffer_tags[i].resize(numFaces()*numInterpolationSlots4(vertexAttribs[i].getStride()));
715 
716     /* cleanup some state for static scenes */
717     /* if (scene_ == nullptr || scene_->isStaticAccel())
718     {
719       vertexCreaseMap.clear();
720       edgeCreaseMap.clear();
721     } */
722 
723     /* clear modified state of all buffers */
724     faceVertices.clearLocalModified();
725     holes.clearLocalModified();
726     for (auto& buffer : vertices) buffer.clearLocalModified();
727     levels.clearLocalModified();
728     edge_creases.clearLocalModified();
729     edge_crease_weights.clearLocalModified();
730     vertex_creases.clearLocalModified();
731     vertex_crease_weights.clearLocalModified();
732 
733     double t1 = getSeconds();
734 
735     /* print statistics in verbose mode */
736     if (device->verbosity(2)) {
737       std::cout << "half edge generation = " << 1000.0*(t1-t0) << "ms, " << 1E-6*double(numHalfEdges)/(t1-t0) << "M/s" << std::endl;
738       printStatistics();
739     }
740   }
741 
verify()742   bool SubdivMesh::verify ()
743   {
744     /*! verify consistent size of vertex arrays */
745     if (vertices.size() == 0) return false;
746     for (const auto& buffer : vertices)
747       if (buffer.size() != numVertices())
748         return false;
749 
750     /*! verify vertex indices */
751     if (!topology[0].verify(numVertices()))
752       return false;
753 
754     for (auto& b : vertexAttribs)
755       if (!topology[b.userData].verify(b.size()))
756         return false;
757 
758     /*! verify vertices */
759     for (const auto& buffer : vertices)
760       for (size_t i=0; i<buffer.size(); i++)
761 	if (!isvalid(buffer[i]))
762 	  return false;
763 
764     return true;
765   }
766 
commit()767   void SubdivMesh::commit ()
768   {
769     initializeHalfEdgeStructures();
770     Geometry::commit();
771   }
772 
getFirstHalfEdge(unsigned int faceID)773   unsigned int SubdivMesh::getFirstHalfEdge(unsigned int faceID)
774   {
775     if (faceID >= numFaces())
776       throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid face");
777 
778     return faceStartEdge[faceID];
779   }
780 
getFace(unsigned int edgeID)781   unsigned int SubdivMesh::getFace(unsigned int edgeID)
782   {
783     if (edgeID >= numHalfEdges)
784       throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid edge");
785 
786     return halfEdgeFace[edgeID];
787   }
788 
getNextHalfEdge(unsigned int edgeID)789   unsigned int SubdivMesh::getNextHalfEdge(unsigned int edgeID)
790   {
791     if (edgeID >= numHalfEdges)
792       throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid half edge");
793 
794     return edgeID + topology[0].halfEdges[edgeID].next_half_edge_ofs;
795   }
796 
getPreviousHalfEdge(unsigned int edgeID)797   unsigned int SubdivMesh::getPreviousHalfEdge(unsigned int edgeID)
798   {
799      if (edgeID >= numHalfEdges)
800       throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid half edge");
801 
802     return edgeID + topology[0].halfEdges[edgeID].prev_half_edge_ofs;
803   }
804 
getOppositeHalfEdge(unsigned int topologyID,unsigned int edgeID)805   unsigned int SubdivMesh::getOppositeHalfEdge(unsigned int topologyID, unsigned int edgeID)
806   {
807     if (topologyID >= topology.size())
808       throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid topology");
809 
810     if (edgeID >= numHalfEdges)
811       throw_RTCError(RTC_ERROR_INVALID_ARGUMENT, "invalid half edge");
812 
813     return edgeID + topology[topologyID].halfEdges[edgeID].opposite_half_edge_ofs;
814   }
815 
816 #endif
817 
818   namespace isa
819   {
createSubdivMesh(Device * device)820     SubdivMesh* createSubdivMesh(Device* device) {
821       return new SubdivMeshISA(device);
822     }
823 
interpolate(const RTCInterpolateArguments * const args)824     void SubdivMeshISA::interpolate(const RTCInterpolateArguments* const args)
825     {
826       unsigned int primID = args->primID;
827       float u = args->u;
828       float v = args->v;
829       RTCBufferType bufferType = args->bufferType;
830       unsigned int bufferSlot = args->bufferSlot;
831       float* P = args->P;
832       float* dPdu = args->dPdu;
833       float* dPdv = args->dPdv;
834       float* ddPdudu = args->ddPdudu;
835       float* ddPdvdv = args->ddPdvdv;
836       float* ddPdudv = args->ddPdudv;
837       unsigned int valueCount = args->valueCount;
838 
839       /* calculate base pointer and stride */
840       assert((bufferType == RTC_BUFFER_TYPE_VERTEX && bufferSlot < RTC_MAX_TIME_STEP_COUNT) ||
841              (bufferType == RTC_BUFFER_TYPE_VERTEX_ATTRIBUTE && bufferSlot < RTC_MAX_USER_VERTEX_BUFFERS));
842       const char* src = nullptr;
843       size_t stride = 0;
844       std::vector<SharedLazyTessellationCache::CacheEntry>* baseEntry = nullptr;
845       Topology* topo = nullptr;
846       if (bufferType == RTC_BUFFER_TYPE_VERTEX_ATTRIBUTE) {
847         assert(bufferSlot < vertexAttribs.size());
848         src    = vertexAttribs[bufferSlot].getPtr();
849         stride = vertexAttribs[bufferSlot].getStride();
850         baseEntry = &vertex_attrib_buffer_tags[bufferSlot];
851         int topologyID = vertexAttribs[bufferSlot].userData;
852         topo = &topology[topologyID];
853       } else {
854         assert(bufferSlot < numTimeSteps);
855         src    = vertices[bufferSlot].getPtr();
856         stride = vertices[bufferSlot].getStride();
857         baseEntry = &vertex_buffer_tags[bufferSlot];
858         topo = &topology[0];
859       }
860 
861       bool has_P = P;
862       bool has_dP = dPdu;     assert(!has_dP  || dPdv);
863       bool has_ddP = ddPdudu; assert(!has_ddP || (ddPdvdv && ddPdudu));
864 
865       for (unsigned int i=0; i<valueCount; i+=4)
866       {
867         vfloat4 Pt, dPdut, dPdvt, ddPdudut, ddPdvdvt, ddPdudvt;
868         isa::PatchEval<vfloat4,vfloat4>(baseEntry->at(interpolationSlot(primID,i/4,stride)),commitCounter,
869                                         topo->getHalfEdge(primID),src+i*sizeof(float),stride,u,v,
870                                         has_P ? &Pt : nullptr,
871                                         has_dP ? &dPdut : nullptr,
872                                         has_dP ? &dPdvt : nullptr,
873                                         has_ddP ? &ddPdudut : nullptr,
874                                         has_ddP ? &ddPdvdvt : nullptr,
875                                         has_ddP ? &ddPdudvt : nullptr);
876 
877         if (has_P) {
878           for (size_t j=i; j<min(i+4,valueCount); j++)
879             P[j] = Pt[j-i];
880         }
881         if (has_dP)
882         {
883           for (size_t j=i; j<min(i+4,valueCount); j++) {
884             dPdu[j] = dPdut[j-i];
885             dPdv[j] = dPdvt[j-i];
886           }
887         }
888         if (has_ddP)
889         {
890           for (size_t j=i; j<min(i+4,valueCount); j++) {
891             ddPdudu[j] = ddPdudut[j-i];
892             ddPdvdv[j] = ddPdvdvt[j-i];
893             ddPdudv[j] = ddPdudvt[j-i];
894           }
895         }
896       }
897     }
898 
interpolateN(const RTCInterpolateNArguments * const args)899     void SubdivMeshISA::interpolateN(const RTCInterpolateNArguments* const args)
900     {
901       const void* valid_i = args->valid;
902       const unsigned* primIDs = args->primIDs;
903       const float* u = args->u;
904       const float* v = args->v;
905       unsigned int N = args->N;
906       RTCBufferType bufferType = args->bufferType;
907       unsigned int bufferSlot = args->bufferSlot;
908       float* P = args->P;
909       float* dPdu = args->dPdu;
910       float* dPdv = args->dPdv;
911       float* ddPdudu = args->ddPdudu;
912       float* ddPdvdv = args->ddPdvdv;
913       float* ddPdudv = args->ddPdudv;
914       unsigned int valueCount = args->valueCount;
915 
916       /* calculate base pointer and stride */
917       assert((bufferType == RTC_BUFFER_TYPE_VERTEX && bufferSlot < RTC_MAX_TIME_STEP_COUNT) ||
918              (bufferType == RTC_BUFFER_TYPE_VERTEX_ATTRIBUTE && bufferSlot < RTC_MAX_USER_VERTEX_BUFFERS));
919       const char* src = nullptr;
920       size_t stride = 0;
921       std::vector<SharedLazyTessellationCache::CacheEntry>* baseEntry = nullptr;
922       Topology* topo = nullptr;
923       if (bufferType == RTC_BUFFER_TYPE_VERTEX_ATTRIBUTE) {
924         assert(bufferSlot < vertexAttribs.size());
925         src    = vertexAttribs[bufferSlot].getPtr();
926         stride = vertexAttribs[bufferSlot].getStride();
927         baseEntry = &vertex_attrib_buffer_tags[bufferSlot];
928         int topologyID = vertexAttribs[bufferSlot].userData;
929         topo = &topology[topologyID];
930       } else {
931         assert(bufferSlot < numTimeSteps);
932         src    = vertices[bufferSlot].getPtr();
933         stride = vertices[bufferSlot].getStride();
934         baseEntry = &vertex_buffer_tags[bufferSlot];
935         topo = &topology[0];
936       }
937 
938       const int* valid = (const int*) valid_i;
939 
940       for (size_t i=0; i<N; i+=4)
941       {
942         vbool4 valid1 = vint4(int(i))+vint4(step) < vint4(int(N));
943         if (valid) valid1 &= vint4::loadu(&valid[i]) == vint4(-1);
944         if (none(valid1)) continue;
945 
946         const vuint4 primID = vuint4::loadu(&primIDs[i]);
947         const vfloat4 uu = vfloat4::loadu(&u[i]);
948         const vfloat4 vv = vfloat4::loadu(&v[i]);
949 
950         foreach_unique(valid1,primID,[&](const vbool4& valid1, const unsigned int primID)
951                        {
952                          for (unsigned int j=0; j<valueCount; j+=4)
953                          {
954                            const size_t M = min(4u,valueCount-j);
955                            isa::PatchEvalSimd<vbool4,vint4,vfloat4,vfloat4>(baseEntry->at(interpolationSlot(primID,j/4,stride)),commitCounter,
956                                                                             topo->getHalfEdge(primID),src+j*sizeof(float),stride,valid1,uu,vv,
957                                                                             P ? P+j*N+i : nullptr,
958                                                                             dPdu ? dPdu+j*N+i : nullptr,
959                                                                             dPdv ? dPdv+j*N+i : nullptr,
960                                                                             ddPdudu ? ddPdudu+j*N+i : nullptr,
961                                                                             ddPdvdv ? ddPdvdv+j*N+i : nullptr,
962                                                                             ddPdudv ? ddPdudv+j*N+i : nullptr,
963                                                                             N,M);
964                          }
965                        });
966       }
967     }
968   }
969 }
970