forked from bartvdbraak/blender
da376e0237
Cycles uses code from some great open source projects, many thanks them: * BVH building and traversal code from NVidia's "Understanding the Efficiency of Ray Traversal on GPUs": http://code.google.com/p/understanding-the-efficiency-of-ray-traversal-on-gpus/ * Open Shading Language for a large part of the shading system: http://code.google.com/p/openshadinglanguage/ * Blender for procedural textures and a few other nodes. * Approximate Catmull Clark subdivision from NVidia Mesh tools: http://code.google.com/p/nvidia-mesh-tools/ * Sobol direction vectors from: http://web.maths.unsw.edu.au/~fkuo/sobol/ * Film response functions from: http://www.cs.columbia.edu/CAVE/software/softlib/dorf.php
667 lines
18 KiB
C++
667 lines
18 KiB
C++
/*
|
|
* Copyright 2006, NVIDIA Corporation Ignacio Castano <icastano@nvidia.com>
|
|
*
|
|
* Modifications copyright (c) 2011, Blender Foundation. All rights reserved.
|
|
*
|
|
* Permission is hereby granted, free of charge, to any person obtaining a copy
|
|
* of this software and associated documentation files (the "Software"), to deal
|
|
* in the Software without restriction, including without limitation the rights
|
|
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|
* copies of the Software, and to permit persons to whom the Software is
|
|
* furnished to do so, subject to the following conditions:
|
|
*
|
|
* The above copyright notice and this permission notice shall be included in
|
|
* all copies or substantial portions of the Software.
|
|
*
|
|
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
|
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
|
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
|
|
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
|
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
|
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
|
|
* THE SOFTWARE.
|
|
*/
|
|
|
|
#include "subd_build.h"
|
|
#include "subd_edge.h"
|
|
#include "subd_face.h"
|
|
#include "subd_ring.h"
|
|
#include "subd_mesh.h"
|
|
#include "subd_patch.h"
|
|
#include "subd_stencil.h"
|
|
#include "subd_vert.h"
|
|
|
|
#include "util_algorithm.h"
|
|
#include "util_debug.h"
|
|
#include "util_math.h"
|
|
#include "util_string.h"
|
|
|
|
CCL_NAMESPACE_BEGIN
|
|
|
|
/* Subd Builder */
|
|
|
|
SubdBuilder *SubdBuilder::create(bool linear)
|
|
{
|
|
if(linear)
|
|
return new SubdLinearBuilder();
|
|
else
|
|
return new SubdAccBuilder();
|
|
}
|
|
|
|
/* Gregory ACC Stencil */
|
|
|
|
class GregoryAccStencil {
|
|
public:
|
|
SubdFaceRing *ring;
|
|
StencilMask stencil[20];
|
|
|
|
GregoryAccStencil(SubdFaceRing *ring_)
|
|
{
|
|
ring = ring_;
|
|
|
|
for(int i = 0; i < 20; i++)
|
|
stencil[i].resize(ring->num_verts());
|
|
}
|
|
|
|
StencilMask& get(int i)
|
|
{
|
|
assert(i < 20);
|
|
return stencil[i];
|
|
}
|
|
|
|
float& get(int i, SubdVert *vert)
|
|
{
|
|
assert(i < 20);
|
|
return stencil[i][ring->vert_index(vert)];
|
|
}
|
|
};
|
|
|
|
static float pseudoValence(SubdVert *vert)
|
|
{
|
|
float valence = (float)vert->valence();
|
|
|
|
if(vert->is_boundary()) {
|
|
/* we treat boundary verts as being half a closed mesh. corners are
|
|
special case. n = 4 for corners and n = 2*(n-1) for boundaries. */
|
|
if(valence == 2) return 4;
|
|
return (valence - 1)*2;
|
|
}
|
|
|
|
return valence;
|
|
}
|
|
|
|
/* Subd ACC Builder */
|
|
|
|
SubdAccBuilder::SubdAccBuilder()
|
|
{
|
|
}
|
|
|
|
SubdAccBuilder::~SubdAccBuilder()
|
|
{
|
|
}
|
|
|
|
Patch *SubdAccBuilder::run(SubdFace *face)
|
|
{
|
|
SubdFaceRing ring(face, face->edge);
|
|
GregoryAccStencil stencil(&ring);
|
|
float3 position[20];
|
|
|
|
computeCornerStencil(&ring, &stencil);
|
|
computeEdgeStencil(&ring, &stencil);
|
|
computeInteriorStencil(&ring, &stencil);
|
|
|
|
ring.evaluate_stencils(position, stencil.stencil, 20);
|
|
|
|
if(face->num_edges() == 3) {
|
|
GregoryTrianglePatch *patch = new GregoryTrianglePatch();
|
|
memcpy(patch->hull, position, sizeof(float3)*20);
|
|
return patch;
|
|
}
|
|
else if(face->num_edges() == 4) {
|
|
GregoryQuadPatch *patch = new GregoryQuadPatch();
|
|
memcpy(patch->hull, position, sizeof(float3)*20);
|
|
return patch;
|
|
}
|
|
|
|
assert(0); /* n-gons should have been split already */
|
|
return NULL;
|
|
}
|
|
|
|
/* Gregory Patch */
|
|
|
|
void SubdAccBuilder::computeCornerStencil(SubdFaceRing *ring, GregoryAccStencil *stencil)
|
|
{
|
|
const int cornerIndices[7] = {8, 11, 19, 16, 6, 9, 12};
|
|
int primitiveOffset = ring->is_quad()? 0: 4;
|
|
|
|
SubdEdge *firstEdge = ring->firstEdge();
|
|
|
|
/* compute corner control points */
|
|
int v = 0;
|
|
|
|
for(SubdFace::EdgeIterator it(firstEdge); !it.isDone(); it.advance(), v++) {
|
|
SubdVert *vert = it.current()->from();
|
|
int valence = vert->valence();
|
|
int cid = cornerIndices[primitiveOffset+v];
|
|
|
|
if(vert->is_boundary()) {
|
|
/* compute vertex limit position */
|
|
SubdEdge *edge0 = vert->edge;
|
|
SubdEdge *edge1 = vert->edge->prev;
|
|
|
|
assert(edge0->face == NULL);
|
|
assert(edge0->to() != vert);
|
|
assert(edge1->face == NULL);
|
|
assert(edge1->from() != vert);
|
|
|
|
stencil->get(cid, vert) = 2.0f/3.0f;
|
|
stencil->get(cid, edge0->to()) = 1.0f/6.0f;
|
|
stencil->get(cid, edge1->from()) = 1.0f/6.0f;
|
|
|
|
assert(stencil->get(cid).is_normalized());
|
|
}
|
|
else {
|
|
stencil->get(cid, vert) = 3.0f*valence*valence;
|
|
|
|
for(SubdVert::EdgeIterator eit(vert->edge); !eit.isDone(); eit.advance()) {
|
|
SubdEdge *edge = eit.current();
|
|
assert(vert->co == edge->from()->co);
|
|
|
|
stencil->get(cid, edge->to()) = 12.0f;
|
|
|
|
if(SubdFaceRing::is_triangle(edge->face)) {
|
|
/* distribute weight to all verts */
|
|
stencil->get(cid, vert) += 1.0f;
|
|
stencil->get(cid, edge->to()) += 1.0f;
|
|
stencil->get(cid, edge->next->to()) += 1.0f;
|
|
}
|
|
else
|
|
stencil->get(cid, edge->next->to()) = 3.0f;
|
|
}
|
|
|
|
/* normalize stencil. */
|
|
stencil->get(cid).normalize();
|
|
}
|
|
}
|
|
}
|
|
|
|
void SubdAccBuilder::computeEdgeStencil(SubdFaceRing *ring, GregoryAccStencil *stencil)
|
|
{
|
|
const int cornerIndices[7] = {8, 11, 19, 16, 6, 9, 12};
|
|
const int edge1Indices[7] = {9, 13, 18, 14, 7, 10, 13};
|
|
const int edge2Indices[7] = {12, 10, 15, 17, 14, 8, 11};
|
|
int primitiveOffset = ring->is_quad()? 0: 4;
|
|
|
|
float tangentScales[14] = {
|
|
0.0f, 0.0f, 0.0f, 0.667791f, 1.0f,
|
|
1.11268f, 1.1284f, 1.10289f, 1.06062f,
|
|
1.01262f, 0.963949f, 0.916926f, 0.872541f, 0.831134f
|
|
};
|
|
|
|
SubdEdge *firstEdge = ring->firstEdge();
|
|
|
|
/* compute corner / edge control points */
|
|
int v = 0;
|
|
|
|
for(SubdFace::EdgeIterator it(firstEdge); !it.isDone(); it.advance(), v++) {
|
|
SubdVert *vert = it.current()->from();
|
|
int valence = vert->valence();
|
|
int cid = cornerIndices[primitiveOffset+v];
|
|
|
|
int i1 = 0, i2 = 0, j = 0;
|
|
|
|
for(SubdVert::EdgeIterator eit(vert->edge); !eit.isDone(); eit.advance(), j++) {
|
|
SubdEdge *edge = eit.current();
|
|
|
|
/* find index of "our" edge for edge control points */
|
|
if(edge == it.current())
|
|
i1 = j;
|
|
if(edge == it.current()->prev->pair)
|
|
i2 = j;
|
|
}
|
|
|
|
if(vert->is_boundary()) {
|
|
int num_verts = ring->num_verts();
|
|
StencilMask r0(num_verts);
|
|
StencilMask r1(num_verts);
|
|
|
|
computeBoundaryTangentStencils(ring, vert, r0, r1);
|
|
|
|
int k = valence - 1;
|
|
float omega = M_PI / k;
|
|
|
|
int eid1 = edge1Indices[primitiveOffset + v];
|
|
int eid2 = edge2Indices[primitiveOffset + v];
|
|
|
|
if(it.current()->is_boundary()) {
|
|
assert(it.current()->from() == vert);
|
|
|
|
stencil->get(eid1, vert) = 2.0f / 3.0f;
|
|
stencil->get(eid1, it.current()->to()) = 1.0f / 3.0f;
|
|
|
|
assert(stencil->get(eid1).is_normalized());
|
|
|
|
if(valence == 2) {
|
|
for(int i = 0; i < num_verts; i++)
|
|
stencil->get(eid1)[i] += r0[i] * 0.0001f;
|
|
}
|
|
}
|
|
else {
|
|
stencil->get(eid1) = stencil->get(cid);
|
|
|
|
/* compute index of it.current() around vert */
|
|
int idx = 0;
|
|
|
|
for(SubdVert::EdgeIterator eit(vert->edges()); !eit.isDone(); eit.advance(), idx++)
|
|
if(eit.current() == it.current())
|
|
break;
|
|
|
|
assert(idx != valence);
|
|
|
|
float c = cosf(idx * omega);
|
|
float s = sinf(idx * omega);
|
|
|
|
for(int i = 0; i < num_verts; i++)
|
|
stencil->get(eid1)[i] += (r0[i] * s + r1[i] * c) / 3.0f;
|
|
}
|
|
|
|
if(it.current()->prev->is_boundary()) {
|
|
assert(it.current()->prev->pair->from() == vert);
|
|
|
|
stencil->get(eid2, vert) = 2.0f / 3.0f;
|
|
stencil->get(eid2, it.current()->prev->pair->to()) = 1.0f / 3.0f;
|
|
|
|
assert(stencil->get(eid2).is_normalized());
|
|
|
|
if(valence == 2) {
|
|
for(int i = 0; i < num_verts; i++)
|
|
stencil->get(eid2)[i] += r0[i] * 0.0001f;
|
|
}
|
|
}
|
|
else {
|
|
stencil->get(eid2) = stencil->get(cid);
|
|
|
|
/* compute index of it.current() around vert */
|
|
int idx = 0;
|
|
|
|
for(SubdVert::EdgeIterator eit(vert->edges()); !eit.isDone(); eit.advance(), idx++)
|
|
if(eit.current() == it.current()->prev->pair)
|
|
break;
|
|
|
|
assert(idx != valence);
|
|
|
|
float c = cosf(idx * omega);
|
|
float s = sinf(idx * omega);
|
|
|
|
for(int i = 0; i < num_verts; i++)
|
|
stencil->get(eid2)[i] += (r0[i] * s + r1[i] * c) / 3;
|
|
}
|
|
}
|
|
else {
|
|
float costerm = cosf(M_PI / valence);
|
|
float sqrtterm = sqrtf(4.0f + costerm*costerm);
|
|
|
|
/* float tangentScale = 1.0f; */
|
|
float tangentScale = tangentScales[min(valence, 13U)];
|
|
|
|
float alpha = (1.0f + costerm / sqrtterm) / (3.0f * valence) * tangentScale;
|
|
float beta = 1.0f / (3.0f * valence * sqrtterm) * tangentScale;
|
|
|
|
|
|
int eid1 = edge1Indices[primitiveOffset + v];
|
|
int eid2 = edge2Indices[primitiveOffset + v];
|
|
|
|
stencil->get(eid1) = stencil->get(cid);
|
|
stencil->get(eid2) = stencil->get(cid);
|
|
|
|
int j = 0;
|
|
for(SubdVert::EdgeIterator eit(vert->edges()); !eit.isDone(); eit.advance(), j++) {
|
|
SubdEdge *edge = eit.current();
|
|
assert(vert->co == edge->from()->co);
|
|
|
|
float costerm1_a = cosf(M_PI * 2 * (j-i1) / valence);
|
|
float costerm1_b = cosf(M_PI * (2 * (j-i1)-1) / valence); /* -1 instead of +1 b/c of edge->next->to() */
|
|
|
|
float costerm2_a = cosf(M_PI * 2 * (j-i2) / valence);
|
|
float costerm2_b = cosf(M_PI * (2 * (j-i2)-1) / valence); /* -1 instead of +1 b/c of edge->next->to() */
|
|
|
|
|
|
stencil->get(eid1, edge->to()) += alpha * costerm1_a;
|
|
stencil->get(eid2, edge->to()) += alpha * costerm2_a;
|
|
|
|
if(SubdFaceRing::is_triangle(edge->face)) {
|
|
/* @@ this probably does not provide watertight results!! (1/3 + 1/3 + 1/3 != 1) */
|
|
|
|
/* distribute weight to all verts */
|
|
stencil->get(eid1, vert) += beta * costerm1_b / 3.0f;
|
|
stencil->get(eid1, edge->to()) += beta * costerm1_b / 3.0f;
|
|
stencil->get(eid1, edge->next->to()) += beta * costerm1_b / 3.0f;
|
|
|
|
stencil->get(eid2, vert) += beta * costerm2_b / 3.0f;
|
|
stencil->get(eid2, edge->to()) += beta * costerm2_b / 3.0f;
|
|
stencil->get(eid2, edge->next->to()) += beta * costerm2_b / 3.0f;
|
|
}
|
|
else {
|
|
stencil->get(eid1, edge->next->to()) += beta * costerm1_b;
|
|
stencil->get(eid2, edge->next->to()) += beta * costerm2_b;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
void SubdAccBuilder::computeInteriorStencil(SubdFaceRing *ring, GregoryAccStencil *stencil)
|
|
{
|
|
static int corner1Indices[7] = {8, 11, 19, 16, 6, 9, 12};
|
|
static int corner2Indices[7] = {11, 19, 16, 8, 9, 12, 6};
|
|
static int edge1Indices[7] = {9, 13, 18, 14, 7, 10, 13};
|
|
static int edge2Indices[7] = {10, 15, 17, 12, 8, 11, 14};
|
|
static int interior1Indices[7] = {1, 3, 6, 4, 1, 3, 5};
|
|
static int interior2Indices[7] = {2, 7, 5, 0, 2, 4, 0};
|
|
|
|
int primitiveOffset = ring->is_quad()? 0: 4;
|
|
|
|
SubdFace * face = ring->face();
|
|
SubdEdge *firstEdge = ring->firstEdge();
|
|
|
|
/* interior control points */
|
|
int v = 0;
|
|
for(SubdFace::EdgeIterator it(firstEdge); !it.isDone(); it.advance(), v++) {
|
|
SubdEdge *edge = it.current();
|
|
|
|
if(edge->is_boundary()) {
|
|
float valence1 = pseudoValence(edge->from());
|
|
float valence2 = pseudoValence(edge->to());
|
|
|
|
float weights1[4];
|
|
float weights2[4];
|
|
|
|
if(ring->is_quad()) {
|
|
weights1[0] = 3 * valence1;
|
|
weights1[1] = 6;
|
|
weights1[2] = 3;
|
|
weights1[3] = 6;
|
|
|
|
weights2[0] = 6;
|
|
weights2[1] = 3 * valence2;
|
|
weights2[2] = 6;
|
|
weights2[3] = 3;
|
|
}
|
|
else {
|
|
assert(ring->is_triangle());
|
|
weights1[0] = 3 * valence1 + 1;
|
|
weights1[1] = 7;
|
|
weights1[2] = 7;
|
|
|
|
weights2[0] = 7;
|
|
weights2[1] = 3 * valence2 + 1;
|
|
weights2[2] = 7;
|
|
}
|
|
|
|
int idx1 = interior1Indices[primitiveOffset+v];
|
|
int idx2 = interior2Indices[primitiveOffset+v];
|
|
|
|
int i = 0;
|
|
for(SubdFace::EdgeIterator it(face->edges(edge)); !it.isDone(); it.advance(), i++) {
|
|
SubdVert *vert = it.current()->from();
|
|
stencil->get(idx1, vert) += weights1[i];
|
|
stencil->get(idx2, vert) += weights2[i];
|
|
}
|
|
|
|
stencil->get(idx1).normalize();
|
|
stencil->get(idx2).normalize();
|
|
}
|
|
else {
|
|
SubdVert *e0 = edge->from();
|
|
float costerm0 = cosf(2.0f * M_PI / pseudoValence(e0));
|
|
|
|
SubdVert *f0 = edge->to();
|
|
float costerm1 = cosf(2.0f * M_PI / pseudoValence(f0));
|
|
|
|
/* p0 +------+ q0
|
|
* | |
|
|
* f0 +======+ e0 <=== current edge
|
|
* | |
|
|
* p1 +------+ q1
|
|
*/
|
|
|
|
SubdVert *q0 = edge->next->to();
|
|
SubdVert *p0 = edge->prev->from();
|
|
|
|
SubdVert *p1 = edge->pair->next->to();
|
|
SubdVert *q1 = edge->pair->prev->from();
|
|
|
|
|
|
StencilMask x(ring->num_verts());
|
|
StencilMask y(ring->num_verts());
|
|
|
|
for(int i = 0; i < ring->num_verts(); i++) {
|
|
x[i] =
|
|
(costerm1 * stencil->get(corner1Indices[primitiveOffset+v])[i] -
|
|
(2*costerm0 + costerm1) * stencil->get(edge1Indices[primitiveOffset+v])[i] +
|
|
2*costerm0 * stencil->get(edge2Indices[primitiveOffset+v])[i]) / 3.0f;
|
|
}
|
|
|
|
/* y = (2*( midedgeA1 - midedgeB1) + 4*(centroidA - centroidB))/18.0f; */
|
|
y[ring->vert_index(p0)] = 1;
|
|
y[ring->vert_index(p1)] = -1;
|
|
|
|
/* add centroidA */
|
|
if(ring->is_triangle()) {
|
|
y[ring->vert_index(p0)] += 4.0f / 3.0f;
|
|
y[ring->vert_index(e0)] += 4.0f / 3.0f;
|
|
y[ring->vert_index(f0)] += 4.0f / 3.0f;
|
|
}
|
|
else {
|
|
y[ring->vert_index(p0)] += 1;
|
|
y[ring->vert_index(q0)] += 1;
|
|
y[ring->vert_index(e0)] += 1;
|
|
y[ring->vert_index(f0)] += 1;
|
|
}
|
|
|
|
/* sub centroidB */
|
|
if(SubdFaceRing::is_triangle(edge->pair->face)) {
|
|
y[ring->vert_index(p1)] -= 4.0f / 3.0f;
|
|
y[ring->vert_index(e0)] -= 4.0f / 3.0f;
|
|
y[ring->vert_index(f0)] -= 4.0f / 3.0f;
|
|
|
|
}
|
|
else {
|
|
y[ring->vert_index(p1)] -= 1;
|
|
y[ring->vert_index(q1)] -= 1;
|
|
y[ring->vert_index(e0)] -= 1;
|
|
y[ring->vert_index(f0)] -= 1;
|
|
}
|
|
|
|
y /= 18.0f;
|
|
|
|
if(ring->is_triangle()) {
|
|
x *= 3.0f / 4.0f;
|
|
y *= 3.0f / 4.0f;
|
|
}
|
|
|
|
/* this change makes the triangle boundaries smoother, but distorts the quads next to them */
|
|
/*if(ring->is_triangle() || SubdFaceRing::is_triangle(edge->pair->face))
|
|
{
|
|
y *= 4.0f / 3.0f;
|
|
}*/
|
|
|
|
stencil->get(interior1Indices[primitiveOffset+v]) = stencil->get(edge1Indices[primitiveOffset+v]);
|
|
stencil->get(interior1Indices[primitiveOffset+v]) += x;
|
|
stencil->get(interior1Indices[primitiveOffset+v]) += y;
|
|
|
|
for(int i = 0; i < ring->num_verts(); i++) {
|
|
x[i] =
|
|
(costerm0 * stencil->get(corner2Indices[primitiveOffset+v])[i] -
|
|
(2*costerm1 + costerm0) * stencil->get(edge2Indices[primitiveOffset+v])[i] +
|
|
2*costerm1 * stencil->get(edge1Indices[primitiveOffset+v])[i]) / 3.0f;
|
|
}
|
|
|
|
/* y = (2*( midedgeA2 - midedgeB2) + 4*(centroidA - centroidB))/18.0f; */
|
|
y = 0.0f;
|
|
|
|
/* (2*( midedgeA2 - midedgeB2) */
|
|
y[ring->vert_index(q0)] = 1;
|
|
y[ring->vert_index(q1)] = -1;
|
|
|
|
/* add centroidA */
|
|
if(ring->is_triangle()) {
|
|
y[ring->vert_index(p0)] += 4.0f / 3.0f;
|
|
y[ring->vert_index(e0)] += 4.0f / 3.0f;
|
|
y[ring->vert_index(f0)] += 4.0f / 3.0f;
|
|
}
|
|
else {
|
|
y[ring->vert_index(p0)] += 1;
|
|
y[ring->vert_index(q0)] += 1;
|
|
y[ring->vert_index(e0)] += 1;
|
|
y[ring->vert_index(f0)] += 1;
|
|
}
|
|
|
|
/* sub centroidB */
|
|
if(SubdFaceRing::is_triangle(edge->pair->face)) {
|
|
y[ring->vert_index(p1)] -= 4.0f / 3.0f;
|
|
y[ring->vert_index(e0)] -= 4.0f / 3.0f;
|
|
y[ring->vert_index(f0)] -= 4.0f / 3.0f;
|
|
|
|
}
|
|
else {
|
|
y[ring->vert_index(p1)] -= 1;
|
|
y[ring->vert_index(q1)] -= 1;
|
|
y[ring->vert_index(e0)] -= 1;
|
|
y[ring->vert_index(f0)] -= 1;
|
|
}
|
|
|
|
y /= 18.0f;
|
|
|
|
if(ring->is_triangle()) {
|
|
x *= 3.0f / 4.0f;
|
|
y *= 3.0f / 4.0f;
|
|
}
|
|
|
|
/* this change makes the triangle boundaries smoother, but distorts the quads next to them. */
|
|
/*if(ring->is_triangle() || SubdFaceRing::is_triangle(edge->pair->face))
|
|
y *= 4.0f / 3.0f;*/
|
|
|
|
stencil->get(interior2Indices[primitiveOffset+v]) = stencil->get(edge2Indices[primitiveOffset+v]);
|
|
stencil->get(interior2Indices[primitiveOffset+v]) += x;
|
|
stencil->get(interior2Indices[primitiveOffset+v]) += y;
|
|
}
|
|
}
|
|
}
|
|
|
|
void SubdAccBuilder::computeBoundaryTangentStencils(SubdFaceRing *ring, SubdVert *vert, StencilMask & r0, StencilMask & r1)
|
|
{
|
|
assert(vert->is_boundary());
|
|
assert(r0.size() == ring->num_verts());
|
|
assert(r1.size() == ring->num_verts());
|
|
|
|
SubdEdge *edge0 = vert->edge;
|
|
assert(edge0->face == NULL);
|
|
assert(edge0->to() != vert);
|
|
|
|
SubdEdge *edgek = vert->edge->prev;
|
|
assert(edgek->face == NULL);
|
|
assert(edgek->from() != vert);
|
|
|
|
int valence = vert->valence();
|
|
|
|
int k = valence - 1;
|
|
float omega = M_PI / k;
|
|
float s = sinf(omega);
|
|
float c = cosf(omega);
|
|
|
|
float factor = 1.0f / (3 * k + c);
|
|
|
|
float gamma = -4 * s * factor;
|
|
r0[ring->vert_index(vert)] = gamma;
|
|
/* r1[ring->vert_index(vert)] = 0; */
|
|
|
|
float salpha0 = -((1 + 2 * c) * sqrtf(1 + c)) * factor / sqrtf(1 - c);
|
|
float calpha0 = 1.0f / 2.0f;
|
|
|
|
r0[ring->vert_index(edge0->to())] = salpha0;
|
|
r1[ring->vert_index(edge0->to())] = calpha0;
|
|
|
|
float salphak = salpha0;
|
|
float calphak = -1.0f / 2.0f;
|
|
|
|
r0[ring->vert_index(edgek->from())] = salphak;
|
|
r1[ring->vert_index(edgek->from())] = calphak;
|
|
|
|
int j = 0;
|
|
for(SubdVert::EdgeIterator it(vert->edges()); !it.isDone(); it.advance(), j++) {
|
|
SubdEdge *edge = it.current();
|
|
|
|
if(j == k) break;
|
|
|
|
SubdVert *p = edge->to();
|
|
SubdVert *q = edge->pair->prev->from();
|
|
|
|
float alphaj = 4 * sinf(j * omega) * factor;
|
|
float betaj = (sinf(j * omega) + sinf((j + 1) * omega)) * factor;
|
|
|
|
if(j != 0)
|
|
r0[ring->vert_index(p)] += alphaj;
|
|
|
|
if(edge->pair->prev->prev->prev == edge->pair) {
|
|
r0[ring->vert_index(vert)] += betaj / 3.0f;
|
|
r0[ring->vert_index(edge->pair->from())] += betaj / 3.0f;
|
|
r0[ring->vert_index(q)] += betaj / 3.0f;
|
|
}
|
|
else
|
|
r0[ring->vert_index(q)] += betaj;
|
|
}
|
|
|
|
if(valence == 2) {
|
|
/* r0 perpendicular to r1 */
|
|
r0[ring->vert_index(vert)] = -4.0f / 3.0f;
|
|
r0[ring->vert_index(edgek->from())] = 1.0f / 2.0f;
|
|
r0[ring->vert_index(edge0->to())] = 1.0f / 2.0f;
|
|
r0[ring->vert_index(edge0->next->to())] = 1.0f / 3.0f;
|
|
}
|
|
}
|
|
|
|
/* Subd Linear Builder */
|
|
|
|
SubdLinearBuilder::SubdLinearBuilder()
|
|
{
|
|
}
|
|
|
|
SubdLinearBuilder::~SubdLinearBuilder()
|
|
{
|
|
}
|
|
|
|
Patch *SubdLinearBuilder::run(SubdFace *face)
|
|
{
|
|
Patch *patch;
|
|
float3 *hull;
|
|
|
|
if(face->num_edges() == 3) {
|
|
LinearTrianglePatch *lpatch = new LinearTrianglePatch();
|
|
hull = lpatch->hull;
|
|
patch = lpatch;
|
|
}
|
|
else if(face->num_edges() == 4) {
|
|
LinearQuadPatch *lpatch = new LinearQuadPatch();
|
|
hull = lpatch->hull;
|
|
patch = lpatch;
|
|
}
|
|
else {
|
|
assert(0); /* n-gons should have been split already */
|
|
return NULL;
|
|
}
|
|
|
|
int i = 0;
|
|
|
|
for(SubdFace::EdgeIterator it(face->edge); !it.isDone(); it.advance())
|
|
hull[i++] = it.current()->from()->co;
|
|
|
|
if(face->num_edges() == 4)
|
|
swap(hull[2], hull[3]);
|
|
|
|
return patch;
|
|
}
|
|
|
|
CCL_NAMESPACE_END
|
|
|