float OrthoCamera::GenerateRay(const Sample &sample,
                               Ray *ray) const {
	// Generate raster and camera samples
	Point Pras(sample.imageX, sample.imageY, 0);
	Point Pcamera;
	RasterToCamera(Pras, &Pcamera);
	ray->o = Pcamera;
	ray->d = Vector(0,0,1);
	// Set ray time value
	ray->time = Lerp(sample.time, ShutterOpen, ShutterClose);
	// Modify ray for depth of field
	if (LensRadius > 0.) {
		// Sample point on lens
		float lensU, lensV;
		ConcentricSampleDisk(sample.lensU, sample.lensV,
		                     &lensU, &lensV);
		lensU *= LensRadius;
		lensV *= LensRadius;
		// Compute point on plane of focus
		float ft = (FocalDistance - ClipHither) / ray->d.z;
		Point Pfocus = (*ray)(ft);
		// Update ray for effect of lens
		ray->o.x += lensU * (FocalDistance - ClipHither) / FocalDistance;
		ray->o.y += lensV * (FocalDistance - ClipHither) / FocalDistance;
		ray->d = Pfocus - ray->o;
	}
	ray->mint = 0.;
	ray->maxt = ClipYon - ClipHither;
	ray->d = Normalize(ray->d);
	CameraToWorld(*ray, ray);
	return 1.f;
}

std::vector<float> ComputeRadialWeights(int rsteps, float minRadius, float maxRadius)
{
	std::vector<float> wnd(rsteps);
	for(int x=0;x<rsteps;x++)
		wnd[x]=Lerp(minRadius, maxRadius, x/(float)rsteps) / (0.5f * (minRadius+maxRadius));
	return wnd;
}

Float Turbulence(const Point3f &p, const Vector3f &dpdx, const Vector3f &dpdy,
                 Float omega, int maxOctaves) {
    // Compute number of octaves for antialiased FBm
    Float len2 = std::max(dpdx.LengthSquared(), dpdy.LengthSquared());
    Float n = Clamp(-1 - .5f * Log2(len2), 0, maxOctaves);
    int nInt = std::floor(n);

    // Compute sum of octaves of noise for turbulence
    Float sum = 0, lambda = 1, o = 1;
    for (int i = 0; i < nInt; ++i) {
        sum += o * std::abs(Noise(lambda * p));
        lambda *= 1.99f;
        o *= omega;
    }

    // Account for contributions of clamped octaves in turbulence
    Float nPartial = n - nInt;
    sum += o * Lerp(SmoothStep(.3f, .7f, nPartial), 0.2,
                    std::abs(Noise(lambda * p)));
    for (int i = nInt; i < maxOctaves; ++i) {
        sum += o * 0.2f;
        o *= omega;
    }
    return sum;
}

template <typename T> void ork::TVector2<T>::Serp( const TVector2<T> & PA, const TVector2<T> & PB, const TVector2<T> & PC, const TVector2<T> & PD, T Par )
{
	TVector2<T> PAB, PCD;
	PAB.Lerp( PA, PB, Par );
	PCD.Lerp( PC, PD, Par );
	Lerp( PAB, PCD, Par );
}

Quaternion Quaternion::Slerp(const Quaternion &a, const Quaternion &b, float t)
{
    Quaternion c;

    float dot = Dot(a, b);

    if(dot < 0)
    {
        dot = -dot;
        c = -b;
    }
    else
    {
        c = b;
    }

    if(dot < 0.95f)
    {
        float angle = acosf(dot);

        Quaternion d = a * sin(angle * (1 - t));
        d += c * sin(angle * t);
        d /= sin(angle);

        return d;
        //return (a * sinf(angle * (1 - t)) + c * sinf(angle * t)) / sinf(angle);
        // TODO: Figure out what is wrong with the + operator.
    }
    else
    {
        return Lerp(a, c, t);
    }
}

