main.cpp

#include <stdint.h>
#include <stdlib.h>
#include <math.h>
#include <algorithm>
#include "dds.h"
#include "csv.h"

uint32_t ReverseBits( uint32_t v )
{
    v = ( ( v >> 1 ) & 0x55555555 ) | ( ( v & 0x55555555 ) << 1 );
    v = ( ( v >> 2 ) & 0x33333333 ) | ( ( v & 0x33333333 ) << 2 );
    v = ( ( v >> 4 ) & 0x0F0F0F0F ) | ( ( v & 0x0F0F0F0F ) << 4 );
    v = ( ( v >> 8 ) & 0x00FF00FF ) | ( ( v & 0x00FF00FF ) << 8 );
    v = (   v >> 16               ) | (   v                << 16 );
    return v;
}

union FP32
{
	unsigned	u;
	float		f;
};

// https://gist.github.com/rygorous/2156668
uint16_t FloatToHalf( float ff )
{
	FP32 f32infty = { 255 << 23 };
	FP32 f16infty = { 31 << 23 };
	FP32 magic = { 15 << 23 };
	unsigned sign_mask = 0x80000000u;
	unsigned round_mask = ~0xfffu; 

	uint16_t o = 0;
	FP32 f;
	f.f = ff;

	unsigned sign = f.u & sign_mask;
	f.u ^= sign;

	// NOTE all the integer compares in this function can be safely
	// compiled into signed compares since all operands are below
	// 0x80000000. Important if you want fast straight SSE2 code
	// (since there's no unsigned PCMPGTD).

	if (f.u >= f32infty.u) // Inf or NaN (all exponent bits set)
		o = (f.u > f32infty.u) ? 0x7e00 : 0x7c00; // NaN->qNaN and Inf->Inf
	else // (De)normalized number or zero
	{
		f.u &= round_mask;
		f.f *= magic.f;
		f.u -= round_mask;
		if (f.u > f16infty.u) f.u = f16infty.u; // Clamp to signed infinity if overflowed

		o = uint16_t( f.u >> 13 ); // Take the bits!
	}

	o |= sign >> 16;
	return o;
}

float GSmith( float roughness, float ndotv, float ndotl )
{
    float const m2   = roughness * roughness;
    float const visV = ndotv + sqrt( ndotv * ( ndotv - ndotv * m2 ) + m2 );
    float const visL = ndotl + sqrt( ndotl * ( ndotl - ndotl * m2 ) + m2 );
    return 1.0f / ( visV * visL );
}

float GSmithCorrelated( float roughness, float ndotv, float ndotl )
{
	float m2	= roughness * roughness;
	float visV	= ndotl * sqrt( ndotv * ( ndotv - ndotv * m2 ) + m2 );
	float visL	= ndotv * sqrt( ndotl * ( ndotl - ndotl * m2 ) + m2 );
	return 0.5f / ( visV + visL );
}

float GSmithShclick( float roughness, float ndotv, float ndotl )
{
	float k = roughness * 0.5f;
	float visV = 0.5f / (ndotv * (1.0f - k) + k);
	float visL = 0.5f / (ndotl * (1.0f - k) + k);
	return visV * visL;
}

int main()
{
    float const MATH_PI         = 3.14159f;
    unsigned const LUT_WIDTH    = 32;
    unsigned const LUT_HEIGHT   = 64;
    unsigned const sampleNum    = 128;

    float lutDataRGBA32F[ LUT_WIDTH * LUT_HEIGHT * 4 ];
	uint16_t lutDataRG16F[ LUT_WIDTH * LUT_HEIGHT * 2 ];

    for ( unsigned y = 0; y < LUT_HEIGHT; ++y )
    {
        float const ndotv = ( y + 0.5f ) / LUT_HEIGHT;

        for ( unsigned x = 0; x < LUT_WIDTH; ++x )
        {
            float const gloss       = ( x + 0.5f ) / LUT_WIDTH;
            float const roughness   = powf( 1.0f - gloss, 2.0f );

            float const vx = sqrtf( 1.0f - ndotv * ndotv );
            float const vy = 0.0f;
            float const vz = ndotv;

            float scale = 0.0f;
            float bias  = 0.0f;

            
            for ( unsigned i = 0; i < sampleNum; ++i )
            {
                float const e1 = (float) i / sampleNum;
                float const e2 = (float) ( (double) ReverseBits( i ) / (double) 0x100000000LL );

                float const phi         = 2.0f * MATH_PI * e1;
                float const cosPhi      = cosf( phi );
                float const sinPhi      = sinf( phi );
                float const cosTheta    = sqrtf( ( 1.0f - e2 ) / ( 1.0f + ( roughness * roughness - 1.0f ) * e2 ) );
                float const sinTheta    = sqrtf( 1.0f - cosTheta * cosTheta );

                float const hx  = sinTheta * cosf( phi );
                float const hy  = sinTheta * sinf( phi );
                float const hz  = cosTheta;

                float const vdh = vx * hx + vy * hy + vz * hz;
                float const lx  = 2.0f * vdh * hx - vx;
                float const ly  = 2.0f * vdh * hy - vy;
                float const lz  = 2.0f * vdh * hz - vz;

                float const ndotl = std::max( lz,   0.0f );
                float const ndoth = std::max( hz,   0.0f );
                float const vdoth = std::max( vdh,  0.0f );

                if ( ndotl > 0.0f )
                {
                    float const gsmith      = GSmithShclick( roughness, ndotv, ndotl );
                    float const ndotlVisPDF = ndotl * gsmith * ( 4.0f * vdoth / ndoth );
                    float const fc          = powf( 1.0f - vdoth, 5.0f );

                    scale   += ndotlVisPDF * ( 1.0f - fc );
                    bias    += ndotlVisPDF * fc;
                }
            }
            scale /= sampleNum;
            bias  /= sampleNum;

            lutDataRGBA32F[ x * 4 + y * LUT_WIDTH * 4 + 0 ] = scale;
            lutDataRGBA32F[ x * 4 + y * LUT_WIDTH * 4 + 1 ] = bias;
            lutDataRGBA32F[ x * 4 + y * LUT_WIDTH * 4 + 2 ] = 0.0f;
            lutDataRGBA32F[ x * 4 + y * LUT_WIDTH * 4 + 3 ] = 0.0f;

			lutDataRG16F[ x * 2 + y * LUT_WIDTH * 2 + 0 ] = FloatToHalf( scale );
			lutDataRG16F[ x * 2 + y * LUT_WIDTH * 2 + 1 ] = FloatToHalf( bias );
        }
    }   

    SaveDDS( "integrateDFG_RGBA32F.dds", DDS_FORMAT_R32G32B32A32_FLOAT, 16, LUT_WIDTH, LUT_HEIGHT, lutDataRGBA32F );
	SaveDDS( "integrateDFG_RG16F.dds", DDS_FORMAT_R16G16_FLOAT, 4, LUT_WIDTH, LUT_HEIGHT, lutDataRG16F );
    SaveCSV( "ndotv.csv", LUT_WIDTH );
    SaveCSV( "gloss.csv", LUT_HEIGHT );
    SaveCSV( "scale.csv", lutDataRGBA32F, LUT_WIDTH, LUT_HEIGHT, 0 );
    SaveCSV( "bias.csv",  lutDataRGBA32F, LUT_WIDTH, LUT_HEIGHT, 1 );

    return 0;
}