Brdf.h

#ifndef _BRDF_H
#define _BRDF_H

#include <cstdlib>
#include <cmath>

#include "Geometry.h"
#include "Sampler.h"

namespace gk {

//! repere local de la brdf.
inline float CosTheta( const Vector& w )
{
    return w.z;
}

inline float SinTheta2( const Vector &w )
{
    return std::max(0.f, 1.f - CosTheta(w) * CosTheta(w));
}

inline float SinTheta( const Vector &w )
{
    return sqrtf(SinTheta2(w));
}

inline float CosPhi( const Vector &w )
{
    float sintheta= SinTheta(w);
    if(sintheta == 0.f) 
        return 1.f;
    return Clamp(w.x / sintheta, -1.f, 1.f);
}

inline float SinPhi( const Vector &w )
{
    float sintheta= SinTheta(w);
    if(sintheta == 0.f) 
        return 0.f;
    return Clamp(w.y / sintheta, -1.f, 1.f);
}

//! renvoie vrai si v et w sont du meme cote de la surface.
inline bool SameHemisphere( const Vector &v, const Vector &w )
{
    return (v.z * w.z > 0.f);
}

//! r= direction miroir de w par rapport a n.
inline void Reflect( const Vector& w, const Vector& n, Vector& r )
{
    r= 2.f * Dot(n, w) * n - w;
}

//! renvoie la direction miroir de w par rapport a n.
inline Vector Reflect( const Vector& w, const Vector& n )
{
    return (2.f * Dot(n, w) * n - w);
}

//! genere une direction w dans l'hemisphere. pdf = constante.
inline float SampleUniform( Sampler& sampler, Vector& w )
{
    const float u1= 2.f * M_PI * sampler.uniformFloat();
    const float u2= sampler.uniformFloat();
    
    w= Vector(
        cosf(u1) * sqrtf(1.f - u2*u2), 
        sinf(u1) * sqrtf(1.f - u2*u2), 
        u2);
    return 1.f / (2.f * M_PI);
}

//! renvoie la probabilite de generer la direction w pour une pdf = constante.
inline float pdfUniform( const Vector& w )
{
    return 1.f / (2.f * M_PI);
}

//! genere une direction w dans l'hemisphere. pdf= cos theta / pi.
inline float SampleCos( Sampler& sampler, Vector& w )
{
    const float u1= 2.f * M_PI * sampler.uniformFloat();
    const float u2= sampler.uniformFloat();
    
    w= Vector(
        cosf(u1) * sqrtf(1.f - u2), 
        sinf(u1) * sqrtf(1.f - u2), 
        sqrtf(u2));
    return CosTheta(w) / M_PI;
}

//! renvoie la probabilite de generer la direction w pour une pdf = cos theta / pi
inline float pdfCos( const Vector& w )
{
    return CosTheta(w) / M_PI;
}

//! genere une direction w dans l'hemisphere. pdf = cos**n theta * (n+1) / (2*pi).
inline float SampleCosPow( Sampler& sampler, const float n, Vector& w )
{
    const float u1= 2.f * M_PI * sampler.uniformFloat();
    const float u2= powf(sampler.uniformFloat(), 1.f / (n + 1.f));
    
    w= Vector(
        cosf(u1) * sqrtf(1.f - u2*u2),
        sinf(u1) * sqrtf(1.f - u2*u2),
        u2);
    return powf(CosTheta(w), n) * (n + 1.f) / (2.f * M_PI);
}

//! renvoie la probabilite de generer la direction w pour une pdf= cos**n theta * (n+1) / (2*pi).
inline float pdfCosPow( const float n, const Vector& w )
{
    return powf(CosTheta(w), n) * (n + 1.f) / (2.f * M_PI);
}


//! representation d'une brdf blinn phong, comportements diffus + lobe glossy
class Brdf
{
    Energy diffuse_color;
    Energy specular_color;
    float kd;
    float ks;
    float n;
    
public:
    Brdf( const float _kd, const float _ks, const float _n )
        :
        diffuse_color(1.f),
        specular_color(1.f),
        kd(_kd),
        ks(_ks),
        n(_n)
    {}
        
    Brdf( const float _kd, const Energy& diffuse= gk::Energy(1.f), const float _ks= 0.f, const Energy& specular= Energy(0.f), const float _n= 0.f )
        :
        diffuse_color(diffuse),
        specular_color(specular),
        kd(_kd),
        ks(_ks),
        n(_n)
    {}
    
    //! renvoie l'evaluation du comportement diffus de la brdf.
    Energy fKd( const Vector& wi, const Vector& wo ) const
    {
        if(CosTheta(wi) < 0.f || CosTheta(wo) < 0.f)
            return 0.f;
        return diffuse_color * kd / M_PI;
    }
    