//-----------------------------------------------------------------------------
// Purpose: Fades lookup weight from CurWeight->0.0 
//-----------------------------------------------------------------------------
void CColorCorrection::FadeOutThink( void )
{
	// Check for conditions where we shouldn't fade out
	if ( m_flFadeOutDuration <= 0 || // not set to fade out
			m_flCurWeight <= 0.0f || // already faded out
					   m_bEnabled || // fade in/out mutex
		   m_flMaxWeight == 0.0f  || // min==max
	 m_flStartFadeOutWeight <= 0.0f )// already at min weight
	{
		SetNextThink ( TICK_NEVER_THINK, s_pFadeOutContextThink );
		return;
	}

	// If we started fading out without fully fading in, use a truncated duration
    float flTimeToFade = m_flFadeOutDuration;
	if ( m_flStartFadeOutWeight < m_flMaxWeight )
	{	
		float flWeightRatio		= m_flStartFadeOutWeight / m_flMaxWeight;
		flWeightRatio = clamp ( flWeightRatio, 0.01f, 1.0f );
		flTimeToFade			= m_flFadeOutDuration * flWeightRatio;
	}	
	
	Assert ( flTimeToFade > 0.0f );
	float flFadeRatio = (gpGlobals->curtime - m_flTimeStartFadeOut) / flTimeToFade;
	flFadeRatio = clamp ( flFadeRatio, 0.0f, 1.0f );
	m_flStartFadeOutWeight = clamp ( m_flStartFadeOutWeight, 0.0f, 1.0f );

	m_flCurWeight = Lerp( 1.0f - flFadeRatio, 0.0f, m_flStartFadeOutWeight );

	SetNextThink( gpGlobals->curtime + COLOR_CORRECTION_ENT_THINK_RATE, s_pFadeOutContextThink );
}

// HighContrastOp Method Definitions
void HighContrastOp::Map(const float *y, int xRes, int yRes,
		float maxDisplayY, float *scale) const {
	// Find minimum and maximum image luminances
	float minY = y[0], maxY = y[0];
	for (int i = 0; i < xRes * yRes; ++i) {
		minY = min(minY, y[i]);
		maxY = max(maxY, y[i]);
	}
	float CYmin = C(minY), CYmax = C(maxY);
	// Build luminance image pyramid
	MIPMap<float> pyramid(xRes, yRes,
	                      y, false,
						  4.f, TEXTURE_CLAMP);
	// Apply high contrast tone mapping operator
	ProgressReporter progress(xRes*yRes, "Tone Mapping"); // NOBOOK
	for (int y = 0; y < yRes; ++y) {
		float yc = (float(y) + .5f) / float(yRes);
		for (int x = 0; x < xRes; ++x) {
			float xc = (float(x) + .5f) / float(xRes);
			// Compute local adaptation luminance at $(x,y)$
			float dwidth = 1.f / float(max(xRes, yRes));
			float maxWidth = 32.f / float(max(xRes, yRes));
			float width = dwidth, prevWidth = 0.f;
			float Yadapt;
			float prevlc = 0.f;
			const float maxLocalContrast = .5f;
			while (1) {
				// Compute local contrast at $(x,y)$
				float b0 = pyramid.Lookup(xc, yc, width,
				                          0.f, 0.f, width);
				float b1 = pyramid.Lookup(xc, yc, 2.f*width,
				                          0.f, 0.f, 2.f*width);
				float lc = fabsf((b0 - b1) / b0);
				// If maximum contrast is exceeded, compute adaptation luminance
				if (lc > maxLocalContrast) {
					float t = (maxLocalContrast - prevlc) / (lc - prevlc);
					float w = Lerp(t, prevWidth, width);
					Yadapt = pyramid.Lookup(xc, yc, w,
					                        0.f, 0.f, w);
					break;
				}
				// Increase search region and prepare to compute contrast again
				prevlc = lc;
				prevWidth = width;
				width += dwidth;
				if (width >= maxWidth) {
					Yadapt = pyramid.Lookup(xc, yc, maxWidth,
					                        0.f, 0.f, maxWidth);
					break;
				}
			}
			// Apply tone mapping based on local adaptation luminance
			scale[x + y*xRes] = T(Yadapt, CYmin, CYmax, maxDisplayY) /
				Yadapt;
		}
		progress.Update(xRes); // NOBOOK
	}
	progress.Done(); // NOBOOK
}

