tuto_bvh2.cpp

Go to the documentation of this file.

 
 
#include <algorithm>
#include <vector>
#include <cfloat>
 
#include "vec.h"
#include "mat.h"
#include "color.h"
#include "image.h"
#include "image_io.h"
#include "orbiter.h"
#include "mesh.h"
#include "wavefront.h"
 
 
// rayon 
struct Ray
{
    Point o;            // origine
    float pad;
    Vector d;           // direction
    float tmax;
    
    Ray( const Point& _o, const Point& _e ) :  o(_o), d(Vector(_o, _e)), tmax(1) {}
    Ray( const Point& _o, const Vector& _d ) :  o(_o), d(_d), tmax(FLT_MAX) {}
    Ray( const Point& _o, const Vector& _d, const float _tmax ) :  o(_o), d(_d), tmax(_tmax) {}
};
 
// 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 instance_id;    // indice de l'instance
    int triangle_id;    // indice du triangle dans le mesh
    
    Hit( ) : t(FLT_MAX), u(), v(), instance_id(-1), triangle_id(-1) {}
    Hit( const Ray& ray ) : t(ray.tmax), u(), v(), instance_id(-1), triangle_id(-1) {}
    Hit( const float _t, const float _u, const float _v, const int _id ) : t(_t), u(_u), v(_v), instance_id(-1), triangle_id(_id) {}
    
    operator bool ( ) { return (triangle_id != -1); }
};
 
// 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; }
};
 
 
// 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(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); }
    
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;
    }
    
    // intersection et parcours simple
    void intersect( const int index, const Ray& ray, const Vector& invd, Hit& hit ) const
    {
        const Node& node= nodes[index];
        if(node.bounds.intersect(ray, invd, hit.t))
        {
            if(node.leaf())
            {
                for(int i= node.leaf_begin(); i < node.leaf_end(); i++)
                    if(Hit h= primitives[i].intersect(ray, hit.t))
                        hit= h;
            }
            else // if(node.internal())
            {
                intersect(node.internal_left(), ray, invd, hit);
                intersect(node.internal_right(), ray, invd, hit);
            }
        }
    }
};
 
 
// 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 id;
    
    Triangle( const TriangleData& data, const int _id ) : p(data.a), e1(Vector(data.a, data.b)), e2(Vector(data.a, data.c)), 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, id);
    }
    
    BBox bounds( ) const 
    {
        BBox box(p);
        return box.insert(p+e1).insert(p+e2);
    }
};
 
typedef BVHT<Triangle> BVH;
typedef BVHT<Triangle> BLAS;
 
 
// instance pour le bvh, cf fonctions bounds() et intersect()
struct Instance
{
    Transform object_transform;
    BBox world_bounds;
    BVH *object_bvh;
    int instance_id;
    
    Instance( const BBox& bounds, const Transform& model, BVH& bvh, const int id ) :  
        object_transform(Inverse(model)), world_bounds(transform(bounds, model)), 
        object_bvh(&bvh),
        instance_id(id)
    {}
    
    BBox bounds( ) const { return world_bounds; }
    
    Hit intersect( const Ray &ray, const float htmax ) const
    {
        // transforme le rayon
        Ray object_ray(object_transform(ray.o), object_transform(ray.d), htmax);
        // et intersection dans le bvh de l'objet instancie...
        
        Hit hit= object_bvh->intersect(object_ray, htmax);
        if(hit)
            // si intersection, stocker aussi l'indice de l'instance, cf retrouver la transformation et la matiere associee au mesh/triangle...
            hit.instance_id= instance_id;
        
        return hit;
    }
    
protected:
    BBox transform( const BBox& bbox, const Transform& m )
    {
        BBox bounds= BBox( m(bbox.pmin) );
        // enumere les sommets de la bbox 
        for(unsigned i= 1; i < 8; i++)
        {
            // chaque sommet de la bbox est soit pmin soit pmax sur chaque axe...
            Point p= bbox.pmin; 
            if(i & 1) p.x= bbox.pmax.x;
            if(i & 2) p.y= bbox.pmax.y;
            if(i & 4) p.z= bbox.pmax.z;
            
            // transforme le sommet de l'englobant 
            bounds.insert( m(p) );
        }
        
        return bounds;
    }
};
 
typedef BVHT<Instance> TLAS;
 
 
int main( int argc, char **argv )
{
    const char *mesh_filename= "data/robot.obj";
    const char *orbiter_filename= nullptr;
    
    if(argc > 1) mesh_filename= argv[1];
    if(argc > 2) orbiter_filename= argv[2];
    
    Mesh mesh= read_mesh(mesh_filename);
    if(mesh.triangle_count() == 0)
        return 1;
    
    BBox mesh_bounds;
    mesh.bounds(mesh_bounds.pmin, mesh_bounds.pmax);
    
    // construit le bvh de l'objet
    BVH bvh;
    {
        // recupere les triangles du mesh
        std::vector<Triangle> triangles;
        for(int i= 0; i <mesh.triangle_count(); i++)
        {
            TriangleData data= mesh.triangle(i);
            triangles.push_back( Triangle(data, triangles.size()) );
        }
        
        // construit le bvh...
        bvh.build(triangles);
    }
    
    // instancie l'objet
    TLAS top_bvh;
    {
        std::vector<Instance> instances;
        for(int i= -2; i <= 2; i++)
        {
            Transform model= Translation(i*10, 0, 0);
            instances.push_back( Instance(mesh_bounds, model, bvh, instances.size()) );
        }
        
        // construit le bvh...
        top_bvh.build(instances);
    }
    
    // regle la camera
    Orbiter camera;
    if(orbiter_filename == nullptr || camera.read_orbiter(orbiter_filename) < 0)
        // les objets sont instancies sur une ligne, agrandir la zone observee par la camera...
        camera.lookat(mesh_bounds.pmin * 5, mesh_bounds.pmax * 5);      
    
    // image resultat
    Image image(1024, 640);
    
    // transformations
    Transform model= Identity();
    Transform view= camera.view();
    Transform projection= camera.projection();
    Transform viewport= Viewport(image.width(), image.height());
    Transform inv= Inverse(viewport * projection * view * model);
    
    // calcule l'image en parallele avec openMP
#pragma omp parallel for 
    for(int y= 0; y < image.height(); y++)
    for(int x= 0; x < image.width(); x++)
    {
        // genere le rayon pour le pixel x,y
        Point o= inv( Point(x, y, 0) ); // origine
        Point e= inv( Point(x, y, 1) ); // extremite
        Ray ray(o, e);
        
        // intersections !
        if(Hit hit= top_bvh.intersect(ray))
            image(x, y)= Red(); // touche ! 
    }
    
    write_image(image, "render.png");
    return 0;
}