Patch for textured area lights in snogglethorpe/snogray [2013-08-16]

tex-light-20130816-0.patch

diff --git a/light/surface-light.cc b/light/surface-light.cc
index ddd3ac2..0270862 100644
--- a/light/surface-light.cc
+++ b/light/surface-light.cc
@@ -1,3 +1,5 @@
+#include "image/image.h"
+#include "util/radical-inverse.h"
 // surface-light.cc -- General-purpose area light
 //
 //  Copyright (C) 2010, 2013  Miles Bader <miles@gnu.org>
@@ -23,17 +25,82 @@
 using namespace snogray;
 
 
+
+// Surface texture sampling
+
+static Hist2d
+sample_texture_over_surface (const TexVal<Color> &tex,
+			     const Surface::Sampler &sampler)
+{
+  unsigned w = 64, h = 64;
+  Hist2d hist (w, h);
+
+  float u_step_size = 1.f / w,  v_step_size = 1.f / h;
+  float u = u_step_size / 2;
+  for (unsigned x = 0; x < w; x++)
+    {
+      float v = v_step_size / 2;
+      for (unsigned y = 0; y < h; y++)
+	{
+	  Surface::Sampler::AreaSample samp = sampler.sample (UV (u, v));
+	  TexCoords tex_coords (samp.pos, samp.tex_coords_uv);
+	  Color val = tex.eval (tex_coords);
+	  hist.add (x, y, val.intensity ());
+	  v += v_step_size;
+	}
+      u += u_step_size;
+    }
+
+  Image dump_img (w, h);
+  for (unsigned x = 0; x < w; x++)
+    for (unsigned y = 0; y < h; y++)
+      dump_img (x, y) = hist (x, y);
+  dump_img.save (",dump.png");
+
+  Image dump_img2 (w, h);
+  Hist2dDist dist (hist);
+  for (unsigned x = 0; x < w; x++)
+    for (unsigned y = 0; y < h; y++)
+      dump_img2 (x, y) = Color (0);
+  for (unsigned i = 0; i < 1000; i++)
+    {
+      UV param (radical_inverse (i+1,2), radical_inverse(i+1,3));
+	
+      float pdf;
+      UV pos = dist.sample (param, pdf);
+
+      unsigned col = min (unsigned (pos.u * w), w-1);
+      unsigned row = min (unsigned (pos.v * h), h-1);
+
+      Color old = dump_img2 (col, row);
+      dump_img2 (col, row) = old + 0.1;
+    }
+  dump_img2.save (",dump2.png");
+
+  return hist;
+}
+
+
+// ctor 
+
 SurfaceLight::SurfaceLight (const Surface &surface,
 			    const TexVal<Color> &_intensity)
-  : sampler (surface.make_sampler ()), intensity (_intensity.default_val)
+  : sampler (surface.make_sampler ()), intensity (_intensity)
 {
   if (! sampler)
     throw std::runtime_error
       ("Surface cannot be used as a light");
 
+  // For textured surfaces, we sample the texture over the surface a
+  // somewhat low resolution.  The resulting 2D intensity map will
+  // then be sampled for lighting.
+  //
   if (_intensity.tex)
-    throw std::runtime_error
-      ("textured intensity not supported by SurfaceLight");
+    {
+      Hist2d sampled_tex = sample_texture_over_surface (_intensity, *sampler);
+
+      intensity_dist.set_histogram (sampled_tex);
+    }
 }
 
 