void Interpolator::Update(uint32 frame_time) {
	if (_ValidMethod() == false) {
		if (VIDEO_DEBUG)
			cerr << "VIDEO WARNING: " << __FUNCTION__ << " was called when an invalid method was set" << endl;
		return;
	}

	// update current time
	_current_time += frame_time;

	if (_current_time > _end_time) {
		_current_time = _end_time;
		_finished = true;
	}

	// Calculate a value from 0.0f to 1.0f that tells how far we are in the interpolation
	float progress;

	if (_end_time == 0) {
		progress = 1.0f;
	}
	else {
		progress = static_cast<float>(_current_time) / static_cast<float>(_end_time);
	}

	if (progress > 1.0f) {
		if (VIDEO_DEBUG)
			cerr << "VIDEO WARNING: " << __FUNCTION__ << " calculated a progress value greater than 1.0" << endl;
		progress = 1.0f;
	}

	// Apply a transformation based on the interpolation method
	switch(_method) {
		case VIDEO_INTERPOLATE_EASE:
			progress = _EaseTransform(progress);
			break;
		case VIDEO_INTERPOLATE_SRCA:
			progress = 0.0f;
			break;
		case VIDEO_INTERPOLATE_SRCB:
			progress = 1.0f;
			break;
		case VIDEO_INTERPOLATE_FAST:
			progress = _FastTransform(progress);
			break;
		case VIDEO_INTERPOLATE_SLOW:
			progress = _SlowTransform(progress);
			break;
		case VIDEO_INTERPOLATE_LINEAR:
			// Nothing to do, just use progress value as it is
			break;
		default:
			if (VIDEO_DEBUG)
				cerr << "VIDEO WARNING: " << __FUNCTION__ << " the current method did not match any supported methods" << endl;
			return;
	};

	_current_value = Lerp(progress, _a, _b);
} // void Interpolator::Update(uint32 frame_time)

//------------------------------------------------------------------------------
// Purpose:
//------------------------------------------------------------------------------
void CAI_Spotlight::UpdateSpotlightEndpoint( void )
{
	if ( !m_hSpotlight )
	{
		CreateSpotlightEntities();
	}

	Vector vecStartPoint, vecEndPoint;
	vecStartPoint = m_hSpotlight->GetAbsStartPos();
	ComputeEndpoint( vecStartPoint, &vecEndPoint );

	// If I'm not facing the spotlight turn it off 
	Vector vecSpotDir;
	VectorSubtract( vecEndPoint, vecStartPoint, vecSpotDir );
	float flBeamLength = VectorNormalize(vecSpotDir);
	
	m_hSpotlightTarget->SetAbsOrigin( vecEndPoint );
	m_hSpotlightTarget->SetAbsVelocity( vec3_origin );
	m_hSpotlightTarget->m_vSpotlightOrg = vecStartPoint;
	m_hSpotlightTarget->m_vSpotlightDir = vecSpotDir;

	// Avoid sudden change in where beam fades out when cross disconinuities
	m_flSpotlightCurLength = Lerp( 0.20f, m_flSpotlightCurLength, flBeamLength );

	// Fade out spotlight end if past max length.  
	if (m_flSpotlightCurLength > 2*m_flSpotlightMaxLength)
	{
		m_hSpotlightTarget->SetRenderColorA( 0 );
		m_hSpotlight->SetFadeLength(m_flSpotlightMaxLength);
	}
	else if (m_flSpotlightCurLength > m_flSpotlightMaxLength)		
	{
		m_hSpotlightTarget->SetRenderColorA( (1-((m_flSpotlightCurLength-m_flSpotlightMaxLength)/m_flSpotlightMaxLength)) );
		m_hSpotlight->SetFadeLength(m_flSpotlightMaxLength);
	}
	else
	{
		m_hSpotlightTarget->SetRenderColorA( 1.0 );
		m_hSpotlight->SetFadeLength(m_flSpotlightCurLength);
	}

	// Adjust end width to keep beam width constant
	float flNewWidth = SPOTLIGHT_WIDTH * ( flBeamLength / m_flSpotlightMaxLength );
	
	flNewWidth = min( 100, flNewWidth );

	m_hSpotlight->SetWidth(flNewWidth);
	m_hSpotlight->SetEndWidth(flNewWidth);

	// Adjust width of light on the end.  
	if ( FBitSet (m_nFlags, AI_SPOTLIGHT_NO_DLIGHTS) )
	{
		m_hSpotlightTarget->m_flLightScale = 0.0;
	}
	else
	{
		m_hSpotlightTarget->m_flLightScale = flNewWidth;
	}
}