    //! renvoie l'evaluation du comportement glossy de la brdf.
    Energy fKs( const Vector& wi, const Vector& wo ) const
    {
        if(CosTheta(wi) < 0.f || CosTheta(wo) < 0.f)
            return 0.f;
        
        const Vector h= Normalize(wi + wo);
        //~ return ks / M_PI * (n + 1.f) / (2.f * M_PI) * powf(CosTheta(h), n);
        return specular_color * ks * (n + 1.f) / (2.f * M_PI) * powf(CosTheta(h), n);
    }
    
    //! renvoie l'evaluation de la brdf pour le couple de directions wi, wo (repere local, orientes tels que Dot(wi, n) > 0 et Dot(wo, n) > 0).
    Energy f( const Vector& wi, const Vector& wo ) const
    {
        if(CosTheta(wi) < 0.f || CosTheta(wo) < 0.f)
            return 0.f;
        
        const Vector h= Normalize(wi + wo);
        //~ return kd / M_PI + ks / M_PI * (n + 1.f) / (2.f * M_PI) * powf(CosTheta(h), n);
        return diffuse_color * kd / M_PI + specular_color * ks * (n + 1.f) / (2.f * M_PI) * powf(CosTheta(h), n);
    }
    
    //! genere une direction wo, dependante de wi.
    float sample( Sampler& sampler, const Vector& wi, Vector& wo ) const
    {
        float u= sampler.uniformFloat();
        if(u < kd)
        {
            // diffus
            float pkd= SampleCos(sampler, wo);
            
            // calcule la probabilite que l'autre comportement ait genere cette direction
            Vector h= Normalize(wi + wo);
            float pks= pdfCosPow(n, h) / (4.f * Dot(wi, h));
            
            return kd * pkd + ks * pks;
        }
        
        else if(u < kd + ks)
        {
            // glossy
            Vector h;
            float ph= SampleCosPow(sampler, n, h);
            wo= Reflect(wi, h);
            float pks= ph / (4.f * Dot(wi, h));
            if(CosTheta(wo) < 0.f)
                // wo est sous la surface, il ne peut pas etre genere, ni pas le comportement diffus, ni par le comportement glossy
                return 0.f;
            
            // calcule la probabilite que l'autre comportement ait genere cette direction
            float pkd= pdfCos(wo);
            
            return kd * pkd + ks * pks;
        }
        
        // absorption 1 - kd - ks
        wo= Vector(0.f, 0.f, 1.f);
        return 0.f;
    }
    
    //! renvoie la probabilite de choisir wo, connaissant wi (repere local, orientes tels que Dot(wi, n) > 0 et Dot(wo, n) > 0).
    float pdf( const Vector& wi, const Vector& wo ) const
    {
        if(CosTheta(wo) < 0.f)
            return 0.f;
        
        float pkd= pdfCos(wo);
        Vector h= Normalize(wi + wo);
        float pks= pdfCosPow(n, h) / (4.f * Dot(wi, h));
        
        return kd * pkd + ks * pks;
    }
    
    //! evalue rho / reflectance, connaissant wo (repere local, oriente tel que Dot(wo, n) > 0).
    float rho( Sampler& sampler, const Vector& wo, const int samples= 100 ) const
    {
        return rhod() + rhos(sampler, wo, samples);
    }

    //! evalue rho / reflectance (incidence normale == valeur max).
    float rho( Sampler& sampler, const int samples= 100 ) const
    {
        return rhod() + rhos(sampler, samples);
    }

    //! renvoie rho diffus / reflectance.
    float rhod( ) const
    {
        // par definition.
        return kd; 
        
    #if 0
        // verification de la definition
        Sampler sampler;
        const int samples= 1000;
        float r= 0.f;
        for(int i= 0; i < samples; i++)
        {
            Vector wo;
            float p= SampleCos(sampler, wo);
            
            r= r + kd / M_PI * CosTheta(wo) / p;
        }
        
        return r / (float) samples;
    #endif
    }

    //! evalue rho / reflectance glossy, pour une direction wo (repere local, oriente tel que Dot(wo, n) > 0).
    float rhos( Sampler& sampler, const Vector& wo, const int samples= 100 ) const
    {
        float r= 0.f;
        for(int i= 0; i < samples; i++)
        {
            Vector h;
            float ph= SampleCosPow(sampler, n, h);
            
            Vector w= Reflect(wo, h);
            if(CosTheta(w) < 0.f)
                continue;
            float p= ph / (4.f * Dot(wo, h));
            
            //~ r= r + ks / M_PI * (n + 1.f) / (2.f * M_PI) * powf(CosTheta(h), n) * CosTheta(w) / p;
            r= r + ks * (n + 1.f) / (2.f * M_PI) * powf(CosTheta(h), n) * CosTheta(w) / p;
        }
        
        return r / (float) samples;
    }

    //! evalue rho / reflectance speculaire max (incidence normale).
    float rhos( Sampler& sampler, const int samples= 100 ) const
    {
        return rhos(sampler, Vector(0.f, 0.f, 1.f), samples);
    }
};

}

#endif

---