alexpete
730 lines
26 KiB
C++
730 lines
26 KiB
C++
/*
|
|
* All or portions of this file Copyright (c) Amazon.com, Inc. or its affiliates or
|
|
* its licensors.
|
|
*
|
|
* For complete copyright and license terms please see the LICENSE at the root of this
|
|
* distribution (the "License"). All use of this software is governed by the License,
|
|
* or, if provided, by the license below or the license accompanying this file. Do not
|
|
* remove or modify any license notices. This file is distributed on an "AS IS" BASIS,
|
|
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
|
*
|
|
*/
|
|
// Original file Copyright Crytek GMBH or its affiliates, used under license.
|
|
|
|
#include "Cry3DEngine_precompiled.h"
|
|
|
|
#define IGNORE_ASSERTS
|
|
#if defined(_DEBUG) && defined(IGNORE_ASSERTS)
|
|
# undef assert
|
|
# define assert(cond) ((void)0)
|
|
#endif
|
|
|
|
#include "SkyLightNishita.h"
|
|
|
|
#include <math.h>
|
|
|
|
// constant definitions (all heights & radii given in km or km^-1 )
|
|
const f64 c_maxAtmosphereHeight(100.0);
|
|
const f64 c_earthRadius(6368.0);
|
|
const f32 c_earthRadiusf(6368.0f);
|
|
const f64 c_avgDensityHeightMieInv(1.0 / 1.2);
|
|
const f64 c_avgDensityHeightRayleighInv(1.0 / 7.994);
|
|
|
|
const f64 c_opticalDepthWhenHittingEarth(1e10);
|
|
|
|
const f64 c_pi(3.1415926535897932384626433832795);
|
|
const f32 c_pif(3.1415926535897932384626433832795f);
|
|
|
|
// Machine epsilon is too small to catch rounding error asserts here. We use a large enough number to prevent rounding errors, but small
|
|
// enough to still catch invalid conditions (10^-6).
|
|
static const float floatDiffFactor = 0.000001f;
|
|
|
|
// constants for optical LUT serialization
|
|
const uint32 c_lutFileTag(0x4C594B53); // "SKYL"
|
|
const uint32 c_lutFileVersion(0x00010002);
|
|
const char c_lutFileName[] = "engineassets/sky/optical.lut";
|
|
|
|
static inline f64 MapSaveExpArg(f64 arg)
|
|
{
|
|
const f64 c_saveExpArgRange((f64)650.0); // -650.0 to 650 range is safe not to introduce fp over-/underflows
|
|
return((arg < -c_saveExpArgRange) ? -c_saveExpArgRange : (arg > c_saveExpArgRange) ? c_saveExpArgRange : arg);
|
|
}
|
|
|
|
|
|
static inline f64 exp_precise(f64 arg)
|
|
{
|
|
return(exp(MapSaveExpArg((f64) arg)));
|
|
}
|
|
|
|
namespace
|
|
{
|
|
union eco
|
|
{
|
|
f64 d;
|
|
struct
|
|
{
|
|
int32 i, j;
|
|
} n;
|
|
};
|
|
}
|
|
|
|
static inline f64 exp_fast(f64 arg)
|
|
{
|
|
const f64 eco_m(1048576L / 0.693147180559945309417232121458177);
|
|
const f64 eco_a(1072693248L - 60801L);
|
|
|
|
#if defined(_CPU_X86) || defined(_CPU_AMD64) || defined(_CPU_ARM)// for little endian (tested on Win32 / Win64)
|
|
eco e;
|
|
# ifdef _DEBUG
|
|
e.d = 1.0;
|
|
assert(e.n.j - 1072693248L || e.n.i == 0); // check IEEE-754 conformance
|
|
# endif
|
|
e.n.j = (int32) (eco_m * MapSaveExpArg(arg) + eco_a);
|
|
return((f64)e.d);
|
|
#elif defined(_CPU_G5)
|
|
eco e;
|
|
# ifdef _DEBUG
|
|
e.d = 1.0;
|
|
assert(e.n.i == 1072693248L || e.n.j == 0); // check IEEE-754 conformance
|
|
# endif
|
|
e.n.i = (int32) (eco_m * MapSaveExpArg(arg) + eco_a);
|
|
return((f64)e.d);
|
|
#else // fall back to default exp_sky() implementation for untested/unsupported target platforms
|
|
# pragma message( "Optimized exp_fast() not available for this platform!" )
|
|
# pragma message( "If your target CPU is IEEE-754 conformant then please specify it in either the little or big endian branch (see SkyLightNishita.cpp::exp_fast())." )
|
|
return(exp(arg));
|
|
#endif
|
|
}
|
|
|
|
static inline f64 OpticalScaleFunction(const f64& height, const f64& avgDensityHeightInv)
|
|
{
|
|
assert(height >= 0.0);
|
|
assert(avgDensityHeightInv > 0.0 && avgDensityHeightInv <= 1.0);
|
|
return(exp_precise(-height * avgDensityHeightInv));
|
|
}
|
|
|
|
|
|
static inline f64 IntegrateOpticalDepthInternal(const Vec3d& start, const f64& startScale,
|
|
const Vec3d& end, const f64& endScale, const f64& avgDensityHeightInv, const f64& error)
|
|
{
|
|
assert(_finite(startScale) && _finite(endScale));
|
|
|
|
Vec3d mid(0.5 * (start + end));
|
|
f64 midScale(OpticalScaleFunction(mid.GetLength() - c_earthRadius, avgDensityHeightInv));
|
|
|
|
if (fabs(startScale - midScale) <= error && fabs(midScale - endScale) <= error)
|
|
{
|
|
// integrate section this via simpson rule and stop recursing
|
|
const f64 c_oneSixth(1.0 / 6.0);
|
|
return((startScale + 4.0 * midScale + endScale) * c_oneSixth * (end - start).GetLength());
|
|
}
|
|
else
|
|
{
|
|
// refine section via recursing down left and right branch
|
|
return(IntegrateOpticalDepthInternal(start, startScale, mid, midScale, avgDensityHeightInv, error) +
|
|
IntegrateOpticalDepthInternal(mid, midScale, end, endScale, avgDensityHeightInv, error));
|
|
}
|
|
}
|
|
|
|
|
|
CSkyLightNishita::CSkyLightNishita()
|
|
: m_opticalDepthLUT()
|
|
, m_opticalScaleLUT()
|
|
, m_phaseLUT()
|
|
, m_Km(0.0f)
|
|
, m_Kr(0.0f)
|
|
, m_sunIntensity(20.0f, 20.0f, 20.0f)
|
|
, m_g(0.0f)
|
|
, m_invRGBWaveLength4(1.0f, 1.0f, 1.0f)
|
|
, m_sunDir(0.0f, 0.707106f, 0.707106f)
|
|
, m_inScatteringStepSize(1)
|
|
{
|
|
SetRGBWaveLengths(Vec3(650.0f, 570.0f, 475.0f));
|
|
SetSunDirection(Vec3(0.0f, 0.707106f, 0.707106f));
|
|
SetAtmosphericConditions(Vec3(20.0f, 20.0f, 20.0f), 0.001f, 0.00025f, -0.99f);
|
|
ILog* pLog(C3DEngine::GetLog());
|
|
if (false == LoadOpticalLUTs())
|
|
{
|
|
if (0 != pLog)
|
|
{
|
|
PrintMessage("Sky light: Optical lookup tables couldn't be loaded off disc. Recomputation needed!");
|
|
}
|
|
ComputeOpticalLUTs();
|
|
}
|
|
else
|
|
{
|
|
if (0 != pLog)
|
|
{
|
|
PrintMessage("Sky light: Optical lookup tables loaded off disc.");
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
CSkyLightNishita::~CSkyLightNishita()
|
|
{
|
|
}
|
|
|
|
|
|
CSkyLightNishita::SOpticalDepthLUTEntry CSkyLightNishita::LookupBilerpedOpticalDepthLUTEntry(
|
|
const SOpticalDepthLUTEntry* const __restrict cpOptDepthLUT,
|
|
uint32 heightIndex, const f32 cosVertAngle) const
|
|
{
|
|
uint32 vertAngleIndex;
|
|
f32 vertAngleIndexFrc;
|
|
f32 saveCosVertAngle(clamp_tpl(cosVertAngle, -1.0f, 1.0f));
|
|
f32 _index((f32) (cOLUT_AngularSteps - 1) * (-saveCosVertAngle * 0.5f + 0.5f));
|
|
vertAngleIndex = (uint32) _index;
|
|
vertAngleIndexFrc = _index - floorf(_index);
|
|
|
|
if (vertAngleIndex >= cOLUT_AngularSteps - 1)
|
|
{
|
|
return(cpOptDepthLUT[ OpticalLUTIndex(heightIndex, vertAngleIndex) ]);
|
|
}
|
|
else
|
|
{
|
|
uint32 index(OpticalLUTIndex(heightIndex, vertAngleIndex));
|
|
const SOpticalDepthLUTEntry& a(cpOptDepthLUT[ index ]);
|
|
const SOpticalDepthLUTEntry& b(cpOptDepthLUT[ index + 1 ]);
|
|
|
|
SOpticalDepthLUTEntry res;
|
|
res.mie = a.mie + vertAngleIndexFrc * (b.mie - a.mie);
|
|
res.rayleigh = a.rayleigh + vertAngleIndexFrc * (b.rayleigh - a.rayleigh);
|
|
return(res);
|
|
}
|
|
}
|
|
|
|
CSkyLightNishita::SPhaseLUTEntry CSkyLightNishita::LookupBilerpedPhaseLUTEntry(const f32 cosPhaseAngle) const
|
|
{
|
|
uint32 index;
|
|
f32 indexFrc;
|
|
MapCosPhaseAngleToIndex(cosPhaseAngle, index, indexFrc);
|
|
|
|
if (index >= cPLUT_AngularSteps - 1)
|
|
{
|
|
return(m_phaseLUT[ cPLUT_AngularSteps - 1 ]);
|
|
}
|
|
else
|
|
{
|
|
const SPhaseLUTEntry& a(m_phaseLUT[ index + 0 ]);
|
|
const SPhaseLUTEntry& b(m_phaseLUT[ index + 1 ]);
|
|
|
|
SPhaseLUTEntry res;
|
|
res.mie = a.mie + indexFrc * (b.mie - a.mie);
|
|
res.rayleigh = a.rayleigh + indexFrc * (b.rayleigh - a.rayleigh);
|
|
return(res);
|
|
}
|
|
}
|
|
|
|
void CSkyLightNishita::SamplePartialInScatteringAtHeight(const SOpticalScaleLUTEntry& osAtHeight,
|
|
const f32 outScatteringConstMie, const Vec3& outScatteringConstRayleigh, const SOpticalDepthLUTEntry& odAtHeightSky,
|
|
const SOpticalDepthLUTEntry& odAtViewerSky, const SOpticalDepthLUTEntry& odAtHeightSun,
|
|
Vec3& partialInScatteringMie, Vec3& partialInScatteringRayleigh) const
|
|
{
|
|
assert(odAtHeightSky.mie >= 0.0 && (odAtHeightSky.mie - floatDiffFactor) <= odAtViewerSky.mie);
|
|
assert(odAtHeightSun.mie >= 0.0);
|
|
assert(odAtHeightSky.rayleigh >= 0.0 && (odAtHeightSky.rayleigh - floatDiffFactor) <= odAtViewerSky.rayleigh);
|
|
assert(odAtHeightSun.rayleigh >= 0.0);
|
|
|
|
// mie out-scattering
|
|
f32 sampleExpArgMie(outScatteringConstMie * (-odAtHeightSun.mie - (odAtViewerSky.mie - odAtHeightSky.mie)));
|
|
|
|
// rayleigh out-scattering
|
|
Vec3 sampleExpArgRayleigh(outScatteringConstRayleigh * (-odAtHeightSun.rayleigh - (odAtViewerSky.rayleigh - odAtHeightSky.rayleigh)));
|
|
|
|
// partial in-scattering sampling result
|
|
Vec3 sampleExpArg(Vec3(sampleExpArgMie, sampleExpArgMie, sampleExpArgMie) + sampleExpArgRayleigh);
|
|
Vec3 sampleRes((float)exp_fast(sampleExpArg.x), (float)exp_fast(sampleExpArg.y), (float)exp_fast(sampleExpArg.z));
|
|
|
|
partialInScatteringMie = osAtHeight.mie * sampleRes;
|
|
partialInScatteringRayleigh = osAtHeight.rayleigh * sampleRes;
|
|
}
|
|
|
|
void CSkyLightNishita::ComputeInScatteringNoPremul(const f32 outScatteringConstMie, const Vec3& outScatteringConstRayleigh, const Vec3& skyDir,
|
|
Vec3& inScatteringMieNoPremul, Vec3& inScatteringRayleighNoPremul) const
|
|
{
|
|
// start integration along the "skyDir" from the viewer's point of view
|
|
const Vec3 c_up(0.0f, 0.0f, 1.0f);
|
|
const Vec3 viewer(c_up * c_earthRadiusf);
|
|
Vec3 curRayPos(viewer);
|
|
|
|
// to be reused by ray-sphere intersection code in loop below
|
|
f32 B(2.0f * viewer.Dot(skyDir));
|
|
f32 Bsq(B * B);
|
|
f32 Cpart(viewer.Dot(viewer));
|
|
|
|
// calculate optical depth at viewer
|
|
const SOpticalDepthLUTEntry* const __restrict cpOptDepthLUT = &m_opticalDepthLUT[0];
|
|
|
|
const Vec3& cSunDir(m_sunDir);
|
|
|
|
SOpticalDepthLUTEntry odAtViewerSky(LookupBilerpedOpticalDepthLUTEntry(cpOptDepthLUT, 0, skyDir.Dot(c_up)));
|
|
SOpticalDepthLUTEntry odAtViewerSun(LookupBilerpedOpticalDepthLUTEntry(cpOptDepthLUT, 0, cSunDir.Dot(c_up)));
|
|
|
|
// sample partial in-scattering term at viewer
|
|
Vec3 curSampleMie, curSampleRayleigh;
|
|
|
|
const SOpticalScaleLUTEntry* const __restrict cpOptScaleLUT = &m_opticalScaleLUT[0];
|
|
|
|
SamplePartialInScatteringAtHeight(cpOptScaleLUT[0], outScatteringConstMie, outScatteringConstRayleigh,
|
|
odAtViewerSky, odAtViewerSky, odAtViewerSun, curSampleMie, curSampleRayleigh);
|
|
|
|
// integrate along "skyDir" over all height segments we've precalculated in the optical lookup table
|
|
inScatteringMieNoPremul = Vec3(0.0f, 0.0f, 0.0f);
|
|
inScatteringRayleighNoPremul = Vec3(0.0f, 0.0f, 0.0f);
|
|
const int32 cInScatteringStepSize(m_inScatteringStepSize);
|
|
for (int a(1); a < cOLUT_HeightSteps; a += cInScatteringStepSize)
|
|
{
|
|
// calculate intersection with current "atmosphere shell"
|
|
const SOpticalScaleLUTEntry& crOpticalScaleLUTEntry = cpOptScaleLUT[a];
|
|
SOpticalScaleLUTEntry osAtHeight(crOpticalScaleLUTEntry);
|
|
|
|
f32 C(Cpart - (c_earthRadiusf + osAtHeight.atmosphereLayerHeight) * (c_earthRadiusf + osAtHeight.atmosphereLayerHeight));
|
|
f32 det(Bsq - 4.0f * C);
|
|
assert(det >= 0.0f && (0.5f * (-B - sqrtf(det)) <= 0.0f) && (((int)(0.5f * (-B + sqrtf(det)))*100.0) / 100.0f >= 0.0f));
|
|
|
|
f32 t(0.5f * (-B + sqrtf(det)));
|
|
|
|
Vec3 newRayPos(viewer + t * skyDir);
|
|
|
|
// calculate optical depth at new position
|
|
// since atmosphere bends we need to determine a new up vector to properly index the optical LUT
|
|
Vec3 newUp(newRayPos.GetNormalized());
|
|
SOpticalDepthLUTEntry odAtHeightSky(LookupBilerpedOpticalDepthLUTEntry(cpOptDepthLUT, a, skyDir.Dot(newUp)));
|
|
SOpticalDepthLUTEntry odAtHeightSun(LookupBilerpedOpticalDepthLUTEntry(cpOptDepthLUT, a, cSunDir.Dot(newUp)));
|
|
|
|
// when optimized in clang, values seem to drift a bit and under certain edge conditions
|
|
// raise asserts in SamplePartialInScatteringAtHeight function
|
|
if (odAtHeightSky.mie > odAtViewerSky.mie)
|
|
{
|
|
odAtHeightSky.mie = odAtViewerSky.mie;
|
|
}
|
|
if (odAtHeightSky.rayleigh > odAtViewerSky.rayleigh)
|
|
{
|
|
odAtHeightSky.rayleigh = odAtViewerSky.rayleigh;
|
|
}
|
|
|
|
// sample partial in-scattering term at new position
|
|
Vec3 newSampleMie, newSampleRayleigh;
|
|
SamplePartialInScatteringAtHeight(osAtHeight, outScatteringConstMie, outScatteringConstRayleigh,
|
|
odAtHeightSky, odAtViewerSky, odAtHeightSun, newSampleMie, newSampleRayleigh);
|
|
|
|
// integrate via trapezoid rule
|
|
f32 weight((newRayPos - curRayPos).GetLength() * 0.5f);
|
|
inScatteringMieNoPremul += (curSampleMie + newSampleMie) * weight;
|
|
inScatteringRayleighNoPremul += (curSampleRayleigh + newSampleRayleigh) * weight;
|
|
|
|
// update sampling data
|
|
curRayPos = newRayPos;
|
|
curSampleMie = newSampleMie;
|
|
curSampleRayleigh = newSampleRayleigh;
|
|
}
|
|
}
|
|
|
|
void CSkyLightNishita::ComputeSkyColor(const Vec3& skyDir, Vec3* pInScattering, Vec3* pInScatteringMieNoPremul,
|
|
Vec3* pInScatteringRayleighNoPremul, Vec3* pInScatteringRayleigh) const
|
|
{
|
|
//// get high precision normalized sky direction
|
|
//Vec3 _skyDir( skyDir );
|
|
//assert( _skyDir.GetLengthSquared() > 0.0 );
|
|
//_skyDir.Normalize();
|
|
|
|
assert(fabsf(skyDir.GetLengthSquared() - 1.0f) < 1e-4f);
|
|
|
|
SPhaseLUTEntry phaseLUTEntry(LookupBilerpedPhaseLUTEntry(-skyDir.Dot(m_sunDir)));
|
|
|
|
// initialize constants for mie scattering
|
|
f32 phaseForPhiGMie(phaseLUTEntry.mie);
|
|
f32 outScatteringConstMie(4.0f * c_pif * m_Km);
|
|
Vec3 inScatteringConstMie(m_sunIntensity * m_Km * phaseForPhiGMie);
|
|
|
|
// initialize constants for rayleigh scattering
|
|
f32 phaseForPhiGRayleigh(phaseLUTEntry.rayleigh);
|
|
Vec3 outScatteringConstRayleigh(4.0f * (float)c_pi * m_Kr * m_invRGBWaveLength4);
|
|
Vec3 inScatteringConstRayleigh((m_sunIntensity * m_Kr * phaseForPhiGRayleigh).CompMul(m_invRGBWaveLength4));
|
|
|
|
// compute in-scattering
|
|
Vec3 inScatteringMieNoPremul, inScatteringRayleighNoPremul;
|
|
|
|
ComputeInScatteringNoPremul(outScatteringConstMie, outScatteringConstRayleigh, skyDir, inScatteringMieNoPremul, inScatteringRayleighNoPremul);
|
|
|
|
assert(inScatteringMieNoPremul.x >= 0.0f && inScatteringMieNoPremul.y >= 0.0f && inScatteringMieNoPremul.z >= 0.0f);
|
|
assert(inScatteringRayleighNoPremul.x >= 0.0f && inScatteringRayleighNoPremul.y >= 0.0f && inScatteringRayleighNoPremul.z >= 0.0f);
|
|
|
|
// return color
|
|
if (pInScattering)
|
|
{
|
|
*pInScattering = Vec3(inScatteringMieNoPremul.CompMul(inScatteringConstMie) + inScatteringRayleighNoPremul.CompMul(inScatteringConstRayleigh));
|
|
}
|
|
|
|
if (pInScatteringMieNoPremul)
|
|
{
|
|
*pInScatteringMieNoPremul = Vec3(inScatteringMieNoPremul);
|
|
}
|
|
|
|
if (pInScatteringRayleighNoPremul)
|
|
{
|
|
*pInScatteringRayleighNoPremul = Vec3(inScatteringRayleighNoPremul);
|
|
}
|
|
|
|
if (pInScatteringRayleigh)
|
|
{
|
|
*pInScatteringRayleigh = Vec3(inScatteringRayleighNoPremul.CompMul(inScatteringConstRayleigh));
|
|
}
|
|
}
|
|
|
|
void CSkyLightNishita::SetInScatteringIntegralStepSize(int32 stepSize)
|
|
{
|
|
stepSize = stepSize < 1 ? 1 : stepSize > 2 ? 2 : stepSize;
|
|
m_inScatteringStepSize = stepSize;
|
|
}
|
|
|
|
int32 CSkyLightNishita::GetInScatteringIntegralStepSize() const
|
|
{
|
|
return(m_inScatteringStepSize);
|
|
}
|
|
|
|
Vec4 CSkyLightNishita::GetPartialMieInScatteringConst() const
|
|
{
|
|
Vec3 res(m_sunIntensity * m_Km);
|
|
return(Vec4(res.x, res.y, res.z, 0.0f));
|
|
}
|
|
|
|
|
|
Vec4 CSkyLightNishita::GetPartialRayleighInScatteringConst() const
|
|
{
|
|
Vec3 res((m_sunIntensity * m_Kr).CompMul(m_invRGBWaveLength4));
|
|
return(Vec4(res.x, res.y, res.z, 0.0f));
|
|
}
|
|
|
|
|
|
Vec3 CSkyLightNishita::GetSunDirection() const
|
|
{
|
|
return(Vec3(m_sunDir.x, m_sunDir.y, m_sunDir.z));
|
|
}
|
|
|
|
|
|
Vec4 CSkyLightNishita::GetPhaseFunctionConsts() const
|
|
{
|
|
//f32 g2( m_g * m_g );
|
|
//f32 miePart( 1.5f * ( 1.0f - g2 ) / ( 2.0f + g2 ) );
|
|
//return( Vec4( m_g, m_g * m_g, miePart, 0.0f ) );
|
|
|
|
f32 g2(m_g * m_g);
|
|
f32 miePart(1.5f * (1.0f - g2) / (2.0f + g2));
|
|
f32 miePartPow(powf(miePart, -2.0f / 3.0f));
|
|
return(Vec4(miePartPow * -2.0f * m_g, miePartPow * (1.0f + g2), 0.0f, 0.0f));
|
|
}
|
|
|
|
|
|
f64 CSkyLightNishita::IntegrateOpticalDepth(const Vec3d& start, const Vec3d& end, const f64& avgDensityHeightInv, const f64& error) const
|
|
{
|
|
f64 startScale(OpticalScaleFunction(start.GetLength() - c_earthRadius, avgDensityHeightInv));
|
|
f64 endScale(OpticalScaleFunction(end.GetLength() - c_earthRadius, avgDensityHeightInv));
|
|
return(IntegrateOpticalDepthInternal(start, startScale, end, endScale, avgDensityHeightInv, error));
|
|
}
|
|
|
|
|
|
bool CSkyLightNishita::ComputeOpticalDepth(const Vec3d& cameraLookDir, const f64& cameraHeight, const f64& avgDensityHeightInv, float& depth) const
|
|
{
|
|
// init camera position
|
|
Vec3d cameraPos(0.0, cameraHeight + c_earthRadius, 0.0);
|
|
|
|
// check if ray hits earth
|
|
// compute B, and C of quadratic function (A=1, as looking direction is normalized)
|
|
f64 B(2.0 * cameraPos.Dot(cameraLookDir));
|
|
f64 Bsq(B * B);
|
|
f64 Cpart(cameraPos.Dot(cameraPos));
|
|
f64 C(Cpart - c_earthRadius * c_earthRadius);
|
|
f64 det(Bsq - 4.0 * C);
|
|
|
|
bool hitsEarth(det >= 0.0 && ((0.5 * (-B - sqrt(det)) > 1e-4) || (0.5 * (-B + sqrt(det)) > 1e-4)));
|
|
if (false != hitsEarth)
|
|
{
|
|
depth = (float)c_opticalDepthWhenHittingEarth;
|
|
return(false);
|
|
}
|
|
|
|
// find intersection with atmosphere top
|
|
C = Cpart - (c_maxAtmosphereHeight + c_earthRadius) * (c_maxAtmosphereHeight + c_earthRadius);
|
|
det = Bsq - 4.0 * C;
|
|
assert(det >= 0.0); // ray defined outside the atmosphere
|
|
f64 t(0.5 * (-B + sqrt(det)));
|
|
assert(t >= -1e-4);
|
|
if (t < 0.0)
|
|
{
|
|
t = 0.0;
|
|
}
|
|
|
|
// integrate depth along ray from camera to atmosphere top
|
|
f64 _depth(0.0);
|
|
|
|
int numInitialSamples((int) t);
|
|
numInitialSamples = (numInitialSamples < 2) ? 2 : numInitialSamples;
|
|
|
|
Vec3d lastCameraPos(cameraPos);
|
|
for (int i(1); i < numInitialSamples; ++i)
|
|
{
|
|
Vec3d curCameraPos(cameraPos + cameraLookDir * (t * ((float) i / (float) numInitialSamples)));
|
|
_depth += IntegrateOpticalDepth(lastCameraPos, curCameraPos, avgDensityHeightInv, 1e-1);
|
|
lastCameraPos = curCameraPos;
|
|
}
|
|
|
|
assert(_depth >= 0.0 && _depth < 1e25);
|
|
assert(0 != _finite(_depth));
|
|
|
|
depth = (float) _depth;
|
|
return(true);
|
|
}
|
|
|
|
|
|
void CSkyLightNishita::ComputeOpticalLUTs()
|
|
{
|
|
LOADING_TIME_PROFILE_SECTION(GetISystem());
|
|
|
|
ILog* pLog(C3DEngine::GetLog());
|
|
if (0 != pLog)
|
|
{
|
|
PrintMessage("Sky light: Computing optical lookup tables (this might take a while)... ");
|
|
}
|
|
|
|
// reset tables
|
|
m_opticalDepthLUT.resize(0);
|
|
m_opticalDepthLUT.reserve(cOLUT_HeightSteps * cOLUT_AngularSteps);
|
|
|
|
m_opticalScaleLUT.resize(0);
|
|
m_opticalScaleLUT.reserve(cOLUT_HeightSteps);
|
|
|
|
// compute LUTs
|
|
for (int a(0); a < cOLUT_HeightSteps; ++a)
|
|
{
|
|
f64 height(MapIndexToHeight(a));
|
|
|
|
// compute optical depth
|
|
for (int i(0); i < cOLUT_AngularSteps; ++i)
|
|
{
|
|
// init looking direction of camera
|
|
f64 cosVertAngle(MapIndexToCosVertAngle(i));
|
|
Vec3d cameraLookDir(sqrt(1.0 - cosVertAngle * cosVertAngle), cosVertAngle, 0.0);
|
|
|
|
// compute optical depth
|
|
SOpticalDepthLUTEntry e;
|
|
bool b0(ComputeOpticalDepth(cameraLookDir, height, c_avgDensityHeightMieInv, e.mie));
|
|
bool b1(ComputeOpticalDepth(cameraLookDir, height, c_avgDensityHeightRayleighInv, e.rayleigh));
|
|
assert(b0 == b1);
|
|
|
|
// blend out previous values once camera ray hits earth
|
|
if (false == b0 && false == b1 && i > 0)
|
|
{
|
|
e = m_opticalDepthLUT.back();
|
|
e.mie = (f32) ((f32)0.5 * (e.mie + c_opticalDepthWhenHittingEarth));
|
|
e.rayleigh = (f32) ((f32)0.5 * (e.rayleigh + c_opticalDepthWhenHittingEarth));
|
|
}
|
|
|
|
// store result
|
|
m_opticalDepthLUT.push_back(e);
|
|
}
|
|
|
|
{
|
|
// compute optical scale
|
|
SOpticalScaleLUTEntry e;
|
|
e.atmosphereLayerHeight = (f32) height;
|
|
e.mie = (f32) OpticalScaleFunction(height, c_avgDensityHeightMieInv);
|
|
e.rayleigh = (f32) OpticalScaleFunction(height, c_avgDensityHeightRayleighInv);
|
|
m_opticalScaleLUT.push_back(e);
|
|
}
|
|
}
|
|
|
|
// save LUTs for next time
|
|
SaveOpticalLUTs();
|
|
if (0 != pLog)
|
|
{
|
|
PrintMessage(" ... done.\n");
|
|
}
|
|
}
|
|
|
|
|
|
void CSkyLightNishita::ComputePhaseLUT()
|
|
{
|
|
//ILog* pLog( C3DEngine::GetLog() );
|
|
//if( 0 != pLog )
|
|
// PrintMessage( "Sky light: Computing phase lookup table... " );
|
|
|
|
// reset tables
|
|
m_phaseLUT.resize(0);
|
|
m_phaseLUT.reserve(cPLUT_AngularSteps);
|
|
|
|
// compute coefficients
|
|
f32 g(m_g);
|
|
f32 g2(g * g);
|
|
f32 miePart(1.5f * (1.0f - g2) / (2.0f + g2));
|
|
|
|
// calculate entries
|
|
for (int i(0); i < cPLUT_AngularSteps; ++i)
|
|
{
|
|
f32 cosine(MapIndexToCosPhaseAngle(i));
|
|
f32 cosine2(cosine * cosine);
|
|
|
|
//f32 t = 1.0f + g2 - 2.0f * g * cosine;
|
|
//if (fabsf(t) < 1e-5f)
|
|
//{
|
|
// PrintMessage( "Sky light: g = %.10f", g );
|
|
// PrintMessage( "Sky light: g2 = %.10f", g2 );
|
|
// PrintMessage( "Sky light: cosine = %.10f", cosine );
|
|
// PrintMessage( "Sky light: cosine2 = %.10f", cosine2 );
|
|
// PrintMessage( "Sky light: t = %.10f", t );
|
|
//}
|
|
|
|
f32 miePhase(miePart * (1.0f + cosine2) / powf(1.0f + g2 - 2.0f * g * cosine, 1.5f));
|
|
f32 rayleighPhase(0.75f * (1.0f + cosine2));
|
|
|
|
SPhaseLUTEntry e;
|
|
e.mie = (float) miePhase;
|
|
e.rayleigh = (float) rayleighPhase;
|
|
m_phaseLUT.push_back(e);
|
|
}
|
|
//if( 0 != pLog )
|
|
// PrintMessage( " ... done.\n" );
|
|
}
|
|
|
|
|
|
f64 CSkyLightNishita::MapIndexToHeight(uint32 index) const
|
|
{
|
|
// a function that maps well to mie and rayleigh at the same time
|
|
// that is, a lot of indices will map below the average density height for mie & rayleigh scattering
|
|
assert(index < cOLUT_HeightSteps);
|
|
f64 x((f64)index / (cOLUT_HeightSteps - 1));
|
|
return(c_maxAtmosphereHeight * exp_precise(10.0 * (x - 1.0)) * x);
|
|
}
|
|
|
|
|
|
f64 CSkyLightNishita::MapIndexToCosVertAngle(uint32 index) const
|
|
{
|
|
assert(index < cOLUT_AngularSteps);
|
|
return(1.0 - 2.0 * ((f64)index / (cOLUT_AngularSteps - 1)));
|
|
}
|
|
|
|
|
|
f32 CSkyLightNishita::MapIndexToCosPhaseAngle(uint32 index) const
|
|
{
|
|
assert(index < cPLUT_AngularSteps);
|
|
return(1.0f - 2.0f * ((f32)index / ((f32)cPLUT_AngularSteps - 1)));
|
|
}
|
|
|
|
|
|
void CSkyLightNishita::MapCosPhaseAngleToIndex(const f32 cosPhaseAngle, uint32& index, f32& indexFrc) const
|
|
{
|
|
//assert( -1 <= cosPhaseAngle && 1 >= cosPhaseAngle );
|
|
f32 saveCosPhaseAngle(clamp_tpl(cosPhaseAngle, -1.0f, 1.0f));
|
|
f32 _index((f32) (cPLUT_AngularSteps - 1) * (-saveCosPhaseAngle * 0.5f + 0.5f));
|
|
index = (uint32) _index;
|
|
indexFrc = _index - floorf(_index);
|
|
}
|
|
|
|
|
|
uint32 CSkyLightNishita::OpticalLUTIndex(uint32 heightIndex, uint32 cosVertAngleIndex) const
|
|
{
|
|
assert(heightIndex < cOLUT_HeightSteps && cosVertAngleIndex < cOLUT_AngularSteps);
|
|
return(heightIndex * cOLUT_AngularSteps + cosVertAngleIndex);
|
|
}
|
|
|
|
|
|
bool CSkyLightNishita::LoadOpticalLUTs()
|
|
{
|
|
auto pPak(C3DEngine::GetPak());
|
|
if (0 != pPak)
|
|
{
|
|
AZ::IO::HandleType fileHandle = pPak->FOpen(c_lutFileName, "rb");
|
|
if (fileHandle != AZ::IO::InvalidHandle)
|
|
{
|
|
size_t itemsRead(0);
|
|
|
|
// read in file tag
|
|
uint32 fileTag(0);
|
|
itemsRead = pPak->FRead(&fileTag, 1, fileHandle);
|
|
if (itemsRead != 1 || fileTag != c_lutFileTag)
|
|
{
|
|
// file tag mismatch
|
|
pPak->FClose(fileHandle);
|
|
return(false);
|
|
}
|
|
|
|
// read in file format version
|
|
uint32 fileVersion(0);
|
|
itemsRead = pPak->FRead(&fileVersion, 1, fileHandle);
|
|
if (itemsRead != 1 || fileVersion != c_lutFileVersion)
|
|
{
|
|
// file version mismatch
|
|
pPak->FClose(fileHandle);
|
|
return(false);
|
|
}
|
|
|
|
// read in optical depth LUT
|
|
m_opticalDepthLUT.resize(cOLUT_HeightSteps * cOLUT_AngularSteps);
|
|
itemsRead = pPak->FRead(&m_opticalDepthLUT[0], m_opticalDepthLUT.size(), fileHandle);
|
|
if (itemsRead != m_opticalDepthLUT.size())
|
|
{
|
|
pPak->FClose(fileHandle);
|
|
return(false);
|
|
}
|
|
|
|
// read in optical scale LUT
|
|
m_opticalScaleLUT.resize(cOLUT_HeightSteps);
|
|
itemsRead = pPak->FRead(&m_opticalScaleLUT[0], m_opticalScaleLUT.size(), fileHandle);
|
|
if (itemsRead != m_opticalScaleLUT.size())
|
|
{
|
|
pPak->FClose(fileHandle);
|
|
return(false);
|
|
}
|
|
|
|
// check if we read entire file
|
|
uint64_t curPos(pPak->FTell(fileHandle));
|
|
pPak->FSeek(fileHandle, 0, SEEK_END);
|
|
uint64_t endPos(pPak->FTell(fileHandle));
|
|
if (curPos != endPos)
|
|
{
|
|
pPak->FClose(fileHandle);
|
|
return(false);
|
|
}
|
|
|
|
// LUT successfully read
|
|
pPak->FClose(fileHandle);
|
|
return(true);
|
|
}
|
|
}
|
|
|
|
return(false);
|
|
}
|
|
|
|
|
|
void CSkyLightNishita::SaveOpticalLUTs() const
|
|
{
|
|
// only save on little endian PCs so the load function can do proper endian swapping
|
|
#if defined(_CPU_X86) || defined(_CPU_AMD64)
|
|
auto pPak(C3DEngine::GetPak());
|
|
if (0 != pPak)
|
|
{
|
|
AZ::IO::HandleType fileHandle = pPak->FOpen(c_lutFileName, "wb");
|
|
if (fileHandle != AZ::IO::InvalidHandle)
|
|
{
|
|
// write out file tag
|
|
pPak->FWrite(&c_lutFileTag, 1, sizeof(c_lutFileTag), fileHandle);
|
|
|
|
// write out file format version
|
|
pPak->FWrite(&c_lutFileVersion, 1, sizeof(c_lutFileVersion), fileHandle);
|
|
|
|
// write out optical depth LUT
|
|
assert(m_opticalDepthLUT.size() == cOLUT_HeightSteps * cOLUT_AngularSteps);
|
|
pPak->FWrite(&m_opticalDepthLUT[0], 1, sizeof(SOpticalDepthLUTEntry) * m_opticalDepthLUT.size(), fileHandle);
|
|
|
|
// write out optical scale LUT
|
|
assert(m_opticalScaleLUT.size() == cOLUT_HeightSteps);
|
|
pPak->FWrite(&m_opticalScaleLUT[0], 1, sizeof(SOpticalScaleLUTEntry) * m_opticalScaleLUT.size(), fileHandle);
|
|
|
|
// close file
|
|
pPak->FClose(fileHandle);
|
|
}
|
|
}
|
|
#endif
|
|
}
|