float MapManager::GetHeightByPosition( float x, float z )
{
	if ( !m_HeightMap )
	{
		return 0.0f;
	}
	
	x /= m_PixelSize;
	z /= m_PixelSize;
	x = static_cast<float>(m_HeightMapWidth) / 2.0f + x;
	z = static_cast<float>(m_HeightMapHeight) / 2.0f + z;
	
	int col = static_cast<int>( std::floor( x ) );
	int row = static_cast<int>( std::floor( z ) );
	
	col = __min( col, m_HeightMapWidth - 1 );
	row = __min( row, m_HeightMapHeight - 1 );
	col = __max( 0, col );
	row = __max( 0, row );

	float leftBottom = GetHeightInMap( col, row + 1 );
	float rightBottom = GetHeightInMap( col + 1, row + 1 );
	float leftTop = GetHeightInMap( col , row );
	float rightTop = GetHeightInMap( col + 1, row );

	float dx = x - col;
	float dz = z - row;
	
	if ( dx < 0 )
	{
		dx = -dx;
	}
	if ( dz < 0 )
	{
		dz = -dz;
	}

	float heightBottom = Lerp( leftBottom, rightBottom, dx );
	float heightTop = Lerp( leftTop, rightTop, dx );
	float height = Lerp( heightTop, heightBottom, dz );
	
	// Log( "(%x, %x) %4f %4f %4f %4f %4f %4f %4f \n",
	//	 col, row, leftBottom, rightBottom, leftTop, rightTop, heightBottom, heightTop, height);

	return height;
}

Point Cylinder::Sample(float u1, float u2, Normal *Ns) const {
    float z = Lerp(u1, zmin, zmax);
    float t = u2 * phiMax;
    Point p = Point(radius * cosf(t), radius * sinf(t), z);
    *Ns = Normalize((*ObjectToWorld)(Normal(p.x, p.y, 0.)));
    if (ReverseOrientation) *Ns *= -1.f;
    return (*ObjectToWorld)(p);
}

	void TimedFadeAction::Update(const Context & context)
	{
		float32 dt = float32(context.mTimeManager->GetSeconds());
		m_CurrentSeconds += dt;
		m_pSpriteComponent->SetColorMultiplier(
			Lerp(m_StartColor, m_EndColor, m_CurrentSeconds / m_Seconds)
			);
	}

Vector UniformSampleCone(float u1, float u2, float costhetamax,
                         const Vector &x, const Vector &y, const Vector &z) {
    float costheta = Lerp(u1, costhetamax, 1.f);
    float sintheta = sqrtf(1.f - costheta*costheta);
    float phi = u2 * 2.f * M_PI;
    return cosf(phi) * sintheta * x + sinf(phi) * sintheta * y +
           costheta * z;
}