@@ -46,11 +113,49 @@ SurfaceLight::SurfaceLight (const Surface &surface,
 Light::Sample
 SurfaceLight::sample (const Intersect &isec, const UV &param) const
 {
-  // Sample the surface, in world-space.
+  // Viewpoint from which we will be sampling.
   //
-  Pos ipos = isec.normal_frame.origin;
-  Surface::Sampler::AngularSample samp
-    = sampler->sample_from_viewpoint (ipos, param);
+  Pos viewpoint = isec.normal_frame.origin;
+
+  // Sample we will create.
+  //
+  Surface::Sampler::AngularSample samp;
+
+  // Sample the surface.  Surfaces with a textured intensity need to
+  // use a more complex method.
+  //
+  if (intensity_dist.empty ())
+    {
+      // Sample a constant-intensity surface, using the possibly more
+      // efficient viewpoint-dependent
+      // Surface::Sampler::sample_from_viewpoint.
+      //
+      samp = sampler->sample_from_viewpoint (viewpoint, param);
+    }
+  else
+    {
+      // Purturb PARAM to follow the intensity distribution, so that a
+      // uniform PARM is more likely to end up in bright areas.
+      //
+      float tex_pdf;
+      UV purturbed_param = intensity_dist.sample (param, tex_pdf);
+
+      // Use the purturbed parameter to sample the surface.
+      //
+      // Because INTENSITY_DIST is viewpoint-independent, we can't use
+      // Surface::Sampler::sample_from_viewpoint, so instead we just
+      // use Surface::Sampler::sample (which might be somewhat less
+      // efficient), and convert the resulting viewpoint-independent
+      // sample to a viewpoint-dependent angular sample.
+      //
+      samp
+	= Surface::Sampler::AngularSample (sampler->sample (purturbed_param),
+					   viewpoint);
+
+      // Reflect the PDF of the intensity distribution in the sample PDF.
+      //
+      samp.pdf *= tex_pdf;
+    }
 
   if (samp.pdf > 0)
     {
@@ -61,7 +166,8 @@ SurfaceLight::sample (const Intersect &isec, const UV &param) const
       // Only process samples which are in front of ISEC.
       //
       if (dir.z > 0)
-	return Sample (intensity, samp.pdf, dir, samp.dist);
+	return Sample (intensity.eval (samp.tex_coords),
+		       samp.pdf, dir, samp.dist);
     }
 
   return Sample ();
@@ -75,9 +181,29 @@ SurfaceLight::sample (const Intersect &isec, const UV &param) const
 Light::FreeSample
 SurfaceLight::sample (const UV &param, const UV &dir_param) const
 {
+  // A copy of PARAM possibly purturbed to follow the light's intensity
+  // distribution.
+  //
+  UV purturbed_param = param;
+
+  // Any PDF modification resulting from said purturbation, or 1 if
+  // there was none.
+  //
+  float tex_pdf = 1;
+
+  // Maybe do extra calculations for textured lights.
+  //
+  if (! intensity_dist.empty ())
+    {
+      // Purturb PARAM to follow the intensity distribution, so that a
+      // uniform PARM is more likely to end up in bright areas.
+      //
+      purturbed_param = intensity_dist.sample (param, tex_pdf);
+    }
+
   // Sample the surface, in world-space.
   //
-  Surface::Sampler::AreaSample samp = sampler->sample (param);
+  Surface::Sampler::AreaSample samp = sampler->sample (purturbed_param);
 
   // Choose a direction in SAMP's normal-frame-of-reference according
   // to DIR_PARAM.
@@ -99,9 +225,10 @@ SurfaceLight::sample (const UV &param, const UV &dir_param) const
   // cosine terms cancel out, and we can just use POS_PDF / pi
   // instead.
   //
-  float pdf = samp.pdf * INV_PIf;
+  float pdf = samp.pdf * INV_PIf * tex_pdf;
 
-  return FreeSample (intensity, pdf, samp.pos, dir);
+  return FreeSample (intensity.eval (TexCoords (samp.pos, samp.tex_coords_uv)),
+		     pdf, samp.pos, dir);
 }
 
 
@@ -123,7 +250,7 @@ SurfaceLight::eval (const Intersect &isec, const Vec &dir) const
     = sampler->eval_from_viewpoint (ipos, wdir);
 
   if (samp.pdf > 0)
-    return Value (intensity, samp.pdf, samp.dist);
+    return Value (intensity.eval (samp.tex_coords), samp.pdf, samp.dist);
 
   return Value ();
 }
diff --git a/light/surface-light.h b/light/surface-light.h
index df906cd..c8b51b8 100644
--- a/light/surface-light.h
+++ b/light/surface-light.h
@@ -16,8 +16,10 @@
 #include "util/unique-ptr.h"
 #include "color/color.h"
 #include "texture/tex.h"
-#include "light.h"
 #include "surface/surface.h"
+#include "geometry/hist-2d-dist.h"
+
+#include "light.h"
 
 
 namespace snogray {
@@ -53,7 +55,11 @@ private:
 
   // Radiant emittance of this light (W / m^2).
   //
-  Color intensity;
+  TexVal<Color> intensity;
+
+  // Intensity distribution over the surface.
+  //
+  Hist2dDist intensity_dist;
 };
 
 
diff --git a/surface/cylinder.cc b/surface/cylinder.cc
index b10f2bd..785c259 100644
--- a/surface/cylinder.cc
+++ b/surface/cylinder.cc
@@ -240,7 +240,8 @@ Cylinder::Sampler::sample (const UV &param) const
   float theta = param.u * 2 * PIf;
   Vec radius (cos (theta), sin (theta), 0);
   Vec norm = cylinder.normal_to_world (radius).unit();
-  return sample_with_approx_area_pdf (PosSampler (cylinder), param, norm);
+  return sample_with_approx_area_pdf (PosSampler (cylinder),
+				      param, norm, param);
 }
 
 // Return a sample of this surface from VIEWPOINT, based on the
