1// Copyright 2009-2021 Intel Corporation
2// SPDX-License-Identifier: Apache-2.0
3
4#pragma once
5
6#include "../common/geometry.h"
7#include "../common/buffer.h"
8#include "half_edge.h"
9#include "catmullclark_coefficients.h"
10
11namespace embree
12{
13 struct __aligned(64) FinalQuad {
14 Vec3fa vtx[4];
15 };
16
17 template<typename Vertex, typename Vertex_t = Vertex>
18 struct __aligned(64) CatmullClark1RingT
19 {
20 ALIGNED_STRUCT_(64);
21
22 int border_index; //!< edge index where border starts
23 unsigned int face_valence; //!< number of adjacent quad faces
24 unsigned int edge_valence; //!< number of adjacent edges (2*face_valence)
25 float vertex_crease_weight; //!< weight of vertex crease (0 if no vertex crease)
26 DynamicStackArray<float,16,MAX_RING_FACE_VALENCE> crease_weight; //!< edge crease weights for each adjacent edge
27 float vertex_level; //!< maximum level of all adjacent edges
28 float edge_level; //!< level of first edge
29 unsigned int eval_start_index; //!< topology dependent index to start evaluation
30 unsigned int eval_unique_identifier; //!< topology dependent unique identifier for this ring
31 Vertex vtx; //!< center vertex
32 DynamicStackArray<Vertex,32,MAX_RING_EDGE_VALENCE> ring; //!< ring of neighboring vertices
33
34 public:
35 CatmullClark1RingT ()
36 : eval_start_index(0), eval_unique_identifier(0) {} // FIXME: default constructor should be empty
37
38 /*! calculates number of bytes required to serialize this structure */
39 __forceinline size_t bytes() const
40 {
41 size_t ofs = 0;
42 ofs += sizeof(border_index);
43 ofs += sizeof(face_valence);
44 assert(2*face_valence == edge_valence);
45 ofs += sizeof(vertex_crease_weight);
46 ofs += face_valence*sizeof(float);
47 ofs += sizeof(vertex_level);
48 ofs += sizeof(edge_level);
49 ofs += sizeof(eval_start_index);
50 ofs += sizeof(eval_unique_identifier);
51 ofs += sizeof(vtx);
52 ofs += edge_valence*sizeof(Vertex);
53 return ofs;
54 }
55
56 template<typename Ty>
57 static __forceinline void store(char* ptr, size_t& ofs, const Ty& v) {
58 *(Ty*)&ptr[ofs] = v; ofs += sizeof(Ty);
59 }
60
61 template<typename Ty>
62 static __forceinline void load(char* ptr, size_t& ofs, Ty& v) {
63 v = *(Ty*)&ptr[ofs]; ofs += sizeof(Ty);
64 }
65
66 /*! serializes the ring to some memory location */
67 __forceinline void serialize(char* ptr, size_t& ofs) const
68 {
69 store(ptr,ofs,border_index);
70 store(ptr,ofs,face_valence);
71 store(ptr,ofs,vertex_crease_weight);
72 for (size_t i=0; i<face_valence; i++)
73 store(ptr,ofs,crease_weight[i]);
74 store(ptr,ofs,vertex_level);
75 store(ptr,ofs,edge_level);
76 store(ptr,ofs,eval_start_index);
77 store(ptr,ofs,eval_unique_identifier);
78 Vertex_t::storeu(&ptr[ofs],vtx); ofs += sizeof(Vertex);
79 for (size_t i=0; i<edge_valence; i++) {
80 Vertex_t::storeu(&ptr[ofs],ring[i]); ofs += sizeof(Vertex);
81 }
82 }
83
84 /*! deserializes the ring from some memory location */
85 __forceinline void deserialize(char* ptr, size_t& ofs)
86 {
87 load(ptr,ofs,border_index);
88 load(ptr,ofs,face_valence);
89 edge_valence = 2*face_valence;
90 load(ptr,ofs,vertex_crease_weight);
91 for (size_t i=0; i<face_valence; i++)
92 load(ptr,ofs,crease_weight[i]);
93 load(ptr,ofs,vertex_level);
94 load(ptr,ofs,edge_level);
95 load(ptr,ofs,eval_start_index);
96 load(ptr,ofs,eval_unique_identifier);
97 vtx = Vertex_t::loadu(&ptr[ofs]); ofs += sizeof(Vertex);
98 for (size_t i=0; i<edge_valence; i++) {
99 ring[i] = Vertex_t::loadu(&ptr[ofs]); ofs += sizeof(Vertex);
100 }
101 }
102
103 __forceinline bool hasBorder() const {
104 return border_index != -1;
105 }
106
107 __forceinline const Vertex& front(size_t i) const {
108 assert(edge_valence>i);
109 return ring[i];
110 }
111
112 __forceinline const Vertex& back(size_t i) const {
113 assert(edge_valence>=i);
114 return ring[edge_valence-i];
115 }
116
117 __forceinline bool has_last_face() const {
118 return (size_t)border_index != (size_t)edge_valence-2;
119 }
120
121 __forceinline bool has_opposite_front(size_t i) const {
122 return (size_t)border_index != 2*i;
123 }
124
125 __forceinline bool has_opposite_back(size_t i) const {
126 return (size_t)border_index != ((size_t)edge_valence-2-2*i);
127 }
128
129 __forceinline BBox3fa bounds() const
130 {
131 BBox3fa bounds ( vtx );
132 for (size_t i = 0; i<edge_valence ; i++)
133 bounds.extend( ring[i] );
134 return bounds;
135 }
136
137 /*! initializes the ring from the half edge structure */
138 __forceinline void init(const HalfEdge* const h, const char* vertices, size_t stride)
139 {
140 border_index = -1;
141 vtx = Vertex_t::loadu(vertices+h->getStartVertexIndex()*stride);
142 vertex_crease_weight = h->vertex_crease_weight;
143
144 HalfEdge* p = (HalfEdge*) h;
145
146 unsigned i=0;
147 unsigned min_vertex_index = (unsigned)-1;
148 unsigned min_vertex_index_face = (unsigned)-1;
149 edge_level = p->edge_level;
150 vertex_level = 0.0f;
151
152 do
153 {
154 vertex_level = max(vertex_level,p->edge_level);
155 crease_weight[i/2] = p->edge_crease_weight;
156 assert(p->hasOpposite() || p->edge_crease_weight == float(inf));
157
158 /* store first two vertices of face */
159 p = p->next();
160 const unsigned index0 = p->getStartVertexIndex();
161 ring[i++] = Vertex_t::loadu(vertices+index0*stride);
162 if (index0 < min_vertex_index) { min_vertex_index = index0; min_vertex_index_face = i>>1; }
163 p = p->next();
164
165 const unsigned index1 = p->getStartVertexIndex();
166 ring[i++] = Vertex_t::loadu(vertices+index1*stride);
167 p = p->next();
168
169 /* continue with next face */
170 if (likely(p->hasOpposite()))
171 p = p->opposite();
172
173 /* if there is no opposite go the long way to the other side of the border */
174 else
175 {
176 /* find minimum start vertex */
177 const unsigned index0 = p->getStartVertexIndex();
178 if (index0 < min_vertex_index) { min_vertex_index = index0; min_vertex_index_face = i>>1; }
179
180 /*! mark first border edge and store dummy vertex for face between the two border edges */
181 border_index = i;
182 crease_weight[i/2] = inf;
183 ring[i++] = Vertex_t::loadu(vertices+index0*stride);
184 ring[i++] = vtx; // dummy vertex
185
186 /*! goto other side of border */
187 p = (HalfEdge*) h;
188 while (p->hasOpposite())
189 p = p->opposite()->next();
190 }
191
192 } while (p != h);
193
194 edge_valence = i;
195 face_valence = i >> 1;
196 eval_unique_identifier = min_vertex_index;
197 eval_start_index = min_vertex_index_face;
198
199 assert( hasValidPositions() );
200 }
201
202 __forceinline void subdivide(CatmullClark1RingT& dest) const
203 {
204 dest.edge_level = 0.5f*edge_level;
205 dest.vertex_level = 0.5f*vertex_level;
206 dest.face_valence = face_valence;
207 dest.edge_valence = edge_valence;
208 dest.border_index = border_index;
209 dest.vertex_crease_weight = max(0.0f,vertex_crease_weight-1.0f);
210 dest.eval_start_index = eval_start_index;
211 dest.eval_unique_identifier = eval_unique_identifier;
212
213 /* calculate face points */
214 Vertex_t S = Vertex_t(0.0f);
215 for (size_t i=0; i<face_valence; i++)
216 {
217 size_t face_index = i + eval_start_index; if (face_index >= face_valence) face_index -= face_valence; assert(face_index < face_valence);
218 size_t index0 = 2*face_index+0; if (index0 >= edge_valence) index0 -= edge_valence; assert(index0 < edge_valence);
219 size_t index1 = 2*face_index+1; if (index1 >= edge_valence) index1 -= edge_valence; assert(index1 < edge_valence);
220 size_t index2 = 2*face_index+2; if (index2 >= edge_valence) index2 -= edge_valence; assert(index2 < edge_valence);
221 S += dest.ring[index1] = ((vtx + ring[index1]) + (ring[index0] + ring[index2])) * 0.25f;
222 }
223
224 /* calculate new edge points */
225 size_t num_creases = 0;
226 array_t<size_t,MAX_RING_FACE_VALENCE> crease_id;
227
228 for (size_t i=0; i<face_valence; i++)
229 {
230 size_t face_index = i + eval_start_index;
231 if (face_index >= face_valence) face_index -= face_valence;
232 const float edge_crease = crease_weight[face_index];
233 dest.crease_weight[face_index] = max(edge_crease-1.0f,0.0f);
234
235 size_t index = 2*face_index;
236 size_t prev_index = face_index == 0 ? edge_valence-1 : 2*face_index-1;
237 size_t next_index = 2*face_index+1;
238
239 const Vertex_t v = vtx + ring[index];
240 const Vertex_t f = dest.ring[prev_index] + dest.ring[next_index];
241 S += ring[index];
242
243 /* fast path for regular edge points */
244 if (likely(edge_crease <= 0.0f)) {
245 dest.ring[index] = (v+f) * 0.25f;
246 }
247
248 /* slower path for hard edge rule */
249 else {
250 crease_id[num_creases++] = face_index;
251 dest.ring[index] = v*0.5f;
252
253 /* even slower path for blended edge rule */
254 if (unlikely(edge_crease < 1.0f)) {
255 dest.ring[index] = lerp((v+f)*0.25f,v*0.5f,edge_crease);
256 }
257 }
258 }
259
260 /* compute new vertex using smooth rule */
261 const float inv_face_valence = 1.0f / (float)face_valence;
262 const Vertex_t v_smooth = (Vertex_t) madd(inv_face_valence,S,(float(face_valence)-2.0f)*vtx)*inv_face_valence;
263 dest.vtx = v_smooth;
264
265 /* compute new vertex using vertex_crease_weight rule */
266 if (unlikely(vertex_crease_weight > 0.0f))
267 {
268 if (vertex_crease_weight >= 1.0f) {
269 dest.vtx = vtx;
270 } else {
271 dest.vtx = lerp(v_smooth,vtx,vertex_crease_weight);
272 }
273 return;
274 }
275
276 /* no edge crease rule and dart rule */
277 if (likely(num_creases <= 1))
278 return;
279
280 /* compute new vertex using crease rule */
281 if (likely(num_creases == 2))
282 {
283 /* update vertex using crease rule */
284 const size_t crease0 = crease_id[0], crease1 = crease_id[1];
285 const Vertex_t v_sharp = (Vertex_t)(ring[2*crease0] + 6.0f*vtx + ring[2*crease1]) * (1.0f / 8.0f);
286 dest.vtx = v_sharp;
287
288 /* update crease_weights using chaikin rule */
289 const float crease_weight0 = crease_weight[crease0], crease_weight1 = crease_weight[crease1];
290 dest.crease_weight[crease0] = max(0.25f*(3.0f*crease_weight0 + crease_weight1)-1.0f,0.0f);
291 dest.crease_weight[crease1] = max(0.25f*(3.0f*crease_weight1 + crease_weight0)-1.0f,0.0f);
292
293 /* interpolate between sharp and smooth rule */
294 const float v_blend = 0.5f*(crease_weight0+crease_weight1);
295 if (unlikely(v_blend < 1.0f)) {
296 dest.vtx = lerp(v_smooth,v_sharp,v_blend);
297 }
298 }
299
300 /* compute new vertex using corner rule */
301 else {
302 dest.vtx = vtx;
303 }
304 }
305
306 __forceinline bool isRegular1() const
307 {
308 if (border_index == -1) {
309 if (face_valence == 4) return true;
310 } else {
311 if (face_valence < 4) return true;
312 }
313 return false;
314 }
315
316 __forceinline size_t numEdgeCreases() const
317 {
318 ssize_t numCreases = 0;
319 for (size_t i=0; i<face_valence; i++) {
320 numCreases += crease_weight[i] > 0.0f;
321 }
322 return numCreases;
323 }
324
325 enum Type {
326 TYPE_NONE = 0, //!< invalid type
327 TYPE_REGULAR = 1, //!< regular patch when ignoring creases
328 TYPE_REGULAR_CREASES = 2, //!< regular patch when considering creases
329 TYPE_GREGORY = 4, //!< gregory patch when ignoring creases
330 TYPE_GREGORY_CREASES = 8, //!< gregory patch when considering creases
331 TYPE_CREASES = 16 //!< patch has crease features
332 };
333
334 __forceinline Type type() const
335 {
336 /* check if there is an edge crease anywhere */
337 const size_t numCreases = numEdgeCreases();
338 const bool noInnerCreases = hasBorder() ? numCreases == 2 : numCreases == 0;
339
340 Type crease_mask = (Type) (TYPE_REGULAR | TYPE_GREGORY);
341 if (noInnerCreases ) crease_mask = (Type) (crease_mask | TYPE_REGULAR_CREASES | TYPE_GREGORY_CREASES);
342 if (numCreases != 0) crease_mask = (Type) (crease_mask | TYPE_CREASES);
343
344 /* calculate if this vertex is regular */
345 bool hasBorder = border_index != -1;
346 if (face_valence == 2 && hasBorder) {
347 if (vertex_crease_weight == 0.0f ) return (Type) (crease_mask & (TYPE_REGULAR | TYPE_REGULAR_CREASES | TYPE_GREGORY | TYPE_GREGORY_CREASES | TYPE_CREASES));
348 else if (vertex_crease_weight == float(inf)) return (Type) (crease_mask & (TYPE_REGULAR | TYPE_REGULAR_CREASES | TYPE_GREGORY | TYPE_GREGORY_CREASES | TYPE_CREASES));
349 else return TYPE_CREASES;
350 }
351 else if (vertex_crease_weight != 0.0f) return TYPE_CREASES;
352 else if (face_valence == 3 && hasBorder) return (Type) (crease_mask & (TYPE_REGULAR | TYPE_REGULAR_CREASES | TYPE_GREGORY | TYPE_GREGORY_CREASES | TYPE_CREASES));
353 else if (face_valence == 4 && !hasBorder) return (Type) (crease_mask & (TYPE_REGULAR | TYPE_REGULAR_CREASES | TYPE_GREGORY | TYPE_GREGORY_CREASES | TYPE_CREASES));
354 else return (Type) (crease_mask & (TYPE_GREGORY | TYPE_GREGORY_CREASES | TYPE_CREASES));
355 }
356
357 __forceinline bool isFinalResolution(float res) const {
358 return vertex_level <= res;
359 }
360
361 /* computes the limit vertex */
362 __forceinline Vertex getLimitVertex() const
363 {
364 /* return hard corner */
365 if (unlikely(std::isinf(vertex_crease_weight)))
366 return vtx;
367
368 /* border vertex rule */
369 if (unlikely(border_index != -1))
370 {
371 const unsigned int second_border_index = border_index+2 >= int(edge_valence) ? 0 : border_index+2;
372 return (4.0f * vtx + (ring[border_index] + ring[second_border_index])) * 1.0f/6.0f;
373 }
374
375 Vertex_t F( 0.0f );
376 Vertex_t E( 0.0f );
377
378 assert(eval_start_index < face_valence);
379
380 for (size_t i=0; i<face_valence; i++) {
381 size_t index = i+eval_start_index;
382 if (index >= face_valence) index -= face_valence;
383 F += ring[2*index+1];
384 E += ring[2*index];
385 }
386
387 const float n = (float)face_valence;
388 return (Vertex_t)(n*n*vtx+4.0f*E+F) / ((n+5.0f)*n);
389 }
390
391 /* gets limit tangent in the direction of edge vtx -> ring[0] */
392 __forceinline Vertex getLimitTangent() const
393 {
394 if (unlikely(std::isinf(vertex_crease_weight)))
395 return ring[0] - vtx;
396
397 /* border vertex rule */
398 if (unlikely(border_index != -1))
399 {
400 if (border_index != (int)edge_valence-2 ) {
401 return ring[0] - vtx;
402 }
403 else
404 {
405 const unsigned int second_border_index = border_index+2 >= int(edge_valence) ? 0 : border_index+2;
406 return (ring[second_border_index] - ring[border_index]) * 0.5f;
407 }
408 }
409
410 Vertex_t alpha( 0.0f );
411 Vertex_t beta ( 0.0f );
412
413 const size_t n = face_valence;
414
415 assert(eval_start_index < face_valence);
416
417 Vertex_t q( 0.0f );
418 for (size_t i=0; i<face_valence; i++)
419 {
420 size_t index = i+eval_start_index;
421 if (index >= face_valence) index -= face_valence;
422 const float a = CatmullClarkPrecomputedCoefficients::table.limittangent_a(index,n);
423 const float b = CatmullClarkPrecomputedCoefficients::table.limittangent_b(index,n);
424 alpha += a * ring[2*index];
425 beta += b * ring[2*index+1];
426 }
427
428 const float sigma = CatmullClarkPrecomputedCoefficients::table.limittangent_c(n);
429 return sigma * (alpha + beta);
430 }
431
432 /* gets limit tangent in the direction of edge vtx -> ring[edge_valence-2] */
433 __forceinline Vertex getSecondLimitTangent() const
434 {
435 if (unlikely(std::isinf(vertex_crease_weight)))
436 return ring[2] - vtx;
437
438 /* border vertex rule */
439 if (unlikely(border_index != -1))
440 {
441 if (border_index != 2) {
442 return ring[2] - vtx;
443 }
444 else {
445 const unsigned int second_border_index = border_index+2 >= int(edge_valence) ? 0 : border_index+2;
446 return (ring[border_index] - ring[second_border_index]) * 0.5f;
447 }
448 }
449
450 Vertex_t alpha( 0.0f );
451 Vertex_t beta ( 0.0f );
452
453 const size_t n = face_valence;
454
455 assert(eval_start_index < face_valence);
456
457 for (size_t i=0; i<face_valence; i++)
458 {
459 size_t index = i+eval_start_index;
460 if (index >= face_valence) index -= face_valence;
461
462 size_t prev_index = index == 0 ? face_valence-1 : index-1; // need to be bit-wise exact in cosf eval
463 const float a = CatmullClarkPrecomputedCoefficients::table.limittangent_a(prev_index,n);
464 const float b = CatmullClarkPrecomputedCoefficients::table.limittangent_b(prev_index,n);
465 alpha += a * ring[2*index];
466 beta += b * ring[2*index+1];
467 }
468
469 const float sigma = CatmullClarkPrecomputedCoefficients::table.limittangent_c(n);
470 return sigma* (alpha + beta);
471 }
472
473 /* gets surface normal */
474 const Vertex getNormal() const {
475 return cross(getLimitTangent(),getSecondLimitTangent());
476 }
477
478 /* returns center of the n-th quad in the 1-ring */
479 __forceinline Vertex getQuadCenter(const size_t index) const
480 {
481 const Vertex_t &p0 = vtx;
482 const Vertex_t &p1 = ring[2*index+0];
483 const Vertex_t &p2 = ring[2*index+1];
484 const Vertex_t &p3 = index == face_valence-1 ? ring[0] : ring[2*index+2];
485 const Vertex p = (p0+p1+p2+p3) * 0.25f;
486 return p;
487 }
488
489 /* returns center of the n-th edge in the 1-ring */
490 __forceinline Vertex getEdgeCenter(const size_t index) const {
491 return (vtx + ring[index*2]) * 0.5f;
492 }
493
494 bool hasValidPositions() const
495 {
496 for (size_t i=0; i<edge_valence; i++) {
497 if (!isvalid(ring[i]))
498 return false;
499 }
500 return true;
501 }
502
503 friend __forceinline embree_ostream operator<<(embree_ostream o, const CatmullClark1RingT &c)
504 {
505 o << "vtx " << c.vtx << " size = " << c.edge_valence << ", " <<
506 "hard_edge = " << c.border_index << ", face_valence " << c.face_valence <<
507 ", edge_level = " << c.edge_level << ", vertex_level = " << c.vertex_level << ", eval_start_index: " << c.eval_start_index << ", ring: " << embree_endl;
508
509 for (unsigned int i=0; i<min(c.edge_valence,(unsigned int)MAX_RING_FACE_VALENCE); i++) {
510 o << i << " -> " << c.ring[i];
511 if (i % 2 == 0) o << " crease = " << c.crease_weight[i/2];
512 o << embree_endl;
513 }
514 return o;
515 }
516 };
517
518 typedef CatmullClark1RingT<Vec3fa,Vec3fa_t> CatmullClark1Ring3fa;
519
520 template<typename Vertex, typename Vertex_t = Vertex>
521 struct __aligned(64) GeneralCatmullClark1RingT
522 {
523 ALIGNED_STRUCT_(64);
524
525 typedef CatmullClark1RingT<Vertex,Vertex_t> CatmullClark1Ring;
526
527 struct Face
528 {
529 __forceinline Face() {}
530 __forceinline Face (int size, float crease_weight)
531 : size(size), crease_weight(crease_weight) {}
532
533 // FIXME: add member that returns total number of vertices
534
535 int size; // number of vertices-2 of nth face in ring
536 float crease_weight;
537 };
538
539 Vertex vtx;
540 DynamicStackArray<Vertex,32,MAX_RING_EDGE_VALENCE> ring;
541 DynamicStackArray<Face,16,MAX_RING_FACE_VALENCE> faces;
542 unsigned int face_valence;
543 unsigned int edge_valence;
544 int border_face;
545 float vertex_crease_weight;
546 float vertex_level; //!< maximum level of adjacent edges
547 float edge_level; // level of first edge
548 bool only_quads; // true if all faces are quads
549 unsigned int eval_start_face_index;
550 unsigned int eval_start_vertex_index;
551 unsigned int eval_unique_identifier;
552
553 public:
554 GeneralCatmullClark1RingT()
555 : eval_start_face_index(0), eval_start_vertex_index(0), eval_unique_identifier(0) {}
556
557 __forceinline bool isRegular() const
558 {
559 if (border_face == -1 && face_valence == 4) return true;
560 return false;
561 }
562
563 __forceinline bool has_last_face() const {
564 return border_face != (int)face_valence-1;
565 }
566
567 __forceinline bool has_second_face() const {
568 return (border_face == -1) || (border_face >= 2);
569 }
570
571 bool hasValidPositions() const
572 {
573 for (size_t i=0; i<edge_valence; i++) {
574 if (!isvalid(ring[i]))
575 return false;
576 }
577 return true;
578 }
579
580 __forceinline void init(const HalfEdge* const h, const char* vertices, size_t stride)
581 {
582 only_quads = true;
583 border_face = -1;
584 vtx = Vertex_t::loadu(vertices+h->getStartVertexIndex()*stride);
585 vertex_crease_weight = h->vertex_crease_weight;
586 HalfEdge* p = (HalfEdge*) h;
587
588 unsigned int e=0, f=0;
589 unsigned min_vertex_index = (unsigned)-1;
590 unsigned min_vertex_index_face = (unsigned)-1;
591 unsigned min_vertex_index_vertex = (unsigned)-1;
592 edge_level = p->edge_level;
593 vertex_level = 0.0f;
594 do
595 {
596 HalfEdge* p_prev = p->prev();
597 HalfEdge* p_next = p->next();
598 const float crease_weight = p->edge_crease_weight;
599 assert(p->hasOpposite() || p->edge_crease_weight == float(inf));
600 vertex_level = max(vertex_level,p->edge_level);
601
602 /* find minimum start vertex */
603 unsigned vertex_index = p_next->getStartVertexIndex();
604 if (vertex_index < min_vertex_index) { min_vertex_index = vertex_index; min_vertex_index_face = f; min_vertex_index_vertex = e; }
605
606 /* store first N-2 vertices of face */
607 unsigned int vn = 0;
608 for (p = p_next; p!=p_prev; p=p->next()) {
609 ring[e++] = Vertex_t::loadu(vertices+p->getStartVertexIndex()*stride);
610 vn++;
611 }
612 faces[f++] = Face(vn,crease_weight);
613 only_quads &= (vn == 2);
614
615 /* continue with next face */
616 if (likely(p->hasOpposite()))
617 p = p->opposite();
618
619 /* if there is no opposite go the long way to the other side of the border */
620 else
621 {
622 /* find minimum start vertex */
623 unsigned vertex_index = p->getStartVertexIndex();
624 if (vertex_index < min_vertex_index) { min_vertex_index = vertex_index; min_vertex_index_face = f; min_vertex_index_vertex = e; }
625
626 /*! mark first border edge and store dummy vertex for face between the two border edges */
627 border_face = f;
628 faces[f++] = Face(2,inf);
629 ring[e++] = Vertex_t::loadu(vertices+p->getStartVertexIndex()*stride);
630 ring[e++] = vtx; // dummy vertex
631
632 /*! goto other side of border */
633 p = (HalfEdge*) h;
634 while (p->hasOpposite())
635 p = p->opposite()->next();
636 }
637
638 } while (p != h);
639
640 edge_valence = e;
641 face_valence = f;
642 eval_unique_identifier = min_vertex_index;
643 eval_start_face_index = min_vertex_index_face;
644 eval_start_vertex_index = min_vertex_index_vertex;
645
646 assert( hasValidPositions() );
647 }
648
649 __forceinline void subdivide(CatmullClark1Ring& dest) const
650 {
651 dest.edge_level = 0.5f*edge_level;
652 dest.vertex_level = 0.5f*vertex_level;
653 dest.face_valence = face_valence;
654 dest.edge_valence = 2*face_valence;
655 dest.border_index = border_face == -1 ? -1 : 2*border_face; // FIXME:
656 dest.vertex_crease_weight = max(0.0f,vertex_crease_weight-1.0f);
657 dest.eval_start_index = eval_start_face_index;
658 dest.eval_unique_identifier = eval_unique_identifier;
659 assert(dest.face_valence <= MAX_RING_FACE_VALENCE);
660
661 /* calculate face points */
662 Vertex_t S = Vertex_t(0.0f);
663 for (size_t face=0, v=eval_start_vertex_index; face<face_valence; face++) {
664 size_t f = (face + eval_start_face_index)%face_valence;
665
666 Vertex_t F = vtx;
667 for (size_t k=v; k<=v+faces[f].size; k++) F += ring[k%edge_valence]; // FIXME: optimize
668 S += dest.ring[2*f+1] = F/float(faces[f].size+2);
669 v+=faces[f].size;
670 v%=edge_valence;
671 }
672
673 /* calculate new edge points */
674 size_t num_creases = 0;
675 array_t<size_t,MAX_RING_FACE_VALENCE> crease_id;
676 Vertex_t C = Vertex_t(0.0f);
677 for (size_t face=0, j=eval_start_vertex_index; face<face_valence; face++)
678 {
679 size_t i = (face + eval_start_face_index)%face_valence;
680
681 const Vertex_t v = vtx + ring[j];
682 Vertex_t f = dest.ring[2*i+1];
683 if (i == 0) f += dest.ring[dest.edge_valence-1];
684 else f += dest.ring[2*i-1];
685 S += ring[j];
686 dest.crease_weight[i] = max(faces[i].crease_weight-1.0f,0.0f);
687
688 /* fast path for regular edge points */
689 if (likely(faces[i].crease_weight <= 0.0f)) {
690 dest.ring[2*i] = (v+f) * 0.25f;
691 }
692
693 /* slower path for hard edge rule */
694 else {
695 C += ring[j]; crease_id[num_creases++] = i;
696 dest.ring[2*i] = v*0.5f;
697
698 /* even slower path for blended edge rule */
699 if (unlikely(faces[i].crease_weight < 1.0f)) {
700 dest.ring[2*i] = lerp((v+f)*0.25f,v*0.5f,faces[i].crease_weight);
701 }
702 }
703 j+=faces[i].size;
704 j%=edge_valence;
705 }
706
707 /* compute new vertex using smooth rule */
708 const float inv_face_valence = 1.0f / (float)face_valence;
709 const Vertex_t v_smooth = (Vertex_t) madd(inv_face_valence,S,(float(face_valence)-2.0f)*vtx)*inv_face_valence;
710 dest.vtx = v_smooth;
711
712 /* compute new vertex using vertex_crease_weight rule */
713 if (unlikely(vertex_crease_weight > 0.0f))
714 {
715 if (vertex_crease_weight >= 1.0f) {
716 dest.vtx = vtx;
717 } else {
718 dest.vtx = lerp(vtx,v_smooth,vertex_crease_weight);
719 }
720 return;
721 }
722
723 if (likely(num_creases <= 1))
724 return;
725
726 /* compute new vertex using crease rule */
727 if (likely(num_creases == 2)) {
728 const Vertex_t v_sharp = (Vertex_t)(C + 6.0f * vtx) * (1.0f / 8.0f);
729 const float crease_weight0 = faces[crease_id[0]].crease_weight;
730 const float crease_weight1 = faces[crease_id[1]].crease_weight;
731 dest.vtx = v_sharp;
732 dest.crease_weight[crease_id[0]] = max(0.25f*(3.0f*crease_weight0 + crease_weight1)-1.0f,0.0f);
733 dest.crease_weight[crease_id[1]] = max(0.25f*(3.0f*crease_weight1 + crease_weight0)-1.0f,0.0f);
734 const float v_blend = 0.5f*(crease_weight0+crease_weight1);
735 if (unlikely(v_blend < 1.0f)) {
736 dest.vtx = lerp(v_sharp,v_smooth,v_blend);
737 }
738 }
739
740 /* compute new vertex using corner rule */
741 else {
742 dest.vtx = vtx;
743 }
744 }
745
746 void convert(CatmullClark1Ring& dst) const
747 {
748 dst.edge_level = edge_level;
749 dst.vertex_level = vertex_level;
750 dst.vtx = vtx;
751 dst.face_valence = face_valence;
752 dst.edge_valence = 2*face_valence;
753 dst.border_index = border_face == -1 ? -1 : 2*border_face;
754 for (size_t i=0; i<face_valence; i++)
755 dst.crease_weight[i] = faces[i].crease_weight;
756 dst.vertex_crease_weight = vertex_crease_weight;
757 for (size_t i=0; i<edge_valence; i++) dst.ring[i] = ring[i];
758
759 dst.eval_start_index = eval_start_face_index;
760 dst.eval_unique_identifier = eval_unique_identifier;
761
762 assert( dst.hasValidPositions() );
763 }
764
765
766 /* gets limit tangent in the direction of edge vtx -> ring[0] */
767 __forceinline Vertex getLimitTangent() const
768 {
769 CatmullClark1Ring cc_vtx;
770
771 /* fast path for quad only rings */
772 if (only_quads)
773 {
774 convert(cc_vtx);
775 return cc_vtx.getLimitTangent();
776 }
777
778 subdivide(cc_vtx);
779 return 2.0f * cc_vtx.getLimitTangent();
780 }
781
782 /* gets limit tangent in the direction of edge vtx -> ring[edge_valence-2] */
783 __forceinline Vertex getSecondLimitTangent() const
784 {
785 CatmullClark1Ring cc_vtx;
786
787 /* fast path for quad only rings */
788 if (only_quads)
789 {
790 convert(cc_vtx);
791 return cc_vtx.getSecondLimitTangent();
792 }
793
794 subdivide(cc_vtx);
795 return 2.0f * cc_vtx.getSecondLimitTangent();
796 }
797
798
799 /* gets limit vertex */
800 __forceinline Vertex getLimitVertex() const
801 {
802 CatmullClark1Ring cc_vtx;
803
804 /* fast path for quad only rings */
805 if (only_quads)
806 convert(cc_vtx);
807 else
808 subdivide(cc_vtx);
809 return cc_vtx.getLimitVertex();
810 }
811
812 friend __forceinline embree_ostream operator<<(embree_ostream o, const GeneralCatmullClark1RingT &c)
813 {
814 o << "vtx " << c.vtx << " size = " << c.edge_valence << ", border_face = " << c.border_face << ", " << " face_valence = " << c.face_valence <<
815 ", edge_level = " << c.edge_level << ", vertex_level = " << c.vertex_level << ", ring: " << embree_endl;
816 for (size_t v=0, f=0; f<c.face_valence; v+=c.faces[f++].size) {
817 for (size_t i=v; i<v+c.faces[f].size; i++) {
818 o << i << " -> " << c.ring[i];
819 if (i == v) o << " crease = " << c.faces[f].crease_weight;
820 o << embree_endl;
821 }
822 }
823 return o;
824 }
825 };
826}
827