tuto_mcdirect.cpp

 
#include <cstdio>
#include <cassert>
#include <vector>
#include <random>
 
#include "vec.h"
#include "mat.h"
#include "orbiter.h"
#include "gltf.h"
 
#include "color.h"
#include "image.h"
#include "image_hdr.h"
 
// rayon 
struct Ray
{
    Point o;            // origine
    float pad;
    Vector d;           // direction
    float tmax;         // tmax= 1 ou \inf, le rayon est un segment ou une demi droite infinie
    
    Ray( const Point& _o, const Point& _e ) :  o(_o), d(Vector(_o, _e)), tmax(1) {} // segment, t entre 0 et 1
    Ray( const Point& _o, const Vector& _d ) :  o(_o), d(_d), tmax(FLT_MAX) {}  // demi droite, t entre 0 et \inf
    Ray( const Point& _o, const Vector& _d, const float _tmax ) :  o(_o), d(_d), tmax(_tmax) {} // explicite
};
 
// intersection avec un triangle
struct Hit
{
    float t;            // p(t)= o + td, position du point d'intersection sur le rayon
    float u, v;         // p(u, v), position du point d'intersection sur le triangle
    int node_id;    // indice de l'instance, pour retrouver la transformation
    int mesh_id;        // indexation globale du triangle dans la scene gltf
    int primitive_id;
    int triangle_id;
    int pad;
    
    Hit( ) : t(FLT_MAX), u(), v(), node_id(-1), mesh_id(-1), primitive_id(-1), triangle_id(-1) {}
    Hit( const Ray& ray ) : t(ray.tmax), u(), v(), node_id(-1), mesh_id(-1), primitive_id(-1), triangle_id(-1) {}
    
    Hit( const float _t, const float _u, const float _v, const int _node_id, const int _mesh_id, const int _primitive_id, const int _id ) : t(_t), u(_u), v(_v), 
        node_id(_node_id), mesh_id(_mesh_id), primitive_id(_primitive_id), triangle_id(_id) {}
    
    operator bool ( ) { return (triangle_id != -1); }   // renvoie vrai si l'intersection est definie / existe
};
 
// intersection avec une boite / un englobant
struct BBoxHit
{
    float tmin, tmax;
    
    BBoxHit() : tmin(FLT_MAX), tmax(-FLT_MAX) {}
    BBoxHit( const float _tmin, const float _tmax ) : tmin(_tmin), tmax(_tmax) {}
    
    operator bool( ) const { return tmin <= tmax; }   // renvoie vrai si l'intersection est definie / existe
};
 
 
// boite englobante
struct BBox
{
    Point pmin, pmax;
    
    BBox( ) : pmin(), pmax() {}
    
    BBox( const Point& p ) : pmin(p), pmax(p) {}
    BBox( const BBox& box ) : pmin(box.pmin), pmax(box.pmax) {}
    BBox( const BBox& a, const BBox& b ) : pmin(min(a.pmin, b.pmin)), pmax(max(a.pmax, b.pmax)) {}
    
    BBox& insert( const Point& p ) { pmin= min(pmin, p); pmax= max(pmax, p); return *this; }
    BBox& insert( const BBox& box ) { pmin= min(pmin, box.pmin); pmax= max(pmax, box.pmax); return *this; }
    
    float centroid( const int axis ) const { return (pmin(axis) + pmax(axis)) / 2; }
    Point centroid( ) const { return (pmin + pmax) / 2; }
    
    BBoxHit intersect( const Ray& ray, const Vector& invd, const float htmax ) const
    {
        Point rmin= pmin;
        Point rmax= pmax;
        if(ray.d.x < 0) std::swap(rmin.x, rmax.x);
        if(ray.d.y < 0) std::swap(rmin.y, rmax.y);
        if(ray.d.z < 0) std::swap(rmin.z, rmax.z);
        Vector dmin= (rmin - ray.o) * invd;
        Vector dmax= (rmax - ray.o) * invd;
        
        float tmin= std::max(dmin.z, std::max(dmin.y, std::max(dmin.x, 0.f)));
        float tmax= std::min(dmax.z, std::min(dmax.y, std::min(dmax.x, htmax)));
        return BBoxHit(tmin, tmax);
    }
};
 
 
// construction de l'arbre / BVH
struct Node
{
    BBox bounds;
    int left;
    int right;
    
    bool internal( ) const { return right > 0; }                        // renvoie vrai si le noeud est un noeud interne
    int internal_left( ) const { assert(internal()); return left; }     // renvoie le fils gauche du noeud interne 
    int internal_right( ) const { assert(internal()); return right; }   // renvoie le fils droit
    
    bool leaf( ) const { return right < 0; }                            // renvoie vrai si le noeud est une feuille
    int leaf_begin( ) const { assert(leaf()); return -left; }           // renvoie le premier objet de la feuille
    int leaf_end( ) const { assert(leaf()); return -right; }            // renvoie le dernier objet
};
 
// creation d'un noeud interne
Node make_node( const BBox& bounds, const int left, const int right )
{
    Node node { bounds, left, right };
    assert(node.internal());    // verifie que c'est bien un noeud...
    return node;
}
 
// creation d'une feuille
Node make_leaf( const BBox& bounds, const int begin, const int end )
{
    Node node { bounds, -begin, -end };
    assert(node.leaf());        // verifie que c'est bien une feuille...
    return node;
}
 
 
// bvh parametre par le type des primitives, cf triangle et instance...
template < typename T >
struct BVHT
{
    // construit un bvh pour l'ensemble de primitives
    int build( const std::vector<T>& _primitives )
    {
        primitives= _primitives;  // copie les primitives pour les trier
        nodes.clear();          // efface les noeuds
        nodes.reserve(primitives.size());
        
        // construit l'arbre... 
        root= build(0, primitives.size());
        return root;
    }
    
    // intersection avec un rayon, entre 0 et htmax
    Hit intersect( const Ray& ray, const float htmax ) const
    {
        Hit hit;
        hit.t= htmax;
        Vector invd= Vector(1 / ray.d.x, 1 / ray.d.y, 1 / ray.d.z);
        intersect_fast(root, ray, invd, hit);
        return hit;
    }
    