diff --git a/surface/ellipse.cc b/surface/ellipse.cc
index 7f86d88..55d9ae6 100644
--- a/surface/ellipse.cc
+++ b/surface/ellipse.cc
@@ -160,11 +160,11 @@ Ellipse::Sampler::sample (const UV &param) const
 {
   dist_t dx, dy;
   disk_sample (dist_t (0.5), param, dx, dy);
-  dx += dist_t (0.5);
-  dy += dist_t (0.5);
-  Pos pos = ellipse.corner + ellipse.edge1 * dx + ellipse.edge2 * dy;
 
-  return AreaSample (pos, ellipse.normal, pdf);
+  UV uv (dx + dist_t (0.5), dy + dist_t (0.5));
+  Pos pos = ellipse.corner + ellipse.edge1 * uv.u + ellipse.edge2 * uv.v;
+
+  return AreaSample (pos, ellipse.normal, pdf, uv);
 }
 
 // If a ray from VIEWPOINT in direction DIR intersects this
@@ -182,7 +182,8 @@ Ellipse::Sampler::eval_from_viewpoint (const Pos &viewpoint, const Vec &dir)
   if (ellipse.intersects (viewpoint, dir, t, u, v))
     {
       Pos pos = viewpoint + t * dir;
-      return AngularSample (AreaSample (pos, ellipse.normal, pdf), viewpoint);
+      return AngularSample (AreaSample (pos, ellipse.normal, pdf, UV (u, v)),
+			    viewpoint);
     }
   return AngularSample ();
 }
diff --git a/surface/sphere.cc b/surface/sphere.cc
index adfbf40..10f58f2 100644
--- a/surface/sphere.cc
+++ b/surface/sphere.cc
@@ -185,7 +185,7 @@ Sphere::Sampler::sample (const UV &param) const
   Pos pos = center + norm * radius;
   float pdf = 1 / (float (radius*radius) * 4 * PIf); // 1 / area
 
-  return AreaSample (pos, norm, pdf);
+  return AreaSample (pos, norm, pdf, sphere.tex_coords (norm));
 }
 
 // If a ray from VIEWPOINT in direction DIR intersects this
@@ -206,7 +206,9 @@ Sphere::Sampler::eval_from_viewpoint (const Pos &viewpoint, const Vec &dir)
       Vec norm = vec.unit ();
       float area_pdf =		// 1 / area
 	1 / (float (sphere.radius * sphere.radius) * 4 * PIf);
-      return AngularSample (AreaSample (pos, norm, area_pdf), viewpoint);
+      return AngularSample (AreaSample (pos, norm, area_pdf,
+					sphere.tex_coords (norm)),
+			    viewpoint);
     }
   return AngularSample ();
 }
diff --git a/surface/sphere2.cc b/surface/sphere2.cc
index 4b8f40a..624f3b5 100644
--- a/surface/sphere2.cc
+++ b/surface/sphere2.cc
@@ -172,8 +172,10 @@ struct PosSampler
 Surface::Sampler::AreaSample
 Sphere2::Sampler::sample (const UV &param) const
 {
-  Vec norm = sphere.normal_to_world (sphere_sample (param)).unit ();
-  return sample_with_approx_area_pdf (PosSampler (sphere), param, norm);
+  Vec onorm = sphere_sample (param);
+  Vec norm = sphere.normal_to_world (onorm).unit ();
+  return sample_with_approx_area_pdf (PosSampler (sphere), param, norm,
+				      sphere.tex_coords (onorm));
 }
 
 // Return a sample of this surface from VIEWPOINT, based on the
