From ee4085b388c52262afd8cc2319cf4d379ce07ab5 Mon Sep 17 00:00:00 2001 From: Robert Beckebans Date: Thu, 8 Apr 2021 12:06:14 +0200 Subject: [PATCH] Added spherical harmonics math --- neo/idlib/Lib.h | 1 + neo/idlib/math/SphericalHarmonics.h | 427 ++++++++++++++++++++++++++++ 2 files changed, 428 insertions(+) create mode 100644 neo/idlib/math/SphericalHarmonics.h diff --git a/neo/idlib/Lib.h b/neo/idlib/Lib.h index 60f3ad62..94a4a17d 100644 --- a/neo/idlib/Lib.h +++ b/neo/idlib/Lib.h @@ -275,6 +275,7 @@ public: #include "math/Curve.h" #include "math/Ode.h" #include "math/Lcp.h" +#include "math/SphericalHarmonics.h" // bounding volumes #include "bv/Sphere.h" diff --git a/neo/idlib/math/SphericalHarmonics.h b/neo/idlib/math/SphericalHarmonics.h new file mode 100644 index 00000000..939e5c60 --- /dev/null +++ b/neo/idlib/math/SphericalHarmonics.h @@ -0,0 +1,427 @@ +/* +The MIT License (MIT) + +Copyright (c) 2015 Yuriy O'Donnell +Copyright (c) 2021 Robert Beckebans (id Tech 4.x integration) + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +*/ + +#ifndef __MATH_SPHERICAL_HARMONICS_H__ +#define __MATH_SPHERICAL_HARMONICS_H__ + +// https://graphics.stanford.edu/papers/envmap/envmap.pdf + +template +struct SphericalHarmonicsT +{ + T data[( L + 1 ) * ( L + 1 )]; + + const T& operator []( size_t i ) const + { + return data[i]; + } + T& operator []( size_t i ) + { + return data[i]; + } + + T& at( int l, int m ) + { + return data[l * l + l + m]; + } + const T& at( int l, int m ) const + { + return data[l * l + l + m]; + } +}; + +typedef SphericalHarmonicsT SphericalHarmonicsL1; +typedef SphericalHarmonicsT SphericalHarmonicsL2; +typedef SphericalHarmonicsT SphericalHarmonicsL1RGB; +typedef SphericalHarmonicsT SphericalHarmonicsL2RGB; + +template +SphericalHarmonicsL1 shEvaluateL1( idVec3 p ); +SphericalHarmonicsL2 shEvaluateL2( idVec3 p ); + +inline size_t shSize( size_t L ) +{ + return ( L + 1 ) * ( L + 1 ); +} + +template +inline void shAddWeighted( SphericalHarmonicsT& accumulatorSh, const SphericalHarmonicsT& sh, const Tw& weight ) +{ + for( size_t i = 0; i < shSize( L ); ++i ) + { + accumulatorSh[i] += sh[i] * weight; + } +} + +template +inline Ta shDot( const SphericalHarmonicsT& shA, const SphericalHarmonicsT& shB ) +{ + Ta result = Ta( 0 ); + for( size_t i = 0; i < shSize( L ); ++i ) + { + result += shA[i] * shB[i]; + } + return result; +} + +template +inline SphericalHarmonicsT shEvaluate( idVec3 p ) +{ + // From Peter-Pike Sloan's Stupid SH Tricks + // http://www.ppsloan.org/publications/StupidSH36.pdf + // https://github.com/dariomanesku/cmft/blob/master/src/cmft/cubemapfilter.cpp#L130 + + static_assert( L <= 4, "Spherical Harmonics above L4 are not supported" ); + + SphericalHarmonicsT result; + + const float x = -p.x; + const float y = -p.y; + const float z = p.z; + + const float x2 = x * x; + const float y2 = y * y; + const float z2 = z * z; + + const float z3 = z2 * z; + + const float x4 = x2 * x2; + const float y4 = y2 * y2; + const float z4 = z2 * z2; + + const float sqrtPi = sqrt( idMath::PI ); + + size_t i = 0; + + result[i++] = 1.0f / ( 2.0f * sqrtPi ); + + if( L >= 1 ) + { + result[i++] = -sqrt( 3.0f / ( 4.0f * idMath::PI ) ) * y; + result[i++] = sqrt( 3.0f / ( 4.0f * idMath::PI ) ) * z; + result[i++] = -sqrt( 3.0f / ( 4.0f * idMath::PI ) ) * x; + } + + if( L >= 2 ) + { + result[i++] = sqrt( 15.0f / ( 4.0f * idMath::PI ) ) * y * x; + result[i++] = -sqrt( 15.0f / ( 4.0f * idMath::PI ) ) * y * z; + result[i++] = sqrt( 5.0f / ( 16.0f * idMath::PI ) ) * ( 3.0f * z2 - 1.0f ); + result[i++] = -sqrt( 15.0f / ( 4.0f * idMath::PI ) ) * x * z; + result[i++] = sqrt( 15.0f / ( 16.0f * idMath::PI ) ) * ( x2 - y2 ); + } + + if( L >= 3 ) + { + result[i++] = -sqrt( 70.0f / ( 64.0f * idMath::PI ) ) * y * ( 3.0f * x2 - y2 ); + result[i++] = sqrt( 105.0f / ( 4.0f * idMath::PI ) ) * y * x * z; + result[i++] = -sqrt( 21.0f / ( 16.0f * idMath::PI ) ) * y * ( -1.0f + 5.0f * z2 ); + result[i++] = sqrt( 7.0f / ( 16.0f * idMath::PI ) ) * ( 5.0f * z3 - 3.0f * z ); + result[i++] = -sqrt( 42.0f / ( 64.0f * idMath::PI ) ) * x * ( -1.0f + 5.0f * z2 ); + result[i++] = sqrt( 105.0f / ( 16.0f * idMath::PI ) ) * ( x2 - y2 ) * z; + result[i++] = -sqrt( 70.0f / ( 64.0f * idMath::PI ) ) * x * ( x2 - 3.0f * y2 ); + } + + if( L >= 4 ) + { + result[i++] = 3.0f * sqrt( 35.0f / ( 16.0f * idMath::PI ) ) * x * y * ( x2 - y2 ); + result[i++] = -3.0f * sqrt( 70.0f / ( 64.0f * idMath::PI ) ) * y * z * ( 3.0f * x2 - y2 ); + result[i++] = 3.0f * sqrt( 5.0f / ( 16.0f * idMath::PI ) ) * y * x * ( -1.0f + 7.0f * z2 ); + result[i++] = -3.0f * sqrt( 10.0f / ( 64.0f * idMath::PI ) ) * y * z * ( -3.0f + 7.0f * z2 ); + result[i++] = ( 105.0f * z4 - 90.0f * z2 + 9.0f ) / ( 16.0f * sqrtPi ); + result[i++] = -3.0f * sqrt( 10.0f / ( 64.0f * idMath::PI ) ) * x * z * ( -3.0f + 7.0f * z2 ); + result[i++] = 3.0f * sqrt( 5.0f / ( 64.0f * idMath::PI ) ) * ( x2 - y2 ) * ( -1.0f + 7.0f * z2 ); + result[i++] = -3.0f * sqrt( 70.0f / ( 64.0f * idMath::PI ) ) * x * z * ( x2 - 3.0f * y2 ); + result[i++] = 3.0f * sqrt( 35.0f / ( 4.0f * ( 64.0f * idMath::PI ) ) ) * ( x4 - 6.0f * y2 * x2 + y4 ); + } + + return result; +} + +inline SphericalHarmonicsL1 shEvaluateL1( idVec3 p ) +{ + return shEvaluate<1>( p ); +} + +inline SphericalHarmonicsL2 shEvaluateL2( idVec3 p ) +{ + return shEvaluate<2>( p ); +} + +inline float shEvaluateDiffuseL1Geomerics( const SphericalHarmonicsL1& sh, const idVec3& n ) +{ + // http://www.geomerics.com/wp-content/uploads/2015/08/CEDEC_Geomerics_ReconstructingDiffuseLighting1.pdf + + float R0 = sh[0]; + + idVec3 R1 = 0.5f * idVec3( sh[3], sh[1], sh[2] ); + float lenR1 = R1.Length(); + + //float q = 0.5f * (1.0f + dot(R1 / lenR1, n)); + float q = 0.5f * ( 1.0f + ( R1 / lenR1 ) * n ); + + float p = 1.0f + 2.0f * lenR1 / R0; + float a = ( 1.0f - lenR1 / R0 ) / ( 1.0f + lenR1 / R0 ); + + return R0 * ( a + ( 1.0f - a ) * ( p + 1.0f ) * pow( q, p ) ); +} + +template +inline SphericalHarmonicsT shConvolveDiffuse( SphericalHarmonicsT& sh ) +{ + SphericalHarmonicsT result; + + // https://cseweb.ucsd.edu/~ravir/papers/envmap/envmap.pdf equation 8 + + const float A[5] = + { + pi, + pi * 2.0f / 3.0f, + pi * 1.0f / 4.0f, + 0.0f, + -pi * 1.0f / 24.0f + }; + + int i = 0; + for( int l = 0; l <= ( int )L; ++l ) + { + for( int m = -l; m <= l; ++m ) + { + result[i] = sh[i] * A[l]; + ++i; + } + } + + return result; +} + +template +inline T shEvaluateDiffuse( const SphericalHarmonicsT& sh, const idVec3& direction ) +{ + static_assert( L <= 4, "Spherical Harmonics above L4 are not supported" ); + + SphericalHarmonicsT directionSh = shEvaluate( direction ); + + // https://cseweb.ucsd.edu/~ravir/papers/envmap/envmap.pdf equation 8 + + const float A[5] = + { + pi, + pi * 2.0f / 3.0f, + pi * 1.0f / 4.0f, + 0.0f, + -pi * 1.0f / 24.0f + }; + + size_t i = 0; + + T result = sh[i] * directionSh[i] * A[0]; + ++i; + + if( L >= 1 ) + { + result += sh[i] * directionSh[i] * A[1]; + ++i; + result += sh[i] * directionSh[i] * A[1]; + ++i; + result += sh[i] * directionSh[i] * A[1]; + ++i; + } + + if( L >= 2 ) + { + result += sh[i] * directionSh[i] * A[2]; + ++i; + result += sh[i] * directionSh[i] * A[2]; + ++i; + result += sh[i] * directionSh[i] * A[2]; + ++i; + result += sh[i] * directionSh[i] * A[2]; + ++i; + result += sh[i] * directionSh[i] * A[2]; + ++i; + } + + // L3 and other odd bands > 1 have 0 factor + + if( L >= 4 ) + { + i = 16; + + result += sh[i] * directionSh[i] * A[4]; + ++i; + result += sh[i] * directionSh[i] * A[4]; + ++i; + result += sh[i] * directionSh[i] * A[4]; + ++i; + result += sh[i] * directionSh[i] * A[4]; + ++i; + result += sh[i] * directionSh[i] * A[4]; + ++i; + result += sh[i] * directionSh[i] * A[4]; + ++i; + result += sh[i] * directionSh[i] * A[4]; + ++i; + result += sh[i] * directionSh[i] * A[4]; + ++i; + result += sh[i] * directionSh[i] * A[4]; + ++i; + } + + return result; +} + + +template +inline T shEvaluateDiffuseL1( const SphericalHarmonicsT& sh, const idVec3& direction ) +{ + return shEvaluateDiffuse( sh, direction ); +} + +template +inline T shEvaluateDiffuseL2( const SphericalHarmonicsT& sh, const idVec3& direction ) +{ + return shEvaluateDiffuse( sh, direction ); +} + +template +float shFindWindowingLambda( const SphericalHarmonicsT& sh, float maxLaplacian ) +{ + // http://www.ppsloan.org/publications/StupidSH36.pdf + // Appendix A7 Solving for Lambda to Reduce the Squared Laplacian + + float tableL[L + 1]; + float tableB[L + 1]; + + tableL[0] = 0.0f; + tableB[0] = 0.0f; + + for( int l = 1; l <= ( int )L; ++l ) + { + tableL[l] = float( sqr( l ) * sqr( l + 1 ) ); + + float B = 0.0f; + for( int m = -1; m <= l; ++m ) + { + B += sqr( sh.at( l, m ) ); + } + tableB[l] = B; + } + + float squaredLaplacian = 0.0f; + + for( int l = 1; l <= ( int )L; ++l ) + { + squaredLaplacian += tableL[l] * tableB[l]; + } + + const float targetSquaredLaplacian = maxLaplacian * maxLaplacian; + if( squaredLaplacian <= targetSquaredLaplacian ) + { + return 0.0f; + } + + float lambda = 0.0f; + + const u32 iterationLimit = 10000000; + for( u32 i = 0; i < iterationLimit; ++i ) + { + float f = 0.0f; + float fd = 0.0f; + + for( int l = 1; l <= ( int )L; ++l ) + { + f += tableL[l] * tableB[l] / sqr( 1.0f + lambda * tableL[l] ); + fd += ( 2.0f * sqr( tableL[l] ) * tableB[l] ) / cube( 1.0f + lambda * tableL[l] ); + } + + f = targetSquaredLaplacian - f; + + float delta = -f / fd; + lambda += delta; + + if( abs( delta ) < 1e-6f ) + { + break; + } + } + + return lambda; +} + +template +void shApplyWindowing( SphericalHarmonicsT& sh, float lambda ) +{ + // From Peter-Pike Sloan's Stupid SH Tricks + // http://www.ppsloan.org/publications/StupidSH36.pdf + + int i = 0; + for( int l = 0; l <= ( int )L; ++l ) + { + float s = 1.0f / ( 1.0f + lambda * l * l * ( l + 1.0f ) * ( l + 1.0f ) ); + for( int m = -l; m <= l; ++m ) + { + sh[i++] *= s; + } + } +} + +#if 0 +template +inline T shMeanSquareError( const SphericalHarmonicsT& sh, const idArray& radianceSamples ) +{ + T errorSquaredSum = T( 0.0f ); + + for( const RadianceSample& sample : radianceSamples ) + { + auto directionSh = shEvaluate( sample.direction ); + auto reconstructedValue = shDot( sh, directionSh ); + auto error = sample.value - reconstructedValue; + errorSquaredSum += error * error; + } + + float sampleWeight = 1.0f / radianceSamples.size(); + return errorSquaredSum * sampleWeight; +} + +template +inline float shMeanSquareErrorScalar( const SphericalHarmonicsT& sh, const idArray& radianceSamples ) +{ + return dot( shMeanSquareError( sh, radianceSamples ), T( 1.0f / 3.0f ) ); +} +#endif + +template +inline SphericalHarmonicsT shLuminance( const SphericalHarmonicsT& sh ) +{ + SphericalHarmonicsT result; + for( size_t i = 0; i < shSize( L ); ++i ) + { + result[i] = rgbLuminance( sh[i] ); + } + return result; +} + +#endif // __MATH_SPHERICAL_HARMONICS_H__