    // intersection avec un rayon, entre 0 et ray.tmax
    Hit intersect( const Ray& ray ) const { return intersect(ray, ray.tmax); }
    
    bool occluded( const Ray& ray ) const
    {
        Vector invd= Vector(1 / ray.d.x, 1 / ray.d.y, 1 / ray.d.z);
        return occluded(root, ray, invd, 1);
    }
    
    bool occluded( const Point& p, const Point& q ) const { Ray ray(p, q); return occluded(ray); }
    
    bool visible( const Point& p, const Point& q ) const { return !occluded(p, q); }
    bool visible( const Ray& ray ) const { return !occluded(ray); }
    
    
protected:
    std::vector<Node> nodes;
    std::vector<T> primitives;
    int root;
    
    int build( const int begin, const int end )
    {
        if(end - begin < 2)
        {
            // inserer une feuille et renvoyer son indice
            int index= nodes.size();
            nodes.push_back( make_leaf( primitive_bounds(begin, end), begin, end ) );
            return index;
        }
        
        // axe le plus etire de l'englobant des centres des englobants des primitives... 
        BBox cbounds= centroid_bounds(begin, end);
        Vector d= Vector(cbounds.pmin, cbounds.pmax);
        int axis;
        if(d.x > d.y && d.x > d.z)  // x plus grand que y et z ?
            axis= 0;
        else if(d.y > d.z)          // y plus grand que z ? (et que x implicitement)
            axis= 1;
        else                        // x et y ne sont pas les plus grands...
            axis= 2;
 
        // coupe l'englobant au milieu
        float cut= cbounds.centroid(axis);
        
        // repartit les primitives 
        T *pm= std::partition(primitives.data() + begin, primitives.data() + end, 
            [axis, cut]( const T& primitive ) 
            { 
                return primitive.bounds().centroid(axis) < cut; 
            }
        );
        int m= std::distance(primitives.data(), pm);
        
        // la repartition peut echouer, et toutes les primitives sont dans la meme moitiee de l'englobant
        // forcer quand meme un decoupage en 2 ensembles 
        if(m == begin || m == end)
            m= (begin + end) / 2;
        assert(m != begin);
        assert(m != end);
        
        // construire le fils gauche, les triangles se trouvent dans [begin .. m)
        int left= build(begin, m);
        
        // on recommence pour le fils droit, les triangles se trouvent dans [m .. end)
        int right= build(m, end);
        
        // construire le noeud et renvoyer son indice
        int index= nodes.size();
        nodes.push_back( make_node( BBox(nodes[left].bounds, nodes[right].bounds), left, right ) );
        return index;
    }
    
    // englobant des primitives
    BBox primitive_bounds( const int begin, const int end )
    {
        BBox bbox= primitives[begin].bounds();
        for(int i= begin +1; i < end; i++)
            bbox.insert(primitives[i].bounds());
        
        return bbox;
    }
    
    // englobant des centres des primitives
    BBox centroid_bounds( const int begin, const int end )
    {
        BBox bbox= primitives[begin].bounds().centroid();
        for(int i= begin +1; i < end; i++)
            bbox.insert(primitives[i].bounds().centroid());
        
        return bbox;
    }
    
    void intersect_fast( const int index, const Ray& ray, const Vector& invd, Hit& hit ) const
    {
        const Node& node= nodes[index];
    
        if(node.leaf())
        {
            //~ if(node.bounds.intersect(ray, invd, hit.t)) // deja teste !!
            {
                for(int i= node.leaf_begin(); i < node.leaf_end(); i++)
                    if(Hit h= primitives[i].intersect(ray, hit.t))
                        hit= h;
            }
        }
        else
        {
            BBoxHit left= nodes[node.internal_left()].bounds.intersect(ray, invd, hit.t);
            BBoxHit right= nodes[node.internal_right()].bounds.intersect(ray, invd, hit.t);
            
            // visiter les 2 fils 
            if(left && right)
            {
                if(left.tmin < right.tmin)  // gauche puis droit
                {
                    intersect_fast(node.internal_left(), ray, invd, hit);
                    intersect_fast(node.internal_right(), ray, invd, hit);
                }
                else    // d'abord droit
                {
                    intersect_fast(node.internal_right(), ray, invd, hit);
                    intersect_fast(node.internal_left(), ray, invd, hit);
                }
            }
            
            // visiter un seul fils
            else if(left)
                intersect_fast(node.internal_left(), ray, invd, hit);
            else if(right)
                intersect_fast(node.internal_right(), ray, invd, hit);
        }
    }    
    
    bool occluded( const int index, const Ray& ray, const Vector& invd, const float htmax ) const
    {
        const Node& node= nodes[index];
        if(node.bounds.intersect(ray, invd, htmax))
        {
            if(node.leaf())
            {
                for(int i= node.leaf_begin(); i < node.leaf_end(); i++)
                    if(primitives[i].intersect(ray, htmax))
                        return true;
            }
            else // if(node.internal())
            {
                if(occluded(node.internal_left(), ray, invd, htmax) || occluded(node.internal_right(), ray, invd, htmax))
                    return true;
            }
        }
        
        return false;
    }
};
 
 
// triangle pour le bvh, cf fonction bounds() et intersect()
struct Triangle
{
    Point p;            // sommet a du triangle
    Vector e1, e2;      // aretes ab, ac du triangle
    int node_id;
    int mesh_id;
    int primitive_id;
    int triangle_id;
    
    Triangle( const vec3& a, const vec3& b, const vec3& c, const int _node_id, const int _mesh_id, const int _primitive_id, const int _id ) : 
        p(a), e1(Vector(a, b)), e2(Vector(a, c)), 
        node_id(_node_id), mesh_id(_mesh_id), primitive_id(_primitive_id), triangle_id(_id) {}
    
