Compositor: Improve Fog Glow glare realism

This patch improves the realism of the Fog Glow mode of the Glare node
based on the Photopic model described in:

  Physically-Based Glare Effects for Digital Images" by G. Spencer, P.
  Shirley, K. Zimmerman, and D. P. Greenberg.

This is a breaking change that can't be versioned, but it is worth it
for the superior realism of the new model.

Pull Request: https://projects.blender.org/blender/blender/pulls/140646

source/blender/blenlib/BLI_math_base.hh

@@ -201,6 +201,11 @@ template<typename T> inline T square(const T &a)
   return a * a;
 }
 
+template<typename T> inline T cube(const T &a)
+{
+  return a * a * a;
+}
+
 template<typename T> inline T exp(const T &x)
 {
   return std::exp(x);

source/blender/compositor/cached_resources/COM_fog_glow_kernel.hh

@@ -9,6 +9,7 @@
 #include <memory>
 
 #include "BLI_map.hh"
+#include "BLI_math_angle_types.hh"
 #include "BLI_math_vector_types.hh"
 
 #include "COM_cached_resource.hh"
@@ -22,8 +23,9 @@ class FogGlowKernelKey {
  public:
   int kernel_size;
   int2 spatial_size;
+  math::AngleRadian field_of_view;
 
-  FogGlowKernelKey(int kernel_size, int2 spatial_size);
+  FogGlowKernelKey(int kernel_size, int2 spatial_size, math::AngleRadian field_of_view);
 
   uint64_t hash() const;
 };
@@ -46,7 +48,7 @@ class FogGlowKernel : public CachedResource {
   std::complex<float> *frequencies_ = nullptr;
 
  public:
-  FogGlowKernel(int kernel_size, int2 spatial_size);
+  FogGlowKernel(int kernel_size, int2 spatial_size, math::AngleRadian field_of_view);
 
   ~FogGlowKernel();
 
@@ -69,7 +71,7 @@ class FogGlowKernelContainer : CachedResourceContainer {
    * container, if one exists, return it, otherwise, return a newly created one and add it to the
    * container. In both cases, tag the cached resource as needed to keep it cached for the next
    * evaluation. */
-  FogGlowKernel &get(int kernel_size, int2 spatial_size);
+  FogGlowKernel &get(int kernel_size, int2 spatial_size, math::AngleRadian field_of_view);
 };
 
 }  // namespace blender::compositor

source/blender/compositor/cached_resources/intern/fog_glow_kernel.cc

@@ -16,7 +16,7 @@
 #include "BLI_index_range.hh"
 #include "BLI_math_base.h"
 #include "BLI_math_base.hh"
-#include "BLI_math_numbers.hh"
+#include "BLI_math_vector.hh"
 #include "BLI_math_vector_types.hh"
 #include "BLI_task.hh"
 
@@ -28,46 +28,50 @@ namespace blender::compositor {
  * Fog Glow Kernel Key.
  */
 
-FogGlowKernelKey::FogGlowKernelKey(int kernel_size, int2 spatial_size)
-    : kernel_size(kernel_size), spatial_size(spatial_size)
+FogGlowKernelKey::FogGlowKernelKey(int kernel_size,
+                                   int2 spatial_size,
+                                   math::AngleRadian field_of_view)
+    : kernel_size(kernel_size), spatial_size(spatial_size), field_of_view(field_of_view)
 {
 }
 
 uint64_t FogGlowKernelKey::hash() const
 {
-  return get_default_hash(kernel_size, spatial_size);
+  return get_default_hash(kernel_size, spatial_size, field_of_view.degree());
 }
 
 bool operator==(const FogGlowKernelKey &a, const FogGlowKernelKey &b)
 {
-  return a.kernel_size == b.kernel_size && a.spatial_size == b.spatial_size;
+  return a.kernel_size == b.kernel_size && a.spatial_size == b.spatial_size &&
+         a.field_of_view == b.field_of_view;
 }
 
 /* --------------------------------------------------------------------
  * Fog Glow Kernel.
  */
 
-/* Given the x and y location in the range from 0 to kernel_size - 1, where kernel_size is odd,
- * compute the fog glow kernel value. The equations are arbitrary and were chosen using visual
- * judgment. The kernel is not normalized and need normalization. */
-[[maybe_unused]] static float compute_fog_glow_kernel_value(int x, int y, int kernel_size)
+/* Given the texel coordinates and the constant field-of-view-per-pixel value, under the assumption
+ * of a relatively small field of view as discussed in Section 3.2, this function computes the
+ * fog glow kernel value. The kernel value is derived from Equation (5) of the following paper:
+ *
+ *   Spencer, Greg, et al. "Physically-Based Glare Effects for Digital Images."
+ *   Proceedings of the 22nd Annual Conference on Computer Graphics and Interactive Techniques,
+ *   1995.
+ */
+
+[[maybe_unused]] static float compute_fog_glow_kernel_value(
+    int2 texel, math::AngleRadian field_of_view_per_pixel)
 {
-  const int half_kernel_size = kernel_size / 2;
-  const float scale = 0.25f * math::sqrt(math::square(kernel_size));
-  const float v = ((y - half_kernel_size) / float(half_kernel_size));
-  const float u = ((x - half_kernel_size) / float(half_kernel_size));
-  const float r = (math::square(u) + math::square(v)) * scale;
-  const float d = -math::sqrt(math::sqrt(math::sqrt(r))) * 9.0f;
-  const float kernel_value = math::exp(d);
+  const float theta_degree = math::length(float2(texel)) * field_of_view_per_pixel.degree();
+  const float f0 = 2.61f * 1e6f * math::exp(-math::square(theta_degree / 0.02f));
+  const float f1 = 20.91f / math::cube(theta_degree + 0.02f);
+  const float f2 = 72.37f / math::square(theta_degree + 0.02f);
+  const float kernel_value = 0.384f * f0 + 0.478f * f1 + 0.138f * f2;
 
-  const float window = (0.5f + 0.5f * math::cos(u * math::numbers::pi)) *
-                       (0.5f + 0.5f * math::cos(v * math::numbers::pi));
-  const float windowed_kernel_value = window * kernel_value;
-
-  return windowed_kernel_value;
+  return kernel_value;
 }
 
-FogGlowKernel::FogGlowKernel(int kernel_size, int2 spatial_size)
+FogGlowKernel::FogGlowKernel(int kernel_size, int2 spatial_size, math::AngleRadian field_of_view)
 {
 #if defined(WITH_FFTW3)
 
@@ -91,27 +95,26 @@ FogGlowKernel::FogGlowKernel(int kernel_size, int2 spatial_size)
   /* Use a double to sum the kernel since floats are not stable with threaded summation. */
   threading::EnumerableThreadSpecific<double> sum_by_thread([]() { return 0.0; });
 
-  /* Compute the kernel while zero padding to match the padded image size. */
+  /* Compute the entire kernel's spatial space using compute_fog_glow_kernel_value. */
   threading::parallel_for(IndexRange(spatial_size.y), 1, [&](const IndexRange sub_y_range) {
     double &sum = sum_by_thread.local();
     for (const int64_t y : sub_y_range) {
       for (const int64_t x : IndexRange(spatial_size.x)) {
+        const int2 texel = int2(x, y);
+        const int2 center_texel = spatial_size / 2;
+        const int2 kernel_texel = texel - center_texel;
+        const math::AngleRadian field_of_view_per_pixel = field_of_view / kernel_size;
+
+        const float kernel_value = compute_fog_glow_kernel_value(kernel_texel,
+                                                                 field_of_view_per_pixel);
+        sum += kernel_value;
+
         /* We offset the computed kernel with wrap around such that it is centered at the zero
          * point, which is the expected format for doing circular convolutions in the frequency
          * domain. */
-        const int half_kernel_size = kernel_size / 2;
-        int64_t output_x = mod_i(x - half_kernel_size, spatial_size.x);
-        int64_t output_y = mod_i(y - half_kernel_size, spatial_size.y);
-
-        const bool is_inside_kernel = x < kernel_size && y < kernel_size;
-        if (is_inside_kernel) {
-          const float kernel_value = compute_fog_glow_kernel_value(x, y, kernel_size);
-          kernel_spatial_domain[output_x + output_y * spatial_size.x] = kernel_value;
-          sum += kernel_value;
-        }
-        else {
-          kernel_spatial_domain[output_x + output_y * spatial_size.x] = 0.0f;
-        }
+        int64_t output_x = mod_i(kernel_texel.x, spatial_size.x);
+        int64_t output_y = mod_i(kernel_texel.y, spatial_size.y);
+        kernel_spatial_domain[output_x + output_y * spatial_size.x] = kernel_value;
       }
     }
   });
@@ -127,7 +130,7 @@ FogGlowKernel::FogGlowKernel(int kernel_size, int2 spatial_size)
    * Fourier transform is linear. */
   normalization_factor_ = float(std::accumulate(sum_by_thread.begin(), sum_by_thread.end(), 0.0));
 #else
-  UNUSED_VARS(kernel_size, spatial_size);
+  UNUSED_VARS(kernel_size, spatial_size, field_of_view);
 #endif
 }
 
@@ -164,12 +167,15 @@ void FogGlowKernelContainer::reset()
   }
 }
 