FloatingPointMatrix3x2<T>
FloatingPointMatrix3x2<T>::Lerp(const FloatingPointMatrix3x2& source1,
    const FloatingPointMatrix3x2& source2, T amount) noexcept
{
    FloatingPointMatrix3x2 result;
    Lerp(source1, source2, amount, result);
    return std::move(result);
}

static real64 PerlinNoise3DFunction(real64 x, real64 y, real64 z) 
{
	// Compute noise cell coordinates and offsets
	int32_t ix = Floor2Int(x);
	int32_t iy = Floor2Int(y);
	int32_t iz = Floor2Int(z);
	real64 dx = x - ix, dy = y - iy, dz = z - iz;
	// Compute gradient weights
	ix &= (NOISE_PERM_SIZE-1);
	iy &= (NOISE_PERM_SIZE-1);
	iz &= (NOISE_PERM_SIZE-1);
	real64 w000 = Grad3d(ix,   iy,   iz,   dx,   dy,   dz);
	real64 w100 = Grad3d(ix+1, iy,   iz,   dx-1, dy,   dz);
	real64 w010 = Grad3d(ix,   iy+1, iz,   dx,   dy-1, dz);
	real64 w110 = Grad3d(ix+1, iy+1, iz,   dx-1, dy-1, dz);
	real64 w001 = Grad3d(ix,   iy,   iz+1, dx,   dy,   dz-1);
	real64 w101 = Grad3d(ix+1, iy,   iz+1, dx-1, dy,   dz-1);
	real64 w011 = Grad3d(ix,   iy+1, iz+1, dx,   dy-1, dz-1);
	real64 w111 = Grad3d(ix+1, iy+1, iz+1, dx-1, dy-1, dz-1);
	// Compute trilinear interpolation of weights
	real64 wx = PerlinFade(dx);
	real64 wy = PerlinFade(dy);
	real64 wz = PerlinFade(dz);
	real64 x00 = Lerp(wx, w000, w100);
	real64 x10 = Lerp(wx, w010, w110);
	real64 x01 = Lerp(wx, w001, w101);
	real64 x11 = Lerp(wx, w011, w111);
	real64 y0 = Lerp(wy, x00, x10);
	real64 y1 = Lerp(wy, x01, x11);
	return Lerp(wz, y0, y1);
}

Float Noise(Float x, Float y, Float z) {
    // Compute noise cell coordinates and offsets
    int ix = std::floor(x), iy = std::floor(y), iz = std::floor(z);
    Float dx = x - ix, dy = y - iy, dz = z - iz;

    // Compute gradient weights
    ix &= NoisePermSize - 1;
    iy &= NoisePermSize - 1;
    iz &= NoisePermSize - 1;
    Float w000 = Grad(ix, iy, iz, dx, dy, dz);
    Float w100 = Grad(ix + 1, iy, iz, dx - 1, dy, dz);
    Float w010 = Grad(ix, iy + 1, iz, dx, dy - 1, dz);
    Float w110 = Grad(ix + 1, iy + 1, iz, dx - 1, dy - 1, dz);
    Float w001 = Grad(ix, iy, iz + 1, dx, dy, dz - 1);
    Float w101 = Grad(ix + 1, iy, iz + 1, dx - 1, dy, dz - 1);
    Float w011 = Grad(ix, iy + 1, iz + 1, dx, dy - 1, dz - 1);
    Float w111 = Grad(ix + 1, iy + 1, iz + 1, dx - 1, dy - 1, dz - 1);

    // Compute trilinear interpolation of weights
    Float wx = NoiseWeight(dx), wy = NoiseWeight(dy), wz = NoiseWeight(dz);
    Float x00 = Lerp(wx, w000, w100);
    Float x10 = Lerp(wx, w010, w110);
    Float x01 = Lerp(wx, w001, w101);
    Float x11 = Lerp(wx, w011, w111);
    Float y0 = Lerp(wy, x00, x10);
    Float y1 = Lerp(wy, x01, x11);
    return Lerp(wz, y0, y1);
}

