godotengine
/
godot
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
godot/thirdparty/embree/kernels/subdiv/gregory_patch.h

894 lines
42 KiB
C++
Raw Normal View History

Upgrade Embree to the latest official release. Since Embree v3.13.0 supports AARCH64, switch back to the official repo instead of using Embree-aarch64. `thirdparty/embree/patches/godot-changes.patch` should now contain an accurate diff of the changes done to the library. (cherry picked from commit 767e374dced69b45db0afb30ca2ccf0bbbeef672)
2021-05-20 10:49:33 +00:00
// Copyright 2009-2021 Intel Corporation
Add Embree thirdparty library
2020-12-19 13:50:20 +00:00
// SPDX-License-Identifier: Apache-2.0
#pragma once
#include "catmullclark_patch.h"
#include "bezier_patch.h"
#include "bezier_curve.h"
#include "catmullclark_coefficients.h"
namespace embree
{
template<typename Vertex, typename Vertex_t = Vertex>
class __aligned(64) GregoryPatchT
{
typedef CatmullClarkPatchT<Vertex,Vertex_t> CatmullClarkPatch;
typedef GeneralCatmullClarkPatchT<Vertex,Vertex_t> GeneralCatmullClarkPatch;
typedef CatmullClark1RingT<Vertex,Vertex_t> CatmullClark1Ring;
typedef BezierCurveT<Vertex> BezierCurve;
public:
Vertex v[4][4];
Vertex f[2][2];
__forceinline GregoryPatchT() {}
__forceinline GregoryPatchT(const CatmullClarkPatch& patch) {
init(patch);
}
__forceinline GregoryPatchT(const CatmullClarkPatch& patch,
const BezierCurve* border0, const BezierCurve* border1, const BezierCurve* border2, const BezierCurve* border3)
{
init_crackfix(patch,border0,border1,border2,border3);
}
__forceinline GregoryPatchT (const HalfEdge* edge, const char* vertices, size_t stride) {
init(CatmullClarkPatch(edge,vertices,stride));
}
__forceinline Vertex& p0() { return v[0][0]; }
__forceinline Vertex& p1() { return v[0][3]; }
__forceinline Vertex& p2() { return v[3][3]; }
__forceinline Vertex& p3() { return v[3][0]; }
__forceinline Vertex& e0_p() { return v[0][1]; }
__forceinline Vertex& e0_m() { return v[1][0]; }
__forceinline Vertex& e1_p() { return v[1][3]; }
__forceinline Vertex& e1_m() { return v[0][2]; }
__forceinline Vertex& e2_p() { return v[3][2]; }
__forceinline Vertex& e2_m() { return v[2][3]; }
__forceinline Vertex& e3_p() { return v[2][0]; }
__forceinline Vertex& e3_m() { return v[3][1]; }
__forceinline Vertex& f0_p() { return v[1][1]; }
__forceinline Vertex& f1_p() { return v[1][2]; }
__forceinline Vertex& f2_p() { return v[2][2]; }
__forceinline Vertex& f3_p() { return v[2][1]; }
__forceinline Vertex& f0_m() { return f[0][0]; }
__forceinline Vertex& f1_m() { return f[0][1]; }
__forceinline Vertex& f2_m() { return f[1][1]; }
__forceinline Vertex& f3_m() { return f[1][0]; }
__forceinline const Vertex& p0() const { return v[0][0]; }
__forceinline const Vertex& p1() const { return v[0][3]; }
__forceinline const Vertex& p2() const { return v[3][3]; }
__forceinline const Vertex& p3() const { return v[3][0]; }
__forceinline const Vertex& e0_p() const { return v[0][1]; }
__forceinline const Vertex& e0_m() const { return v[1][0]; }
__forceinline const Vertex& e1_p() const { return v[1][3]; }
__forceinline const Vertex& e1_m() const { return v[0][2]; }
__forceinline const Vertex& e2_p() const { return v[3][2]; }
__forceinline const Vertex& e2_m() const { return v[2][3]; }
__forceinline const Vertex& e3_p() const { return v[2][0]; }
__forceinline const Vertex& e3_m() const { return v[3][1]; }
__forceinline const Vertex& f0_p() const { return v[1][1]; }
__forceinline const Vertex& f1_p() const { return v[1][2]; }
__forceinline const Vertex& f2_p() const { return v[2][2]; }
__forceinline const Vertex& f3_p() const { return v[2][1]; }
__forceinline const Vertex& f0_m() const { return f[0][0]; }
__forceinline const Vertex& f1_m() const { return f[0][1]; }
__forceinline const Vertex& f2_m() const { return f[1][1]; }
__forceinline const Vertex& f3_m() const { return f[1][0]; }
__forceinline Vertex initCornerVertex(const CatmullClarkPatch& irreg_patch, const size_t index) {
return irreg_patch.ring[index].getLimitVertex();
}
__forceinline Vertex initPositiveEdgeVertex(const CatmullClarkPatch& irreg_patch, const size_t index, const Vertex& p_vtx) {
return madd(1.0f/3.0f,irreg_patch.ring[index].getLimitTangent(),p_vtx);
}
__forceinline Vertex initNegativeEdgeVertex(const CatmullClarkPatch& irreg_patch, const size_t index, const Vertex& p_vtx) {
return madd(1.0f/3.0f,irreg_patch.ring[index].getSecondLimitTangent(),p_vtx);
}
__forceinline Vertex initPositiveEdgeVertex2(const CatmullClarkPatch& irreg_patch, const size_t index, const Vertex& p_vtx)
{
CatmullClark1Ring3fa r0,r1,r2;
irreg_patch.ring[index].subdivide(r0);
r0.subdivide(r1);
r1.subdivide(r2);
return madd(8.0f/3.0f,r2.getLimitTangent(),p_vtx);
}
__forceinline Vertex initNegativeEdgeVertex2(const CatmullClarkPatch& irreg_patch, const size_t index, const Vertex& p_vtx)
{
CatmullClark1Ring3fa r0,r1,r2;
irreg_patch.ring[index].subdivide(r0);
r0.subdivide(r1);
r1.subdivide(r2);
return madd(8.0f/3.0f,r2.getSecondLimitTangent(),p_vtx);
}
void initFaceVertex(const CatmullClarkPatch& irreg_patch,
const size_t index,
const Vertex& p_vtx,
const Vertex& e0_p_vtx,
const Vertex& e1_m_vtx,
const unsigned int face_valence_p1,
const Vertex& e0_m_vtx,
const Vertex& e3_p_vtx,
const unsigned int face_valence_p3,
Vertex& f_p_vtx,
Vertex& f_m_vtx)
{
const unsigned int face_valence = irreg_patch.ring[index].face_valence;
const unsigned int edge_valence = irreg_patch.ring[index].edge_valence;
const unsigned int border_index = irreg_patch.ring[index].border_index;
const Vertex& vtx = irreg_patch.ring[index].vtx;
const Vertex e_i = irreg_patch.ring[index].getEdgeCenter(0);
const Vertex c_i_m_1 = irreg_patch.ring[index].getQuadCenter(0);
const Vertex e_i_m_1 = irreg_patch.ring[index].getEdgeCenter(1);
Vertex c_i, e_i_p_1;
const bool hasHardEdge0 =
std::isinf(irreg_patch.ring[index].vertex_crease_weight) &&
std::isinf(irreg_patch.ring[index].crease_weight[0]);
if (unlikely((border_index == edge_valence-2) || hasHardEdge0))
{
/* mirror quad center and edge mid-point */
c_i = madd(2.0f, e_i - c_i_m_1, c_i_m_1);
e_i_p_1 = madd(2.0f, vtx - e_i_m_1, e_i_m_1);
}
else
{
c_i = irreg_patch.ring[index].getQuadCenter( face_valence-1 );
e_i_p_1 = irreg_patch.ring[index].getEdgeCenter( face_valence-1 );
}
Vertex c_i_m_2, e_i_m_2;
const bool hasHardEdge1 =
std::isinf(irreg_patch.ring[index].vertex_crease_weight) &&
std::isinf(irreg_patch.ring[index].crease_weight[1]);
if (unlikely(border_index == 2 || hasHardEdge1))
{
/* mirror quad center and edge mid-point */
c_i_m_2 = madd(2.0f, e_i_m_1 - c_i_m_1, c_i_m_1);
e_i_m_2 = madd(2.0f, vtx - e_i, + e_i);
}
else
{
c_i_m_2 = irreg_patch.ring[index].getQuadCenter( 1 );
e_i_m_2 = irreg_patch.ring[index].getEdgeCenter( 2 );
}
const float d = 3.0f;
//const float c = cosf(2.0f*M_PI/(float)face_valence);
//const float c_e_p = cosf(2.0f*M_PI/(float)face_valence_p1);
//const float c_e_m = cosf(2.0f*M_PI/(float)face_valence_p3);
const float c = CatmullClarkPrecomputedCoefficients::table.cos_2PI_div_n(face_valence);
const float c_e_p = CatmullClarkPrecomputedCoefficients::table.cos_2PI_div_n(face_valence_p1);
const float c_e_m = CatmullClarkPrecomputedCoefficients::table.cos_2PI_div_n(face_valence_p3);
const Vertex r_e_p = 1.0f/3.0f * (e_i_m_1 - e_i_p_1) + 2.0f/3.0f * (c_i_m_1 - c_i);
const Vertex r_e_m = 1.0f/3.0f * (e_i - e_i_m_2) + 2.0f/3.0f * (c_i_m_1 - c_i_m_2);
f_p_vtx = 1.0f / d * (c_e_p * p_vtx + (d - 2.0f*c - c_e_p) * e0_p_vtx + 2.0f*c* e1_m_vtx + r_e_p);
f_m_vtx = 1.0f / d * (c_e_m * p_vtx + (d - 2.0f*c - c_e_m) * e0_m_vtx + 2.0f*c* e3_p_vtx + r_e_m);
}
__noinline void init(const CatmullClarkPatch& patch)
{
assert( patch.ring[0].hasValidPositions() );
assert( patch.ring[1].hasValidPositions() );
assert( patch.ring[2].hasValidPositions() );
assert( patch.ring[3].hasValidPositions() );
p0() = initCornerVertex(patch,0);
p1() = initCornerVertex(patch,1);
p2() = initCornerVertex(patch,2);
p3() = initCornerVertex(patch,3);
e0_p() = initPositiveEdgeVertex(patch,0, p0());
e1_p() = initPositiveEdgeVertex(patch,1, p1());
e2_p() = initPositiveEdgeVertex(patch,2, p2());
e3_p() = initPositiveEdgeVertex(patch,3, p3());
e0_m() = initNegativeEdgeVertex(patch,0, p0());
e1_m() = initNegativeEdgeVertex(patch,1, p1());
e2_m() = initNegativeEdgeVertex(patch,2, p2());
e3_m() = initNegativeEdgeVertex(patch,3, p3());
const unsigned int face_valence_p0 = patch.ring[0].face_valence;
const unsigned int face_valence_p1 = patch.ring[1].face_valence;
const unsigned int face_valence_p2 = patch.ring[2].face_valence;
const unsigned int face_valence_p3 = patch.ring[3].face_valence;
initFaceVertex(patch,0,p0(),e0_p(),e1_m(),face_valence_p1,e0_m(),e3_p(),face_valence_p3,f0_p(),f0_m() );
initFaceVertex(patch,1,p1(),e1_p(),e2_m(),face_valence_p2,e1_m(),e0_p(),face_valence_p0,f1_p(),f1_m() );
initFaceVertex(patch,2,p2(),e2_p(),e3_m(),face_valence_p3,e2_m(),e1_p(),face_valence_p1,f2_p(),f2_m() );
initFaceVertex(patch,3,p3(),e3_p(),e0_m(),face_valence_p0,e3_m(),e2_p(),face_valence_p3,f3_p(),f3_m() );
}
__noinline void init_crackfix(const CatmullClarkPatch& patch,
const BezierCurve* border0,
const BezierCurve* border1,
const BezierCurve* border2,
const BezierCurve* border3)
{
assert( patch.ring[0].hasValidPositions() );
assert( patch.ring[1].hasValidPositions() );
assert( patch.ring[2].hasValidPositions() );
assert( patch.ring[3].hasValidPositions() );
p0() = initCornerVertex(patch,0);
p1() = initCornerVertex(patch,1);
p2() = initCornerVertex(patch,2);
p3() = initCornerVertex(patch,3);
e0_p() = initPositiveEdgeVertex(patch,0, p0());
e1_p() = initPositiveEdgeVertex(patch,1, p1());
e2_p() = initPositiveEdgeVertex(patch,2, p2());
e3_p() = initPositiveEdgeVertex(patch,3, p3());
e0_m() = initNegativeEdgeVertex(patch,0, p0());
e1_m() = initNegativeEdgeVertex(patch,1, p1());
e2_m() = initNegativeEdgeVertex(patch,2, p2());
e3_m() = initNegativeEdgeVertex(patch,3, p3());
if (unlikely(border0 != nullptr))
{
p0() = border0->v0;
e0_p() = border0->v1;
e1_m() = border0->v2;
p1() = border0->v3;
}
if (unlikely(border1 != nullptr))
{
p1() = border1->v0;
e1_p() = border1->v1;
e2_m() = border1->v2;
p2() = border1->v3;
}
if (unlikely(border2 != nullptr))
{
p2() = border2->v0;
e2_p() = border2->v1;
e3_m() = border2->v2;
p3() = border2->v3;
}
if (unlikely(border3 != nullptr))
{
p3() = border3->v0;
e3_p() = border3->v1;
e0_m() = border3->v2;
p0() = border3->v3;
}
const unsigned int face_valence_p0 = patch.ring[0].face_valence;
const unsigned int face_valence_p1 = patch.ring[1].face_valence;
const unsigned int face_valence_p2 = patch.ring[2].face_valence;
const unsigned int face_valence_p3 = patch.ring[3].face_valence;
initFaceVertex(patch,0,p0(),e0_p(),e1_m(),face_valence_p1,e0_m(),e3_p(),face_valence_p3,f0_p(),f0_m() );
initFaceVertex(patch,1,p1(),e1_p(),e2_m(),face_valence_p2,e1_m(),e0_p(),face_valence_p0,f1_p(),f1_m() );
initFaceVertex(patch,2,p2(),e2_p(),e3_m(),face_valence_p3,e2_m(),e1_p(),face_valence_p1,f2_p(),f2_m() );
initFaceVertex(patch,3,p3(),e3_p(),e0_m(),face_valence_p0,e3_m(),e2_p(),face_valence_p3,f3_p(),f3_m() );
}
void computeGregoryPatchFacePoints(const unsigned int face_valence,
const Vertex& r_e_p,
const Vertex& r_e_m,
const Vertex& p_vtx,
const Vertex& e0_p_vtx,
const Vertex& e1_m_vtx,
const unsigned int face_valence_p1,
const Vertex& e0_m_vtx,
const Vertex& e3_p_vtx,
const unsigned int face_valence_p3,
Vertex& f_p_vtx,
Vertex& f_m_vtx,
const float d = 3.0f)
{
//const float c = cosf(2.0*M_PI/(float)face_valence);
//const float c_e_p = cosf(2.0*M_PI/(float)face_valence_p1);
//const float c_e_m = cosf(2.0*M_PI/(float)face_valence_p3);
const float c = CatmullClarkPrecomputedCoefficients::table.cos_2PI_div_n(face_valence);
const float c_e_p = CatmullClarkPrecomputedCoefficients::table.cos_2PI_div_n(face_valence_p1);
const float c_e_m = CatmullClarkPrecomputedCoefficients::table.cos_2PI_div_n(face_valence_p3);
f_p_vtx = 1.0f / d * (c_e_p * p_vtx + (d - 2.0f*c - c_e_p) * e0_p_vtx + 2.0f*c* e1_m_vtx + r_e_p);
f_m_vtx = 1.0f / d * (c_e_m * p_vtx + (d - 2.0f*c - c_e_m) * e0_m_vtx + 2.0f*c* e3_p_vtx + r_e_m);
f_p_vtx = 1.0f / d * (c_e_p * p_vtx + (d - 2.0f*c - c_e_p) * e0_p_vtx + 2.0f*c* e1_m_vtx + r_e_p);
f_m_vtx = 1.0f / d * (c_e_m * p_vtx + (d - 2.0f*c - c_e_m) * e0_m_vtx + 2.0f*c* e3_p_vtx + r_e_m);
}
__noinline void init(const GeneralCatmullClarkPatch& patch)
{
assert(patch.size() == 4);
#if 0
CatmullClarkPatch qpatch; patch.init(qpatch);
init(qpatch);
#else
const float face_valence_p0 = patch.ring[0].face_valence;
const float face_valence_p1 = patch.ring[1].face_valence;
const float face_valence_p2 = patch.ring[2].face_valence;
const float face_valence_p3 = patch.ring[3].face_valence;
Vertex p0_r_p, p0_r_m;
patch.ring[0].computeGregoryPatchEdgePoints( p0(), e0_p(), e0_m(), p0_r_p, p0_r_m );
Vertex p1_r_p, p1_r_m;
patch.ring[1].computeGregoryPatchEdgePoints( p1(), e1_p(), e1_m(), p1_r_p, p1_r_m );
Vertex p2_r_p, p2_r_m;
patch.ring[2].computeGregoryPatchEdgePoints( p2(), e2_p(), e2_m(), p2_r_p, p2_r_m );
Vertex p3_r_p, p3_r_m;
patch.ring[3].computeGregoryPatchEdgePoints( p3(), e3_p(), e3_m(), p3_r_p, p3_r_m );
computeGregoryPatchFacePoints(face_valence_p0, p0_r_p, p0_r_m, p0(), e0_p(), e1_m(), face_valence_p1, e0_m(), e3_p(), face_valence_p3, f0_p(), f0_m() );
computeGregoryPatchFacePoints(face_valence_p1, p1_r_p, p1_r_m, p1(), e1_p(), e2_m(), face_valence_p2, e1_m(), e0_p(), face_valence_p0, f1_p(), f1_m() );
computeGregoryPatchFacePoints(face_valence_p2, p2_r_p, p2_r_m, p2(), e2_p(), e3_m(), face_valence_p3, e2_m(), e1_p(), face_valence_p1, f2_p(), f2_m() );
computeGregoryPatchFacePoints(face_valence_p3, p3_r_p, p3_r_m, p3(), e3_p(), e0_m(), face_valence_p0, e3_m(), e2_p(), face_valence_p3, f3_p(), f3_m() );
#endif
}
__forceinline void convert_to_bezier()
{
f0_p() = (f0_p() + f0_m()) * 0.5f;
f1_p() = (f1_p() + f1_m()) * 0.5f;
f2_p() = (f2_p() + f2_m()) * 0.5f;
f3_p() = (f3_p() + f3_m()) * 0.5f;
f0_m() = Vertex( zero );
f1_m() = Vertex( zero );
f2_m() = Vertex( zero );
f3_m() = Vertex( zero );
}
static __forceinline void computeInnerVertices(const Vertex matrix[4][4], const Vertex f_m[2][2], const float uu, const float vv,
Vertex_t& matrix_11, Vertex_t& matrix_12, Vertex_t& matrix_22, Vertex_t& matrix_21)
{
if (unlikely(uu == 0.0f || uu == 1.0f || vv == 0.0f || vv == 1.0f))
{
matrix_11 = matrix[1][1];
matrix_12 = matrix[1][2];
matrix_22 = matrix[2][2];
matrix_21 = matrix[2][1];
}
else
{
const Vertex_t f0_p = matrix[1][1];
const Vertex_t f1_p = matrix[1][2];
const Vertex_t f2_p = matrix[2][2];
const Vertex_t f3_p = matrix[2][1];
const Vertex_t f0_m = f_m[0][0];
const Vertex_t f1_m = f_m[0][1];
const Vertex_t f2_m = f_m[1][1];
const Vertex_t f3_m = f_m[1][0];
matrix_11 = ( uu * f0_p + vv * f0_m)*rcp(uu+vv);
matrix_12 = ((1.0f-uu) * f1_m + vv * f1_p)*rcp(1.0f-uu+vv);
matrix_22 = ((1.0f-uu) * f2_p + (1.0f-vv) * f2_m)*rcp(2.0f-uu-vv);
matrix_21 = ( uu * f3_m + (1.0f-vv) * f3_p)*rcp(1.0f+uu-vv);
}
}
template<typename vfloat>
static __forceinline void computeInnerVertices(const Vertex v[4][4], const Vertex f[2][2],
size_t i, const vfloat& uu, const vfloat& vv, vfloat& matrix_11, vfloat& matrix_12, vfloat& matrix_22, vfloat& matrix_21)
{
const auto m_border = (uu == 0.0f) | (uu == 1.0f) | (vv == 0.0f) | (vv == 1.0f);
const vfloat f0_p = v[1][1][i];
const vfloat f1_p = v[1][2][i];
const vfloat f2_p = v[2][2][i];
const vfloat f3_p = v[2][1][i];
const vfloat f0_m = f[0][0][i];
const vfloat f1_m = f[0][1][i];
const vfloat f2_m = f[1][1][i];
const vfloat f3_m = f[1][0][i];
const vfloat one_minus_uu = vfloat(1.0f) - uu;
const vfloat one_minus_vv = vfloat(1.0f) - vv;
const vfloat f0_i = ( uu * f0_p + vv * f0_m) * rcp(uu+vv);
const vfloat f1_i = (one_minus_uu * f1_m + vv * f1_p) * rcp(one_minus_uu+vv);
const vfloat f2_i = (one_minus_uu * f2_p + one_minus_vv * f2_m) * rcp(one_minus_uu+one_minus_vv);
const vfloat f3_i = ( uu * f3_m + one_minus_vv * f3_p) * rcp(uu+one_minus_vv);
matrix_11 = select(m_border,f0_p,f0_i);
matrix_12 = select(m_border,f1_p,f1_i);
matrix_22 = select(m_border,f2_p,f2_i);
matrix_21 = select(m_border,f3_p,f3_i);
}
static __forceinline Vertex eval(const Vertex matrix[4][4], const Vertex f[2][2], const float& uu, const float& vv)
{
Vertex_t v_11, v_12, v_22, v_21;
computeInnerVertices(matrix,f,uu,vv,v_11, v_12, v_22, v_21);
const Vec4<float> Bu = BezierBasis::eval(uu);
const Vec4<float> Bv = BezierBasis::eval(vv);
return madd(Bv.x,madd(Bu.x,matrix[0][0],madd(Bu.y,matrix[0][1],madd(Bu.z,matrix[0][2],Bu.w * matrix[0][3]))),
madd(Bv.y,madd(Bu.x,matrix[1][0],madd(Bu.y,v_11 ,madd(Bu.z,v_12 ,Bu.w * matrix[1][3]))),
madd(Bv.z,madd(Bu.x,matrix[2][0],madd(Bu.y,v_21 ,madd(Bu.z,v_22 ,Bu.w * matrix[2][3]))),
Bv.w*madd(Bu.x,matrix[3][0],madd(Bu.y,matrix[3][1],madd(Bu.z,matrix[3][2],Bu.w * matrix[3][3]))))));
}
static __forceinline Vertex eval_du(const Vertex matrix[4][4], const Vertex f[2][2], const float uu, const float vv) // approximative derivative
{
Vertex_t v_11, v_12, v_22, v_21;
computeInnerVertices(matrix,f,uu,vv,v_11, v_12, v_22, v_21);
const Vec4<float> Bu = BezierBasis::derivative(uu);
const Vec4<float> Bv = BezierBasis::eval(vv);
return madd(Bv.x,madd(Bu.x,matrix[0][0],madd(Bu.y,matrix[0][1],madd(Bu.z,matrix[0][2],Bu.w * matrix[0][3]))),
madd(Bv.y,madd(Bu.x,matrix[1][0],madd(Bu.y,v_11 ,madd(Bu.z,v_12 ,Bu.w * matrix[1][3]))),
madd(Bv.z,madd(Bu.x,matrix[2][0],madd(Bu.y,v_21 ,madd(Bu.z,v_22 ,Bu.w * matrix[2][3]))),
Bv.w*madd(Bu.x,matrix[3][0],madd(Bu.y,matrix[3][1],madd(Bu.z,matrix[3][2],Bu.w * matrix[3][3]))))));
}
static __forceinline Vertex eval_dv(const Vertex matrix[4][4], const Vertex f[2][2], const float uu, const float vv) // approximative derivative
{
Vertex_t v_11, v_12, v_22, v_21;
computeInnerVertices(matrix,f,uu,vv,v_11, v_12, v_22, v_21);
const Vec4<float> Bu = BezierBasis::eval(uu);
const Vec4<float> Bv = BezierBasis::derivative(vv);
return madd(Bv.x,madd(Bu.x,matrix[0][0],madd(Bu.y,matrix[0][1],madd(Bu.z,matrix[0][2],Bu.w * matrix[0][3]))),
madd(Bv.y,madd(Bu.x,matrix[1][0],madd(Bu.y,v_11 ,madd(Bu.z,v_12 ,Bu.w * matrix[1][3]))),
madd(Bv.z,madd(Bu.x,matrix[2][0],madd(Bu.y,v_21 ,madd(Bu.z,v_22 ,Bu.w * matrix[2][3]))),
Bv.w*madd(Bu.x,matrix[3][0],madd(Bu.y,matrix[3][1],madd(Bu.z,matrix[3][2],Bu.w * matrix[3][3]))))));
}
static __forceinline Vertex eval_dudu(const Vertex matrix[4][4], const Vertex f[2][2], const float uu, const float vv) // approximative derivative
{
Vertex_t v_11, v_12, v_22, v_21;
computeInnerVertices(matrix,f,uu,vv,v_11, v_12, v_22, v_21);
const Vec4<float> Bu = BezierBasis::derivative2(uu);
const Vec4<float> Bv = BezierBasis::eval(vv);
return madd(Bv.x,madd(Bu.x,matrix[0][0],madd(Bu.y,matrix[0][1],madd(Bu.z,matrix[0][2],Bu.w * matrix[0][3]))),
madd(Bv.y,madd(Bu.x,matrix[1][0],madd(Bu.y,v_11 ,madd(Bu.z,v_12 ,Bu.w * matrix[1][3]))),
madd(Bv.z,madd(Bu.x,matrix[2][0],madd(Bu.y,v_21 ,madd(Bu.z,v_22 ,Bu.w * matrix[2][3]))),
Bv.w*madd(Bu.x,matrix[3][0],madd(Bu.y,matrix[3][1],madd(Bu.z,matrix[3][2],Bu.w * matrix[3][3]))))));
}
static __forceinline Vertex eval_dvdv(const Vertex matrix[4][4], const Vertex f[2][2], const float uu, const float vv) // approximative derivative
{
Vertex_t v_11, v_12, v_22, v_21;
computeInnerVertices(matrix,f,uu,vv,v_11, v_12, v_22, v_21);
const Vec4<float> Bu = BezierBasis::eval(uu);
const Vec4<float> Bv = BezierBasis::derivative2(vv);
return madd(Bv.x,madd(Bu.x,matrix[0][0],madd(Bu.y,matrix[0][1],madd(Bu.z,matrix[0][2],Bu.w * matrix[0][3]))),
madd(Bv.y,madd(Bu.x,matrix[1][0],madd(Bu.y,v_11 ,madd(Bu.z,v_12 ,Bu.w * matrix[1][3]))),
madd(Bv.z,madd(Bu.x,matrix[2][0],madd(Bu.y,v_21 ,madd(Bu.z,v_22 ,Bu.w * matrix[2][3]))),
Bv.w*madd(Bu.x,matrix[3][0],madd(Bu.y,matrix[3][1],madd(Bu.z,matrix[3][2],Bu.w * matrix[3][3]))))));
}
static __forceinline Vertex eval_dudv(const Vertex matrix[4][4], const Vertex f[2][2], const float uu, const float vv) // approximative derivative
{
Vertex_t v_11, v_12, v_22, v_21;
computeInnerVertices(matrix,f,uu,vv,v_11, v_12, v_22, v_21);
const Vec4<float> Bu = BezierBasis::derivative(uu);
const Vec4<float> Bv = BezierBasis::derivative(vv);
return madd(Bv.x,madd(Bu.x,matrix[0][0],madd(Bu.y,matrix[0][1],madd(Bu.z,matrix[0][2],Bu.w * matrix[0][3]))),
madd(Bv.y,madd(Bu.x,matrix[1][0],madd(Bu.y,v_11 ,madd(Bu.z,v_12 ,Bu.w * matrix[1][3]))),
madd(Bv.z,madd(Bu.x,matrix[2][0],madd(Bu.y,v_21 ,madd(Bu.z,v_22 ,Bu.w * matrix[2][3]))),
Bv.w*madd(Bu.x,matrix[3][0],madd(Bu.y,matrix[3][1],madd(Bu.z,matrix[3][2],Bu.w * matrix[3][3]))))));
}
__forceinline Vertex eval(const float uu, const float vv) const {
return eval(v,f,uu,vv);
}
__forceinline Vertex eval_du( const float uu, const float vv) const {
return eval_du(v,f,uu,vv);
}
__forceinline Vertex eval_dv( const float uu, const float vv) const {
return eval_dv(v,f,uu,vv);
}
__forceinline Vertex eval_dudu( const float uu, const float vv) const {
return eval_dudu(v,f,uu,vv);
}
__forceinline Vertex eval_dvdv( const float uu, const float vv) const {
return eval_dvdv(v,f,uu,vv);
}
__forceinline Vertex eval_dudv( const float uu, const float vv) const {
return eval_dudv(v,f,uu,vv);
}
static __forceinline Vertex normal(const Vertex matrix[4][4], const Vertex f_m[2][2], const float uu, const float vv) // FIXME: why not using basis functions
{
/* interpolate inner vertices */
Vertex_t matrix_11, matrix_12, matrix_22, matrix_21;
computeInnerVertices(matrix,f_m,uu,vv,matrix_11, matrix_12, matrix_22, matrix_21);
/* tangentU */
const Vertex_t col0 = deCasteljau(vv, (Vertex_t)matrix[0][0], (Vertex_t)matrix[1][0], (Vertex_t)matrix[2][0], (Vertex_t)matrix[3][0]);
const Vertex_t col1 = deCasteljau(vv, (Vertex_t)matrix[0][1], (Vertex_t)matrix_11 , (Vertex_t)matrix_21 , (Vertex_t)matrix[3][1]);
const Vertex_t col2 = deCasteljau(vv, (Vertex_t)matrix[0][2], (Vertex_t)matrix_12 , (Vertex_t)matrix_22 , (Vertex_t)matrix[3][2]);
const Vertex_t col3 = deCasteljau(vv, (Vertex_t)matrix[0][3], (Vertex_t)matrix[1][3], (Vertex_t)matrix[2][3], (Vertex_t)matrix[3][3]);
const Vertex_t tangentU = deCasteljau_tangent(uu, col0, col1, col2, col3);
/* tangentV */
const Vertex_t row0 = deCasteljau(uu, (Vertex_t)matrix[0][0], (Vertex_t)matrix[0][1], (Vertex_t)matrix[0][2], (Vertex_t)matrix[0][3]);
const Vertex_t row1 = deCasteljau(uu, (Vertex_t)matrix[1][0], (Vertex_t)matrix_11 , (Vertex_t)matrix_12 , (Vertex_t)matrix[1][3]);
const Vertex_t row2 = deCasteljau(uu, (Vertex_t)matrix[2][0], (Vertex_t)matrix_21 , (Vertex_t)matrix_22 , (Vertex_t)matrix[2][3]);
const Vertex_t row3 = deCasteljau(uu, (Vertex_t)matrix[3][0], (Vertex_t)matrix[3][1], (Vertex_t)matrix[3][2], (Vertex_t)matrix[3][3]);
const Vertex_t tangentV = deCasteljau_tangent(vv, row0, row1, row2, row3);
/* normal = tangentU x tangentV */
const Vertex_t n = cross(tangentU,tangentV);
return n;
}
__forceinline Vertex normal( const float uu, const float vv) const {
return normal(v,f,uu,vv);
}
__forceinline void eval(const float u, const float v,
Vertex* P, Vertex* dPdu, Vertex* dPdv,
Vertex* ddPdudu, Vertex* ddPdvdv, Vertex* ddPdudv,
const float dscale = 1.0f) const
{
if (P) {
*P = eval(u,v);
}
if (dPdu) {
assert(dPdu); *dPdu = eval_du(u,v)*dscale;
assert(dPdv); *dPdv = eval_dv(u,v)*dscale;
}
if (ddPdudu) {
assert(ddPdudu); *ddPdudu = eval_dudu(u,v)*sqr(dscale);
assert(ddPdvdv); *ddPdvdv = eval_dvdv(u,v)*sqr(dscale);
assert(ddPdudv); *ddPdudv = eval_dudv(u,v)*sqr(dscale);
}
}
template<class vfloat>
static __forceinline vfloat eval(const Vertex v[4][4], const Vertex f[2][2],
const size_t i, const vfloat& uu, const vfloat& vv, const Vec4<vfloat>& u_n, const Vec4<vfloat>& v_n,
vfloat& matrix_11, vfloat& matrix_12, vfloat& matrix_22, vfloat& matrix_21)
{
const vfloat curve0_x = madd(v_n[0],vfloat(v[0][0][i]),madd(v_n[1],vfloat(v[1][0][i]),madd(v_n[2],vfloat(v[2][0][i]),v_n[3] * vfloat(v[3][0][i]))));
const vfloat curve1_x = madd(v_n[0],vfloat(v[0][1][i]),madd(v_n[1],vfloat(matrix_11 ),madd(v_n[2],vfloat(matrix_21 ),v_n[3] * vfloat(v[3][1][i]))));
const vfloat curve2_x = madd(v_n[0],vfloat(v[0][2][i]),madd(v_n[1],vfloat(matrix_12 ),madd(v_n[2],vfloat(matrix_22 ),v_n[3] * vfloat(v[3][2][i]))));
const vfloat curve3_x = madd(v_n[0],vfloat(v[0][3][i]),madd(v_n[1],vfloat(v[1][3][i]),madd(v_n[2],vfloat(v[2][3][i]),v_n[3] * vfloat(v[3][3][i]))));
return madd(u_n[0],curve0_x,madd(u_n[1],curve1_x,madd(u_n[2],curve2_x,u_n[3] * curve3_x)));
}
template<typename vbool, typename vfloat>
static __forceinline void eval(const Vertex v[4][4], const Vertex f[2][2],
const vbool& valid, const vfloat& uu, const vfloat& vv,
float* P, float* dPdu, float* dPdv, float* ddPdudu, float* ddPdvdv, float* ddPdudv,
const float dscale, const size_t dstride, const size_t N)
{
if (P) {
const Vec4<vfloat> u_n = BezierBasis::eval(uu);
const Vec4<vfloat> v_n = BezierBasis::eval(vv);
for (size_t i=0; i<N; i++) {
vfloat matrix_11, matrix_12, matrix_22, matrix_21;
computeInnerVertices(v,f,i,uu,vv,matrix_11,matrix_12,matrix_22,matrix_21); // FIXME: calculated multiple times
vfloat::store(valid,P+i*dstride,eval(v,f,i,uu,vv,u_n,v_n,matrix_11,matrix_12,matrix_22,matrix_21));
}
}
if (dPdu)
{
{
assert(dPdu);
const Vec4<vfloat> u_n = BezierBasis::derivative(uu);
const Vec4<vfloat> v_n = BezierBasis::eval(vv);
for (size_t i=0; i<N; i++) {
vfloat matrix_11, matrix_12, matrix_22, matrix_21;
computeInnerVertices(v,f,i,uu,vv,matrix_11,matrix_12,matrix_22,matrix_21); // FIXME: calculated multiple times
vfloat::store(valid,dPdu+i*dstride,eval(v,f,i,uu,vv,u_n,v_n,matrix_11,matrix_12,matrix_22,matrix_21)*dscale);
}
}
{
assert(dPdv);
const Vec4<vfloat> u_n = BezierBasis::eval(uu);
const Vec4<vfloat> v_n = BezierBasis::derivative(vv);
for (size_t i=0; i<N; i++) {
vfloat matrix_11, matrix_12, matrix_22, matrix_21;
computeInnerVertices(v,f,i,uu,vv,matrix_11,matrix_12,matrix_22,matrix_21); // FIXME: calculated multiple times
vfloat::store(valid,dPdv+i*dstride,eval(v,f,i,uu,vv,u_n,v_n,matrix_11,matrix_12,matrix_22,matrix_21)*dscale);
}
}
}
if (ddPdudu)
{
{
assert(ddPdudu);
const Vec4<vfloat> u_n = BezierBasis::derivative2(uu);
const Vec4<vfloat> v_n = BezierBasis::eval(vv);
for (size_t i=0; i<N; i++) {
vfloat matrix_11, matrix_12, matrix_22, matrix_21;
computeInnerVertices(v,f,i,uu,vv,matrix_11,matrix_12,matrix_22,matrix_21); // FIXME: calculated multiple times
vfloat::store(valid,ddPdudu+i*dstride,eval(v,f,i,uu,vv,u_n,v_n,matrix_11,matrix_12,matrix_22,matrix_21)*sqr(dscale));
}
}
{
assert(ddPdvdv);
const Vec4<vfloat> u_n = BezierBasis::eval(uu);
const Vec4<vfloat> v_n = BezierBasis::derivative2(vv);
for (size_t i=0; i<N; i++) {
vfloat matrix_11, matrix_12, matrix_22, matrix_21;
computeInnerVertices(v,f,i,uu,vv,matrix_11,matrix_12,matrix_22,matrix_21); // FIXME: calculated multiple times
vfloat::store(valid,ddPdvdv+i*dstride,eval(v,f,i,uu,vv,u_n,v_n,matrix_11,matrix_12,matrix_22,matrix_21)*sqr(dscale));
}
}
{
assert(ddPdudv);
const Vec4<vfloat> u_n = BezierBasis::derivative(uu);
const Vec4<vfloat> v_n = BezierBasis::derivative(vv);
for (size_t i=0; i<N; i++) {
vfloat matrix_11, matrix_12, matrix_22, matrix_21;
computeInnerVertices(v,f,i,uu,vv,matrix_11,matrix_12,matrix_22,matrix_21); // FIXME: calculated multiple times
vfloat::store(valid,ddPdudv+i*dstride,eval(v,f,i,uu,vv,u_n,v_n,matrix_11,matrix_12,matrix_22,matrix_21)*sqr(dscale));
}
}
}
}
template<typename vbool, typename vfloat>
__forceinline void eval(const vbool& valid, const vfloat& uu, const vfloat& vv,
float* P, float* dPdu, float* dPdv, float* ddPdudu, float* ddPdvdv, float* ddPdudv,
const float dscale, const size_t dstride, const size_t N) const {
eval(v,f,valid,uu,vv,P,dPdu,dPdv,ddPdudu,ddPdvdv,ddPdudv,dscale,dstride,N);
}
template<class T>
static __forceinline Vec3<T> eval_t(const Vertex matrix[4][4], const Vec3<T> f[2][2], const T& uu, const T& vv)
{
typedef typename T::Bool M;
const M m_border = (uu == 0.0f) | (uu == 1.0f) | (vv == 0.0f) | (vv == 1.0f);
const Vec3<T> f0_p = Vec3<T>(matrix[1][1].x,matrix[1][1].y,matrix[1][1].z);
const Vec3<T> f1_p = Vec3<T>(matrix[1][2].x,matrix[1][2].y,matrix[1][2].z);
const Vec3<T> f2_p = Vec3<T>(matrix[2][2].x,matrix[2][2].y,matrix[2][2].z);
const Vec3<T> f3_p = Vec3<T>(matrix[2][1].x,matrix[2][1].y,matrix[2][1].z);
const Vec3<T> f0_m = f[0][0];
const Vec3<T> f1_m = f[0][1];
const Vec3<T> f2_m = f[1][1];
const Vec3<T> f3_m = f[1][0];
const T one_minus_uu = T(1.0f) - uu;
const T one_minus_vv = T(1.0f) - vv;
const Vec3<T> f0_i = ( uu * f0_p + vv * f0_m) * rcp(uu+vv);
const Vec3<T> f1_i = (one_minus_uu * f1_m + vv * f1_p) * rcp(one_minus_uu+vv);
const Vec3<T> f2_i = (one_minus_uu * f2_p + one_minus_vv * f2_m) * rcp(one_minus_uu+one_minus_vv);
const Vec3<T> f3_i = ( uu * f3_m + one_minus_vv * f3_p) * rcp(uu+one_minus_vv);
const Vec3<T> F0( select(m_border,f0_p.x,f0_i.x), select(m_border,f0_p.y,f0_i.y), select(m_border,f0_p.z,f0_i.z) );
const Vec3<T> F1( select(m_border,f1_p.x,f1_i.x), select(m_border,f1_p.y,f1_i.y), select(m_border,f1_p.z,f1_i.z) );
const Vec3<T> F2( select(m_border,f2_p.x,f2_i.x), select(m_border,f2_p.y,f2_i.y), select(m_border,f2_p.z,f2_i.z) );
const Vec3<T> F3( select(m_border,f3_p.x,f3_i.x), select(m_border,f3_p.y,f3_i.y), select(m_border,f3_p.z,f3_i.z) );
const T B0_u = one_minus_uu * one_minus_uu * one_minus_uu;
const T B0_v = one_minus_vv * one_minus_vv * one_minus_vv;
const T B1_u = 3.0f * (one_minus_uu * uu * one_minus_uu);
const T B1_v = 3.0f * (one_minus_vv * vv * one_minus_vv);
const T B2_u = 3.0f * (uu * one_minus_uu * uu);
const T B2_v = 3.0f * (vv * one_minus_vv * vv);
const T B3_u = uu * uu * uu;
const T B3_v = vv * vv * vv;
const T x = madd(B0_v,madd(B0_u,matrix[0][0].x,madd(B1_u,matrix[0][1].x,madd(B2_u,matrix[0][2].x,B3_u * matrix[0][3].x))),
madd(B1_v,madd(B0_u,matrix[1][0].x,madd(B1_u,F0.x ,madd(B2_u,F1.x ,B3_u * matrix[1][3].x))),
madd(B2_v,madd(B0_u,matrix[2][0].x,madd(B1_u,F3.x ,madd(B2_u,F2.x ,B3_u * matrix[2][3].x))),
B3_v*madd(B0_u,matrix[3][0].x,madd(B1_u,matrix[3][1].x,madd(B2_u,matrix[3][2].x,B3_u * matrix[3][3].x))))));
const T y = madd(B0_v,madd(B0_u,matrix[0][0].y,madd(B1_u,matrix[0][1].y,madd(B2_u,matrix[0][2].y,B3_u * matrix[0][3].y))),
madd(B1_v,madd(B0_u,matrix[1][0].y,madd(B1_u,F0.y ,madd(B2_u,F1.y ,B3_u * matrix[1][3].y))),
madd(B2_v,madd(B0_u,matrix[2][0].y,madd(B1_u,F3.y ,madd(B2_u,F2.y ,B3_u * matrix[2][3].y))),
B3_v*madd(B0_u,matrix[3][0].y,madd(B1_u,matrix[3][1].y,madd(B2_u,matrix[3][2].y,B3_u * matrix[3][3].y))))));
const T z = madd(B0_v,madd(B0_u,matrix[0][0].z,madd(B1_u,matrix[0][1].z,madd(B2_u,matrix[0][2].z,B3_u * matrix[0][3].z))),
madd(B1_v,madd(B0_u,matrix[1][0].z,madd(B1_u,F0.z ,madd(B2_u,F1.z ,B3_u * matrix[1][3].z))),
madd(B2_v,madd(B0_u,matrix[2][0].z,madd(B1_u,F3.z ,madd(B2_u,F2.z ,B3_u * matrix[2][3].z))),
B3_v*madd(B0_u,matrix[3][0].z,madd(B1_u,matrix[3][1].z,madd(B2_u,matrix[3][2].z,B3_u * matrix[3][3].z))))));
return Vec3<T>(x,y,z);
}
template<class T>
__forceinline Vec3<T> eval(const T& uu, const T& vv) const
{
Vec3<T> ff[2][2];
ff[0][0] = Vec3<T>(f[0][0]);
ff[0][1] = Vec3<T>(f[0][1]);
ff[1][1] = Vec3<T>(f[1][1]);
ff[1][0] = Vec3<T>(f[1][0]);
return eval_t(v,ff,uu,vv);
}
template<class T>
static __forceinline Vec3<T> normal_t(const Vertex matrix[4][4], const Vec3<T> f[2][2], const T& uu, const T& vv)
{
typedef typename T::Bool M;
const Vec3<T> f0_p = Vec3<T>(matrix[1][1].x,matrix[1][1].y,matrix[1][1].z);
const Vec3<T> f1_p = Vec3<T>(matrix[1][2].x,matrix[1][2].y,matrix[1][2].z);
const Vec3<T> f2_p = Vec3<T>(matrix[2][2].x,matrix[2][2].y,matrix[2][2].z);
const Vec3<T> f3_p = Vec3<T>(matrix[2][1].x,matrix[2][1].y,matrix[2][1].z);
const Vec3<T> f0_m = f[0][0];
const Vec3<T> f1_m = f[0][1];
const Vec3<T> f2_m = f[1][1];
const Vec3<T> f3_m = f[1][0];
const T one_minus_uu = T(1.0f) - uu;
const T one_minus_vv = T(1.0f) - vv;
const Vec3<T> f0_i = ( uu * f0_p + vv * f0_m) * rcp(uu+vv);
const Vec3<T> f1_i = (one_minus_uu * f1_m + vv * f1_p) * rcp(one_minus_uu+vv);
const Vec3<T> f2_i = (one_minus_uu * f2_p + one_minus_vv * f2_m) * rcp(one_minus_uu+one_minus_vv);
const Vec3<T> f3_i = ( uu * f3_m + one_minus_vv * f3_p) * rcp(uu+one_minus_vv);
#if 1
const M m_corner0 = (uu == 0.0f) & (vv == 0.0f);
const M m_corner1 = (uu == 1.0f) & (vv == 0.0f);
const M m_corner2 = (uu == 1.0f) & (vv == 1.0f);
const M m_corner3 = (uu == 0.0f) & (vv == 1.0f);
const Vec3<T> matrix_11( select(m_corner0,f0_p.x,f0_i.x), select(m_corner0,f0_p.y,f0_i.y), select(m_corner0,f0_p.z,f0_i.z) );
const Vec3<T> matrix_12( select(m_corner1,f1_p.x,f1_i.x), select(m_corner1,f1_p.y,f1_i.y), select(m_corner1,f1_p.z,f1_i.z) );
const Vec3<T> matrix_22( select(m_corner2,f2_p.x,f2_i.x), select(m_corner2,f2_p.y,f2_i.y), select(m_corner2,f2_p.z,f2_i.z) );
const Vec3<T> matrix_21( select(m_corner3,f3_p.x,f3_i.x), select(m_corner3,f3_p.y,f3_i.y), select(m_corner3,f3_p.z,f3_i.z) );
#else
const M m_border = (uu == 0.0f) | (uu == 1.0f) | (vv == 0.0f) | (vv == 1.0f);
const Vec3<T> matrix_11( select(m_border,f0_p.x,f0_i.x), select(m_border,f0_p.y,f0_i.y), select(m_border,f0_p.z,f0_i.z) );
const Vec3<T> matrix_12( select(m_border,f1_p.x,f1_i.x), select(m_border,f1_p.y,f1_i.y), select(m_border,f1_p.z,f1_i.z) );
const Vec3<T> matrix_22( select(m_border,f2_p.x,f2_i.x), select(m_border,f2_p.y,f2_i.y), select(m_border,f2_p.z,f2_i.z) );
const Vec3<T> matrix_21( select(m_border,f3_p.x,f3_i.x), select(m_border,f3_p.y,f3_i.y), select(m_border,f3_p.z,f3_i.z) );
#endif
const Vec3<T> matrix_00 = Vec3<T>(matrix[0][0].x,matrix[0][0].y,matrix[0][0].z);
const Vec3<T> matrix_10 = Vec3<T>(matrix[1][0].x,matrix[1][0].y,matrix[1][0].z);
const Vec3<T> matrix_20 = Vec3<T>(matrix[2][0].x,matrix[2][0].y,matrix[2][0].z);
const Vec3<T> matrix_30 = Vec3<T>(matrix[3][0].x,matrix[3][0].y,matrix[3][0].z);
const Vec3<T> matrix_01 = Vec3<T>(matrix[0][1].x,matrix[0][1].y,matrix[0][1].z);
const Vec3<T> matrix_02 = Vec3<T>(matrix[0][2].x,matrix[0][2].y,matrix[0][2].z);
const Vec3<T> matrix_03 = Vec3<T>(matrix[0][3].x,matrix[0][3].y,matrix[0][3].z);
const Vec3<T> matrix_31 = Vec3<T>(matrix[3][1].x,matrix[3][1].y,matrix[3][1].z);
const Vec3<T> matrix_32 = Vec3<T>(matrix[3][2].x,matrix[3][2].y,matrix[3][2].z);
const Vec3<T> matrix_33 = Vec3<T>(matrix[3][3].x,matrix[3][3].y,matrix[3][3].z);
const Vec3<T> matrix_13 = Vec3<T>(matrix[1][3].x,matrix[1][3].y,matrix[1][3].z);
const Vec3<T> matrix_23 = Vec3<T>(matrix[2][3].x,matrix[2][3].y,matrix[2][3].z);
/* tangentU */
const Vec3<T> col0 = deCasteljau(vv, matrix_00, matrix_10, matrix_20, matrix_30);
const Vec3<T> col1 = deCasteljau(vv, matrix_01, matrix_11, matrix_21, matrix_31);
const Vec3<T> col2 = deCasteljau(vv, matrix_02, matrix_12, matrix_22, matrix_32);
const Vec3<T> col3 = deCasteljau(vv, matrix_03, matrix_13, matrix_23, matrix_33);
const Vec3<T> tangentU = deCasteljau_tangent(uu, col0, col1, col2, col3);
/* tangentV */
const Vec3<T> row0 = deCasteljau(uu, matrix_00, matrix_01, matrix_02, matrix_03);
const Vec3<T> row1 = deCasteljau(uu, matrix_10, matrix_11, matrix_12, matrix_13);
const Vec3<T> row2 = deCasteljau(uu, matrix_20, matrix_21, matrix_22, matrix_23);
const Vec3<T> row3 = deCasteljau(uu, matrix_30, matrix_31, matrix_32, matrix_33);
const Vec3<T> tangentV = deCasteljau_tangent(vv, row0, row1, row2, row3);
/* normal = tangentU x tangentV */
const Vec3<T> n = cross(tangentU,tangentV);
return n;
}
template<class T>
__forceinline Vec3<T> normal(const T& uu, const T& vv) const
{
Vec3<T> ff[2][2];
ff[0][0] = Vec3<T>(f[0][0]);
ff[0][1] = Vec3<T>(f[0][1]);
ff[1][1] = Vec3<T>(f[1][1]);
ff[1][0] = Vec3<T>(f[1][0]);
return normal_t(v,ff,uu,vv);
}
__forceinline BBox<Vertex> bounds() const
{
const Vertex *const cv = &v[0][0];
BBox<Vertex> bounds (cv[0]);
for (size_t i=1; i<16; i++)
bounds.extend( cv[i] );
bounds.extend(f[0][0]);
bounds.extend(f[1][0]);
bounds.extend(f[1][1]);
bounds.extend(f[1][1]);
return bounds;
}
friend embree_ostream operator<<(embree_ostream o, const GregoryPatchT& p)
{
for (size_t y=0; y<4; y++)
for (size_t x=0; x<4; x++)
o << "v[" << y << "][" << x << "] " << p.v[y][x] << embree_endl;
for (size_t y=0; y<2; y++)
for (size_t x=0; x<2; x++)
o << "f[" << y << "][" << x << "] " << p.f[y][x] << embree_endl;
return o;
}
};
typedef GregoryPatchT<Vec3fa,Vec3fa_t> GregoryPatch3fa;
template<typename Vertex, typename Vertex_t>
__forceinline BezierPatchT<Vertex,Vertex_t>::BezierPatchT (const HalfEdge* edge, const char* vertices, size_t stride)
{
CatmullClarkPatchT<Vertex,Vertex_t> patch(edge,vertices,stride);
GregoryPatchT<Vertex,Vertex_t> gpatch(patch);
gpatch.convert_to_bezier();
for (size_t y=0; y<4; y++)
for (size_t x=0; x<4; x++)
matrix[y][x] = (Vertex_t)gpatch.v[y][x];
}
template<typename Vertex, typename Vertex_t>
__forceinline BezierPatchT<Vertex,Vertex_t>::BezierPatchT(const CatmullClarkPatchT<Vertex,Vertex_t>& patch)
{
GregoryPatchT<Vertex,Vertex_t> gpatch(patch);
gpatch.convert_to_bezier();
for (size_t y=0; y<4; y++)
for (size_t x=0; x<4; x++)
matrix[y][x] = (Vertex_t)gpatch.v[y][x];
}
template<typename Vertex, typename Vertex_t>
__forceinline BezierPatchT<Vertex,Vertex_t>::BezierPatchT(const CatmullClarkPatchT<Vertex,Vertex_t>& patch,
const BezierCurveT<Vertex>* border0,
const BezierCurveT<Vertex>* border1,
const BezierCurveT<Vertex>* border2,
const BezierCurveT<Vertex>* border3)
{
GregoryPatchT<Vertex,Vertex_t> gpatch(patch,border0,border1,border2,border3);
gpatch.convert_to_bezier();
for (size_t y=0; y<4; y++)
for (size_t x=0; x<4; x++)
matrix[y][x] = (Vertex_t)gpatch.v[y][x];
}
}