-FogGlowKernel &FogGlowKernelContainer::get(int kernel_size, int2 spatial_size)
+FogGlowKernel &FogGlowKernelContainer::get(int kernel_size,
+                                           int2 spatial_size,
+                                           math::AngleRadian field_of_view)
 {
-  const FogGlowKernelKey key(kernel_size, spatial_size);
+  const FogGlowKernelKey key(kernel_size, spatial_size, field_of_view);
 
-  auto &kernel = *map_.lookup_or_add_cb(
-      key, [&]() { return std::make_unique<FogGlowKernel>(kernel_size, spatial_size); });
+  auto &kernel = *map_.lookup_or_add_cb(key, [&]() {
+    return std::make_unique<FogGlowKernel>(kernel_size, spatial_size, field_of_view);
+  });
 
   kernel.needed = true;
   return kernel;

source/blender/nodes/composite/nodes/node_composite_glare.cc

@@ -21,6 +21,7 @@
 #include "BLI_assert.h"
 #include "BLI_fftw.hh"
 #include "BLI_index_range.hh"
+#include "BLI_math_angle_types.hh"
 #include "BLI_math_base.hh"
 #include "BLI_math_color.h"
 #include "BLI_math_vector.hh"
@@ -2031,14 +2032,14 @@ class GlareOperation : public NodeOperation {
   {
 #if defined(WITH_FFTW3)
 
-    const int kernel_size = compute_fog_glow_kernel_size(highlights);
+    const int kernel_size = int(math::reduce_max(highlights.domain().size));
 
-    /* Since we will be doing a circular convolution, we need to zero pad our input image by half
+    /* Since we will be doing a circular convolution, we need to zero pad our input image by
      * the kernel size to avoid the kernel affecting the pixels at the other side of image.
      * Therefore, zero boundary is assumed. */
-    const int needed_padding_amount = kernel_size / 2;
+    const int needed_padding_amount = kernel_size;
     const int2 image_size = highlights.domain().size;
-    const int2 needed_spatial_size = image_size + needed_padding_amount;
+    const int2 needed_spatial_size = image_size + needed_padding_amount - 1;
     const int2 spatial_size = fftw::optimal_size_for_real_transform(needed_spatial_size);
 
     /* The FFTW real to complex transforms utilizes the hermitian symmetry of real transforms and
@@ -2108,7 +2109,7 @@ class GlareOperation : public NodeOperation {
     });
 
     const FogGlowKernel &fog_glow_kernel = context().cache_manager().fog_glow_kernels.get(
-        kernel_size, spatial_size);
+        kernel_size, spatial_size, this->compute_fog_glow_field_of_view());
 
     /* Multiply the kernel and the image in the frequency domain to perform the convolution. The
      * FFT is not normalized, meaning the result of the FFT followed by an inverse FFT will result
@@ -2197,23 +2198,21 @@ class GlareOperation : public NodeOperation {
     return fog_glow_result;
   }
 
-  /* Computes the size of the fog glow kernel that will be convolved with the image, which is
-   * essentially the extent of the glare in pixels. */
-  int compute_fog_glow_kernel_size(const Result &highlights)
+  /* Computes the field of view of the glare based on the give size as per:
+   *
+   *   Spencer, Greg, et al. "Physically-Based Glare Effects for Digital Images."
+   *   Proceedings of the 22nd Annual Conference on Computer Graphics and Interactive Techniques,
+   *   1995.
+   *
+   * We choose a minimum field of view of 10 degrees using visual judgement on typical setups,
+   * otherwise, a too small field of view would make the evaluation domain of the glare lie almost
+   * entirely in the central Gaussian of the function, losing the exponential characteristic of the
+   * function. Additionally, we take the power of the size with 1/3 to adjust the rate of change of
+   * the size to make the apparent size of the glare more linear with respect to the size input. */
+  math::AngleRadian compute_fog_glow_field_of_view()
   {
-    /* The input size is relative to the larger dimension of the image. */
-    const int size = int(math::reduce_max(highlights.domain().size) * this->get_size());
-
-    /* Make sure size is at least 3 pixels for implicitly since code deals with half kernel sizes
-     * which will be zero if less than 3, causing zero division. */
-    const int safe_size = math::max(3, size);
-
-    /* Make sure the kernel size is odd since an even one will typically introduce a tiny offset as
-     * it has no exact center value. */
-    const bool is_even = safe_size % 2 == 0;
-    const int odd_size = safe_size + (is_even ? 1 : 0);
-
-    return odd_size;
+    return math::AngleRadian::from_degree(
+        math::interpolate(180.0f, 10.0f, math::pow(this->get_size(), 1.0f / 3.0f)));
   }
 
   /* ----------