Merge pull request #960 from davidgyu/mtl_patch_normals

Metal patch shader changes for degenerate normals
This commit is contained in:
barfowl 2018-01-16 15:31:11 -08:00 committed by GitHub
commit dcd82f4137
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
2 changed files with 191 additions and 201 deletions

View File

@ -57,6 +57,9 @@ static void initShapes() {
g_defaultShapes.push_back( ShapeDesc("catmark_chaikin0", catmark_chaikin0, kCatmark ) );
g_defaultShapes.push_back( ShapeDesc("catmark_chaikin1", catmark_chaikin1, kCatmark ) );
g_defaultShapes.push_back( ShapeDesc("catmark_chaikin2", catmark_chaikin2, kCatmark ) );
g_defaultShapes.push_back( ShapeDesc("catmark_single_crease", catmark_single_crease, kCatmark ) );
g_defaultShapes.push_back( ShapeDesc("catmark_inf_crease0", catmark_inf_crease0, kCatmark ) );
g_defaultShapes.push_back( ShapeDesc("catmark_inf_crease1", catmark_inf_crease1, kCatmark ) );
g_defaultShapes.push_back( ShapeDesc("catmark_fan", catmark_fan, kCatmark ) );
g_defaultShapes.push_back( ShapeDesc("catmark_flap", catmark_flap, kCatmark ) );
g_defaultShapes.push_back( ShapeDesc("catmark_flap2", catmark_flap2, kCatmark ) );
@ -91,7 +94,6 @@ static void initShapes() {
g_defaultShapes.push_back( ShapeDesc("catmark_tent", catmark_tent, kCatmark ) );
g_defaultShapes.push_back( ShapeDesc("catmark_torus", catmark_torus, kCatmark ) );
g_defaultShapes.push_back( ShapeDesc("catmark_torus_creases0", catmark_torus_creases0, kCatmark ) );
g_defaultShapes.push_back( ShapeDesc("catmark_single_crease", catmark_single_crease, kCatmark ) );
g_defaultShapes.push_back( ShapeDesc("catmark_smoothtris0", catmark_smoothtris0, kCatmark ) );
g_defaultShapes.push_back( ShapeDesc("catmark_smoothtris1", catmark_smoothtris1, kCatmark ) );
// g_defaultShapes.push_back( ShapeDesc("catmark_square_hedit0", catmark_square_hedit0, kCatmark ) );
@ -120,3 +122,4 @@ static void initShapes() {
g_defaultShapes.push_back( ShapeDesc("loop_pole64", loop_pole64, kLoop ) );
g_defaultShapes.push_back( ShapeDesc("loop_pole360", loop_pole360, kLoop ) );
}
//------------------------------------------------------------------------------

View File

@ -590,6 +590,14 @@ OsdComputeBSplineBoundaryPoints(thread VertexType* cpt, int3 patchParam)
}
}
template<typename PerPatchVertexBezier>
void
OsdEvalPatchBezier(int3 patchParam, float2 UV,
PerPatchVertexBezier cv,
thread float3& P, thread float3& dPu, thread float3& dPv,
thread float3& N, thread float3& dNu, thread float3& dNv,
thread float2& vSegments);
void
OsdEvalPatchGregory(int3 patchParam, float2 UV, thread float3* cv,
thread float3& P, thread float3& dPu, thread float3& dPv,
@ -624,113 +632,28 @@ OsdEvalPatchGregory(int3 patchParam, float2 UV, thread float3* cv,
float d21 = u+V;
float d22 = U+V;
float3 q[16];
OsdPerPatchVertexBezier bezcv[16];
float2 vSegments;
q[ 5] = (d11 == 0.0) ? cv[3] : (u*cv[3] + v*cv[4])/d11;
q[ 6] = (d12 == 0.0) ? cv[8] : (U*cv[9] + v*cv[8])/d12;
q[ 9] = (d21 == 0.0) ? cv[18] : (u*cv[19] + V*cv[18])/d21;
q[10] = (d22 == 0.0) ? cv[13] : (U*cv[13] + V*cv[14])/d22;
bezcv[ 5].P = (d11 == 0.0) ? cv[3] : (u*cv[3] + v*cv[4])/d11;
bezcv[ 6].P = (d12 == 0.0) ? cv[8] : (U*cv[9] + v*cv[8])/d12;
bezcv[ 9].P = (d21 == 0.0) ? cv[18] : (u*cv[19] + V*cv[18])/d21;
bezcv[10].P = (d22 == 0.0) ? cv[13] : (U*cv[13] + V*cv[14])/d22;
q[ 0] = cv[0];
q[ 1] = cv[1];
q[ 2] = cv[7];
q[ 3] = cv[5];
q[ 4] = cv[2];
q[ 7] = cv[6];
q[ 8] = cv[16];
q[11] = cv[12];
q[12] = cv[15];
q[13] = cv[17];
q[14] = cv[11];
q[15] = cv[10];
bezcv[ 0].P = cv[0];
bezcv[ 1].P = cv[1];
bezcv[ 2].P = cv[7];
bezcv[ 3].P = cv[5];
bezcv[ 4].P = cv[2];
bezcv[ 7].P = cv[6];
bezcv[ 8].P = cv[16];
bezcv[11].P = cv[12];
bezcv[12].P = cv[15];
bezcv[13].P = cv[17];
bezcv[14].P = cv[11];
bezcv[15].P = cv[10];
P = float3(0,0,0);
dPu = float3(0,0,0);
dPv = float3(0,0,0);
#if OSD_COMPUTE_NORMAL_DERIVATIVES
float B[4], D[4], C[4];
float3 BUCP[4] = {float3(0,0,0),float3(0,0,0),float3(0,0,0),float3(0,0,0)},
DUCP[4] = {float3(0,0,0),float3(0,0,0),float3(0,0,0),float3(0,0,0)},
CUCP[4] = {float3(0,0,0),float3(0,0,0),float3(0,0,0),float3(0,0,0)};
float3 dUU = float3(0,0,0);
float3 dVV = float3(0,0,0);
float3 dUV = float3(0,0,0);
OsdUnivar4x4(UV.x, B, D, C);
for (int i=0; i<4; ++i) {
for (int j=0; j<4; ++j) {
float3 A = q[4*i + j];
BUCP[i] += A * B[j];
DUCP[i] += A * D[j];
CUCP[i] += A * C[j];
}
}
OsdUnivar4x4(UV.y, B, D, C);
for (int i=0; i<4; ++i) {
P += B[i] * BUCP[i];
dPu += B[i] * DUCP[i];
dPv += D[i] * BUCP[i];
dUU += B[i] * CUCP[i];
dVV += C[i] * BUCP[i];
dUV += D[i] * DUCP[i];
}
int level = OsdGetPatchFaceLevel(patchParam);
dPu *= 3 * level;
dPv *= 3 * level;
dUU *= 6 * level;
dVV *= 6 * level;
dUV *= 9 * level;
float3 n = cross(dPu, dPv);
N = normalize(n);
float E = dot(dPu, dPu);
float F = dot(dPu, dPv);
float G = dot(dPv, dPv);
float e = dot(N, dUU);
float f = dot(N, dUV);
float g = dot(N, dVV);
dNu = (f*F-e*G)/(E*G-F*F) * dPu + (e*F-f*E)/(E*G-F*F) * dPv;
dNv = (g*F-f*G)/(E*G-F*F) * dPu + (f*F-g*E)/(E*G-F*F) * dPv;
dNu = dNu/length(n) - n * (dot(dNu,n)/powr(dot(n,n), 1.5));
dNv = dNv/length(n) - n * (dot(dNv,n)/powr(dot(n,n), 1.5));
#else //OSD_COMPUTE_NORMAL_DERIVATIVES
float B[4], D[4];
float3 BUCP[4] = {float3(0,0,0),float3(0,0,0),float3(0,0,0),float3(0,0,0)},
DUCP[4] = {float3(0,0,0),float3(0,0,0),float3(0,0,0),float3(0,0,0)};
OsdUnivar4x4(UV.x, B, D);
for (int i=0; i<4; ++i) {
for (int j=0; j<4; ++j) {
float3 A = q[4*i + j];
BUCP[i] += A * B[j];
DUCP[i] += A * D[j];
}
}
OsdUnivar4x4(UV.y, B, D);
for (int i=0; i<4; ++i) {
P += B[i] * BUCP[i];
dPu += B[i] * DUCP[i];
dPv += D[i] * BUCP[i];
}
int level = OsdGetPatchFaceLevel(patchParam);
dPu *= 3 * level;
dPv *= 3 * level;
N = normalize(cross(dPu, dPv));
dNu = float3(0,0,0);
dNv = float3(0,0,0);
#endif //OSD_COMPUTE_NORMAL_DERIVATIVES
OsdEvalPatchBezier(patchParam, UV, bezcv, P, dPu, dPv, N, dNu, dNv, vSegments);
}
// ----------------------------------------------------------------------------
@ -1559,20 +1482,30 @@ OsdEvalPatchBezier(int3 patchParam, float2 UV,
thread float3& N, thread float3& dNu, thread float3& dNv,
thread float2& vSegments)
{
#if OSD_COMPUTE_NORMAL_DERIVATIVES
float B[4], D[4], C[4];
float3 BUCP[4] = {float3(0,0,0),float3(0,0,0),float3(0,0,0),float3(0,0,0)},
DUCP[4] = {float3(0,0,0),float3(0,0,0),float3(0,0,0),float3(0,0,0)},
CUCP[4] = {float3(0,0,0),float3(0,0,0),float3(0,0,0),float3(0,0,0)};
OsdUnivar4x4(UV.x, B, D, C);
#else
float B[4], D[4];
float3 BUCP[4] = {float3(0,0,0),float3(0,0,0),float3(0,0,0),float3(0,0,0)},
DUCP[4] = {float3(0,0,0),float3(0,0,0),float3(0,0,0),float3(0,0,0)};
OsdUnivar4x4(UV.x, B, D);
#endif
//
// Use the recursive nature of the basis functions to compute a 2x2 set
// of intermediate points (via repeated linear interpolation). These
// points define a bilinear surface tangent to the desired surface at P
// and so containing dPu and dPv. The cost of computing P, dPu and dPv
// this way is comparable to that of typical tensor product evaluation
// (if not faster).
//
// If N = dPu X dPv degenerates, it often results from an edge of the
// 2x2 bilinear hull collapsing or two adjacent edges colinear. In both
// cases, the expected non-planar quad degenerates into a triangle, and
// the tangent plane of that triangle provides the desired normal N.
//
// ----------------------------------------------------------------
// Reduce 4x4 points to 2x4 -- two levels of linear interpolation in U
// and so 3 original rows contributing to each of the 2 resulting rows:
float u = UV.x;
float uinv = 1.0f - u;
float u0 = uinv * uinv;
float u1 = u * uinv * 2.0f;
float u2 = u * u;
float3 LROW[4], RROW[4];
#if OSD_PATCH_ENABLE_SINGLE_CREASE
#if USE_PTVS_SHARPNESS
float sharpness = OsdGetPatchSharpness(patchParam);
@ -1580,112 +1513,166 @@ OsdEvalPatchBezier(int3 patchParam, float2 UV,
float Sc = ceil(sharpness);
float s0 = 1 - exp2(-Sf);
float s1 = 1 - exp2(-Sc);
vSegments = float2(s0, s1);
#else //USE_PTVS_SHARPNESS
#else // USE_PTVS_SHARPNESS
vSegments = cv[0].vSegments;
#endif //USE_PTVS_SHARPNESS
#endif // USE_PTVS_SHARPNESS
float s = OsdGetPatchSingleCreaseSegmentParameter(patchParam, UV);
for (int i=0; i<4; ++i) {
for (int j=0; j<4; ++j) {
int k = 4*i + j;
float3 A = (s <= vSegments.x) ? cv[k].P
: ((s <= vSegments.y) ? cv[k].P1
: cv[k].P2);
BUCP[i] += A * B[j];
DUCP[i] += A * D[j];
#if OSD_COMPUTE_NORMAL_DERIVATIVES
CUCP[i] += A * C[j];
#endif //OSD_COMPUTE_NORMAL_DERIVATIVES
for (int i = 0; i < 4; ++i) {
int j = i*4;
if (s <= vSegments.x) {
LROW[i] = u0 * cv[ j ].P + u1 * cv[j+1].P + u2 * cv[j+2].P;
RROW[i] = u0 * cv[j+1].P + u1 * cv[j+2].P + u2 * cv[j+3].P;
} else if (s <= vSegments.y) {
LROW[i] = u0 * cv[ j ].P1 + u1 * cv[j+1].P1 + u2 * cv[j+2].P1;
RROW[i] = u0 * cv[j+1].P1 + u1 * cv[j+2].P1 + u2 * cv[j+3].P1;
} else {
LROW[i] = u0 * cv[ j ].P2 + u1 * cv[j+1].P2 + u2 * cv[j+2].P2;
RROW[i] = u0 * cv[j+1].P2 + u1 * cv[j+2].P2 + u2 * cv[j+3].P2;
}
}
#else //OSD_PATCH_ENABLE_SINGLE_CREASE
// ----------------------------------------------------------------
for (int i=0; i<4; ++i) {
for (int j=0; j<4; ++j) {
float3 A = cv[4*i + j].P;
BUCP[i] += A * B[j];
DUCP[i] += A * D[j];
#if OSD_COMPUTE_NORMAL_DERIVATIVES
CUCP[i] += A * C[j];
#endif //OSD_COMPUTE_NORMAL_DERIVATIVES
}
}
#endif //OSD_PATCH_ENABLE_SINGLE_CREASE
// ----------------------------------------------------------------
#else
LROW[0] = u0 * cv[ 0].P + u1 * cv[ 1].P + u2 * cv[ 2].P;
LROW[1] = u0 * cv[ 4].P + u1 * cv[ 5].P + u2 * cv[ 6].P;
LROW[2] = u0 * cv[ 8].P + u1 * cv[ 9].P + u2 * cv[10].P;
LROW[3] = u0 * cv[12].P + u1 * cv[13].P + u2 * cv[14].P;
#if OSD_COMPUTE_NORMAL_DERIVATIVES
// used for weingarten term
OsdUnivar4x4(UV.y, B, D, C);
RROW[0] = u0 * cv[ 1].P + u1 * cv[ 2].P + u2 * cv[ 3].P;
RROW[1] = u0 * cv[ 5].P + u1 * cv[ 6].P + u2 * cv[ 7].P;
RROW[2] = u0 * cv[ 9].P + u1 * cv[10].P + u2 * cv[11].P;
RROW[3] = u0 * cv[13].P + u1 * cv[14].P + u2 * cv[15].P;
#endif
P = B[0] * BUCP[0];
dPu = B[0] * DUCP[0];
dPv = D[0] * BUCP[0];
// Reduce 2x4 points to 2x2 -- two levels of linear interpolation in V
// and so 3 original pairs contributing to each of the 2 resulting:
float v = UV.y;
float vinv = 1.0f - v;
float3 dUU = B[0] * CUCP[0];
float3 dVV = C[0] * BUCP[0];
float3 dUV = D[0] * DUCP[0];
float v0 = vinv * vinv;
float v1 = v * vinv * 2.0f;
float v2 = v * v;
for (int k=1; k<4; ++k) {
P += B[k] * BUCP[k];
dPu += B[k] * DUCP[k];
dPv += D[k] * BUCP[k];
float3 LPAIR[2], RPAIR[2];
LPAIR[0] = v0 * LROW[0] + v1 * LROW[1] + v2 * LROW[2];
RPAIR[0] = v0 * RROW[0] + v1 * RROW[1] + v2 * RROW[2];
dUU += B[k] * CUCP[k];
dVV += C[k] * BUCP[k];
dUV += D[k] * DUCP[k];
}
LPAIR[1] = v0 * LROW[1] + v1 * LROW[2] + v2 * LROW[3];
RPAIR[1] = v0 * RROW[1] + v1 * RROW[2] + v2 * RROW[3];
// Interpolate points on the edges of the 2x2 bilinear hull from which
// both position and partials are trivially determined:
float3 DU0 = vinv * LPAIR[0] + v * LPAIR[1];
float3 DU1 = vinv * RPAIR[0] + v * RPAIR[1];
float3 DV0 = uinv * LPAIR[0] + u * RPAIR[0];
float3 DV1 = uinv * LPAIR[1] + u * RPAIR[1];
int level = OsdGetPatchFaceLevel(patchParam);
dPu *= 3 * level;
dPv *= 3 * level;
dPu = (DU1 - DU0) * 3 * level;
dPv = (DV1 - DV0) * 3 * level;
P = u * DU1 + uinv * DU0;
// Compute the normal and test for degeneracy:
//
// We need a geometric measure of the size of the patch for a suitable
// tolerance. Magnitudes of the partials are generally proportional to
// that size -- the sum of the partials is readily available, cheap to
// compute, and has proved effective in most cases (though not perfect).
// The size of the bounding box of the patch, or some approximation to
// it, would be better but more costly to compute.
//
float proportionalNormalTolerance = 0.00001f;
float nEpsilon = (length(dPu) + length(dPv)) * proportionalNormalTolerance;
N = cross(dPu, dPv);
float nLength = length(N);
if (nLength > nEpsilon) {
N = N / nLength;
} else {
float3 diagCross = cross(RPAIR[1] - LPAIR[0], LPAIR[1] - RPAIR[0]);
float diagCrossLength = length(diagCross);
if (diagCrossLength > nEpsilon) {
N = diagCross / diagCrossLength;
}
}
#ifndef OSD_COMPUTE_NORMAL_DERIVATIVES
dNu = float3(0,0,0);
dNv = float3(0,0,0);
#else
//
// Compute 2nd order partials of P(u,v) in order to compute 1st order partials
// for the un-normalized n(u,v) = dPu X dPv, then project into the tangent
// plane of normalized N. With resulting dNu and dNv we can make another
// attempt to resolve a still-degenerate normal.
//
// We don't use the Weingarten equations here as they require N != 0 and also
// are a little less numerically stable/accurate in single precision.
//
float B0u[4], B1u[4], B2u[4];
float B0v[4], B1v[4], B2v[4];
OsdUnivar4x4(UV.x, B0u, B1u, B2u);
OsdUnivar4x4(UV.y, B0v, B1v, B2v);
float3 dUU = float3(0,0,0);
float3 dVV = float3(0,0,0);
float3 dUV = float3(0,0,0);
for (int i=0; i<4; ++i) {
for (int j=0; j<4; ++j) {
#if OSD_PATCH_ENABLE_SINGLE_CREASE
int k = 4*i + j;
float3 CV = (s <= vSegments.x) ? cv[k].P
: ((s <= vSegments.y) ? cv[k].P1
: cv[k].P2);
#else
float3 CV = cv[4*i + j].P;
#endif
dUU += (B0v[i] * B2u[j]) * CV;
dVV += (B2v[i] * B0u[j]) * CV;
dUV += (B1v[i] * B1u[j]) * CV;
}
}
dUU *= 6 * level;
dVV *= 6 * level;
dUV *= 9 * level;
float3 n = cross(dPu, dPv);
float ln = 1.0 / length(n);
N = ln * n;
dNu = cross(dUU, dPv) + cross(dPu, dUV);
dNv = cross(dUV, dPv) + cross(dPu, dVV);
float E = dot(dPu, dPu);
float F = dot(dPu, dPv);
float G = dot(dPv, dPv);
float e = dot(N, dUU);
float f = dot(N, dUV);
float g = dot(N, dVV);
float EGFF = 1.0 / (E*G - F*F);
float nLengthInv = 1.0;
if (nLength > nEpsilon) {
nLengthInv = 1.0 / nLength;
} else {
// N may have been resolved above if degenerate, but if N was resolved
// we don't have an accurate length for its un-normalized value, and that
// length is needed to project the un-normalized dNu and dNv into the
// tangent plane of N.
//
// So compute N more accurately with available second derivatives, i.e.
// with a 1st order Taylor approximation to un-normalized N(u,v).
dNu = (f*F-e*G) * EGFF * dPu + (e*F-f*E) * EGFF * dPv;
dNv = (g*F-f*G) * EGFF * dPu + (f*F-g*E) * EGFF * dPv;
float DU = (UV.x == 1.0f) ? -1.0f : 1.0f;
float DV = (UV.y == 1.0f) ? -1.0f : 1.0f;
float powrn = 1.0 / powr(dot(n,n), 1.5);
N = DU * dNu + DV * dNv;
dNu = dNu * ln - n * (dot(dNu,n) * powrn);
dNv = dNv * ln - n * (dot(dNv,n) * powrn);
#else //OSD_COMPUTE_NORMAL_DERIVATIVES
OsdUnivar4x4(UV.y, B, D);
P = B[0] * BUCP[0];
dPu = B[0] * DUCP[0];
dPv = D[0] * BUCP[0];
for (int k=1; k<4; ++k) {
P += B[k] * BUCP[k];
dPu += B[k] * DUCP[k];
dPv += D[k] * BUCP[k];
nLength = length(N);
if (nLength > nEpsilon) {
nLengthInv = 1.0f / nLength;
N = N * nLengthInv;
}
}
int level = OsdGetPatchFaceLevel(patchParam);
dPu *= 3 * level;
dPv *= 3 * level;
N = normalize(cross(dPu, dPv));
dNu = float3(0,0,0);
dNv = float3(0,0,0);
#endif //OSD_COMPUTE_NORMAL_DERIVATIVES
// Project derivatives of non-unit normals into tangent plane of N:
dNu = (dNu - dot(dNu,N) * N) * nLengthInv;
dNv = (dNv - dot(dNv,N) * N) * nLengthInv;
#endif
}
// compute single-crease patch matrix