Merge pull request #948 from barfowl/glsl_patch_normals

Improved GLSL patch shaders to compute normals in common degenerate cases
This commit is contained in:
David G Yu 2017-12-12 08:50:30 -08:00 committed by GitHub
commit 226f1a27fd
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23

View File

@ -1230,113 +1230,187 @@ OsdEvalPatchBezier(ivec3 patchParam, vec2 UV,
out vec3 P, out vec3 dPu, out vec3 dPv, out vec3 P, out vec3 dPu, out vec3 dPv,
out vec3 N, out vec3 dNu, out vec3 dNv) out vec3 N, out vec3 dNu, out vec3 dNv)
{ {
#ifdef OSD_COMPUTE_NORMAL_DERIVATIVES //
float B[4], D[4], C[4]; // Use the recursive nature of the basis functions to compute a 2x2 set
vec3 BUCP[4] = vec3[4](vec3(0), vec3(0), vec3(0), vec3(0)), // of intermediate points (via repeated linear interpolation). These
DUCP[4] = vec3[4](vec3(0), vec3(0), vec3(0), vec3(0)), // points define a bilinear surface tangent to the desired surface at P
CUCP[4] = vec3[4](vec3(0), vec3(0), vec3(0), vec3(0)); // and so containing dPu and dPv. The cost of computing P, dPu and dPv
OsdUnivar4x4(UV.x, B, D, C); // this way is comparable to that of typical tensor product evaluation
#else // (if not faster).
float B[4], D[4]; //
vec3 BUCP[4] = vec3[4](vec3(0), vec3(0), vec3(0), vec3(0)), // If N = dPu X dPv degenerates, it often results from an edge of the
DUCP[4] = vec3[4](vec3(0), vec3(0), vec3(0), vec3(0)); // 2x2 bilinear hull collapsing or two adjacent edges colinear. In both
OsdUnivar4x4(UV.x, B, D); // cases, the expected non-planar quad degenerates into a triangle, and
#endif // the tangent plane of that triangle provides the desired normal N.
//
// ---------------------------------------------------------------- // Reduce 4x4 points to 2x4 -- two levels of linear interpolation in U
#if defined OSD_PATCH_ENABLE_SINGLE_CREASE // 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;
vec3 LROW[4], RROW[4];
#ifndef OSD_PATCH_ENABLE_SINGLE_CREASE
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;
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;
#else
vec2 vSegments = cv[0].vSegments; vec2 vSegments = cv[0].vSegments;
float s = OsdGetPatchSingleCreaseSegmentParameter(patchParam, UV); float s = OsdGetPatchSingleCreaseSegmentParameter(patchParam, UV);
for (int i=0; i<4; ++i) { for (int i = 0; i < 4; ++i) {
for (int j=0; j<4; ++j) { int j = i*4;
int k = 4*i + j; if (s <= vSegments.x) {
LROW[i] = u0 * cv[ j ].P + u1 * cv[j+1].P + u2 * cv[j+2].P;
vec3 A = (s <= vSegments.x) ? cv[k].P RROW[i] = u0 * cv[j+1].P + u1 * cv[j+2].P + u2 * cv[j+3].P;
: ((s <= vSegments.y) ? cv[k].P1 } else if (s <= vSegments.y) {
: cv[k].P2); 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;
BUCP[i] += A * B[j]; } else {
DUCP[i] += A * D[j]; LROW[i] = u0 * cv[ j ].P2 + u1 * cv[j+1].P2 + u2 * cv[j+2].P2;
#ifdef OSD_COMPUTE_NORMAL_DERIVATIVES RROW[i] = u0 * cv[j+1].P2 + u1 * cv[j+2].P2 + u2 * cv[j+3].P2;
CUCP[i] += A * C[j];
#endif
} }
} }
#endif
// 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;
float v0 = vinv * vinv;
float v1 = v * vinv * 2.0f;
float v2 = v * v;
vec3 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];
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:
vec3 DU0 = vinv * LPAIR[0] + v * LPAIR[1];
vec3 DU1 = vinv * RPAIR[0] + v * RPAIR[1];
vec3 DV0 = uinv * LPAIR[0] + u * RPAIR[0];
vec3 DV1 = uinv * LPAIR[1] + u * RPAIR[1];
int level = OsdGetPatchFaceLevel(patchParam);
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 {
vec3 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 = vec3(0);
dNv = vec3(0);
#else #else
// ---------------------------------------------------------------- //
for (int i=0; i<4; ++i) { // Compute 2nd order partials of P(u,v) in order to compute 1st order partials
for (int j=0; j<4; ++j) { // for the un-normalized n(u,v) = dPu X dPv, then project into the tangent
vec3 A = cv[4*i + j].P; // plane of normalized N. With resulting dNu and dNv we can make another
BUCP[i] += A * B[j]; // attempt to resolve a still-degenerate normal.
DUCP[i] += A * D[j]; //
#ifdef OSD_COMPUTE_NORMAL_DERIVATIVES // We don't use the Weingarten equations here as they require N != 0 and also
CUCP[i] += A * C[j]; // are a little less numerically stable/accurate in single precision.
#endif //
} float B0u[4], B1u[4], B2u[4];
} float B0v[4], B1v[4], B2v[4];
#endif
// ----------------------------------------------------------------
P = vec3(0); OsdUnivar4x4(UV.x, B0u, B1u, B2u);
dPu = vec3(0); OsdUnivar4x4(UV.y, B0v, B1v, B2v);
dPv = vec3(0);
#ifdef OSD_COMPUTE_NORMAL_DERIVATIVES
// used for weingarten term
OsdUnivar4x4(UV.y, B, D, C);
vec3 dUU = vec3(0); vec3 dUU = vec3(0);
vec3 dVV = vec3(0); vec3 dVV = vec3(0);
vec3 dUV = vec3(0); vec3 dUV = vec3(0);
for (int k=0; k<4; ++k) { for (int i=0; i<4; ++i) {
P += B[k] * BUCP[k]; for (int j=0; j<4; ++j) {
dPu += B[k] * DUCP[k]; #ifdef OSD_PATCH_ENABLE_SINGLE_CREASE
dPv += D[k] * BUCP[k]; int k = 4*i + j;
vec3 CV = (s <= vSegments.x) ? cv[k].P
dUU += B[k] * CUCP[k]; : ((s <= vSegments.y) ? cv[k].P1
dVV += C[k] * BUCP[k]; : cv[k].P2);
dUV += D[k] * DUCP[k]; #else
vec3 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;
}
} }
int level = OsdGetPatchFaceLevel(patchParam);
dPu *= 3 * level;
dPv *= 3 * level;
dUU *= 6 * level; dUU *= 6 * level;
dVV *= 6 * level; dVV *= 6 * level;
dUV *= 9 * level; dUV *= 9 * level;
vec3 n = cross(dPu, dPv); dNu = cross(dUU, dPv) + cross(dPu, dUV);
N = normalize(n); dNv = cross(dUV, dPv) + cross(dPu, dVV);
float E = dot(dPu, dPu); float nLengthInv = 1.0;
float F = dot(dPu, dPv); if (nLength > nEpsilon) {
float G = dot(dPv, dPv); nLengthInv = 1.0 / nLength;
float e = dot(N, dUU); } else {
float f = dot(N, dUV); // N may have been resolved above if degenerate, but if N was resolved
float g = dot(N, dVV); // 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)/(E*G-F*F) * dPu + (e*F-f*E)/(E*G-F*F) * dPv; float DU = (UV.x == 1.0f) ? -1.0f : 1.0f;
dNv = (g*F-f*G)/(E*G-F*F) * dPu + (f*F-g*E)/(E*G-F*F) * dPv; float DV = (UV.y == 1.0f) ? -1.0f : 1.0f;
dNu = dNu/length(n) - n * (dot(dNu,n)/pow(dot(n,n), 1.5)); N = DU * dNu + DV * dNv;
dNv = dNv/length(n) - n * (dot(dNv,n)/pow(dot(n,n), 1.5));
#else
OsdUnivar4x4(UV.y, B, D);
for (int k=0; k<4; ++k) { nLength = length(N);
P += B[k] * BUCP[k]; if (nLength > nEpsilon) {
dPu += B[k] * DUCP[k]; nLengthInv = 1.0f / nLength;
dPv += D[k] * BUCP[k]; N = N * nLengthInv;
}
} }
int level = OsdGetPatchFaceLevel(patchParam);
dPu *= 3 * level;
dPv *= 3 * level;
N = normalize(cross(dPu, dPv)); // Project derivatives of non-unit normals into tangent plane of N:
dNu = vec3(0); dNu = (dNu - dot(dNu,N) * N) * nLengthInv;
dNv = vec3(0); dNv = (dNv - dot(dNv,N) * N) * nLengthInv;
#endif #endif
} }
@ -1391,113 +1465,27 @@ OsdEvalPatchGregory(ivec3 patchParam, vec2 UV, vec3 cv[20],
float d21 = u+V; float d21 = u+V;
float d22 = U+V; float d22 = U+V;
vec3 q[16]; OsdPerPatchVertexBezier bezcv[16];
q[ 5] = (d11 == 0.0) ? cv[3] : (u*cv[3] + v*cv[4])/d11; bezcv[ 5].P = (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; bezcv[ 6].P = (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; bezcv[ 9].P = (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[10].P = (d22 == 0.0) ? cv[13] : (U*cv[13] + V*cv[14])/d22;
q[ 0] = cv[0]; bezcv[ 0].P = cv[0];
q[ 1] = cv[1]; bezcv[ 1].P = cv[1];
q[ 2] = cv[7]; bezcv[ 2].P = cv[7];
q[ 3] = cv[5]; bezcv[ 3].P = cv[5];
q[ 4] = cv[2]; bezcv[ 4].P = cv[2];
q[ 7] = cv[6]; bezcv[ 7].P = cv[6];
q[ 8] = cv[16]; bezcv[ 8].P = cv[16];
q[11] = cv[12]; bezcv[11].P = cv[12];
q[12] = cv[15]; bezcv[12].P = cv[15];
q[13] = cv[17]; bezcv[13].P = cv[17];
q[14] = cv[11]; bezcv[14].P = cv[11];
q[15] = cv[10]; bezcv[15].P = cv[10];
P = vec3(0); OsdEvalPatchBezier(patchParam, UV, bezcv, P, dPu, dPv, N, dNu, dNv);
dPu = vec3(0);
dPv = vec3(0);
#ifdef OSD_COMPUTE_NORMAL_DERIVATIVES
float B[4], D[4], C[4];
vec3 BUCP[4] = vec3[4](vec3(0), vec3(0), vec3(0), vec3(0)),
DUCP[4] = vec3[4](vec3(0), vec3(0), vec3(0), vec3(0)),
CUCP[4] = vec3[4](vec3(0), vec3(0), vec3(0), vec3(0));
vec3 dUU = vec3(0);
vec3 dVV = vec3(0);
vec3 dUV = vec3(0);
OsdUnivar4x4(UV.x, B, D, C);
for (int i=0; i<4; ++i) {
for (int j=0; j<4; ++j) {
vec3 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;
vec3 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)/pow(dot(n,n), 1.5));
dNv = dNv/length(n) - n * (dot(dNv,n)/pow(dot(n,n), 1.5));
#else
float B[4], D[4];
vec3 BUCP[4] = vec3[4](vec3(0), vec3(0), vec3(0), vec3(0)),
DUCP[4] = vec3[4](vec3(0), vec3(0), vec3(0), vec3(0));
OsdUnivar4x4(UV.x, B, D);
for (int i=0; i<4; ++i) {
for (int j=0; j<4; ++j) {
vec3 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 = vec3(0);
dNv = vec3(0);
#endif
} }
// ---------------------------------------------------------------------------- // ----------------------------------------------------------------------------