    /* calcule l'intersection ray/triangle
        cf "fast, minimum storage ray-triangle intersection" 
        
        renvoie faux s'il n'y a pas d'intersection valide (une intersection peut exister mais peut ne pas se trouver dans l'intervalle [0 tmax] du rayon.)
        renvoie vrai + les coordonnees barycentriques (u, v) du point d'intersection + sa position le long du rayon (t).
        convention barycentrique : p(u, v)= (1 - u - v) * a + u * b + v * c
    */
    Hit intersect( const Ray &ray, const float htmax ) const
    {
        Vector pvec= cross(ray.d, e2);
        float det= dot(e1, pvec);
        
        float inv_det= 1 / det;
        Vector tvec(p, ray.o);
        
        float u= dot(tvec, pvec) * inv_det;
        if(u < 0 || u > 1) return Hit();
        
        Vector qvec= cross(tvec, e1);
        float v= dot(ray.d, qvec) * inv_det;
        if(v < 0 || u + v > 1) return Hit();
        
        float t= dot(e2, qvec) * inv_det;
        if(t < 0 || t > htmax) return Hit();
        
        return Hit(t, u, v, node_id, mesh_id, primitive_id, triangle_id);
    }
    
    BBox bounds( ) const 
    {
        BBox box(p);
        return box.insert(p+e1).insert(p+e2);
    }
};
 
typedef BVHT<Triangle> BVH;
 
 
// renvoie la positions sur le rayon
Point hit_position( const Hit& hit, const Ray& ray ) 
{ 
    assert(hit.triangle_id != -1); 
    return ray.o + hit.t * ray.d; 
}
 
 
//
static GLTFMaterial default_material= GLTFMaterial();
 
// renvoie la matiere du point d'intersection
const GLTFMaterial& hit_material( const Hit& hit, const GLTFScene& scene )
{
    assert(hit.mesh_id != -1);
    assert(hit.primitive_id != -1);
    
    const GLTFMesh& mesh= scene.meshes[hit.mesh_id];
    const GLTFPrimitives& primitives= mesh.primitives[hit.primitive_id];
    
    if(primitives.material_index == -1)
        return default_material;
    
    assert(primitives.material_index < int(scene.materials.size()));
    return scene.materials[primitives.material_index];
}
 
Color material_diffuse( const GLTFMaterial& material )
{
    return (1 - material.metallic) * material.color / float(M_PI);
}
 
 
// renvoie la normale interpolee au point d'intersection dans le repere de la scene
Vector hit_normal( const Hit& hit, const GLTFScene& scene )
{
    assert(hit.node_id != -1);
    assert(hit.mesh_id != -1);
    assert(hit.primitive_id != -1);
    assert(hit.triangle_id != -1);
    
    const GLTFMesh& mesh= scene.meshes[hit.mesh_id];
    const GLTFPrimitives& primitives= mesh.primitives[hit.primitive_id];
    
    // indice des sommets
    int a= primitives.indices[hit.triangle_id];
    int b= primitives.indices[hit.triangle_id+1];
    int c= primitives.indices[hit.triangle_id+2];
 
    // normales des sommets
    assert(primitives.normals.size());
    Vector na= primitives.normals[a];
    Vector nb= primitives.normals[b];
    Vector nc= primitives.normals[c];
    
    // interpole la normale au point d'intersection
    // attention : il faut utiliser la meme convetion barycentrique que la fonction d'intersection rayon/triangle !!
    Vector n= (1 - hit.u - hit.v) * na + hit.u * nb + hit.v * nc;
    
    // transforme la normale dans le repere de la scene
    const GLTFNode& node= scene.nodes[hit.node_id];
    // les normales ne se transforment pas exactement comme les positions... 
    // les sommets sont transformes par node.model, comment transformer la normale pour quelle reste perpendiculaire au triangle
    Transform T= node.model.normal();   // todo !! a precalculer, inverser la matrice a chaque fois est un peu trop lent...
    
    return normalize( T(n) );
}
 
 
// construit un repere tbn, a partir d'un seul vecteur.
// cf "generating a consistently oriented tangent space" 
// http://people.compute.dtu.dk/jerf/papers/abstracts/onb.html
// et
// cf "Building an Orthonormal Basis, Revisited", Pixar, 2017
// http://jcgt.org/published/0006/01/01/
// pour corriger un probleme numerique...
struct World
{
    World( ) : t(), b(), n() {}
    World( const Vector& _n ) : n(_n) 
    {
        float sign= std::copysign(1.0f, n.z);
        float a= -1.0f / (sign + n.z);
        float d= n.x * n.y * a;
        t= Vector(1.0f + sign * n.x * n.x * a, sign * d, -sign * n.x);
        b= Vector(d, sign + n.y * n.y * a, -n.y);        
    }
    
    // transforme le vecteur du repere local vers le repere du monde
    Vector operator( ) ( const Vector& local )  const { return local.x * t + local.y * b + local.z * n; }
    
    // transforme le vecteur du repere du monde vers le repere local
    Vector local( const Vector& global ) const { return Vector(dot(global, t), dot(global, b), dot(global, n)); }
    
    Vector t;
    Vector b;
    Vector n;
};
 
 
struct Sampler
{
    std::uniform_real_distribution<float> u01;
    std::default_random_engine rng;
    
    //~ Sampler( const unsigned _seed ) : u01(0, 1 - 2*std::numeric_limits<float>::epsilon()), rng(_seed) {}
    Sampler( const unsigned _seed ) : u01(), rng(_seed) {}
    void seed( const unsigned _seed ) { rng.seed(_seed); }
    
    float sample( ) { return u01(rng); }    // renvoie un reel [0, 1)
    unsigned sample_binary( ) { return rng(); }     // renvoie un entier 32bits
    int sample_range( const int n ) { return int(sample() * n); }   // renvoie un entier [0 n)
};
 
 
struct Source
{
    Point a, b, c;
    Vector n;
    Color emission;
    float area;
    
    Source( const Point& _a, const Point& _b, const Point& _c, const Color& _emission ) : a(_a), b(_b), c(_c), n(), emission(_emission), area()
    {
        Vector ng= cross( Vector(a, b), Vector(a, c) );
        area= length(ng) / 2;
        n= normalize(ng);
        
        assert(area * emission.max() > 0);
    }
    
