Subversion Repositories gelsvn

#ifndef __PATHTRACER_CORE__HPP__
#define __PATHTRACER_CORE__HPP__

#include "CGLA/Vec2f.h"
#include "CGLA/Vec3f.h"

#include <limits>

struct ray
{
        ray(void) : origin(0.f), direction(0.f), distance(0.f), depth(0) {}

        CGLA::Vec3f origin;
        CGLA::Vec3f direction;

        float distance;
        int depth;
};

struct hit_info
{
        hit_info(void) : 
                distance(0.f),
                position(0.f), geometric_normal(0.f), shading_normal(0.f),
                inside(false),
                emitted(0.f),
                diffuse(0.f), glossy(0.f), shininess(0.f),
                reflection(0.f), refraction(0.f),
                ior(0.f), extinction(0.f) {}

        float distance;
        CGLA::Vec3f position;
        CGLA::Vec3f geometric_normal;
        CGLA::Vec3f shading_normal;
        
        //The tangent, bitangent, and shading_normal form an orthonormal basis
        CGLA::Vec3f tangent;
        CGLA::Vec3f bitangent;

        //texture coordinates
        CGLA::Vec2f texcoords;

        //flag to indicate if the inside of the object (wrt the normal) was hit
        bool inside;

        //emitted radiance
        CGLA::Vec3f emitted;

        //diffuse reflectivity (rho)
        CGLA::Vec3f diffuse;
        
        //parameters for the modified Phong reflectance model
        CGLA::Vec3f glossy;
        float shininess;

        //parameters for perfect specular reflection/refraction
        CGLA::Vec3f reflection;
        CGLA::Vec3f refraction;

        //index of refraction + complex part
        float ior;
        float extinction;
};

//RGB to intensity
inline float intensity(const CGLA::Vec3f& v)
{
        return v[0] * 0.27f + v[1] * 0.67f + v[2] * 0.06f;
}

//standard fresnel formula (see eg. Jensen p. 23)
inline float fresnel_dielectric(float cosi, float cost, float eta)
{
        assert(cosi>0.f && cost>0.f);

        float etai = 1.f;
        float etat = eta;

        float Fs = (etai*cosi - etat*cost) / (etai*cosi + etat*cost);
        float Fp = (etat*cosi - etai*cost) / (etat*cosi + etai*cost);
        
        float F = 0.5f * (Fs*Fs + Fp*Fp);
        assert(F>=0.f && F<=1.f);
        return F; 
}

//standard fresnel formula (see eg. Pharr et all p. 422)
inline float fresnel_conductor(float cosi, float eta, float extinction)
{
        assert(cosi >= 0.f);

        float z = eta*eta + extinction*extinction;

        float Fp = (z * cosi*cosi - 2.f * eta * cosi + 1.f)     /
                (z + 2.f * eta * cosi + 1.f);

        float Fs = (z - 2.f * eta * cosi + cosi * cosi) /
                (z + 2.f * eta * cosi + cosi * cosi);

        float F = (Fp + Fs) * 0.5f;
        assert(F>=0.f && F<=1.f);
        return F;
}

//reflect vector around normal. n should be normalized.
inline CGLA::Vec3f reflect(const CGLA::Vec3f& n, const CGLA::Vec3f& r)
{
        return 2.f * n * dot(n, r) - r;
}

//eta = eta_from / eta_to
inline bool refract(const CGLA::Vec3f & n, const CGLA::Vec3f & i, 
                                        float eta, CGLA::Vec3f & t)
{
        float c1 = dot(i,n);
        float c2 = 1.f - eta*eta * (1.f - c1*c1);

        if (c2 < 0.f)
                return false;

        float c3 = std::sqrt(c2);
        t = -eta*i + (eta*c1 - c3) * n;

        return true;
}

//helper function, clamps v to range [i; a]
template <class T>
T clamp(const T& v, const T& i, const T& a)
{
        return std::min(std::max(v, i), a);
}

//mersenne twister
extern "C" double genrand_real2(void);
inline float random(void)
{
        return float(genrand_real2());
}

//evaluation of the diffuse brdf
inline CGLA::Vec3f lambertian_brdf(const hit_info& hi,
                                                                   const CGLA::Vec3f& wi,
                                                                   const CGLA::Vec3f& wo)
{
        return hi.diffuse / float(M_PI);
}

//sampling of diffuse brdf, wi = sampled direction, should return probability
inline float sample_lambertian(const hit_info& hi,
                                                        const CGLA::Vec3f& wo, CGLA::Vec3f& wi)
{
        CGLA::Vec3f x, y, z = hi.shading_normal;
        orthogonal(z, x, y);
        float e0 = random();
        float e1 = random();
        float cost = std::sqrt(e0);
        float sint = std::sqrt(1.f - e0);
        float cosp = std::cos(e1 * 2.f * float(M_PI));
        float sinp = std::sin(e1 * 2.f * float(M_PI));
        CGLA::Vec3f dir(sint*cosp, sint*sinp, cost);

        wi = dir(0) * x + dir(1) * y + dir(2) * z;
        return 1.f / float(M_PI) * cost;
}

//evaluation of the phong brdf 
inline CGLA::Vec3f phong_brdf(const hit_info& hi,
                                                           const CGLA::Vec3f& wi,const CGLA::Vec3f& wo)
{
        float dwi = dot(wi, hi.geometric_normal);
        float dwo = dot(wo, hi.geometric_normal);
        float same_hemisphere = dwi * dwo > 0.f;

        if (same_hemisphere)
        {
                //glossy contribution
                if (intensity(hi.glossy) > 0.f)
                {
                        float cost = std::max(dot(reflect(hi.shading_normal, wi), wo),
                                0.f);
                        float cosn = std::pow(cost, hi.shininess);
                        return hi.glossy * (hi.shininess+2.f)/(2.f * float(M_PI))*cosn;
                }
        }
        return CGLA::Vec3f(0.f);
}

//sampling of phong brdf, wi = sampled direction, should return probability
inline float sample_phong(const hit_info& hi,
                                                  const CGLA::Vec3f& wo, CGLA::Vec3f& wi)
{
        float e0 = random();
        float e1 = random();
        
        float cost = std::pow(e0, 1.f/(hi.shininess + 1.f));
        float sint = std::sqrt(1.f - std::pow(e0, 2.f/(hi.shininess + 1.f)));
        float cosp = std::cos(e1 * 2.f * float(M_PI));
        float sinp = std::sin(e1 * 2.f * float(M_PI));
        
        CGLA::Vec3f dir(sint*cosp, sint*sinp, cost);

        CGLA::Vec3f x, y, z = reflect(hi.shading_normal, wo);
        orthogonal(z, x, y);

        float cosn = std::pow(dir(2), hi.shininess);
        wi = dir(0) * x + dir(1) * y + dir(2) * z;

        return (hi.shininess+1.f) / (2.f * float(M_PI)) * cosn;
}

//convenience function to evaluate the non-specular parts of the bsdf
inline CGLA::Vec3f bsdf_evaluate(const hit_info& hi, 
                                                                 const CGLA::Vec3f& wi, const CGLA::Vec3f& wo)
{
        CGLA::Vec3f fs(0.f);

        fs += lambertian_brdf(hi, wi, wo);
        fs += phong_brdf(hi, wi, wo);

        return fs;
}

//epsilon for shadow testing
static const float epsilon = 1e-5f;

static const float step = 1e-1f;

#endif

//02566 framework, Anders Wang Kristensen, awk@imm.dtu.dk, 2007