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 11 namespace 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: CatmullClark1RingTCatmullClark1RingT35 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 */ bytesCatmullClark1RingT39 __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> storeCatmullClark1RingT57 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> loadCatmullClark1RingT62 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 */ serializeCatmullClark1RingT67 __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 */ deserializeCatmullClark1RingT85 __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 hasBorderCatmullClark1RingT103 __forceinline bool hasBorder() const { 104 return border_index != -1; 105 } 106 frontCatmullClark1RingT107 __forceinline const Vertex& front(size_t i) const { 108 assert(edge_valence>i); 109 return ring[i]; 110 } 111 backCatmullClark1RingT112 __forceinline const Vertex& back(size_t i) const { 113 assert(edge_valence>=i); 114 return ring[edge_valence-i]; 115 } 116 has_last_faceCatmullClark1RingT117 __forceinline bool has_last_face() const { 118 return (size_t)border_index != (size_t)edge_valence-2; 119 } 120 has_opposite_frontCatmullClark1RingT121 __forceinline bool has_opposite_front(size_t i) const { 122 return (size_t)border_index != 2*i; 123 } 124 has_opposite_backCatmullClark1RingT125 __forceinline bool has_opposite_back(size_t i) const { 126 return (size_t)border_index != ((size_t)edge_valence-2-2*i); 127 } 128 boundsCatmullClark1RingT129 __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 */ initCatmullClark1RingT138 __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 subdivideCatmullClark1RingT202 __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 isRegular1CatmullClark1RingT306 __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 numEdgeCreasesCatmullClark1RingT316 __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 typeCatmullClark1RingT334 __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 isFinalResolutionCatmullClark1RingT357 __forceinline bool isFinalResolution(float res) const { 358 return vertex_level <= res; 359 } 360 361 /* computes the limit vertex */ getLimitVertexCatmullClark1RingT362 __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 egde vtx -> ring[0] */ getLimitTangentCatmullClark1RingT392 __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 egde vtx -> ring[edge_valence-2] */ getSecondLimitTangentCatmullClark1RingT433 __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 */ getNormalCatmullClark1RingT474 const Vertex getNormal() const { 475 return cross(getLimitTangent(),getSecondLimitTangent()); 476 } 477 478 /* returns center of the n-th quad in the 1-ring */ getQuadCenterCatmullClark1RingT479 __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 */ getEdgeCenterCatmullClark1RingT490 __forceinline Vertex getEdgeCenter(const size_t index) const { 491 return (vtx + ring[index*2]) * 0.5f; 492 } 493 hasValidPositionsCatmullClark1RingT494 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 { FaceGeneralCatmullClark1RingT::Face529 __forceinline Face() {} FaceGeneralCatmullClark1RingT::Face530 __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: GeneralCatmullClark1RingTGeneralCatmullClark1RingT554 GeneralCatmullClark1RingT() 555 : eval_start_face_index(0), eval_start_vertex_index(0), eval_unique_identifier(0) {} 556 isRegularGeneralCatmullClark1RingT557 __forceinline bool isRegular() const 558 { 559 if (border_face == -1 && face_valence == 4) return true; 560 return false; 561 } 562 has_last_faceGeneralCatmullClark1RingT563 __forceinline bool has_last_face() const { 564 return border_face != (int)face_valence-1; 565 } 566 has_second_faceGeneralCatmullClark1RingT567 __forceinline bool has_second_face() const { 568 return (border_face == -1) || (border_face >= 2); 569 } 570 hasValidPositionsGeneralCatmullClark1RingT571 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 initGeneralCatmullClark1RingT580 __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 subdivideGeneralCatmullClark1RingT649 __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 convertGeneralCatmullClark1RingT746 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 egde vtx -> ring[0] */ getLimitTangentGeneralCatmullClark1RingT767 __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 egde vtx -> ring[edge_valence-2] */ getSecondLimitTangentGeneralCatmullClark1RingT783 __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 */ getLimitVertexGeneralCatmullClark1RingT800 __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