--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Framework/ColorSpaces.cpp	Fri Oct 18 13:08:55 2024 +0200
@@ -0,0 +1,274 @@
+/**
+ * Orthanc - A Lightweight, RESTful DICOM Store
+ * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics
+ * Department, University Hospital of Liege, Belgium
+ * Copyright (C) 2017-2023 Osimis S.A., Belgium
+ * Copyright (C) 2024-2024 Orthanc Team SRL, Belgium
+ * Copyright (C) 2021-2024 Sebastien Jodogne, ICTEAM UCLouvain, Belgium
+ *
+ * This program is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU Affero General Public License
+ * as published by the Free Software Foundation, either version 3 of
+ * the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Affero General Public License for more details.
+ * 
+ * You should have received a copy of the GNU Affero General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ **/
+
+
+#include "ColorSpaces.h"
+
+#include <SerializationToolbox.h>
+#include <Toolbox.h>
+
+#include <boost/math/special_functions/round.hpp>
+
+
+namespace OrthancWSI
+{
+  RGBColor::RGBColor(const sRGBColor& srgb)
+  {
+    if (srgb.GetR() < 0)
+    {
+      r_ = 0;
+    }
+    else if (srgb.GetR() >= 1)
+    {
+      r_ = 255;
+    }
+    else
+    {
+      r_ = boost::math::iround(srgb.GetR() * 255.0f);
+    }
+
+    if (srgb.GetG() < 0)
+    {
+      g_ = 0;
+    }
+    else if (srgb.GetG() >= 1)
+    {
+      g_ = 255;
+    }
+    else
+    {
+      g_ = boost::math::iround(srgb.GetG() * 255.0f);
+    }
+
+    if (srgb.GetB() < 0)
+    {
+      b_ = 0;
+    }
+    else if (srgb.GetB() >= 1)
+    {
+      b_ = 255;
+    }
+    else
+    {
+      b_ = boost::math::iround(srgb.GetB() * 255.0f);
+    }
+  }
+
+
+  sRGBColor::sRGBColor(const RGBColor& rgb)
+  {
+    r_ = static_cast<float>(rgb.GetR()) / 255.0f;
+    g_ = static_cast<float>(rgb.GetG()) / 255.0f;
+    b_ = static_cast<float>(rgb.GetB()) / 255.0f;
+  }
+
+
+  static float ApplyGammaXYZ(float value)
+  {
+    // https://www.image-engineering.de/library/technotes/958-how-to-convert-between-srgb-and-ciexyz
+    if (value <= 0.0031308f)
+    {
+      return value * 12.92f;
+    }
+    else
+    {
+      return 1.055f * powf(value, 1.0f / 2.4f) - 0.055f;
+    }
+  }
+
+
+  sRGBColor::sRGBColor(const XYZColor& xyz)
+  {
+    // https://en.wikipedia.org/wiki/SRGB#From_CIE_XYZ_to_sRGB
+    // https://www.image-engineering.de/library/technotes/958-how-to-convert-between-srgb-and-ciexyz
+    const float sr =  3.2404542f * xyz.GetX() - 1.5371385f * xyz.GetY() - 0.4985314f * xyz.GetZ();
+    const float sg = -0.9692660f * xyz.GetX() + 1.8760108f * xyz.GetY() + 0.0415560f * xyz.GetZ();
+    const float sb =  0.0556434f * xyz.GetX() - 0.2040259f * xyz.GetY() + 1.0572252f * xyz.GetZ();
+
+    r_ = ApplyGammaXYZ(sr);
+    g_ = ApplyGammaXYZ(sg);
+    b_ = ApplyGammaXYZ(sb);
+  }
+
+
+  static float LinearizeXYZ(float value)
+  {
+    // https://www.image-engineering.de/library/technotes/958-how-to-convert-between-srgb-and-ciexyz
+    if (value <= 0.04045f)
+    {
+      return value / 12.92f;
+    }
+    else
+    {
+      return powf((value + 0.055f) / 1.055f, 2.4f);
+    }
+  }
+
+
+  XYZColor::XYZColor(const sRGBColor& srgb)
+  {
+    // https://en.wikipedia.org/wiki/SRGB#From_sRGB_to_CIE_XYZ
+    // https://www.image-engineering.de/library/technotes/958-how-to-convert-between-srgb-and-ciexyz
+    const float linearizedR = LinearizeXYZ(srgb.GetR());
+    const float linearizedG = LinearizeXYZ(srgb.GetG());
+    const float linearizedB = LinearizeXYZ(srgb.GetB());
+
+    x_ = 0.4124564f * linearizedR + 0.3575761f * linearizedG + 0.1804375f * linearizedB;
+    y_ = 0.2126729f * linearizedR + 0.7151522f * linearizedG + 0.0721750f * linearizedB;
+    z_ = 0.0193339f * linearizedR + 0.1191920f * linearizedG + 0.9503041f * linearizedB;
+  }
+
+
+  static const float LAB_DELTA = 6.0f / 29.0f;
+
+  static float LABf_inv(float t)
+  {
+    if (t > LAB_DELTA)
+    {
+      return powf(t, 3.0f);
+    }
+    else
+    {
+      return 3.0f * LAB_DELTA * LAB_DELTA * (t - 4.0f / 29.0f);
+    }
+  }
+
+
+  // Those correspond to Standard Illuminant D65
+  // https://en.wikipedia.org/wiki/CIELAB_color_space#From_CIEXYZ_to_CIELAB
+  static const float X_N = 95.0489f;
+  static const float Y_N = 100.0f;
+  static const float Z_N = 108.8840f;
+
+  XYZColor::XYZColor(const LABColor& lab)
+  {
+    // https://en.wikipedia.org/wiki/CIELAB_color_space#From_CIELAB_to_CIEXYZ
+    const float shared = (lab.GetL() + 16.0f) / 116.0f;
+
+    x_ = X_N * LABf_inv(shared + lab.GetA() / 500.0f) / 100.0f;
+    y_ = Y_N * LABf_inv(shared) / 100.0f;
+    z_ = Z_N * LABf_inv(shared - lab.GetB() / 200.0f) / 100.0f;
+  }
+
+
+  static float LABf(float t)
+  {
+    if (t > powf(LAB_DELTA, 3.0f))
+    {
+      return powf(t, 1.0f / 3.0f);
+    }
+    else
+    {
+      return t / 3.0f * powf(LAB_DELTA, -2.0f) + 4.0f / 29.0f;
+    }
+  }
+
+
+  LABColor::LABColor(const XYZColor& xyz)
+  {
+    // https://en.wikipedia.org/wiki/CIELAB_color_space#From_CIEXYZ_to_CIELAB
+
+    const float fx = LABf(xyz.GetX() * 100.0f / X_N);
+    const float fy = LABf(xyz.GetY() * 100.0f / Y_N);
+    const float fz = LABf(xyz.GetZ() * 100.0f / Z_N);
+
+    l_ = 116.0f * fy - 16.0f;
+    a_ = 500.0f * (fx - fy);
+    b_ = 200.0f * (fy - fz);
+  }
+
+
+  static uint16_t EncodeUint16(float value,
+                               float minValue,
+                               float maxValue)
+  {
+    if (value <= minValue)
+    {
+      return 0;
+    }
+    else if (value >= maxValue)
+    {
+      return 0xffff;
+    }
+    else
+    {
+      float lambda = (value - minValue) / (maxValue - minValue);
+      assert(lambda >= 0 && lambda <= 1);
+      return static_cast<uint16_t>(boost::math::iround(lambda * static_cast<float>(0xffff)));
+    }
+  }
+
+
+  void LABColor::EncodeDicomRecommendedAbsentPixelCIELab(uint16_t target[3]) const
+  {
+    /**
+     * "An L value linearly scaled to 16 bits, such that 0x0000
+     * corresponds to an L of 0.0, and 0xFFFF corresponds to an L of
+     * 100.0."
+     **/
+    target[0] = EncodeUint16(GetL(), 0.0f, 100.0f);
+
+    /**
+     * "An a* then a b* value, each linearly scaled to 16 bits and
+     * offset to an unsigned range, such that 0x0000 corresponds to an
+     * a* or b* of -128.0, 0x8080 corresponds to an a* or b* of 0.0
+     * and 0xFFFF corresponds to an a* or b* of 127.0"
+     **/
+    target[1] = EncodeUint16(GetA(), -128.0f, 127.0f);
+    target[2] = EncodeUint16(GetB(), -128.0f, 127.0f);
+  }
+
+
+  LABColor LABColor::DecodeDicomRecommendedAbsentPixelCIELab(uint16_t l,
+                                                             uint16_t a,
+                                                             uint16_t b)
+  {
+    return LABColor(static_cast<float>(l) / static_cast<float>(0xffff) * 100.0f,
+                    -128.0f + static_cast<float>(a) / static_cast<float>(0xffff) * 255.0f,
+                    -128.0f + static_cast<float>(b) / static_cast<float>(0xffff) * 255.0f);
+  }
+
+
+  bool LABColor::DecodeDicomRecommendedAbsentPixelCIELab(LABColor& target,
+                                                         const std::string& tag)
+  {
+    std::vector<std::string> channels;
+    Orthanc::Toolbox::TokenizeString(channels, tag, '\\');
+
+    unsigned int l, a, b;
+    if (channels.size() == 3 &&
+        Orthanc::SerializationToolbox::ParseUnsignedInteger32(l, channels[0]) &&
+        Orthanc::SerializationToolbox::ParseUnsignedInteger32(a, channels[1]) &&
+        Orthanc::SerializationToolbox::ParseUnsignedInteger32(b, channels[2]) &&
+        l <= 0xffffu &&
+        a <= 0xffffu &&
+        b <= 0xffffu)
+    {
+      target = LABColor::DecodeDicomRecommendedAbsentPixelCIELab(l, a, b);
+      return true;
+    }
+    else
+    {
+      return false;
+    }
+  }
+}