    Point sample( const float u1, const float u2 ) const
    {
    #if 0
        // cf GI compendium eq 18, https://people.cs.kuleuven.be/~philip.dutre/GI/TotalCompendium.pdf
        // ou PBRT http://www.pbr-book.org/3ed-2018/Monte_Carlo_Integration/2D_Sampling_with_Multidimensional_Transformations.html#SamplingaTriangle
        float r1= std::sqrt(u1);
        float alpha= 1 - r1;
        float beta= (1 - u2) * r1;
        float gamma= u2 * r1;
        
        return alpha*a + beta*b + gamma*c;
        
    #else
        // cf A Low-Distortion Map Between Triangle and Square
        // https://hal.archives-ouvertes.fr/hal-02073696v2/document
        float t1= u1 / 2;
        float t2= u2 / 2;
        float offset= t2 - t1;
        if(offset > 0)
            t2+= offset;
        else
            t1-= offset;
            
        return t1*a + t2*b + (1 - t1 - t2)*c;
    #endif
    }
    
    float pdf( const Point& p ) const 
    {
        return 1 / area;
    }
};
 
 
float blinn( const float roughness )
{
    if(roughness < 0.001)
        return 1000;
    
    return 2 / (roughness*roughness*roughness*roughness) - 2;
}
 
// matiere : melange diffuse + reflet
// modele blinn-phong normalise, cf http://perso.univ-lyon1.fr/jean-claude.iehl/Public/educ/M1IMAGE/html/group__brdf.html
struct Brdf
{
    World world;
    Color diffuse;
    float kd;
    float ns;
    
    Brdf( ) : world(), diffuse(), kd(), ns() {}
    Brdf( const Vector& ng, const Color& color, const float plastic, const float n ) : world(ng), diffuse(color), kd(1 - plastic), ns(n) {}
    
    Color f( const Vector& l, const Vector& o ) const
    {
        float cos_theta= dot(world.n, l);
        if(cos_theta <= 0)
            return Black();
        
        Vector h= normalize(o + l);
        float cos_theta_h= dot(world.n, h);
        if(cos_theta_h <= 0)
            return Black();
        
        // evalue la brdf : mixture d'un terme diffus et d'un reflet,
        Color d= 1 / float(M_PI) * diffuse;
        Color s= (ns + 8) / float(8*M_PI) * std::pow(cos_theta_h, ns) * White();
        return kd * d + (1 - kd) * s;
    }
    
    Vector sample( const float u1, const float u2, const float u3, const Vector& o ) const
    {
        // echantillonne la mixture diffus + reflet
        if(u1 <= kd)
        {
            // terme diffus
            // genere une direction cos theta / pi, cf GI compendium, eq 35
            float phi= float(2*M_PI) * u3;
            float cos_theta= std::sqrt(u2);
            float sin_theta= std::sqrt(1 - cos_theta*cos_theta);
            
            return world(Vector(std::cos(phi) * sin_theta, std::sin(phi) * sin_theta, cos_theta));
        }
        else
        {
            // terme reflechissant
            // genere une direction h
            // genere une direction cos^n theta / pi, cf GI compendium, eq 35+
            float phi= float(2*M_PI) * u3;
            float cos_theta= std::pow(u2, 1 / float(ns +1));
            assert(cos_theta >= 0);
            float sin_theta= std::sqrt(1 - cos_theta*cos_theta);
            Vector h= world(Vector(std::cos(phi) * sin_theta, std::sin(phi) * sin_theta, cos_theta));
            // genere une direction reflechie, l= reflect(o | h)
            return -o + 2 * dot(h, o) * h;    
            // attention : peut etre sous le plan tangent...
        }
    }
    
    float pdf( const Vector& l, const Vector& o ) const
    {
        float cos_theta= dot(world.n, l);
        if(cos_theta <= 0)
            return 0;
        
        // pdf du terme diffus 
        float d= cos_theta / float(M_PI);
        
        // pdf du terme reflechissant
        Vector h= normalize(o + l);
        if(dot(world.n, h) <= 0)
            return 0;
        if(dot(l, h) <= 0)
            return 0;
        
        float cos_theta_h= std::max(float(0), dot(world.n, h));
        float s= (ns + 1) / float(2*M_PI) * std::pow(cos_theta_h, ns) / (4 * dot(l, h));
        // le terme 1 / dot(l, h) est introduit par le changement de variable : direction h vers direction reflechie 
        
        // pdf de la mixture
        return kd*d + (1 - kd)*s;
    }
};
 
 
 
 
 
struct MCColor
{
    Color color;
    Color M;
    Color M2;
    int n;
};
 
const float epsilon= 0.001;
    