@@ -200,7 +202,8 @@ Sphere2::Sampler::sample_from_viewpoint (const Pos &viewpoint, const UV &param)
   //
   Vec norm = sphere.normal_to_world (onorm).unit ();
   Surface::Sampler::AreaSample area_sample =
-    sample_with_approx_area_pdf (PosSampler (sphere), samp_param, norm);
+    sample_with_approx_area_pdf (PosSampler (sphere), samp_param, norm,
+				 sphere.tex_coords (onorm));
 
   // Because we mirror samples to always point towards VIEWPOINT, double
   // the PDF, as the same number of samples is concentrated into half
diff --git a/surface/surface-sampler.cc b/surface/surface-sampler.cc
index ec535f5..7bc4653 100644
--- a/surface/surface-sampler.cc
+++ b/surface/surface-sampler.cc
@@ -25,7 +25,8 @@ using namespace snogray;
 //
 Surface::Sampler::AngularSample::AngularSample (const AreaSample &area_sample,
 						const Pos &viewpoint)
-  : normal (area_sample.normal), pdf (0)
+  : normal (area_sample.normal), pdf (0),
+    tex_coords (area_sample.pos, area_sample.tex_coords_uv)
 {
   Vec view_vec = area_sample.pos - viewpoint;
 
diff --git a/surface/surface-sampler.h b/surface/surface-sampler.h
index ec33b3c..da656d1 100644
--- a/surface/surface-sampler.h
+++ b/surface/surface-sampler.h
@@ -30,8 +30,9 @@ public:
   //
   struct AreaSample
   {
-    AreaSample (const Pos &_pos, const Vec &_norm, float _pdf)
-      : pos (_pos), normal (_norm), pdf (_pdf)
+    AreaSample (const Pos &_pos, const Vec &_norm, float _pdf,
+		const UV &tex_coords_uv)
+      : pos (_pos), normal (_norm), pdf (_pdf), tex_coords_uv (tex_coords_uv)
     { }
     AreaSample () : pdf (0) { }
 
@@ -48,6 +49,10 @@ public:
     // surface.
     //
     float pdf;
+
+    // 2D-texture coordinates of the point on the surface.
+    //
+    UV tex_coords_uv;
   };
 
   // A sample of the surface area from a particular viewpoint.
@@ -58,8 +63,9 @@ public:
   struct AngularSample
   {
     AngularSample (const Vec &_dir, const Vec &_norm,
-		   float _pdf, dist_t _dist)
-      : dir (_dir), normal (_norm), pdf (_pdf), dist (_dist)
+		   float _pdf, dist_t _dist, const TexCoords &_tex_coords)
+      : dir (_dir), normal (_norm), pdf (_pdf), dist (_dist),
+	tex_coords (_tex_coords)
     { }
     AngularSample () : pdf (0) { }
 
@@ -85,6 +91,10 @@ public:
     // The distance from the viewpoint to the sample.
     //
     dist_t dist;
+
+    // Texture coordinates of the point on the surface.
+    //
+    TexCoords tex_coords;
   };
 
   // Return a sample of this surface.
@@ -139,7 +149,8 @@ protected:
   template<typename PosSampleFun>
   static Surface::Sampler::AreaSample
   sample_with_approx_area_pdf (const PosSampleFun &pos_sample_fun,
-			       const UV &param, const Vec &norm)
+			       const UV &param, const Vec &norm,
+			       const UV &tex_coords_uv)
   {
     float delta = 0.0001f;  // this value seems to work well
     Pos pos = pos_sample_fun (param);
@@ -148,7 +159,7 @@ protected:
     float sample_area = cross (pos_du - pos, pos_dv - pos).length ();
     float param_area = delta * delta;
     float pdf = sample_area == 0 ? 0 : param_area / sample_area;
-    return AreaSample (pos, norm, pdf);
+    return AreaSample (pos, norm, pdf, tex_coords_uv);
   }
 };
 
diff --git a/surface/tripar.cc b/surface/tripar.cc
index 7c07ea6..d7bc994 100644
--- a/surface/tripar.cc
+++ b/surface/tripar.cc
@@ -182,7 +182,7 @@ Tripar::Sampler::sample (const UV &param) const
   // pdf
   float pdf = 1 / area;
 
-  return AreaSample (pos, norm, pdf);
+  return AreaSample (pos, norm, pdf, UV (u, v));
 }
 
 // If a ray from VIEWPOINT in direction DIR intersects this