void MtlBlinnSW3D::Transmit(ShadeContext3D &sc)
{
	BlinnBlock &block = (BlinnBlock &) sc.GetMaterialBlock(this);
	const Vector4f &vertexColor = IsSet(block.VertexColor) ? *block.VertexColor : WHITE4F;

	Color4f diffuseMtl;

	if (DiffuseMtl)
	{
		DiffuseMtl->Shade(sc);
		diffuseMtl = *block.Color;
	}
	else
		diffuseMtl = WHITE4F;

	Color4f Kd = Mul(PreDiffuse, Mul(diffuseMtl, vertexColor));

	if (fabs(Kd.A) < 0.0001f)
		Kd = Color4f(1.0f, 1.0f, 1.0f, 0.0f);
	else
	{
		Kd.R /= Kd.A;
		Kd.G /= Kd.A;
		Kd.B /= Kd.A;
	}

	float32 y = Kd.Luminance();
	Color4f satColor = Lerp(Color4f(y, y, y, 1.0f), Kd, Saturation);
	satColor.A = 1.0f;

	Color4f textureDetail = Lerp(Color4f(0.0f, 0.0f, 0.0f, 1.0f), satColor, ColorDetail);
	textureDetail.A = 1.0f;

	float32 alphaDetail = Lerp(0.0f, 1.0f - Kd.A, AlphaDetail);
	Color4f outTrans = (1.0f - alphaDetail) * textureDetail + Color4f(alphaDetail);			// compute 1 - (1 - ad) * (1 - td)

	outTrans = Color4f(1.0f) - Mul(Color4f(1.0f) - outTrans, Color4f(1.0f) - Transmittance);

	sc.Transmittance.R = max(min(outTrans.R, 1.0f), 0.0f);
	sc.Transmittance.G = max(min(outTrans.G, 1.0f), 0.0f);
	sc.Transmittance.B = max(min(outTrans.B, 1.0f), 0.0f);
	sc.Transmittance.A = outTrans.A;

	FuASSERT(sc.Transmittance.A >= 0.0f && sc.Transmittance.A <= 1.002f, (""));
}

//-----------------------------------------------------------------------------
// Purpose:
//-----------------------------------------------------------------------------
void C_EntityParticleTrail::AddParticle( float flInitialDeltaTime, const Vector &vecMins, const Vector &vecMaxs, const matrix3x4_t &boxToWorld )
{
    // Select a random point somewhere in the hitboxes of the entity.
    Vector vecLocalPosition, vecWorldPosition;
    vecLocalPosition.x			= Lerp( random->RandomFloat( 0.0f, 1.0f ), vecMins.x, vecMaxs.x );
    vecLocalPosition.y			= Lerp( random->RandomFloat( 0.0f, 1.0f ), vecMins.y, vecMaxs.y );
    vecLocalPosition.z			= Lerp( random->RandomFloat( 0.0f, 1.0f ), vecMins.z, vecMaxs.z );
    VectorTransform( vecLocalPosition, boxToWorld, vecWorldPosition );

    // Don't emit the particle unless it's inside the model
    if ( m_hConstraintEntity.Get() )
    {
        Ray_t ray;
        trace_t tr;
        ray.Init( vecWorldPosition, vecWorldPosition );
        enginetrace->ClipRayToEntity( ray, MASK_ALL, m_hConstraintEntity, &tr );

        if ( !tr.startsolid )
            return;
    }

    // Make a new particle
    SimpleParticle *pParticle = (SimpleParticle *)m_ParticleEffect.AddParticle( sizeof(SimpleParticle), m_hMaterial );
    if ( pParticle == NULL )
        return;

    pParticle->m_Pos			= vecWorldPosition;
    pParticle->m_flRoll			= Helper_RandomInt( 0, 360 );
    pParticle->m_flRollDelta	= Helper_RandomFloat( -2.0f, 2.0f );

    pParticle->m_flLifetime		= flInitialDeltaTime;
    pParticle->m_flDieTime		= m_Info.m_flLifetime;

    pParticle->m_uchColor[0]	= 64;
    pParticle->m_uchColor[1]	= 140;
    pParticle->m_uchColor[2]	= 225;
    pParticle->m_uchStartAlpha	= Helper_RandomInt( 64, 64 );
    pParticle->m_uchEndAlpha	= 0;

    pParticle->m_uchStartSize	= m_Info.m_flStartSize;
    pParticle->m_uchEndSize		= m_Info.m_flEndSize;

    pParticle->m_vecVelocity	= vec3_origin;
    VectorMA( pParticle->m_Pos, flInitialDeltaTime, pParticle->m_vecVelocity, pParticle->m_Pos );
}