Image direct( const Orbiter& camera, const GLTFScene& scene, const BVH& bvh, const int samples, const std::vector<Source>& sources )
{
    Image image(camera.width(), camera.height());
    
    const int AA= 1;
    const int N= samples;
    
    printf("render  %dx%d %d samples...\n", image.width(), image.height(), samples);
    
    // normalisation pdf source
    float total= 0;
    for(unsigned i= 0; i < sources.size(); i++)
        total+= sources[i].area;
    
    Transform model= Identity();
    Transform view= camera.view();
    Transform projection= camera.projection();
    Transform viewport= camera.viewport();
    Transform inv= Inverse(viewport *projection * view);        // image vers scene
    
    #pragma omp parallel for schedule(dynamic, 1)
    for(int py= 0; py < image.height(); py++)
    {
        std::random_device hwseed;
        Sampler rng( hwseed() );        
        
        for(int px= 0; px < image.width(); px++)
        {
            Color M;
            Color M2;
            int n= 0;
            
            Color emission;
            for(int pa= 0; pa < AA; pa++)
            {
                //~ float ax= rng.sample();
                //~ float ay= rng.sample();
                float ax= 0;
                float ay= 0;
                
                Point o= inv( Point(px +ax, py +ay, 0) );
                Point e= inv( Point(px +ax, py +ay, 1) );
                Vector d= Vector(o, e);
                
                Ray ray(o, d);
                if(Hit hit= bvh.intersect(ray))
                {
                    Point p= hit_position(hit, ray);
                    Vector pn= hit_normal(hit, scene);
                    Color diffuse= material_diffuse( hit_material(hit, scene) );
                    
                    if(dot(ray.d, pn) > 0)      // retourne la normale si elle n'est pas orientee vers l'origine du rayon...
                        pn= -pn;
                    
                    emission= emission + hit_material(hit, scene).emission;
                    
                    World world(pn);
                    for(int i= 0; i < N; i++)
                    {
                        //~ int id= -1;
                        int id= int(sources.size()) -1;
                        {
                            float u= rng.sample();
                            // construit et parcours la fonction de repartition
                            float cdf= 0;
                            for(unsigned i= 0; i < sources.size(); i++)
                            {
                                cdf+= sources[i].area / total;
                                if(u < cdf)
                                {
                                    id= i;
                                    break;
                                }
                            }
                            
                            if(u > cdf) 
                            { 
                                //~ printf("cdf %.9f %.9f / %.9f\n", u, cdf, total); 
                                if(id == -1) id= int(sources.size()) -1; 
                            }
                            
                            assert(id != -1);
                            assert(id < int(sources.size()));
                        }
                        
                        // construit le point q sur la source selectionnee
                        Point q= sources[id].sample(rng.sample(), rng.sample());
                        float pdf= 1 / total;
                        
                        Color sample;
                        if(bvh.visible( p+epsilon*pn, q+epsilon*sources[id].n ))
                        {
                            Vector l= Vector(p, q);
                            float cos_theta= std::max( float(0), dot( normalize(pn), normalize(l) ) );
                            float cos_theta_q= std::max( float(0), -dot( normalize(sources[id].n), normalize(l) ) );
                            
                            sample= diffuse * sources[id].emission * cos_theta * cos_theta_q / length2(l) / pdf;
                        }
                        
                        // accumule le sample,
                        // cf https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
                        n++;
                        Color Mn1= M;
                        M= Mn1 + (sample - Mn1) / float(n);
                        M2= M2 + (sample - Mn1) * (sample - M);
                    }
                }
            }
            
            emission= emission / float(AA);
            image(px, py)= Color(M + emission, 1);
            //~ image(px, py)= Color(M2 / float(n), 1);
        }
    }
    
    return image;
}
    
Image direct( const Orbiter& camera, const GLTFScene& scene, const BVH& bvh, const int samples )
{
    Image image(camera.width(), camera.height());
    
    const int AA= 1;
    const int N= samples;
    
    printf("render %dx%d %d samples...\n", image.width(), image.height(), samples);
    
    Transform model= Identity();
    Transform view= camera.view();
    Transform projection= camera.projection();
    Transform viewport= camera.viewport();
    Transform inv= Inverse(viewport *projection * view);        // image vers scene
    
    #pragma omp parallel for schedule(dynamic, 1)
    for(int py= 0; py < image.height(); py++)
    {
        std::random_device hwseed;
        Sampler rng( hwseed() );        
        
        for(int px= 0; px < image.width(); px++)
        {
            Color M;
            Color M2;
            int n= 0;
            
            Color emission;
            for(int pa= 0; pa < AA; pa++)
            {
                //~ float ax= rng.sample();
                //~ float ay= rng.sample();
                float ax= 0;
                float ay= 0;
                
                Point o= inv( Point(px +ax, py +ay, 0) );
                Point e= inv( Point(px +ax, py +ay, 1) );
                Vector d= Vector(o, e);
                
                Ray ray(o, d);
                if(Hit hit= bvh.intersect(ray))
                {
                    Point p= hit_position(hit, ray);
                    Vector pn= hit_normal(hit, scene);
                    Color diffuse= material_diffuse( hit_material(hit, scene) );
                    
                    if(dot(ray.d, pn) > 0)      // retourne la normale si elle n'est pas orientee vers l'origine du rayon...
                        pn= -pn;
                    
                    emission= emission + hit_material(hit, scene).emission;
                    
                    World world(pn);                    
                    for(int i= 0; i < N; i++)
                    {
                        Vector l;
                        float pdf;
                        #if 0
                        {
                            // genere une direction p(l)= 1 / (2pi)
                            float cos_theta= rng.sample();
                            float sin_theta= std::sqrt(1 - cos_theta * cos_theta);
                            float phi= float(2*M_PI) * rng.sample();
                            
                            // construit les composantes de la direction l et transforme dans le repere de la scene
                            l= world( Vector(std::cos(phi) * sin_theta, std::sin(phi) * sin_theta, cos_theta) );
                            
                            // evalue la pdf
                            pdf= 1 / float(2*M_PI);
                        }
                        #else
                        {
                            // genere une direction p(l)= cos / (pi)
                            float cos_theta= std::sqrt( rng.sample() );
                            float sin_theta= std::sqrt(1 - cos_theta * cos_theta);
                            float phi= float(2*M_PI) * rng.sample();
                            
                            // construit les composantes de la direction l et transforme dans le repere de la scene
                            l= world( Vector(std::cos(phi) * sin_theta, std::sin(phi) * sin_theta, cos_theta) );
                            
                            // evalue la pdf
                            pdf= cos_theta / float(M_PI);
                        }
                        #endif
                        
                        Color sample;
                        if(Hit source= bvh.intersect( Ray(p+epsilon*pn, l) ))
                        {
                            Color Le= hit_material(source, scene).emission;
                            float cos_theta= std::max( float(0), dot( normalize(pn), normalize(l) ) );
                            
                            sample= diffuse * Le * cos_theta / pdf; 
                        }
                        
                        // accumule le sample,
                        // cf https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
                        n++;
                        Color Mn1= M;
                        M= Mn1 + (sample - Mn1) / float(n);
                        M2= M2 + (sample - Mn1) * (sample - M);
                    }
                }
            }
            
            emission= emission / float(AA);
            image(px, py)= Color(M + emission, 1);
            //~ image(px, py)= Color(M2 / float(n), 1);
        }
    }
    
    return image;
}
 
 
Color direct( const Point& p, const World& world, const GLTFScene& scene, const BVH& bvh, Sampler& rng, const std::vector<Source>& sources, const float sources_area )
{
    int id= int(sources.size()) -1;
    {
        float u= rng.sample();
        // construit et parcours la fonction de repartition
        float cdf= 0;
        for(unsigned i= 0; i < sources.size(); i++)
        {
            cdf+= sources[i].area / sources_area;
            if(u < cdf)
            {
                id= i;
                break;
            }
        }
    }
    
    // construit le point q sur la source selectionnee
    Point q= sources[id].sample(rng.sample(), rng.sample());
    float pdf= 1 / sources_area;
    
    Vector l= Vector(p, q);
    if(bvh.visible( Ray(p+epsilon*world.n, q+epsilon*sources[id].n) ))
    {
        float cos_theta= std::max(float(0), dot( normalize(world.n), normalize(l) ));
        float cos_theta_q= std::max(float(0), -dot( normalize(sources[id].n), normalize(l) ));
        
        return sources[id].emission * cos_theta * cos_theta_q / length2(l) / pdf; 
    }
    
    return Black();
}
 
