From 917fc36dde746910ac29fab7fb1892962269f483 Mon Sep 17 00:00:00 2001 From: David G Yu Date: Tue, 21 Nov 2017 18:01:57 -0800 Subject: [PATCH] HLSL patch shader changes for degenerate normals Updated HLSL patch shaders to resolve degenerate normals. This fix was ported from the GLSL patch shader source. Also, added missing inf sharp test cases to dxViewer. --- examples/dxViewer/init_shapes.h | 7 +- opensubdiv/osd/hlslPatchCommon.hlsl | 360 ++++++++++++++-------------- 2 files changed, 179 insertions(+), 188 deletions(-) diff --git a/examples/dxViewer/init_shapes.h b/examples/dxViewer/init_shapes.h index dc80d707..5dd4a350 100644 --- a/examples/dxViewer/init_shapes.h +++ b/examples/dxViewer/init_shapes.h @@ -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)); @@ -71,6 +74,7 @@ static void initShapes() { g_defaultShapes.push_back(ShapeDesc("catmark_gregory_test5", catmark_gregory_test5, kCatmark)); g_defaultShapes.push_back(ShapeDesc("catmark_gregory_test6", catmark_gregory_test6, kCatmark)); g_defaultShapes.push_back(ShapeDesc("catmark_gregory_test7", catmark_gregory_test7, kCatmark)); + g_defaultShapes.push_back( ShapeDesc("catmark_gregory_test8", catmark_gregory_test8, kCatmark ) ); g_defaultShapes.push_back(ShapeDesc("catmark_hole_test1", catmark_hole_test1, kCatmark)); g_defaultShapes.push_back(ShapeDesc("catmark_hole_test2", catmark_hole_test2, kCatmark)); g_defaultShapes.push_back(ShapeDesc("catmark_hole_test3", catmark_hole_test3, kCatmark)); @@ -90,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 ) ); @@ -118,5 +121,5 @@ static void initShapes() { g_defaultShapes.push_back(ShapeDesc("loop_pole8", loop_pole8, kLoop)); g_defaultShapes.push_back(ShapeDesc("loop_pole64", loop_pole64, kLoop)); g_defaultShapes.push_back(ShapeDesc("loop_pole360", loop_pole360, kLoop)); - } +//------------------------------------------------------------------------------ diff --git a/opensubdiv/osd/hlslPatchCommon.hlsl b/opensubdiv/osd/hlslPatchCommon.hlsl index e306f9f9..72727840 100644 --- a/opensubdiv/osd/hlslPatchCommon.hlsl +++ b/opensubdiv/osd/hlslPatchCommon.hlsl @@ -1105,113 +1105,187 @@ OsdEvalPatchBezier(int3 patchParam, float2 UV, out float3 P, out float3 dPu, out float3 dPv, out float3 N, out float3 dNu, out float3 dNv) { -#ifdef 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. + // - // ---------------------------------------------------------------- -#if defined OSD_PATCH_ENABLE_SINGLE_CREASE + // 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]; +#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 float2 vSegments = cv[0].vSegments; 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]; -#ifdef OSD_COMPUTE_NORMAL_DERIVATIVES - CUCP[i] += A * C[j]; -#endif + 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; } } +#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; + + 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]; + + 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 = (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 - // ---------------------------------------------------------------- - 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]; -#ifdef OSD_COMPUTE_NORMAL_DERIVATIVES - CUCP[i] += A * C[j]; -#endif - } - } -#endif - // ---------------------------------------------------------------- + // + // 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]; - P = float3(0,0,0); - dPu = float3(0,0,0); - dPv = float3(0,0,0); - -#ifdef OSD_COMPUTE_NORMAL_DERIVATIVES - // used for weingarten term - OsdUnivar4x4(UV.y, B, D, C); + 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 k=0; k<4; ++k) { - P += B[k] * BUCP[k]; - dPu += B[k] * DUCP[k]; - dPv += D[k] * BUCP[k]; - - dUU += B[k] * CUCP[k]; - dVV += C[k] * BUCP[k]; - dUV += D[k] * DUCP[k]; + for (int i=0; i<4; ++i) { + for (int j=0; j<4; ++j) { +#ifdef 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; + } } - 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); + 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 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)/(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; + float DU = (UV.x == 1.0f) ? -1.0f : 1.0f; + float DV = (UV.y == 1.0f) ? -1.0f : 1.0f; - 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 - OsdUnivar4x4(UV.y, B, D); + N = DU * dNu + DV * dNv; - for (int k=0; 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); + // 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 } @@ -1266,113 +1340,27 @@ OsdEvalPatchGregory(int3 patchParam, float2 UV, float3 cv[20], float d21 = u+V; float d22 = U+V; - float3 q[16]; + OsdPerPatchVertexBezier bezcv[16]; - 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); - -#ifdef 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)/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]; - 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 + OsdEvalPatchBezier(patchParam, UV, bezcv, P, dPu, dPv, N, dNu, dNv); } // ----------------------------------------------------------------------------