// David Eberly, Geometric Tools, Redmond WA 98052
// Copyright (c) 1998-2020
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt
// https://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
// Version: 4.0.2019.08.13
#pragma once
#include <Mathematics/NURBSSurface.h>
#include <Mathematics/Vector3.h>
#include <functional>
// The algorithm for representing a circle as a NURBS curve or a sphere as a
// NURBS surface is described in
// https://www.geometrictools.com/Documentation/NURBSCircleSphere.pdf
// The implementations are related to the documents as shown next.
// NURBSEighthSphereDegree4 implements Section 3.1.2 (triangular domain)
// NURBSHalfSphereDegree3 implements Section 3.2 (rectangular domain)
// NURBSFullSphereDegree3 implements Section 2.3 (rectangular domain)
// TODO: The class NURBSSurface currently assumes a rectangular domain.
// Once support is added for triangular domains, make that new class a
// base class of the sphere-representing NURBS. This will allow sharing
// of the NURBS basis functions and evaluation framework.
namespace WwiseGTE
{
template <typename Real>
class NURBSEighthSphereDegree4
{
public:
// Construction. The eigth sphere is x^2 + y^2 + z^2 = 1 for x >= 0,
// y >= 0 and z >= 0.
NURBSEighthSphereDegree4()
{
Real const sqrt2 = std::sqrt((Real)2);
Real const sqrt3 = std::sqrt((Real)3);
Real const a0 = (sqrt3 - (Real)1) / sqrt3;
Real const a1 = (sqrt3 + (Real)1) / ((Real)2 * sqrt3);
Real const a2 = (Real)1 - ((Real)5 - sqrt2) * ((Real)7 - sqrt3) / (Real)46;
Real const b0 = (Real)4 * sqrt3 * (sqrt3 - (Real)1);
Real const b1 = (Real)3 * sqrt2;
Real const b2 = (Real)4;
Real const b3 = sqrt2 * ((Real)3 + (Real)2 * sqrt2 - sqrt3) / sqrt3;
mControls[0][0] = { (Real)0, (Real)0, (Real)1 }; // P004
mControls[0][1] = { (Real)0, a0, (Real)1 }; // P013
mControls[0][2] = { (Real)0, a1, a1 }; // P022
mControls[0][3] = { (Real)0, (Real)1, a0 }; // P031
mControls[0][4] = { (Real)0, (Real)1, (Real)0 }; // P040
mControls[1][0] = { a0, (Real)0, (Real)1 }; // P103
mControls[1][1] = { a2, a2, (Real)1 }; // P112
mControls[1][2] = { a2, (Real)1, a2 }; // P121
mControls[1][3] = { a0, (Real)1, (Real)0 }; // P130
mControls[1][4] = { (Real)0, (Real)0, (Real)0 }; // unused
mControls[2][0] = { a1, (Real)0, a1 }; // P202
mControls[2][1] = { (Real)1, a2, a2 }; // P211
mControls[2][2] = { a1, a1, (Real)0 }; // P220
mControls[2][3] = { (Real)0, (Real)0, (Real)0 }; // unused
mControls[2][4] = { (Real)0, (Real)0, (Real)0 }; // unused
mControls[3][0] = { (Real)1, (Real)0, a0 }; // P301
mControls[3][1] = { (Real)1, a0, (Real)0 }; // P310
mControls[3][2] = { (Real)0, (Real)0, (Real)0 }; // unused
mControls[3][3] = { (Real)0, (Real)0, (Real)0 }; // unused
mControls[3][4] = { (Real)0, (Real)0, (Real)0 }; // unused
mControls[4][0] = { (Real)1, (Real)0, (Real)0 }; // P400
mControls[4][1] = { (Real)0, (Real)0, (Real)0 }; // unused
mControls[4][2] = { (Real)0, (Real)0, (Real)0 }; // unused
mControls[4][3] = { (Real)0, (Real)0, (Real)0 }; // unused
mControls[4][4] = { (Real)0, (Real)0, (Real)0 }; // unused
mWeights[0][0] = b0; // w004
mWeights[0][1] = b1; // w013
mWeights[0][2] = b2; // w022
mWeights[0][3] = b1; // w031
mWeights[0][4] = b0; // w040
mWeights[1][0] = b1; // w103
mWeights[1][1] = b3; // w112
mWeights[1][2] = b3; // w121
mWeights[1][3] = b1; // w130
mWeights[1][4] = (Real)0; // unused
mWeights[2][0] = b2; // w202
mWeights[2][1] = b3; // w211
mWeights[2][2] = b2; // w220
mWeights[2][3] = (Real)0; // unused
mWeights[2][4] = (Real)0; // unused
mWeights[3][0] = b1; // w301
mWeights[3][1] = b1; // w310
mWeights[3][2] = (Real)0; // unused
mWeights[3][3] = (Real)0; // unused
mWeights[3][4] = (Real)0; // unused
mWeights[4][0] = b0; // w400
mWeights[4][1] = (Real)0; // unused
mWeights[4][2] = (Real)0; // unused
mWeights[4][3] = (Real)0; // unused
mWeights[4][4] = (Real)0; // unused
}
// Evaluation of the surface. The function supports derivative
// calculation through order 2; that is, maxOrder <= 2 is required.
// If you want only the position, pass in maxOrder of 0. If you want
// the position and first-order derivatives, pass in maxOrder of 1,
// and so on. The output 'values' are ordered as: position X;
// first-order derivatives dX/du, dX/dv; second-order derivatives
// d2X/du2, d2X/dudv, d2X/dv2.
void Evaluate(Real u, Real v, unsigned int maxOrder, Vector<3, Real> values[6]) const
{
// TODO: Some of the polynomials are used in other polynomials.
// Optimize the code by eliminating the redundant computations.
Real w = (Real)1 - u - v;
Real uu = u * u, uv = u * v, uw = u * w, vv = v * v, vw = v * w, ww = w * w;
// Compute the order-0 polynomials. Only the elements to be used
// are filled in. The other terms are uninitialized but never
// used.
Real B[5][5];
B[0][0] = ww * ww;
B[0][1] = (Real)4 * vw * ww;
B[0][2] = (Real)6 * vv * ww;
B[0][3] = (Real)4 * vv * vw;
B[0][4] = vv * vv;
B[1][0] = (Real)4 * uw * ww;
B[1][1] = (Real)12 * uv * ww;
B[1][2] = (Real)12 * uv * vw;
B[1][3] = (Real)4 * uv * vv;
B[2][0] = (Real)6 * uu * ww;
B[2][1] = (Real)12 * uu * vw;
B[2][2] = (Real)6 * uu * vv;
B[3][0] = (Real)4 * uu * uw;
B[3][1] = (Real)4 * uu * uv;
B[4][0] = uu * uu;
// Compute the NURBS position.
Vector<3, Real> N{ (Real)0, (Real)0, (Real)0 };
Real D(0);
for (int j1 = 0; j1 <= 4; ++j1)
{
for (int j0 = 0; j0 <= 4 - j1; ++j0)
{
Real product = mWeights[j1][j0] * B[j1][j0];
N += product * mControls[j1][j0];
D += product;
}
}
values[0] = N / D;
if (maxOrder >= 1)
{
// Compute the order-1 polynomials. Only the elements to be
// used are filled in. The other terms are uninitialized but
// never used.
Real WmU = w - u;
Real WmTwoU = WmU - u;
Real WmThreeU = WmTwoU - u;
Real TwoWmU = w + WmU;
Real ThreeWmU = w + TwoWmU;
Real WmV = w - v;
Real WmTwoV = WmV - v;
Real WmThreeV = WmTwoV - v;
Real TwoWmV = w + WmV;
Real ThreeWmV = w + TwoWmV;
Real Dsqr = D * D;
Real Bu[5][5];
Bu[0][0] = (Real)-4 * ww * w;
Bu[0][1] = (Real)-12 * v * ww;
Bu[0][2] = (Real)-12 * vv * w;
Bu[0][3] = (Real)-4 * v * vv;
Bu[0][4] = (Real)0;
Bu[1][0] = (Real)4 * ww * WmThreeU;
Bu[1][1] = (Real)12 * vw * WmTwoU;
Bu[1][2] = (Real)12 * vv * WmU;
Bu[1][3] = (Real)4 * vv;
Bu[2][0] = (Real)12 * uw * WmU;
Bu[2][1] = (Real)12 * uv * TwoWmU;
Bu[2][2] = (Real)12 * u * vv;
Bu[3][0] = (Real)4 * uu * ThreeWmU;
Bu[3][1] = (Real)12 * uu * v;
Bu[4][0] = (Real)4 * uu * u;
Real Bv[5][5];
Bv[0][0] = (Real)-4 * ww * w;
Bv[0][1] = (Real)4 * ww * WmThreeV;
Bv[0][2] = (Real)12 * vw * WmV;
Bv[0][3] = (Real)4 * vv * ThreeWmV;
Bv[0][4] = (Real)4 * vv * v;
Bv[1][0] = (Real)-12 * u * ww;
Bv[1][1] = (Real)12 * uw * WmTwoV;
Bv[1][2] = (Real)12 * uv * TwoWmV;
Bv[1][3] = (Real)12 * u * vv;
Bv[2][0] = (Real)-12 * uu * w;
Bv[2][1] = (Real)12 * uu * WmV;
Bv[2][2] = (Real)12 * uu * v;
Bv[3][0] = (Real)-4 * uu * u;
Bv[3][1] = (Real)4 * uu * u;
Bv[4][0] = (Real)0;
Vector<3, Real> Nu{ (Real)0, (Real)0, (Real)0 };
Vector<3, Real> Nv{ (Real)0, (Real)0, (Real)0 };
Real Du(0), Dv(0);
for (int j1 = 0; j1 <= 4; ++j1)
{
for (int j0 = 0; j0 <= 4 - j1; ++j0)
{
Real product = mWeights[j1][j0] * Bu[j1][j0];
Nu += product * mControls[j1][j0];
Du += product;
product = mWeights[j1][j0] * Bv[j1][j0];
Nv += product * mControls[j1][j0];
Dv += product;
}
}
Vector<3, Real> numerDU = D * Nu - Du * N;
Vector<3, Real> numerDV = D * Nv - Dv * N;
values[1] = numerDU / Dsqr;
values[2] = numerDV / Dsqr;
if (maxOrder >= 2)
{
// Compute the order-2 polynomials. Only the elements to
// be used are filled in. The other terms are
// uninitialized but never used.
Real Dcub = Dsqr * D;
Real Buu[5][5];
Buu[0][0] = (Real)12 * ww;
Buu[0][1] = (Real)24 * vw;
Buu[0][2] = (Real)12 * vv;
Buu[0][3] = (Real)0;
Buu[0][4] = (Real)0;
Buu[1][0] = (Real)-24 * w * WmU;
Buu[1][1] = (Real)-24 * v * TwoWmU;
Buu[1][2] = (Real)-24 * vv;
Buu[1][3] = (Real)0;
Buu[2][0] = (Real)12 * (ww - (Real)4 * uw + uu);
Buu[2][1] = (Real)24 * v * WmTwoU;
Buu[2][2] = (Real)12 * vv;
Buu[3][0] = (Real)24 * u * WmU;
Buu[3][1] = (Real)24 * uv;
Buu[4][0] = (Real)12 * uu;
Real Buv[5][5];
Buv[0][0] = (Real)12 * ww;
Buv[0][1] = (Real)-12 * w * WmTwoV;
Buv[0][2] = (Real)-12 * v * TwoWmV;
Buv[0][3] = (Real)-12 * vv;
Buv[0][4] = (Real)0;
Buv[1][0] = (Real)-12 * w * WmTwoU;
Buv[1][1] = (Real)12 * (ww + (Real)2 * (uv - uw - vw));
Buv[1][2] = (Real)12 * v * ((Real)2 * WmU - v);
Buv[1][3] = (Real)12 * vv;
Buv[2][0] = (Real)-12 * u * TwoWmU;
Buv[2][1] = (Real)12 * u * ((Real)2 * WmV - u);
Buv[2][2] = (Real)24 * uv;
Buv[3][0] = (Real)-12 * uu;
Buv[3][1] = (Real)12 * uu;
Buv[4][0] = (Real)0;
Real Bvv[5][5];
Bvv[0][0] = (Real)12 * ww;
Bvv[0][1] = (Real)-24 * w * WmV;
Bvv[0][2] = (Real)12 * (ww - (Real)4 * vw + vv);
Bvv[0][3] = (Real)24 * v * WmV;
Bvv[0][4] = (Real)12 * vv;
Bvv[1][0] = (Real)24 * uw;
Bvv[1][1] = (Real)-24 * u * TwoWmV;
Bvv[1][2] = (Real)24 * u * WmTwoV;
Bvv[1][3] = (Real)24 * uv;
Bvv[2][0] = (Real)12 * uu;
Bvv[2][1] = (Real)-24 * uu;
Bvv[2][2] = (Real)12 * uu;
Bvv[3][0] = (Real)0;
Bvv[3][1] = (Real)0;
Bvv[4][0] = (Real)0;
Vector<3, Real> Nuu{ (Real)0, (Real)0, (Real)0 };
Vector<3, Real> Nuv{ (Real)0, (Real)0, (Real)0 };
Vector<3, Real> Nvv{ (Real)0, (Real)0, (Real)0 };
Real Duu(0), Duv(0), Dvv(0);
for (int j1 = 0; j1 <= 4; ++j1)
{
for (int j0 = 0; j0 <= 4 - j1; ++j0)
{
Real product = mWeights[j1][j0] * Buu[j1][j0];
Nuu += product * mControls[j1][j0];
Duu += product;
product = mWeights[j1][j0] * Buv[j1][j0];
Nuv += product * mControls[j1][j0];
Duv += product;
product = mWeights[j1][j0] * Bvv[j1][j0];
Nvv += product * mControls[j1][j0];
Dvv += product;
}
}
Vector<3, Real> termDuu = D * (D * Nuu - Duu * N);
Vector<3, Real> termDuv = D * (D * Nuv - Duv * N - Du * Nv - Dv * Nu);
Vector<3, Real> termDvv = D * (D * Nvv - Dvv * N);
values[3] = (D * termDuu - (Real)2 * Du * numerDU) / Dcub;
values[4] = (D * termDuv + (Real)2 * Du * Dv * N) / Dcub;
values[5] = (D * termDvv - (Real)2 * Dv * numerDV) / Dcub;
}
}
}
private:
// For simplicity of the implementation, 2-dimensional arrays
// of size 5-by-5 are used. Only array[r][c] is used where
// 0 <= r <= 4 and 0 <= c < 4 - r.
Vector3<Real> mControls[5][5];
Real mWeights[5][5];
};
template <typename Real>
class NURBSHalfSphereDegree3 : public NURBSSurface<3, Real>
{
public:
NURBSHalfSphereDegree3()
:
NURBSSurface<3, Real>(BasisFunctionInput<Real>(4, 3),
BasisFunctionInput<Real>(4, 3), nullptr, nullptr)
{
// weight[j][i] is mWeights[i + 4 * j], 0 <= i < 4, 0 <= j < 4
Real const oneThird = (Real)1 / (Real)3;
Real const oneNinth = (Real)1 / (Real)9;
this->mWeights[0] = (Real)1;
this->mWeights[1] = oneThird;
this->mWeights[2] = oneThird;
this->mWeights[3] = (Real)1;
this->mWeights[4] = oneThird;
this->mWeights[5] = oneNinth;
this->mWeights[6] = oneNinth;
this->mWeights[7] = oneThird;
this->mWeights[8] = oneThird;
this->mWeights[9] = oneNinth;
this->mWeights[10] = oneNinth;
this->mWeights[11] = oneThird;
this->mWeights[12] = (Real)1;
this->mWeights[13] = oneThird;
this->mWeights[14] = oneThird;
this->mWeights[15] = (Real)1;
// control[j][i] is mControls[i + 4 * j], 0 <= i < 4, 0 <= j < 4
this->mControls[0] = { (Real)0, (Real)0, (Real)1 };
this->mControls[1] = { (Real)0, (Real)0, (Real)1 };
this->mControls[2] = { (Real)0, (Real)0, (Real)1 };
this->mControls[3] = { (Real)0, (Real)0, (Real)1 };
this->mControls[4] = { (Real)2, (Real)0, (Real)1 };
this->mControls[5] = { (Real)2, (Real)4, (Real)1 };
this->mControls[6] = { (Real)-2, (Real)4, (Real)1 };
this->mControls[7] = { (Real)-2, (Real)0, (Real)1 };
this->mControls[8] = { (Real)2, (Real)0, (Real)-1 };
this->mControls[9] = { (Real)2, (Real)4, (Real)-1 };
this->mControls[10] = { (Real)-2, (Real)4, (Real)-1 };
this->mControls[11] = { (Real)-2, (Real)0, (Real)-1 };
this->mControls[12] = { (Real)0, (Real)0, (Real)-1 };
this->mControls[13] = { (Real)0, (Real)0, (Real)-1 };
this->mControls[14] = { (Real)0, (Real)0, (Real)-1 };
this->mControls[15] = { (Real)0, (Real)0, (Real)-1 };
}
};
template <typename Real>
class NURBSFullSphereDegree3 : public NURBSSurface<3, Real>
{
public:
NURBSFullSphereDegree3()
:
NURBSSurface<3, Real>(BasisFunctionInput<Real>(4, 3),
CreateBasisFunctionInputV(), nullptr, nullptr)
{
// weight[j][i] is mWeights[i + 4 * j], 0 <= i < 4, 0 <= j < 7
Real const oneThird = (Real)1 / (Real)3;
Real const oneNinth = (Real)1 / (Real)9;
this->mWeights[0] = (Real)1;
this->mWeights[4] = oneThird;
this->mWeights[8] = oneThird;
this->mWeights[12] = (Real)1;
this->mWeights[16] = oneThird;
this->mWeights[20] = oneThird;
this->mWeights[24] = (Real)1;
this->mWeights[1] = oneThird;
this->mWeights[5] = oneNinth;
this->mWeights[9] = oneNinth;
this->mWeights[13] = oneThird;
this->mWeights[17] = oneNinth;
this->mWeights[21] = oneNinth;
this->mWeights[25] = oneThird;
this->mWeights[2] = oneThird;
this->mWeights[6] = oneNinth;
this->mWeights[10] = oneNinth;
this->mWeights[14] = oneThird;
this->mWeights[18] = oneNinth;
this->mWeights[22] = oneNinth;
this->mWeights[26] = oneThird;
this->mWeights[3] = (Real)1;
this->mWeights[7] = oneThird;
this->mWeights[11] = oneThird;
this->mWeights[15] = (Real)1;
this->mWeights[19] = oneThird;
this->mWeights[23] = oneThird;
this->mWeights[27] = (Real)1;
// control[j][i] is mControls[i + 4 * j], 0 <= i < 4, 0 <= j < 7
this->mControls[0] = { (Real)0, (Real)0, (Real)1 };
this->mControls[4] = { (Real)0, (Real)0, (Real)1 };
this->mControls[8] = { (Real)0, (Real)0, (Real)1 };
this->mControls[12] = { (Real)0, (Real)0, (Real)1 };
this->mControls[16] = { (Real)0, (Real)0, (Real)1 };
this->mControls[20] = { (Real)0, (Real)0, (Real)1 };
this->mControls[24] = { (Real)0, (Real)0, (Real)1 };
this->mControls[1] = { (Real)2, (Real)0, (Real)1 };
this->mControls[5] = { (Real)2, (Real)4, (Real)1 };
this->mControls[9] = { (Real)-2, (Real)4, (Real)1 };
this->mControls[13] = { (Real)-2, (Real)0, (Real)1 };
this->mControls[17] = { (Real)-2, (Real)-4, (Real)1 };
this->mControls[21] = { (Real)2, (Real)-4, (Real)1 };
this->mControls[25] = { (Real)2, (Real)0, (Real)1 };
this->mControls[2] = { (Real)2, (Real)0, (Real)-1 };
this->mControls[6] = { (Real)2, (Real)4, (Real)-1 };
this->mControls[10] = { (Real)-2, (Real)4, (Real)-1 };
this->mControls[14] = { (Real)-2, (Real)0, (Real)-1 };
this->mControls[18] = { (Real)-2, (Real)-4, (Real)-1 };
this->mControls[22] = { (Real)2, (Real)-4, (Real)-1 };
this->mControls[26] = { (Real)2, (Real)0, (Real)-1 };
this->mControls[3] = { (Real)0, (Real)0, (Real)-1 };
this->mControls[7] = { (Real)0, (Real)0, (Real)-1 };
this->mControls[11] = { (Real)0, (Real)0, (Real)-1 };
this->mControls[15] = { (Real)0, (Real)0, (Real)-1 };
this->mControls[19] = { (Real)0, (Real)0, (Real)-1 };
this->mControls[23] = { (Real)0, (Real)0, (Real)-1 };
this->mControls[27] = { (Real)0, (Real)0, (Real)-1 };
}
private:
static BasisFunctionInput<Real> CreateBasisFunctionInputV()
{
BasisFunctionInput<Real> input;
input.numControls = 7;
input.degree = 3;
input.uniform = true;
input.periodic = false;
input.numUniqueKnots = 3;
input.uniqueKnots.resize(input.numUniqueKnots);
input.uniqueKnots[0] = { (Real)0, 4 };
input.uniqueKnots[1] = { (Real)0.5, 3 };
input.uniqueKnots[2] = { (Real)1, 4 };
return input;
}
};
}