Color indirect( const Point& p, const World& world, const GLTFScene& scene, const BVH& bvh, Sampler& rng, const std::vector<Source>& sources, const float sources_area )
{
    // genere une direction p(l)= cos / pi
    Vector l;
    float pdf;
    {
        float cos_theta= std::sqrt( rng.sample() );
        float sin_theta= std::sqrt(1 - cos_theta * cos_theta);
        float phi= float(2*M_PI) * rng.sample();
        
        // construit les composantes de la direction l et transforme dans le repere de la scene
        l= world( Vector(std::cos(phi) * sin_theta, std::sin(phi) * sin_theta, cos_theta) );
        
        // evalue la pdf
        pdf= cos_theta / float(M_PI);
    }
    
    // intersection
    Ray ray(p+epsilon*world.n, l);
    if(Hit qhit= bvh.intersect(ray))
    {
        // recuperer les infos sur l'intersection
        Point q= hit_position(qhit, ray);
        Vector qn= hit_normal(qhit, scene);
        
        if(dot(ray.d, qn) > 0)      // retourne la normale si elle n'est pas orientee vers l'origine du rayon...
            qn= -qn;
        World qworld(qn);
        
        float cos_theta= std::max(float(0), dot(world.n, l));
        Color qdiffuse= material_diffuse( hit_material(qhit, scene) );
        return qdiffuse * direct(q, qworld, scene, bvh, rng, sources, sources_area) * cos_theta / pdf;
    }
    
    return Black();
}
 
 
Image indirect( const Orbiter& camera, const GLTFScene& scene, const BVH& bvh, const int samples, const std::vector<Source>& sources )
{
    Image image(camera.width(), camera.height());
    
    const int AA= 1;
    const int N= samples;
    
    printf("render  %dx%d %d samples...\n", image.width(), image.height(), samples);
    
    // normalisation pdf source
    float total= 0;
    for(unsigned i= 0; i < sources.size(); i++)
        total+= sources[i].area;
    
    Transform model= Identity();
    Transform view= camera.view();
    Transform projection= camera.projection();
    Transform viewport= camera.viewport();
    Transform inv= Inverse(viewport *projection * view);        // image vers scene
    
    #pragma omp parallel for schedule(dynamic, 1)
    for(int py= 0; py < image.height(); py++)
    {
        std::random_device hwseed;
        Sampler rng( hwseed() );        
        
        for(int px= 0; px < image.width(); px++)
        {
            Color M;
            Color M2;
            int n= 0;
            
            Color emission;
            for(int pa= 0; pa < AA; pa++)
            {
                //~ float ax= rng.sample();
                //~ float ay= rng.sample();
                float ax= 0;
                float ay= 0;
                
                Point o= inv( Point(px +ax, py +ay, 0) );
                Point e= inv( Point(px +ax, py +ay, 1) );
                Vector d= Vector(o, e);
                
                Ray ray(o, d);
                if(Hit hit= bvh.intersect(ray))
                {
                    Point p= hit_position(hit, ray);
                    Vector pn= hit_normal(hit, scene);
                    Color diffuse= material_diffuse( hit_material(hit, scene) );
                    
                    emission= emission + hit_material(hit, scene).emission;
                    
                    if(dot(ray.d, pn) > 0)      // retourne la normale si elle n'est pas orientee vers l'origine du rayon...
                        pn= -pn;
                    
                    World world(pn);
                    for(int i= 0; i < N; i++)
                    {
                        Color d= diffuse * direct(p, world, scene, bvh, rng, sources, total);
                        Color i1= diffuse * indirect(p, world, scene, bvh, rng, sources, total);
                        Color sample= d + i1;
                        //~ Color sample= d;
                        
                        // accumule le sample,
                        // cf https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
                        n++;
                        Color Mn1= M;
                        M= Mn1 + (sample - Mn1) / float(n);
                        M2= M2 + (sample - Mn1) * (sample - M);
                    }
                }
            }
            
            emission= emission / float(AA);
            
            //~ color= color / float(N);
            //~ color= color / float(AA);
            
            //~ image(px, py)= Color(color + emission, 1);
            //~ image(px, py)= Color(M2 / float(n), 1);
            image(px, py)= Color(M + emission, 1);
        }
    }
    
    return image;
}
 
 
Color direct( const Point& p, const World& world, const GLTFScene& scene, const BVH& bvh, Sampler& rng )
{
    // genere une direction p(l)= cos / pi
    Vector l;
    float pdf;
    {
        float cos_theta= std::sqrt( rng.sample() );
        float sin_theta= std::sqrt(1 - cos_theta * cos_theta);
        float phi= float(2*M_PI) * rng.sample();
        
        // construit les composantes de la direction l et transforme dans le repere de la scene
        l= world( Vector(std::cos(phi) * sin_theta, std::sin(phi) * sin_theta, cos_theta) );
        
        // evalue la pdf
        pdf= cos_theta / float(M_PI);
    }
    
    // intersection
    Ray ray(p+epsilon*world.n, l);
    if(Hit qhit= bvh.intersect(ray))
    {
        // recuperer les infos sur l'intersection
        Vector qn= hit_normal(qhit, scene);
        
        Color qemission= Black();
        if(dot(l, qn) < 0)
            qemission= hit_material(qhit, scene).emission;        
        
        float cos_theta= std::max(float(0), dot(world.n, l));
        return qemission * cos_theta / pdf;
    }
    
    return Black();
}
 
