Revision control
++ // Copyright 2009-2021 Intel Corporation
++ // 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];
++ }
++ }