diff --git a/code/projects/path_sobol.cpp b/code/projects/path_sobol.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..41cb35b7c564ac82175b29e7402b6396df3dc5b7
--- /dev/null
+++ b/code/projects/path_sobol.cpp
@@ -0,0 +1,380 @@
+
+
+#include <omp.h>
+#include <chrono>
+#include <random>
+#include <stdint.h>
+
+#include "options.h"
+
+#include "vec.h"
+#include "mat.h"
+
+#include "color.h"
+#include "framebuffer.h"
+#include "image.h"
+#include "image_io.h"
+
+#include "gltf/gltf.h"
+#include "gltf/scene.h"
+#include "orbiter.h"
+
+#include "sampler.h"
+#include "sampler_sobol.h"
+
+// utilise sobol' global
+typedef Sobol Sampler;
+
+Options options;
+
+
+Color direct( const Point& o, const Point& p, const Brdf& pbrdf, const Scene& scene, Sampler& rng )
+{
+ Color color;
+
+ /* strategie 1 : genere un point sur une source
+ */
+ const Source& source= scene.sources.sample( rng.sample() );
+ Point s= source.sample( rng.sample(), rng.sample() );
+ Vector sn= source.n;
+
+ Vector l= normalize( Vector(p, s) );
+ float cos_theta= std::max(float(0), dot(pbrdf.world.n, l));
+ float cos_theta_s= std::max(float(0), dot(sn, -l));
+ float G= cos_theta_s / distance2(p, s);
+ float pdf= source.pdf(s) * scene.sources.pdf(source);
+
+ if(pdf * cos_theta * cos_theta_s > 0
+ && scene.visible( offset_point(p, pbrdf.world.n), offset_point(s, sn) ))
+ {
+ Color fr= pbrdf.f(l, normalize(o - p));
+ color= color + source.emission * G * fr * cos_theta / pdf;
+ }
+
+ /* strategie 2 : genere une direction en fonction de la matiere
+
+ Vector l= pbrdf.sample( rng.sample(), rng.sample(), rng.sample(), normalize(o - p));
+ float pdf= pbrdf.pdf( l, normalize(o - p) );
+ float cos_theta= std::max(float(0), dot(pbrdf.world.n, l));
+ if(pdf * cos_theta > 0
+ && scene.visible(offset_point(p, pbrdf.world.n), l))
+ {
+ Color fr= pbrdf.f(l, normalize(o - p));
+ color= color + Color(0.02) * fr * cos_theta / pdf;
+ }
+ */
+
+ return color;
+}
+
+// mis multi-sample, evalue les 2 integrales + ponderation de chaque echantillon
+Color direct_mis( const Point& o, const Point& p, const Brdf& pbrdf, Hit& qhit, const Scene& scene, Sampler& rng )
+{
+ Color color;
+
+ // strategie 1 : selectionne une source + un point sur la source
+ {
+ const Source& source= scene.sources.sample( rng.sample() );
+ Point s= source.sample( rng.sample(), rng.sample() );
+
+ Vector sn= source.n;
+ Vector l= normalize( Vector(p, s) );
+ float cos_theta= std::max(float(0), dot(pbrdf.world.n, l));
+ float cos_theta_s= std::max(float(0), dot(sn, -l));
+ float G= cos_theta_s / distance2(p, s);
+ if(G * cos_theta > 0)
+ {
+ float pdf1= source.pdf(s) * scene.sources.pdf(source);
+ if(pdf1 > 0)
+ {
+ if(scene.visible( offset_point(p, pbrdf.world.n), offset_point(s, sn) ) )
+ {
+ float pdf2= pbrdf.pdf(l, normalize( Vector(p, o) )) * G;
+ float w= (pdf1*pdf1) / (pdf1*pdf1 + pdf2*pdf2); // pdf1 et pdf2 sur les aires, power heuristic
+
+ Color fr= pbrdf.f(l, normalize(o - p));
+ color= color + w * source.emission * fr * G * cos_theta / pdf1;
+ }
+ }
+ }
+ }
+
+ // strategie 2 : selectionne un terme de la brdf et genere une direcion
+ {
+ Vector l= pbrdf.sample( rng.sample(), rng.sample(), rng.sample(), normalize(o - p) );
+ float pdf2= pbrdf.pdf(l, normalize(o - p) );
+ if(pdf2 > 0)
+ {
+ Ray ray( offset_point(p, pbrdf.world.n), l );
+ if(Hit hit= scene.intersect(ray))
+ {
+ Color emission= scene.material(hit).emission;
+ if(emission.max() > 0)
+ {
+ Point s= ray.o + ray.d * hit.t;
+ Vector sn= scene.normal(hit);
+ float cos_theta= std::max(float(0), dot(pbrdf.world.n, l));
+ float cos_theta_s= std::max(float(0), dot(sn, -l));
+ float G= cos_theta_s / distance2(p, s);
+ if(G * cos_theta > 0)
+ {
+ pdf2= pdf2 * G;
+ float pdf1= 1 / scene.sources.area;
+ float w= (pdf2*pdf2) / (pdf1*pdf1 + pdf2*pdf2); // pdf1 et pdf2 sur les aires, power heuristic
+
+ Color fr= pbrdf.f(l, normalize(o - p));
+ color= color + w * emission * fr * G * cos_theta / pdf2;
+ }
+ }
+ else
+ // conserve le point pour prolonger le chemin...
+ qhit= hit;
+ }
+ }
+ }
+
+ return color;
+}
+
+
+Color path( const int depth, const Point& o, const Point& p, const Brdf& pbrdf, const Scene& scene, Sampler& rng )
+{
+ Hit qhit;
+
+ //~ Color color= direct(o, p, pbrdf, scene, rng );
+ Color color= direct_mis(o, p, pbrdf, qhit, scene, rng );
+ if(depth == 0)
+ return color;
+
+ // recuperer l'intersection mis du direct... cf direct_mis, elimine un rayon par rebond
+ if(qhit)
+ {
+ Point q= scene.position(qhit);
+ Vector l= normalize(q - p);
+ float pdf= pbrdf.pdf(l, normalize(o - p));
+ if(pdf > 0) // ne devrait pas arriver...
+ {
+ Vector qn= scene.normal(qhit);
+ Brdf qbrdf= scene.brdf(qhit, qn);
+
+ float cos_theta= std::max(float(0), dot(pbrdf.world.n, l));
+ Color fr= pbrdf.f(l, normalize(o - p));
+ return color + path(depth -1, p, q, qbrdf, scene, rng) * fr * cos_theta / pdf;
+ }
+ }
+
+ /* sinon genere une direction pour prolonger le chemin...
+ Vector l= pbrdf.sample( rng.sample(), rng.sample(), rng.sample(), normalize(o - p) );
+ float pdf= pbrdf.pdf(l, normalize(Vector(p, o)));
+ if(pdf > 0)
+ {
+ Ray ray( offset_point(p, pbrdf.world.n), l );
+ if(Hit hit= scene.intersect(ray))
+ {
+ Point q= scene.position(hit, ray);
+ Vector qn= scene.normal(hit);
+ Brdf qbrdf= scene.brdf(hit, qn);
+
+ float cos_theta= std::max(float(0), dot(pbrdf.world.n, l));
+ Color fr= pbrdf.f(l, normalize(o - p));
+ return color + path(depth -1, p, q, qbrdf, scene, rng) * fr * cos_theta / pdf;
+ }
+ }
+ */
+
+ return color;
+}
+
+
+
+struct alpha_filter
+{
+ const Scene& scene;
+ FCRNG& rng;
+
+ alpha_filter( const Scene& _scene, FCRNG& _rng ) : scene(_scene), rng(_rng) {}
+
+ bool operator() ( const Hit& hit ) const
+ {
+ const GLTFMaterial& material= scene.material(hit);
+ if(material.color_texture == -1)
+ return true;
+ if(scene.textures[material.color_texture].opaque)
+ return true;
+
+ vec3 uv_lod= scene.texcoords_lod(hit);
+ Color diffuse= scene.textures[material.color_texture].sample(uv_lod.x, uv_lod.y, uv_lod.z);
+
+ return (rng.sample() <= diffuse.a);
+ }
+};
+
+
+
+int main( const int argc, const char **argv )
+{
+ options= Options(argv, argv + argc);
+
+ Scene scene= read_scene(options.filename);
+ if(scene.gltf.nodes.empty())
+ return 1; // pas de geometrie, pas d'image
+
+ // transformations
+ Transform view;
+ Transform projection;
+ float aspect= 1;
+
+ if(options.orbiter_filename)
+ {
+ Orbiter camera;
+ if(camera.read_orbiter(options.orbiter_filename) == false)
+ return 1; // pas de camera, pas d'image
+
+ view= camera.view();
+ projection= camera.projection();
+ aspect= camera.aspect();
+ }
+ else if(scene.gltf.cameras.size())
+ {
+ view= scene.camera.view;
+ projection= scene.camera.projection;
+ aspect= scene.camera.aspect;
+ }
+ else
+ {
+ printf("[error] no glTF camera...\n");
+ return 1; // pas de camera, pas d'image
+ }
+
+ if(scene.sources.size() == 0)
+ {
+ printf("[error] no emissive triangles / light sources...\n");
+ return 1;
+ }
+
+
+ // proportions largeur/hauteur de la camera gltf
+ int samples= options.samples;
+ int width= 1024;
+ int height= width / aspect;
+ Framebuffer image(width, height);
+ Image aov_normal(width, height);
+ Image aov_albedo(width, height);
+
+ // charge les textures + parametres de filtrage en fonction de l'image
+ read_textures(options.filename, scene, width, height);
+ //~ read_textures(options.filename, scene, 256, 256);
+
+
+ //
+ Transform viewport= Viewport(image.width(), image.height());
+ Transform inv= Inverse(viewport * projection * view);
+
+ unsigned maxdims= 0;
+ unsigned lines= 0;
+ auto cpu_start= std::chrono::high_resolution_clock::now();
+
+ // determine les dimensions d'une image carre / puissance de 2
+ int pow2= 1;
+ while(pow2 < width && pow2 < height)
+ pow2*= 2;
+
+ printf("image %dx%d\n", width, height);
+ printf("square image %dx%d\n", pow2, pow2);
+
+ // c'est parti ! parcours tous samples de la sequence de sobol
+ const size_t n= width*height*samples;
+
+#pragma omp parallel for schedule(dynamic, 4096)
+ for(size_t i= 0; i < n; i++)
+ {
+ Sobol rng;
+ rng.index(i);
+ // les dimensions 0, 1 de la sequence de Sobol' choisissent le pixel
+ float x= rng.sample();
+ float y= rng.sample();
+
+ // rejette les echantillons en dehors de l'image
+ int py= y * pow2;
+ int px= x * pow2;
+
+ if(py >= height)
+ continue;
+ if(px >= width)
+ continue;
+
+ // generer le rayon
+ x= x * pow2;
+ y= y * pow2;
+ Point origine= inv(Point(x, y, 0));
+ Point extremite= inv(Point(x, y, 1));
+ Ray ray(origine, extremite);
+
+ Color color;
+ Vector normal;
+ Color albedo;
+
+ // todo : pour initialiser le rng d'alpha_filter, il faut recalculer le seed du pixel (px, py), cf path.cpp
+
+ // calculer les intersections avec tous les triangles
+ if(Hit hit= scene.bvh.intersect(ray))
+ //~ if(Hit hit= scene.bvh.intersect(ray, alpha_filter(scene, alpha_rng)))
+ {
+ Point p= scene.position(hit, ray);
+ Vector pn= scene.normal(hit);
+
+ // matiere emissive, si la source est orientee vers la camera
+ if(dot(pn, ray.d) < 0)
+ {
+ const GLTFMaterial& material= scene.material(hit);
+ color= color + material.emission;
+ // suppose que les normales des sources sont correctes...
+ }
+
+ if(dot(pn, ray.d) > 0)
+ // mais ce n'est pas toujours le cas pour le reste de la geometrie...
+ pn= -pn;
+
+ Brdf brdf= scene.brdf(hit, pn);
+ //~ color= color + direct(origine, p, brdf, scene, rng);
+ color= color + path(options.depth, origine, p, brdf, scene, rng);
+
+ normal= normal + pn;
+ albedo= albedo + brdf.f(pn, normalize(ray.o - p));
+ }
+
+ image.splat(px, py, color); // accumule l'echantillon
+
+ aov_normal(px, py)= Color(normal / float(samples), 1); // force une couleur opaque...
+ aov_albedo(px, py)= Color(albedo / float(samples), 1); // force une couleur opaque...
+
+ // #pragma omp compare // ne marche pas sur visual...
+ maxdims= std::max(rng.d, maxdims);
+ }
+
+ printf("sampler dimensions %u\n", maxdims);
+
+ auto cpu_stop= std::chrono::high_resolution_clock::now();
+ auto cpu_time= std::chrono::duration_cast<std::chrono::milliseconds>(cpu_stop - cpu_start).count();
+ printf("cpu %ds %03dms\n", int(cpu_time / 1000), int(cpu_time % 1000));
+
+ char tmp[1024];
+ //~ sprintf(tmp, "%s-%05u-%04u.png", options.prefix, unsigned(cpu_time), samples);
+ sprintf(tmp, "%s-%04u.png", options.prefix, samples);
+ printf("writing '%s'...\n", tmp);
+ write_image_preview(image.flatten(), tmp);
+
+ //~ sprintf(tmp, "%s-%05u-%04u.pfm", options.prefix, unsigned(cpu_time), samples);
+ sprintf(tmp, "%s-%04u.pfm", options.prefix, samples);
+ printf("writing '%s'...\n", tmp);
+ write_image_pfm(image.flatten(), tmp);
+
+ sprintf(tmp, "%s-normals-%04u.pfm", options.prefix, samples);
+ write_image_pfm(aov_normal, tmp);
+
+ sprintf(tmp, "%s-albedo-%04u.pfm", options.prefix, samples);
+ write_image_pfm(aov_albedo, tmp);
+ printf("\n");
+
+ return 0;
+}