Color indirect( const Point& p, const World& world, const GLTFScene& scene, const BVH& bvh, Sampler& rng )
{
    // genere une direction p(l)= cos / pi
    Vector l;
    float pdf;
    {
        float cos_theta= std::sqrt( rng.sample() );
        float sin_theta= std::sqrt(1 - cos_theta * cos_theta);
        float phi= float(2*M_PI) * rng.sample();
        
        // construit les composantes de la direction l et transforme dans le repere de la scene
        l= world( Vector(std::cos(phi) * sin_theta, std::sin(phi) * sin_theta, cos_theta) );
        
        // evalue la pdf
        pdf= cos_theta / float(M_PI);
    }
    
    // intersection
    Ray ray(p+epsilon*world.n, l);
    if(Hit qhit= bvh.intersect(ray))
    {
        // recuperer les infos sur l'intersection
        const GLTFMaterial& qmaterial= hit_material(qhit, scene);
        if(qmaterial.emission.max() > 0)
            return Black();
        
        Point q= hit_position(qhit, ray);
        Vector qn= hit_normal(qhit, scene);
        
        if(dot(ray.d, qn) > 0)      // retourne la normale si elle n'est pas orientee vers l'origine du rayon...
            qn= -qn;
        
        World qworld(qn);
        Color qdiffuse= material_diffuse( qmaterial );
        float cos_theta= std::max(float(0), dot(world.n, l));
        return qdiffuse * direct(q, qworld, scene, bvh, rng) * cos_theta / pdf;
    }
    
    return Black();
}
 
Image indirect( const Orbiter& camera, const GLTFScene& scene, const BVH& bvh, const int samples )
{
    Image image(camera.width(), camera.height());
    
    const int AA= 1;
    const int N= samples;
    
    printf("render %dx%d %d samples...\n", image.width(), image.height(), samples);
    
    Transform model= Identity();
    Transform view= camera.view();
    Transform projection= camera.projection();
    Transform viewport= camera.viewport();
    Transform inv= Inverse(viewport *projection * view);        // image vers scene
    
    #pragma omp parallel for schedule(dynamic, 1)
    for(int py= 0; py < image.height(); py++)
    {
        std::random_device hwseed;
        Sampler rng( hwseed() );        
        
        for(int px= 0; px < image.width(); px++)
        {
            Color color;
            Color emission;
            for(int pa= 0; pa < AA; pa++)
            {
                //~ float ax= rng.sample();
                //~ float ay= rng.sample();
                float ax= 0;
                float ay= 0;
                
                Point o= inv( Point(px +ax, py +ay, 0) );
                Point e= inv( Point(px +ax, py +ay, 1) );
                Vector d= Vector(o, e);
                
                Ray ray(o, d);
                if(Hit hit= bvh.intersect(ray))
                {
                    Point p= hit_position(hit, ray);
                    Vector pn= hit_normal(hit, scene);
                    Color diffuse= material_diffuse( hit_material(hit, scene) );
                    
                    if(dot(d, pn) < 0)
                        emission= emission + hit_material(hit, scene).emission;
                    
                    if(dot(ray.d, pn) > 0)      // retourne la normale si elle n'est pas orientee vers l'origine du rayon...
                        pn= -pn;
                    
                    World world(pn);
                    for(int i= 0; i < N; i++)
                    {
                        Color d= diffuse * direct(p, world, scene, bvh, rng);
                        Color i1= diffuse * indirect(p, world, scene, bvh, rng);
                        color= color + d + i1;
                    }
                }
            }
            
            emission= emission / float(AA);
            color= color / float(AA*N);
            
            image(px, py)= Color(color + emission, 1);
        }
    }
    
    return image;
}
 
Image ao( const Orbiter& camera, const GLTFScene& scene, const BVH& bvh, const int samples, const std::vector<Source>& sources )
{
    Image image(camera.width(), camera.height());
    
    const int AA= 1;
    const int N= samples;
    
    printf("render  %dx%d %d samples...\n", image.width(), image.height(), samples);
    
    // normalisation pdf source
    float total= 0;
    for(unsigned i= 0; i < sources.size(); i++)
        total+= sources[i].area;
    
    Transform model= Identity();
    Transform view= camera.view();
    Transform projection= camera.projection();
    Transform viewport= camera.viewport();
    Transform inv= Inverse(viewport *projection * view);        // image vers scene
    
    #pragma omp parallel for schedule(dynamic, 1)
    for(int py= 0; py < image.height(); py++)
    {
        std::random_device hwseed;
        Sampler rng( hwseed() );        
        
        for(int px= 0; px < image.width(); px++)
        {
            Color M;
            Color M2;
            int n= 0;
            
            Color emission;
            for(int pa= 0; pa < AA; pa++)
            {
                //~ float ax= rng.sample();
                //~ float ay= rng.sample();
                float ax= 0;
                float ay= 0;
                
                Point o= inv( Point(px +ax, py +ay, 0) );
                Point e= inv( Point(px +ax, py +ay, 1) );
                Vector d= Vector(o, e);
                
                Ray ray(o, d);
                if(Hit hit= bvh.intersect(ray))
                {
                    Point p= hit_position(hit, ray);
                    Vector pn= hit_normal(hit, scene);
                    Color diffuse= material_diffuse( hit_material(hit, scene) );
                    
                    if(dot(ray.d, pn) > 0)      // retourne la normale si elle n'est pas orientee vers l'origine du rayon...
                        pn= -pn;
                    
                    emission= emission + hit_material(hit, scene).emission;
                    
                    World world(pn);
                    for(int i= 0; i < N; i++)
                    {
                        Vector l;
                        float pdf;
                        #if 1
                        {
                            // genere une direction p(l)= 1 / (2pi)
                            float cos_theta= rng.sample();
                            float sin_theta= std::sqrt(1 - cos_theta * cos_theta);
                            float phi= float(2*M_PI) * rng.sample();
                            
                            // construit les composantes de la direction l et transforme dans le repere de la scene
                            l= world( Vector(std::cos(phi) * sin_theta, std::sin(phi) * sin_theta, cos_theta) );
                            
                            // evalue la pdf
                            pdf= 1 / float(2*M_PI);
                        }
                        #else
                        {
                            // genere une direction p(l)= cos / (pi)
                            float cos_theta= std::sqrt( rng.sample() );
                            float sin_theta= std::sqrt(1 - cos_theta * cos_theta);
                            float phi= float(2*M_PI) * rng.sample();
                            
                            // construit les composantes de la direction l et transforme dans le repere de la scene
                            l= world( Vector(std::cos(phi) * sin_theta, std::sin(phi) * sin_theta, cos_theta) );
                            
                            // evalue la pdf
                            pdf= cos_theta / float(M_PI);
                        }
                        #endif
                        
                        // verifie la visibilité de la lumière / du ciel dans la direction l
                        Color sample;
                        if(!bvh.intersect( Ray(p + epsilon*pn, l)))
                        {
                            // V(p, l) = 1
                            // il ne reste plus qu'à évaluer les autres termes...
                            
                            Color Li= Color(1);
                            float cos_theta= std::max(float(0), dot( normalize(pn), normalize(l) ));
                            sample= diffuse * Li * cos_theta / pdf;
                        }
                        
                        // accumule le sample,
                        // cf https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
                        n++;
                        Color Mn1= M;
                        M= Mn1 + (sample - Mn1) / float(n);
                        M2= M2 + (sample - Mn1) * (sample - M);                        
                    }
                }
            }
            
            emission= emission / float(AA);
            
            //~ color= color / float(N);
            //~ color= color / float(AA);
            //~ image(px, py)= Color(color + emission, 1);
            
            image(px, py)= Color(M2 / float(n), 1);
            //~ image(px, py)= Color(M + emission, 1);            
        }
    }
    
    return image;
}
 
