Divide by value in Ddf::sample()

src/gui.h

@@ -21,7 +21,7 @@ class Gui: public RenderPlane, public Camera {
     SimpleCamera camera;
     cimg_library::CImg<float> processed;
     void updateDisplay();
-    float glare_cutoff = 1.01f;     // +0.01 to prevent preview render from glowing
+    float glare_cutoff = 1e+9f;     // +0.01 to prevent preview render from glowing
 public:
     Gui(const Camera& _camera);
     void work();

src/libddf/check_ddf.cpp

@@ -46,7 +46,7 @@ void mc_integral_and_max(const Ddf& ddf, float& ddf_integral, float& ddf_max){
 
         vec3 dir = vec3();
         while(dir == vec3())
-            dir = sph.sample();
+            dir = normalize(sph.sample());
 
         float value = ddf.value(dir);
         if(value > ddf_max)
@@ -129,7 +129,7 @@ bool check_ddf(const Ddf& ddf, bool strict_integral, size_t size_alpha, size_t s
     for(size_t i=0; i<N; ++i){
 
         ++total_tries;
-        vec3 vec = ddf.sample();
+        vec3 vec = normalize(ddf.sample());
         if(vec == vec3()){
             --i;
             continue;

src/libddf/ddf.cpp

@@ -63,7 +63,7 @@ vec3 SphericalDdf::sample() const {
     float alpha = acos(u1);
     float phi = 2*M_PI*u2;
     float r = sin(alpha);
-    return vec3(r*cos(phi), r*sin(phi), u1);
+    return vec3(r*cos(phi), r*sin(phi), u1)/float(0.25/M_PI);
 }
 float SphericalDdf::value( vec3 ) const {
     // TODO assert length = 1?
@@ -78,7 +78,7 @@ vec3 UpperHalfDdf::sample() const {
     float alpha = acos(u1);
     float phi = 2*M_PI*u2;
     float r = sin(alpha);
-    return vec3(r*cos(phi), r*sin(phi), u1);
+    return vec3(r*cos(phi), r*sin(phi), u1)/float(0.5/M_PI);
 }
 float UpperHalfDdf::value( vec3 arg ) const {
     // TODO assert length = 1?
@@ -97,7 +97,7 @@ vec3 CosineDdf::sample() const {
     float alpha = acos(cos_alpha);
     float phi = 2*M_PI*u2;
     float r = sin(alpha);
-    return vec3(r*cos(phi), r*sin(phi), cos_alpha);
+    return vec3(r*cos(phi), r*sin(phi), cos_alpha)/float(cos_alpha/M_PI);
 }
 float CosineDdf::value( vec3 arg ) const {
     // TODO assert length = 1?
@@ -141,15 +141,39 @@ vec3 UnionDdf::sample() const {
         return vec3();
     float r = randf();
     float acc = 0.0f;
+    size_t winning_i = 0;
     // TODO What's best used here as i?
     for(size_t i=0; i<components.size(); ++i){
         acc += weights[i];
         if(r<acc){
-            res = components[i]->sample();
+            assert(weights[i]!=0.0f);
+            res = components[i]->sample()/weights[i];
+            assert(isfinite(res.x) && isfinite(res.y) && isfinite(res.z));
+            winning_i = i;
             break;
         }
     }// for
     assert(r<acc);
+
+    // Case with 0,0,0 is used e.g. if no lights are accessible from given point
+    // TODO Try to eliminate this
+    if(res == vec3())
+        return res;
+
+    // now adjust value by considering other components' values
+    float value = 0.0f;
+    for(size_t i=0; i<components.size(); ++i){
+        // use generic value() for others
+        if(i!=winning_i)
+            value += weights[i] * components[i]->value(res);
+        // and use actually returned value for winning_i
+        else
+            value += 1.0f/length(res);
+    }// for
+    assert(isfinite(value) && value != 0.0f);
+
+    res = normalize(res) / value;
+
     return res;
 }
 

src/libddf/ddf.h

@@ -14,7 +14,7 @@ struct Ddf {
     // Sample a signle vec3 value. Always succeedes
     virtual glm::vec3 sample() const = 0;
 
-    // Get value in direction of arg. Will return NaN if singular
+    // Get value in direction of arg. For singulars returns 0
     virtual float value( glm::vec3 arg ) const = 0;
 
     //static std::atomic_size_t object_counter;

src/libddf/ddf_detail.h

@@ -63,8 +63,8 @@ class MirrorDdf: public Ddf {
     virtual glm::vec3 sample() const override {
         return glm::vec3(0.0f, 0.0f, 1.0f);
     }
-    virtual float value( glm::vec3 arg ) const override {
-        return std::numeric_limits<float>::quiet_NaN();
+    virtual float value( glm::vec3 ) const override {
+        return 0.0f;
     }
 };
 

src/lighting/lighting.cpp

@@ -55,7 +55,11 @@ glm::vec3 DdfFromLight::sample() const {
     if(cosinus < 1e-5f)               // if facing back
         return vec3();
     Lighting::last_sample = inter;  // HACK TODO how to make it accessible in a cleaner way?
-    return dir;
+
+    float decay = dot(inter.position-origin, inter.position-origin);
+    float value = decay/cosinus/light->area;
+
+    return dir/value;
 }
 
 float DdfFromLight::value( glm::vec3 direction ) const {

src/main.cpp

@@ -22,6 +22,14 @@
 using namespace glm;
 using namespace std;
 
+size_t n_rays = 16;
+size_t depth_max = 4;
+
+class StatsNode;
+float ray_power_preview(const Geometry& geometry, const Lighting& lighting, vec3 origin, vec3 direction, size_t, int, StatsNode* stats);
+float ray_power_recursive(const Geometry& geometry, const Lighting& lighting, vec3 origin, vec3 direction, size_t depth, int n_rays, StatsNode* stats);
+auto ray_power = ray_power_recursive;
+
 // note: root will manage memory for all children
 class StatsNode {
 public:
@@ -48,10 +56,6 @@ class StatsNode {
 
 boost::pool<> StatsNode::pool(sizeof(StatsNode));
 
-float ray_power_preview(const Geometry& geometry, const Lighting& lighting, vec3 origin, vec3 direction, size_t, int, StatsNode* stats);
-float ray_power_recursive(const Geometry& geometry, const Lighting& lighting, vec3 origin, vec3 direction, size_t depth, int n_rays, StatsNode* stats);
-auto ray_power = ray_power_recursive;
-
 float ray_power_preview(const Geometry& geometry, const Lighting& lighting, vec3 origin, vec3 direction, size_t, int, StatsNode* stats){
 
     // Stats: origin
@@ -91,9 +95,6 @@ float ray_power_preview(const Geometry& geometry, const Lighting& lighting, vec3
     return res;
 }
 
-size_t n_rays = 16;
-size_t depth_max = 4;
-
 // TODO light_hint is not very good solution!
 float ray_power_recursive(const Geometry& geometry, const Lighting& lighting, vec3 origin, vec3 direction, size_t depth, int n_rays, StatsNode* stats){
 
@@ -136,10 +137,11 @@ float ray_power_recursive(const Geometry& geometry, const Lighting& lighting, ve
     //DEBUG for geometry debugging
     //return si->position.y+1.0f;
 
+    unique_ptr<Ddf> light_ddf = lighting.distributionInPoint(si->position);
+
     // TODO better solution?
     Ddf* sdf_tmp = si->sdf.get();
 
-    unique_ptr<Ddf> light_ddf = lighting.distributionInPoint(si->position);
     unique_ptr<Ddf> mix_ddf = unite(move(light_ddf), 1.0f, move(si->sdf), 1.0f);
 
     vec3 new_direction;
@@ -149,18 +151,12 @@ float ray_power_recursive(const Geometry& geometry, const Lighting& lighting, ve
     for(size_t i=0; i<n_rays; ++i){
 
         new_direction = mix_ddf->sample();
-        float mix_val = mix_ddf->value(new_direction);
-
-        // if(i<n_rays/2)
-        //     new_direction = si->sdf->sample();
-        // else
-        //     new_direction = light_ddf->sample();
-
-        // float mix_val = 0.5f*si->sdf->value( new_direction ) + 0.5f*light_ddf->value( new_direction );
+        float mix_val = 1.0f/length(new_direction);
 
-        if( new_direction == vec3() ){
+        if(new_direction == vec3())
             continue;
-        }
+
+        new_direction = normalize(new_direction);
 
         // Stats: false (miss)
         StatsNode* child_stats = new StatsNode();
@@ -247,6 +243,7 @@ int main(int argc, char** argv){
     GridRenderPlane r_plane(640, 640);
 
     //Scene scene = make_scene_lit_corner();
+    //Scene scene = make_scene_square_lit_by_square();
     Scene scene = make_scene_box();
     auto gui = make_shared<Gui>(*scene.camera);
     scene.camera = gui;     // replace camera with more interactive one