lib/jxl/enc_photon_noise.cc

   1 // Copyright (c) the JPEG XL Project Authors. All rights reserved.
   2 //
   3 // Use of this source code is governed by a BSD-style
   4 // license that can be found in the LICENSE file.
   5
   6 #include "lib/jxl/enc_photon_noise.h"
   7
   8 #include <algorithm>
   9
  10 #include "lib/jxl/cms/opsin_params.h"
  11
  12 namespace jxl {
  13
  14 namespace {
  15
  16 // Assumes a daylight-like spectrum.
  17 // https://www.strollswithmydog.com/effective-quantum-efficiency-of-sensor/#:~:text=11%2C260%20photons/um%5E2/lx-s
  18 constexpr float kPhotonsPerLxSPerUm2 = 11260;
  19
  20 // Order of magnitude for cameras in the 2010-2020 decade, taking the CFA into
  21 // account.
  22 constexpr float kEffectiveQuantumEfficiency = 0.20;
  23
  24 // TODO(sboukortt): reevaluate whether these are good defaults, notably whether
  25 // it would be worth making read noise higher at lower ISO settings.
  26 constexpr float kPhotoResponseNonUniformity = 0.005;
  27 constexpr float kInputReferredReadNoise = 3;
  28
  29 // Assumes a 35mm sensor.
  30 constexpr float kSensorAreaUm2 = 36000.f * 24000;
  31
  32 template <typename T>
  33 inline constexpr T Square(const T x) {
  34   return x * x;
  35 }
  36 template <typename T>
  37 inline constexpr T Cube(const T x) {
  38   return x * x * x;
  39 }
  40
  41 }  // namespace
  42
  43 NoiseParams SimulatePhotonNoise(const size_t xsize, const size_t ysize,
  44                                 const float iso) {
  45   const float kOpsinAbsorbanceBiasCbrt =
  46       std::cbrt(jxl::cms::kOpsinAbsorbanceBias[1]);
  47
  48   // Focal plane exposure for 18% of kDefaultIntensityTarget, in lx·s.
  49   // (ISO = 10 lx·s ÷ H)
  50   const float h_18 = 10 / iso;
  51
  52   const float pixel_area_um2 = kSensorAreaUm2 / (xsize * ysize);
  53
  54   const float electrons_per_pixel_18 = kEffectiveQuantumEfficiency *
  55                                        kPhotonsPerLxSPerUm2 * h_18 *
  56                                        pixel_area_um2;
  57
  58   NoiseParams params;
  59
  60   for (size_t i = 0; i < NoiseParams::kNumNoisePoints; ++i) {
  61     const float scaled_index = i / (NoiseParams::kNumNoisePoints - 2.f);
  62     // scaled_index is used for XYB = (0, 2·scaled_index, 2·scaled_index)
  63     const float y = 2 * scaled_index;
  64     // 1 = default intensity target
  65     const float linear = std::max(0.f, Cube(y - kOpsinAbsorbanceBiasCbrt) +
  66                                            jxl::cms::kOpsinAbsorbanceBias[1]);
  67     const float electrons_per_pixel = electrons_per_pixel_18 * (linear / 0.18f);
  68     // Quadrature sum of read noise, photon shot noise (sqrt(S) so simply not
  69     // squared here) and photo response non-uniformity.
  70     // https://doi.org/10.1117/3.725073
  71     // Units are electrons rms.
  72     const float noise =
  73         std::sqrt(Square(kInputReferredReadNoise) + electrons_per_pixel +
  74                   Square(kPhotoResponseNonUniformity * electrons_per_pixel));
  75     const float linear_noise = noise * (0.18f / electrons_per_pixel_18);
  76     const float opsin_derivative =
  77         (1.f / 3) /
  78         Square(std::cbrt(linear - jxl::cms::kOpsinAbsorbanceBias[1]));
  79     const float opsin_noise = linear_noise * opsin_derivative;
  80
  81     // TODO(sboukortt): verify more thoroughly whether the denominator is
  82     // correct.
  83     params.lut[i] =
  84         Clamp1(opsin_noise /
  85                    (0.22f             // norm_const
  86                     * std::sqrt(2.f)  // red_noise + green_noise
  87                     * 1.13f  // standard deviation of a plane of generated noise
  88                     ),
  89                0.f, 1.f);
  90   }
  91
  92   return params;
  93 }
  94
  95 }  // namespace jxl