int main( int argc, char **argv )
{
    const char *filename= "cornell.gltf";
    const char *orbiter_filename= "orbiter.txt";
    if(argc > 1) filename= argv[1];
    if(argc > 2) orbiter_filename= argv[2];
    
    Orbiter camera;
    if(camera.read_orbiter(orbiter_filename) < 0)
        return 1;   // pas de camera, pas d'image
    
    GLTFScene scene= read_gltf_scene(filename);
    if(scene.meshes.size() == 0)
        return 1;   // pas de geometrie, pas d'image
    
    //
    BVH bvh;
    std::vector<Source> sources;
    {
        std::vector<Triangle> triangles;
        
        for(unsigned node_id= 0; node_id < scene.nodes.size(); node_id++)
        {
            const GLTFNode& node= scene.nodes[node_id];
            
            const Transform& model= node.model;
            int mesh_id= node.mesh_index;
            
            const GLTFMesh& mesh= scene.meshes[mesh_id];
            for(unsigned primitive_id= 0; primitive_id < mesh.primitives.size(); primitive_id++)
            {
                const GLTFPrimitives& primitives= mesh.primitives[primitive_id];
                
                assert(primitives.material_index != -1);
                const GLTFMaterial& material= scene.materials[primitives.material_index];
                
                bool emission= (material.emission.max() > 0);
                for(unsigned i= 0; i +2 < primitives.indices.size(); i+= 3)
                {
                    // transforme les sommets dans le repere de la scene
                    Point a= model( Point(primitives.positions[primitives.indices[i]]) );
                    Point b= model( Point(primitives.positions[primitives.indices[i+1]]) );
                    Point c= model( Point(primitives.positions[primitives.indices[i+2]]) );
                    
                    // verifie que le triangle n'est pas degenere... oui ca arrive... sur robot.gltf, par exemple...
                    float area= length( cross( Vector(a, b), Vector(a, c) ) );
                    if(area > 0)
                    {
                        // indice du premier sommet du triangle dans la primitive
                        triangles.push_back( Triangle(a, b, c, node_id, mesh_id, primitive_id, i) );
                        
                        // si le triangle emet de la lumiere, ajouter une source
                        if(emission)
                            sources.push_back( Source(a, b, c, material.emission) );
                    }
                }
            }
        }
        
        //~ assert(sources.size()); // pas de lumiere, pas d'image
        assert(triangles.size());   // pas de triangles, pas d'image
        
        // construit le bvh avec les triangles
        bvh.build(triangles);
    }
    
    
    // c'est parti !!
    char tmp[1024];
    //~ for(int n : { 1, 4, 16, 64 })
    //~ for(int n : { 1, 4, 16, 64, 256, 1024, 4096 })
    //~ for(int n : { 4096 })
    for(int n : { 16384 })
    //~ for(int n : { 256, 1024 })
    {
        //~ Image image= indirect(camera, scene, bvh, n, sources);
        Image image= indirect(camera, scene, bvh, n);
        //~ Image image= direct(camera, scene, bvh, n, sources);
        //~ Image image= direct(camera, scene, bvh, n);
        //~ Image image= ao(camera, scene, bvh, n, sources);
        
        sprintf(tmp, "cornell_direct-%04d.hdr", n); // direct explicite
        //~ sprintf(tmp, "cornell_indirect1-%04d.hdr", n);  // direct explicite
        //~ sprintf(tmp, "cornell_indirect0-%04d.hdr", n);  // direct implicite
        
        //~ sprintf(tmp, "cornell_sources_var-%04d.hdr", n);
        //~ sprintf(tmp, "emission_sources-%04d.hdr", n);
        //~ sprintf(tmp, "ao_uniform-%04d.hdr", n);
        //~ sprintf(tmp, "ao_uniform_var-%04d.hdr", n);
        //~ sprintf(tmp, "ao_cos_var-%04d.hdr", n);
        //~ sprintf(tmp, "directions-%04d.hdr", n);
        //~ sprintf(tmp, "emission_directions_uniform_var-%04d.hdr", n);
        //~ sprintf(tmp, "emission_directions_cos_var-%04d.hdr", n);
        //~ sprintf(tmp, "emission_sources_var-%04d.hdr", n);
        //~ sprintf(tmp, "indirectM2-%04d.hdr", n);
        //~ sprintf(tmp, "indirectM21-%04d.hdr", n);
        write_image_hdr(image, tmp);
    }
    
    return 0;
}