float Noise(float x, float y, float z) {
    // Compute noise cell coordinates and offsets
    int ix = Floor2Int(x), iy = Floor2Int(y), iz = Floor2Int(z);
    float dx = x - ix, dy = y - iy, dz = z - iz;

    // Compute gradient weights
    ix &= (NOISE_PERM_SIZE-1);
    iy &= (NOISE_PERM_SIZE-1);
    iz &= (NOISE_PERM_SIZE-1);
    float w000 = Grad(ix,   iy,   iz,   dx,   dy,   dz);
    float w100 = Grad(ix+1, iy,   iz,   dx-1, dy,   dz);
    float w010 = Grad(ix,   iy+1, iz,   dx,   dy-1, dz);
    float w110 = Grad(ix+1, iy+1, iz,   dx-1, dy-1, dz);
    float w001 = Grad(ix,   iy,   iz+1, dx,   dy,   dz-1);
    float w101 = Grad(ix+1, iy,   iz+1, dx-1, dy,   dz-1);
    float w011 = Grad(ix,   iy+1, iz+1, dx,   dy-1, dz-1);
    float w111 = Grad(ix+1, iy+1, iz+1, dx-1, dy-1, dz-1);

    // Compute trilinear interpolation of weights
    float wx = NoiseWeight(dx), wy = NoiseWeight(dy), wz = NoiseWeight(dz);
    float x00 = Lerp(wx, w000, w100);
    float x10 = Lerp(wx, w010, w110);
    float x01 = Lerp(wx, w001, w101);
    float x11 = Lerp(wx, w011, w111);
    float y0 = Lerp(wy, x00, x10);
    float y1 = Lerp(wy, x01, x11);
    return Lerp(wz, y0, y1);
}

		double noise(double x, double y, double z) const
		{
			const std::int32_t X = static_cast<std::int32_t>(std::floor(x)) & 255;
			const std::int32_t Y = static_cast<std::int32_t>(std::floor(y)) & 255;
			const std::int32_t Z = static_cast<std::int32_t>(std::floor(z)) & 255;

			x -= std::floor(x);
			y -= std::floor(y);
			z -= std::floor(z);

			const double u = Fade(x);
			const double v = Fade(y);
			const double w = Fade(z);

			const std::int32_t A = p[X] + Y, AA = p[A] + Z, AB = p[A + 1] + Z;
			const std::int32_t B = p[X + 1] + Y, BA = p[B] + Z, BB = p[B + 1] + Z;

			return Lerp(w, Lerp(v, Lerp(u, Grad(p[AA], x, y, z),
				Grad(p[BA], x - 1, y, z)),
				Lerp(u, Grad(p[AB], x, y - 1, z),
				Grad(p[BB], x - 1, y - 1, z))),
				Lerp(v, Lerp(u, Grad(p[AA + 1], x, y, z - 1),
				Grad(p[BA + 1], x - 1, y, z - 1)),
				Lerp(u, Grad(p[AB + 1], x, y - 1, z - 1),
				Grad(p[BB + 1], x - 1, y - 1, z - 1))));
		}