/**
 * 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 Lesser 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
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this program. If not, see
 * <http://www.gnu.org/licenses/>.
 **/


#include "../PrecompiledHeaders.h"
#include "ImageProcessing.h"

#include "Image.h"
#include "ImageTraits.h"
#include "PixelTraits.h"
#include "../OrthancException.h"

#ifdef __EMSCRIPTEN__
/* 
   Avoid this error:
   -----------------
   .../boost/math/special_functions/round.hpp:118:12: warning: implicit conversion from 'std::__2::numeric_limits<long long>::type' (aka 'long long') to 'float' changes value from 9223372036854775807 to 9223372036854775808 [-Wimplicit-int-float-conversion]
   .../mnt/c/osi/dev/orthanc/Core/Images/ImageProcessing.cpp:333:28: note: in instantiation of function template specialization 'boost::math::llround<float>' requested here
   .../mnt/c/osi/dev/orthanc/Core/Images/ImageProcessing.cpp:1006:9: note: in instantiation of function template specialization 'Orthanc::MultiplyConstantInternal<unsigned char, true>' requested here
*/
#pragma GCC diagnostic ignored "-Wimplicit-int-float-conversion"
#endif 

#include <boost/math/special_functions/round.hpp>

#include <algorithm>
#include <cassert>
#include <limits>
#include <list>
#include <map>
#include <stdint.h>
#include <string.h>

namespace Orthanc
{
  ImageProcessing::ImagePoint::ImagePoint(int32_t x,
                                          int32_t y) :
    x_(x),
    y_(y)
  {
  }

  int32_t ImageProcessing::ImagePoint::GetX() const
  {
    return x_;
  }

  int32_t ImageProcessing::ImagePoint::GetY() const
  {
    return y_;
  }

  void ImageProcessing::ImagePoint::Set(int32_t x, int32_t y)
  {
    x_ = x;
    y_ = y;
  }

  void ImageProcessing::ImagePoint::ClipTo(int32_t minX, int32_t maxX, int32_t minY, int32_t maxY)
  {
    x_ = std::max(minX, std::min(maxX, x_));
    y_ = std::max(minY, std::min(maxY, y_));
  }

  double ImageProcessing::ImagePoint::GetDistanceTo(const ImagePoint& other) const
  {
    double dx = (double)(other.GetX() - GetX());
    double dy = (double)(other.GetY() - GetY());
    return sqrt(dx * dx + dy * dy);
  }

  double ImageProcessing::ImagePoint::GetDistanceToLine(double a, double b, double c) const // where ax + by + c = 0 is the equation of the line
  {
    return std::abs(a * static_cast<double>(GetX()) + b * static_cast<double>(GetY()) + c) / pow(a * a + b * b, 0.5);
  }

  template <typename TargetType, typename SourceType>
  static void ConvertInternal(ImageAccessor& target,
                              const ImageAccessor& source)
  {
    // WARNING - "::min()" should be replaced by "::lowest()" if
    // dealing with float or double (which is not the case so far)
    assert(sizeof(TargetType) <= 2);  // Safeguard to remember about "float/double"
    const TargetType minValue = std::numeric_limits<TargetType>::min();
    const TargetType maxValue = std::numeric_limits<TargetType>::max();

    const unsigned int width = source.GetWidth();
    const unsigned int height = source.GetHeight();
    
    for (unsigned int y = 0; y < height; y++)
    {
      TargetType* t = reinterpret_cast<TargetType*>(target.GetRow(y));
      const SourceType* s = reinterpret_cast<const SourceType*>(source.GetConstRow(y));

      for (unsigned int x = 0; x < width; x++, t++, s++)
      {
        if (static_cast<int32_t>(*s) < static_cast<int32_t>(minValue))
        {
          *t = minValue;
        }
        else if (static_cast<int32_t>(*s) > static_cast<int32_t>(maxValue))
        {
          *t = maxValue;
        }
        else
        {
          *t = static_cast<TargetType>(*s);
        }
      }
    }
  }


  template <typename SourceType>
  static void ConvertGrayscaleToFloat(ImageAccessor& target,
                                      const ImageAccessor& source)
  {
    assert(sizeof(float) == 4);

    const unsigned int width = source.GetWidth();
    const unsigned int height = source.GetHeight();
    
    for (unsigned int y = 0; y < height; y++)
    {
      float* t = reinterpret_cast<float*>(target.GetRow(y));
      const SourceType* s = reinterpret_cast<const SourceType*>(source.GetConstRow(y));

      for (unsigned int x = 0; x < width; x++, t++, s++)
      {
        *t = static_cast<float>(*s);
      }
    }
  }


  template <PixelFormat TargetFormat>
  static void ConvertFloatToGrayscale(ImageAccessor& target,
                                      const ImageAccessor& source)
  {
    typedef typename PixelTraits<TargetFormat>::PixelType  TargetType;
    
    assert(sizeof(float) == 4);

    const unsigned int width = source.GetWidth();
    const unsigned int height = source.GetHeight();

    for (unsigned int y = 0; y < height; y++)
    {
      TargetType* q = reinterpret_cast<TargetType*>(target.GetRow(y));
      const float* p = reinterpret_cast<const float*>(source.GetConstRow(y));

      for (unsigned int x = 0; x < width; x++, p++, q++)
      {
        PixelTraits<TargetFormat>::FloatToPixel(*q, *p);
      }
    }
  }


  template <typename TargetType>
  static void ConvertColorToGrayscale(ImageAccessor& target,
                                      const ImageAccessor& source)
  {
    assert(source.GetFormat() == PixelFormat_RGB24);

    // WARNING - "::min()" should be replaced by "::lowest()" if
    // dealing with float or double (which is not the case so far)
    assert(sizeof(TargetType) <= 2);  // Safeguard to remember about "float/double"
    const TargetType minValue = std::numeric_limits<TargetType>::min();
    const TargetType maxValue = std::numeric_limits<TargetType>::max();

    const unsigned int width = source.GetWidth();
    const unsigned int height = source.GetHeight();
    
    for (unsigned int y = 0; y < height; y++)
    {
      TargetType* t = reinterpret_cast<TargetType*>(target.GetRow(y));
      const uint8_t* s = reinterpret_cast<const uint8_t*>(source.GetConstRow(y));

      for (unsigned int x = 0; x < width; x++, t++, s += 3)
      {
        // Y = 0.2126 R + 0.7152 G + 0.0722 B
        int32_t v = (2126 * static_cast<int32_t>(s[0]) +
                     7152 * static_cast<int32_t>(s[1]) +
                     0722 * static_cast<int32_t>(s[2])) / 10000;
        
        if (static_cast<int32_t>(v) < static_cast<int32_t>(minValue))
        {
          *t = minValue;
        }
        else if (static_cast<int32_t>(v) > static_cast<int32_t>(maxValue))
        {
          *t = maxValue;
        }
        else
        {
          *t = static_cast<TargetType>(v);
        }
      }
    }
  }


  static void MemsetZeroInternal(ImageAccessor& image)
  {
    const unsigned int height = image.GetHeight();
    const size_t lineSize = image.GetBytesPerPixel() * image.GetWidth();
    const size_t pitch = image.GetPitch();

    uint8_t *p = reinterpret_cast<uint8_t*>(image.GetBuffer());
    
    for (unsigned int y = 0; y < height; y++)
    {
      memset(p, 0, lineSize);
      p += pitch;
    }
  }


  template <typename PixelType>
  static void SetInternal(ImageAccessor& image,
                          int64_t constant)
  {
    if (constant == 0 &&
        (image.GetFormat() == PixelFormat_Grayscale8 ||
         image.GetFormat() == PixelFormat_Grayscale16 ||
         image.GetFormat() == PixelFormat_Grayscale32 ||
         image.GetFormat() == PixelFormat_Grayscale64 ||
         image.GetFormat() == PixelFormat_SignedGrayscale16))
    {
      MemsetZeroInternal(image);
    }
    else
    {
      const unsigned int width = image.GetWidth();
      const unsigned int height = image.GetHeight();

      for (unsigned int y = 0; y < height; y++)
      {
        PixelType* p = reinterpret_cast<PixelType*>(image.GetRow(y));

        for (unsigned int x = 0; x < width; x++, p++)
        {
          *p = static_cast<PixelType>(constant);
        }
      }
    }
  }


  template <typename PixelType>
  static void GetMinMaxValueInternal(PixelType& minValue,
                                     PixelType& maxValue,
                                     const ImageAccessor& source,
                                     const PixelType LowestValue = std::numeric_limits<PixelType>::min())
  {
    // Deal with the special case of empty image
    if (source.GetWidth() == 0 ||
        source.GetHeight() == 0)
    {
      minValue = 0;
      maxValue = 0;
      return;
    }

    minValue = std::numeric_limits<PixelType>::max();
    maxValue = LowestValue;

    const unsigned int height = source.GetHeight();
    const unsigned int width = source.GetWidth();

    for (unsigned int y = 0; y < height; y++)
    {
      const PixelType* p = reinterpret_cast<const PixelType*>(source.GetConstRow(y));

      for (unsigned int x = 0; x < width; x++, p++)
      {
        if (*p < minValue)
        {
          minValue = *p;
        }

        if (*p > maxValue)
        {
          maxValue = *p;
        }
      }
    }
  }



  template <typename PixelType>
  static void AddConstantInternal(ImageAccessor& image,
                                  int64_t constant)
  {
    if (constant == 0)
    {
      return;
    }

    // WARNING - "::min()" should be replaced by "::lowest()" if
    // dealing with float or double (which is not the case so far)
    assert(sizeof(PixelType) <= 2);  // Safeguard to remember about "float/double"
    const int64_t minValue = std::numeric_limits<PixelType>::min();
    const int64_t maxValue = std::numeric_limits<PixelType>::max();

    const unsigned int width = image.GetWidth();
    const unsigned int height = image.GetHeight();
    
    for (unsigned int y = 0; y < height; y++)
    {
      PixelType* p = reinterpret_cast<PixelType*>(image.GetRow(y));

      for (unsigned int x = 0; x < width; x++, p++)
      {
        int64_t v = static_cast<int64_t>(*p) + constant;

        if (v > maxValue)
        {
          *p = std::numeric_limits<PixelType>::max();
        }
        else if (v < minValue)
        {
          *p = std::numeric_limits<PixelType>::min();
        }
        else
        {
          *p = static_cast<PixelType>(v);
        }
      }
    }
  }



  template <typename PixelType,
            bool UseRound>
  static void MultiplyConstantInternal(ImageAccessor& image,
                                       float factor)
  {
    if (std::abs(factor - 1.0f) <= std::numeric_limits<float>::epsilon())
    {
      return;
    }

    // WARNING - "::min()" should be replaced by "::lowest()" if
    // dealing with float or double (which is not the case so far)
    assert(sizeof(PixelType) <= 2);  // Safeguard to remember about "float/double"
    const int64_t minValue = std::numeric_limits<PixelType>::min();
    const int64_t maxValue = std::numeric_limits<PixelType>::max();

    const unsigned int width = image.GetWidth();
    const unsigned int height = image.GetHeight();

    for (unsigned int y = 0; y < height; y++)
    {
      PixelType* p = reinterpret_cast<PixelType*>(image.GetRow(y));

      for (unsigned int x = 0; x < width; x++, p++)
      {
        int64_t v;
        if (UseRound)
        {
          assert(sizeof(long long) == sizeof(int64_t));
          // The "round" operation is very costly
          v = boost::math::llround(static_cast<float>(*p) * factor);
        }
        else
        {
          v = static_cast<int64_t>(static_cast<float>(*p) * factor);
        }

        if (v > maxValue)
        {
          *p = std::numeric_limits<PixelType>::max();
        }
        else if (v < minValue)
        {
          *p = std::numeric_limits<PixelType>::min();
        }
        else
        {
          *p = static_cast<PixelType>(v);
        }
      }
    }
  }


  // Computes "a * x + b" at each pixel => Note that this is not the
  // same convention as in "ShiftScale()", but it is the convention of
  // "ShiftScale2()"
  template <typename TargetType,
            typename SourceType,
            bool UseRound,
            bool Invert>
  static void ShiftScaleIntegerInternal(ImageAccessor& target,
                                        const ImageAccessor& source,
                                        float a,
                                        float b)
  // This function can be applied inplace (source == target)
  {
    assert(target.GetFormat() != PixelFormat_Float32);
    
    if (source.GetWidth() != target.GetWidth() ||
        source.GetHeight() != target.GetHeight())
    {
      throw OrthancException(ErrorCode_IncompatibleImageSize);
    }

    if (&source == &target &&
        source.GetFormat() != target.GetFormat())
    {
      throw OrthancException(ErrorCode_IncompatibleImageFormat);
    }
    
    const TargetType minPixelValue = std::numeric_limits<TargetType>::min();
    const TargetType maxPixelValue = std::numeric_limits<TargetType>::max();
    const float minFloatValue = static_cast<float>(minPixelValue);
    const float maxFloatValue = static_cast<float>(maxPixelValue);

    const unsigned int height = target.GetHeight();
    const unsigned int width = target.GetWidth();
    
    for (unsigned int y = 0; y < height; y++)
    {
      TargetType* p = reinterpret_cast<TargetType*>(target.GetRow(y));
      const SourceType* q = reinterpret_cast<const SourceType*>(source.GetConstRow(y));

      for (unsigned int x = 0; x < width; x++, p++, q++)
      {
        float v = a * static_cast<float>(*q) + b;

        if (v >= maxFloatValue)
        {
          *p = maxPixelValue;
        }
        else if (v <= minFloatValue)
        {
          *p = minPixelValue;
        }
        else if (UseRound)
        {         
          // The "round" operation is very costly
          assert(sizeof(TargetType) < sizeof(int));
          *p = static_cast<TargetType>(boost::math::iround(v));
        }
        else
        {
          *p = static_cast<TargetType>(std::floor(v));
        }

        if (Invert)
        {
          *p = maxPixelValue - *p;
        }
      }
    }
  }


  template <typename SourceType>
  static void ShiftScaleFloatInternal(ImageAccessor& target,
                                      const ImageAccessor& source,
                                      float a,
                                      float b)
  // This function can be applied inplace (source == target)
  {
    assert(target.GetFormat() == PixelFormat_Float32);
    
    if (source.GetWidth() != target.GetWidth() ||
        source.GetHeight() != target.GetHeight())
    {
      throw OrthancException(ErrorCode_IncompatibleImageSize);
    }

    if (&source == &target &&
        source.GetFormat() != target.GetFormat())
    {
      throw OrthancException(ErrorCode_IncompatibleImageFormat);
    }
    
    const unsigned int height = target.GetHeight();
    const unsigned int width = target.GetWidth();
    
    for (unsigned int y = 0; y < height; y++)
    {
      float* p = reinterpret_cast<float*>(target.GetRow(y));
      const SourceType* q = reinterpret_cast<const SourceType*>(source.GetConstRow(y));

      for (unsigned int x = 0; x < width; x++, p++, q++)
      {
        *p = a * static_cast<float>(*q) + b;
      }
    }
  }
  

  template <typename PixelType>
  static void ShiftRightInternal(ImageAccessor& image,
                                 unsigned int shift)
  {
    const unsigned int height = image.GetHeight();
    const unsigned int width = image.GetWidth();

    for (unsigned int y = 0; y < height; y++)
    {
      PixelType* p = reinterpret_cast<PixelType*>(image.GetRow(y));

      for (unsigned int x = 0; x < width; x++, p++)
      {
        *p = *p >> shift;
      }
    }
  }

  template <typename PixelType>
  static void ShiftLeftInternal(ImageAccessor& image,
                                unsigned int shift)
  {
    const unsigned int height = image.GetHeight();
    const unsigned int width = image.GetWidth();

    for (unsigned int y = 0; y < height; y++)
    {
      PixelType* p = reinterpret_cast<PixelType*>(image.GetRow(y));

      for (unsigned int x = 0; x < width; x++, p++)
      {
        *p = *p << shift;
      }
    }
  }

  void ImageProcessing::Copy(ImageAccessor& target,
                             const ImageAccessor& source)
  {
    if (target.GetWidth() != source.GetWidth() ||
        target.GetHeight() != source.GetHeight())
    {
      throw OrthancException(ErrorCode_IncompatibleImageSize);
    }

    if (target.GetFormat() != source.GetFormat())
    {
      throw OrthancException(ErrorCode_IncompatibleImageFormat);
    }

    const unsigned int lineSize = source.GetBytesPerPixel() * source.GetWidth();
    assert(source.GetPitch() >= lineSize && target.GetPitch() >= lineSize);

    const unsigned int height = source.GetHeight();
    for (unsigned int y = 0; y < height; y++)
    {
      memcpy(target.GetRow(y), source.GetConstRow(y), lineSize);
    }
  }

  template <typename TargetType, typename SourceType>
  static void ApplyWindowingInternal(ImageAccessor& target,
                                     const ImageAccessor& source,
                                     float windowCenter,
                                     float windowWidth,
                                     float rescaleSlope,
                                     float rescaleIntercept,
                                     bool invert)
  {
    assert(sizeof(SourceType) == source.GetBytesPerPixel() &&
           sizeof(TargetType) == target.GetBytesPerPixel());
    
    const TargetType maxTargetValue = std::numeric_limits<TargetType>::max();
    const float maxFloatValue = static_cast<float>(maxTargetValue);
    
    const float windowIntercept = windowCenter - windowWidth / 2.0f;
    const float windowSlope = (maxFloatValue + 1.0f) / windowWidth;

    const float a = rescaleSlope * windowSlope;
    const float b = (rescaleIntercept - windowIntercept) * windowSlope;

    if (invert)
    {
      ShiftScaleIntegerInternal<TargetType, SourceType, false, true>(target, source, a, b);
    }
    else
    {
      ShiftScaleIntegerInternal<TargetType, SourceType, false, false>(target, source, a, b);
    }
  }

  void ImageProcessing::ApplyWindowing_Deprecated(ImageAccessor& target,
                                                  const ImageAccessor& source,
                                                  float windowCenter,
                                                  float windowWidth,
                                                  float rescaleSlope,
                                                  float rescaleIntercept,
                                                  bool invert)
  {
    if (target.GetWidth() != source.GetWidth() ||
        target.GetHeight() != source.GetHeight())
    {
      throw OrthancException(ErrorCode_IncompatibleImageSize);
    }

    switch (source.GetFormat())
    {
      case PixelFormat_Float32:
      {
        switch (target.GetFormat())
        {
          case PixelFormat_Grayscale8:
            ApplyWindowingInternal<uint8_t, float>(target, source, windowCenter, windowWidth, rescaleSlope, rescaleIntercept, invert);
            break;
          case PixelFormat_Grayscale16:
            ApplyWindowingInternal<uint16_t, float>(target, source, windowCenter, windowWidth, rescaleSlope, rescaleIntercept, invert);
            break;
          default:
            throw OrthancException(ErrorCode_NotImplemented);
        }
      };break;
      case PixelFormat_Grayscale8:
      {
        switch (target.GetFormat())
        {
          case PixelFormat_Grayscale8:
            ApplyWindowingInternal<uint8_t, uint8_t>(target, source, windowCenter, windowWidth, rescaleSlope, rescaleIntercept, invert);
            break;
          case PixelFormat_Grayscale16:
            ApplyWindowingInternal<uint16_t, uint8_t>(target, source, windowCenter, windowWidth, rescaleSlope, rescaleIntercept, invert);
            break;
          default:
            throw OrthancException(ErrorCode_NotImplemented);
        }
      };break;
      case PixelFormat_Grayscale16:
      {
        switch (target.GetFormat())
        {
          case PixelFormat_Grayscale8:
            ApplyWindowingInternal<uint8_t, uint16_t>(target, source, windowCenter, windowWidth, rescaleSlope, rescaleIntercept, invert);
            break;
          case PixelFormat_Grayscale16:
            ApplyWindowingInternal<uint16_t, uint16_t>(target, source, windowCenter, windowWidth, rescaleSlope, rescaleIntercept, invert);
            break;
          default:
            throw OrthancException(ErrorCode_NotImplemented);
        }
      };break;
      default:
        throw OrthancException(ErrorCode_NotImplemented);
    }
  }


  void ImageProcessing::Convert(ImageAccessor& target,
                                const ImageAccessor& source)
  {
    if (target.GetWidth() != source.GetWidth() ||
        target.GetHeight() != source.GetHeight())
    {
      throw OrthancException(ErrorCode_IncompatibleImageSize);
    }

    const unsigned int width = source.GetWidth();
    const unsigned int height = source.GetHeight();

    if (source.GetFormat() == target.GetFormat())
    {
      Copy(target, source);
      return;
    }

    if (target.GetFormat() == PixelFormat_Grayscale16 &&
        source.GetFormat() == PixelFormat_Grayscale8)
    {
      ConvertInternal<uint16_t, uint8_t>(target, source);
      return;
    }

    if (target.GetFormat() == PixelFormat_SignedGrayscale16 &&
        source.GetFormat() == PixelFormat_Grayscale8)
    {
      ConvertInternal<int16_t, uint8_t>(target, source);
      return;
    }

    if (target.GetFormat() == PixelFormat_Grayscale8 &&
        source.GetFormat() == PixelFormat_Grayscale16)
    {
      ConvertInternal<uint8_t, uint16_t>(target, source);
      return;
    }

    if (target.GetFormat() == PixelFormat_SignedGrayscale16 &&
        source.GetFormat() == PixelFormat_Grayscale16)
    {
      ConvertInternal<int16_t, uint16_t>(target, source);
      return;
    }

    if (target.GetFormat() == PixelFormat_Grayscale8 &&
        source.GetFormat() == PixelFormat_SignedGrayscale16)
    {
      ConvertInternal<uint8_t, int16_t>(target, source);
      return;
    }

    if (target.GetFormat() == PixelFormat_Grayscale16 &&
        source.GetFormat() == PixelFormat_SignedGrayscale16)
    {
      ConvertInternal<uint16_t, int16_t>(target, source);
      return;
    }

    if (target.GetFormat() == PixelFormat_Grayscale8 &&
        source.GetFormat() == PixelFormat_RGB24)
    {
      ConvertColorToGrayscale<uint8_t>(target, source);
      return;
    }

    if (target.GetFormat() == PixelFormat_Grayscale16 &&
        source.GetFormat() == PixelFormat_RGB24)
    {
      ConvertColorToGrayscale<uint16_t>(target, source);
      return;
    }

    if (target.GetFormat() == PixelFormat_SignedGrayscale16 &&
        source.GetFormat() == PixelFormat_RGB24)
    {
      ConvertColorToGrayscale<int16_t>(target, source);
      return;
    }

    if (target.GetFormat() == PixelFormat_Float32 &&
        source.GetFormat() == PixelFormat_Grayscale8)
    {
      ConvertGrayscaleToFloat<uint8_t>(target, source);
      return;
    }

    if (target.GetFormat() == PixelFormat_Float32 &&
        source.GetFormat() == PixelFormat_Grayscale16)
    {
      ConvertGrayscaleToFloat<uint16_t>(target, source);
      return;
    }

    if (target.GetFormat() == PixelFormat_Float32 &&
        source.GetFormat() == PixelFormat_Grayscale32)
    {
      ConvertGrayscaleToFloat<uint32_t>(target, source);
      return;
    }

    if (target.GetFormat() == PixelFormat_Float32 &&
        source.GetFormat() == PixelFormat_SignedGrayscale16)
    {
      ConvertGrayscaleToFloat<int16_t>(target, source);
      return;
    }

    
    if (target.GetFormat() == PixelFormat_Grayscale8 &&
        source.GetFormat() == PixelFormat_RGBA32)
    {
      for (unsigned int y = 0; y < height; y++)
      {
        const uint8_t* p = reinterpret_cast<const uint8_t*>(source.GetConstRow(y));
        uint8_t* q = reinterpret_cast<uint8_t*>(target.GetRow(y));
        for (unsigned int x = 0; x < width; x++, q++)
        {
          *q = static_cast<uint8_t>((2126 * static_cast<uint32_t>(p[0]) +
                                     7152 * static_cast<uint32_t>(p[1]) +
                                     0722 * static_cast<uint32_t>(p[2])) / 10000);
          p += 4;
        }
      }

      return;
    }

    if (target.GetFormat() == PixelFormat_Grayscale8 &&
        source.GetFormat() == PixelFormat_BGRA32)
    {
      for (unsigned int y = 0; y < height; y++)
      {
        const uint8_t* p = reinterpret_cast<const uint8_t*>(source.GetConstRow(y));
        uint8_t* q = reinterpret_cast<uint8_t*>(target.GetRow(y));
        for (unsigned int x = 0; x < width; x++, q++)
        {
          *q = static_cast<uint8_t>((2126 * static_cast<uint32_t>(p[2]) +
                                     7152 * static_cast<uint32_t>(p[1]) +
                                     0722 * static_cast<uint32_t>(p[0])) / 10000);
          p += 4;
        }
      }

      return;
    }

    if (target.GetFormat() == PixelFormat_RGB24 &&
        source.GetFormat() == PixelFormat_RGBA32)
    {
      for (unsigned int y = 0; y < height; y++)
      {
        const uint8_t* p = reinterpret_cast<const uint8_t*>(source.GetConstRow(y));
        uint8_t* q = reinterpret_cast<uint8_t*>(target.GetRow(y));
        for (unsigned int x = 0; x < width; x++)
        {
          q[0] = p[0];
          q[1] = p[1];
          q[2] = p[2];
          p += 4;
          q += 3;
        }
      }

      return;
    }

    if (target.GetFormat() == PixelFormat_RGB24 &&
        source.GetFormat() == PixelFormat_BGRA32)
    {
      for (unsigned int y = 0; y < height; y++)
      {
        const uint8_t* p = reinterpret_cast<const uint8_t*>(source.GetConstRow(y));
        uint8_t* q = reinterpret_cast<uint8_t*>(target.GetRow(y));
        for (unsigned int x = 0; x < width; x++)
        {
          q[0] = p[2];
          q[1] = p[1];
          q[2] = p[0];
          p += 4;
          q += 3;
        }
      }

      return;
    }

    if (target.GetFormat() == PixelFormat_RGBA32 &&
        source.GetFormat() == PixelFormat_RGB24)
    {
      for (unsigned int y = 0; y < height; y++)
      {
        const uint8_t* p = reinterpret_cast<const uint8_t*>(source.GetConstRow(y));
        uint8_t* q = reinterpret_cast<uint8_t*>(target.GetRow(y));
        for (unsigned int x = 0; x < width; x++)
        {
          q[0] = p[0];
          q[1] = p[1];
          q[2] = p[2];
          q[3] = 255;   // Set the alpha channel to full opacity
          p += 3;
          q += 4;
        }
      }

      return;
    }

    if (target.GetFormat() == PixelFormat_RGB24 &&
        source.GetFormat() == PixelFormat_Grayscale8)
    {
      for (unsigned int y = 0; y < height; y++)
      {
        const uint8_t* p = reinterpret_cast<const uint8_t*>(source.GetConstRow(y));
        uint8_t* q = reinterpret_cast<uint8_t*>(target.GetRow(y));
        for (unsigned int x = 0; x < width; x++)
        {
          q[0] = *p;
          q[1] = *p;
          q[2] = *p;
          p += 1;
          q += 3;
        }
      }

      return;
    }

    if ((target.GetFormat() == PixelFormat_RGBA32 ||
         target.GetFormat() == PixelFormat_BGRA32) &&
        source.GetFormat() == PixelFormat_Grayscale8)
    {
      for (unsigned int y = 0; y < height; y++)
      {
        const uint8_t* p = reinterpret_cast<const uint8_t*>(source.GetConstRow(y));
        uint8_t* q = reinterpret_cast<uint8_t*>(target.GetRow(y));
        for (unsigned int x = 0; x < width; x++)
        {
          q[0] = *p;
          q[1] = *p;
          q[2] = *p;
          q[3] = 255;
          p += 1;
          q += 4;
        }
      }

      return;
    }

    if (target.GetFormat() == PixelFormat_BGRA32 &&
        source.GetFormat() == PixelFormat_Grayscale16)
    {
      for (unsigned int y = 0; y < height; y++)
      {
        const uint16_t* p = reinterpret_cast<const uint16_t*>(source.GetConstRow(y));
        uint8_t* q = reinterpret_cast<uint8_t*>(target.GetRow(y));
        for (unsigned int x = 0; x < width; x++)
        {
          uint8_t value = (*p < 256 ? *p : 255);
          q[0] = value;
          q[1] = value;
          q[2] = value;
          q[3] = 255;
          p += 1;
          q += 4;
        }
      }

      return;
    }

    if (target.GetFormat() == PixelFormat_BGRA32 &&
        source.GetFormat() == PixelFormat_SignedGrayscale16)
    {
      for (unsigned int y = 0; y < height; y++)
      {
        const int16_t* p = reinterpret_cast<const int16_t*>(source.GetConstRow(y));
        uint8_t* q = reinterpret_cast<uint8_t*>(target.GetRow(y));
        for (unsigned int x = 0; x < width; x++)
        {
          uint8_t value;
          if (*p < 0)
          {
            value = 0;
          }
          else if (*p > 255)
          {
            value = 255;
          }
          else
          {
            value = static_cast<uint8_t>(*p);
          }

          q[0] = value;
          q[1] = value;
          q[2] = value;
          q[3] = 255;
          p += 1;
          q += 4;
        }
      }

      return;
    }

    if (target.GetFormat() == PixelFormat_BGRA32 &&
        source.GetFormat() == PixelFormat_RGB24)
    {
      for (unsigned int y = 0; y < height; y++)
      {
        const uint8_t* p = reinterpret_cast<const uint8_t*>(source.GetConstRow(y));
        uint8_t* q = reinterpret_cast<uint8_t*>(target.GetRow(y));
        for (unsigned int x = 0; x < width; x++)
        {
          q[0] = p[2];
          q[1] = p[1];
          q[2] = p[0];
          q[3] = 255;
          p += 3;
          q += 4;
        }
      }

      return;
    }

    if ((target.GetFormat() == PixelFormat_BGRA32 &&
         source.GetFormat() == PixelFormat_RGBA32)
        || (target.GetFormat() == PixelFormat_RGBA32 &&
            source.GetFormat() == PixelFormat_BGRA32))
    {
      for (unsigned int y = 0; y < height; y++)
      {
        const uint8_t* p = reinterpret_cast<const uint8_t*>(source.GetConstRow(y));
        uint8_t* q = reinterpret_cast<uint8_t*>(target.GetRow(y));
        for (unsigned int x = 0; x < width; x++)
        {
          q[0] = p[2];
          q[1] = p[1];
          q[2] = p[0];
          q[3] = p[3];
          p += 4;
          q += 4;
        }
      }

      return;
    }

    if (target.GetFormat() == PixelFormat_RGB24 &&
        source.GetFormat() == PixelFormat_RGB48)
    {
      for (unsigned int y = 0; y < height; y++)
      {
        const uint16_t* p = reinterpret_cast<const uint16_t*>(source.GetConstRow(y));
        uint8_t* q = reinterpret_cast<uint8_t*>(target.GetRow(y));
        for (unsigned int x = 0; x < width; x++)
        {
          q[0] = p[0] >> 8;
          q[1] = p[1] >> 8;
          q[2] = p[2] >> 8;
          p += 3;
          q += 3;
        }
      }

      return;
    }

    if (target.GetFormat() == PixelFormat_Grayscale16 &&
        source.GetFormat() == PixelFormat_Float32)
    {
      ConvertFloatToGrayscale<PixelFormat_Grayscale16>(target, source);
      return;
    }

    if (target.GetFormat() == PixelFormat_Grayscale8 &&
        source.GetFormat() == PixelFormat_Float32)
    {
      ConvertFloatToGrayscale<PixelFormat_Grayscale8>(target, source);
      return;
    }

    if (target.GetFormat() == PixelFormat_RGB24 &&
        source.GetFormat() == PixelFormat_Float32)
    {
      ConvertFloatToGrayscale<PixelFormat_RGB24>(target, source);
      return;
    }

    throw OrthancException(ErrorCode_NotImplemented);
  }



  void ImageProcessing::Set(ImageAccessor& image,
                            int64_t value)
  {
    switch (image.GetFormat())
    {
      case PixelFormat_Grayscale8:
        SetInternal<uint8_t>(image, value);
        return;

      case PixelFormat_Grayscale16:
        SetInternal<uint16_t>(image, value);
        return;

      case PixelFormat_Grayscale32:
        SetInternal<uint32_t>(image, value);
        return;

      case PixelFormat_Grayscale64:
        SetInternal<uint64_t>(image, value);
        return;

      case PixelFormat_SignedGrayscale16:
        SetInternal<int16_t>(image, value);
        return;

      case PixelFormat_Float32:
        assert(sizeof(float) == 4);
        SetInternal<float>(image, value);
        return;

      case PixelFormat_RGBA32:
      case PixelFormat_BGRA32:
      case PixelFormat_RGB24:
      {
        uint8_t v = static_cast<uint8_t>(value);
        Set(image, v, v, v, v);  // Use the color version
        return;
      }

      default:
        throw OrthancException(ErrorCode_NotImplemented);
    }
  }


  void ImageProcessing::Set(ImageAccessor& image,
                            uint8_t red,
                            uint8_t green,
                            uint8_t blue,
                            uint8_t alpha)
  {
    uint8_t p[4];
    unsigned int size;

    switch (image.GetFormat())
    {
      case PixelFormat_Grayscale8:
      {
        // New in Orthanc 1.9.0
        uint8_t grayscale = (2126 * static_cast<uint16_t>(red) + 
                             7152 * static_cast<uint16_t>(green) +
                             0722 * static_cast<uint16_t>(blue)) / 10000;
        Orthanc::ImageProcessing::Set(image, grayscale);
        return;
      }
      
      case PixelFormat_RGBA32:
        p[0] = red;
        p[1] = green;
        p[2] = blue;
        p[3] = alpha;
        size = 4;
        break;

      case PixelFormat_BGRA32:
        p[0] = blue;
        p[1] = green;
        p[2] = red;
        p[3] = alpha;
        size = 4;
        break;

      case PixelFormat_RGB24:
        p[0] = red;
        p[1] = green;
        p[2] = blue;
        size = 3;
        break;

      default:
        throw OrthancException(ErrorCode_NotImplemented);
    }

    const unsigned int width = image.GetWidth();
    const unsigned int height = image.GetHeight();

    for (unsigned int y = 0; y < height; y++)
    {
      uint8_t* q = reinterpret_cast<uint8_t*>(image.GetRow(y));

      for (unsigned int x = 0; x < width; x++)
      {
        for (unsigned int i = 0; i < size; i++)
        {
          q[i] = p[i];
        }

        q += size;
      }
    }
  }

  void ImageProcessing::Set(ImageAccessor& image,
                            uint8_t red,
                            uint8_t green,
                            uint8_t blue,
                            ImageAccessor& alpha)
  {
    uint8_t p[4];

    if (alpha.GetWidth() != image.GetWidth() || alpha.GetHeight() != image.GetHeight())
    {
      throw OrthancException(ErrorCode_IncompatibleImageSize);
    }

    if (alpha.GetFormat() != PixelFormat_Grayscale8)
    {
      throw OrthancException(ErrorCode_NotImplemented);
    }

    switch (image.GetFormat())
    {
      case PixelFormat_RGBA32:
        p[0] = red;
        p[1] = green;
        p[2] = blue;
        break;

      case PixelFormat_BGRA32:
        p[0] = blue;
        p[1] = green;
        p[2] = red;
        break;

      default:
        throw OrthancException(ErrorCode_NotImplemented);
    }

    const unsigned int width = image.GetWidth();
    const unsigned int height = image.GetHeight();

    for (unsigned int y = 0; y < height; y++)
    {
      uint8_t* q = reinterpret_cast<uint8_t*>(image.GetRow(y));
      uint8_t* a = reinterpret_cast<uint8_t*>(alpha.GetRow(y));

      for (unsigned int x = 0; x < width; x++)
      {
        for (unsigned int i = 0; i < 3; i++)
        {
          q[i] = p[i];
        }
        q[3] = *a;
        q += 4;
        ++a;
      }
    }
  }


  void ImageProcessing::ShiftRight(ImageAccessor& image,
                                   unsigned int shift)
  {
    if (image.GetWidth() == 0 ||
        image.GetHeight() == 0 ||
        shift == 0)
    {
      // Nothing to do
      return;
    }

    switch (image.GetFormat())
    {
      case PixelFormat_Grayscale8:
      {
        ShiftRightInternal<uint8_t>(image, shift);
        break;
      }

      case PixelFormat_Grayscale16:
      {
        ShiftRightInternal<uint16_t>(image, shift);
        break;
      }
      default:
        throw OrthancException(ErrorCode_NotImplemented);
    }
  }

  void ImageProcessing::ShiftLeft(ImageAccessor& image,
                                  unsigned int shift)
  {
    if (image.GetWidth() == 0 ||
        image.GetHeight() == 0 ||
        shift == 0)
    {
      // Nothing to do
      return;
    }

    switch (image.GetFormat())
    {
      case PixelFormat_Grayscale8:
      {
        ShiftLeftInternal<uint8_t>(image, shift);
        break;
      }

      case PixelFormat_Grayscale16:
      {
        ShiftLeftInternal<uint16_t>(image, shift);
        break;
      }
      default:
        throw OrthancException(ErrorCode_NotImplemented);
    }
  }

  void ImageProcessing::GetMinMaxIntegerValue(int64_t& minValue,
                                              int64_t& maxValue,
                                              const ImageAccessor& image)
  {
    switch (image.GetFormat())
    {
      case PixelFormat_Grayscale8:
      {
        uint8_t a, b;
        GetMinMaxValueInternal<uint8_t>(a, b, image);
        minValue = a;
        maxValue = b;
        break;
      }

      case PixelFormat_Grayscale16:
      {
        uint16_t a, b;
        GetMinMaxValueInternal<uint16_t>(a, b, image);
        minValue = a;
        maxValue = b;
        break;
      }

      case PixelFormat_Grayscale32:
      {
        uint32_t a, b;
        GetMinMaxValueInternal<uint32_t>(a, b, image);
        minValue = a;
        maxValue = b;
        break;
      }

      case PixelFormat_SignedGrayscale16:
      {
        int16_t a, b;
        GetMinMaxValueInternal<int16_t>(a, b, image);
        minValue = a;
        maxValue = b;
        break;
      }

      default:
        throw OrthancException(ErrorCode_NotImplemented);
    }
  }


  void ImageProcessing::GetMinMaxFloatValue(float& minValue,
                                            float& maxValue,
                                            const ImageAccessor& image)
  {
    switch (image.GetFormat())
    {
      case PixelFormat_Float32:
      {
        assert(sizeof(float) == 4);
        float a, b;

        /**
         * WARNING - On floating-point types, the minimal value is
         * "-FLT_MAX" (as implemented by "::lowest()"), not "FLT_MIN"
         * (as implemented by "::min()")
         * https://en.cppreference.com/w/cpp/types/numeric_limits
         **/
        GetMinMaxValueInternal<float>(a, b, image, -std::numeric_limits<float>::max());
        minValue = a;
        maxValue = b;
        break;
      }

      default:
        throw OrthancException(ErrorCode_NotImplemented);
    }
  }



  void ImageProcessing::AddConstant(ImageAccessor& image,
                                    int64_t value)
  {
    switch (image.GetFormat())
    {
      case PixelFormat_Grayscale8:
        AddConstantInternal<uint8_t>(image, value);
        return;

      case PixelFormat_Grayscale16:
        AddConstantInternal<uint16_t>(image, value);
        return;

      case PixelFormat_SignedGrayscale16:
        AddConstantInternal<int16_t>(image, value);
        return;

      default:
        throw OrthancException(ErrorCode_NotImplemented);
    }
  }


  void ImageProcessing::MultiplyConstant(ImageAccessor& image,
                                         float factor,
                                         bool useRound)
  {
    switch (image.GetFormat())
    {
      case PixelFormat_Grayscale8:
        if (useRound)
        {
          MultiplyConstantInternal<uint8_t, true>(image, factor);
        }
        else
        {
          MultiplyConstantInternal<uint8_t, false>(image, factor);
        }
        return;

      case PixelFormat_Grayscale16:
        if (useRound)
        {
          MultiplyConstantInternal<uint16_t, true>(image, factor);
        }
        else
        {
          MultiplyConstantInternal<uint16_t, false>(image, factor);
        }
        return;

      case PixelFormat_SignedGrayscale16:
        if (useRound)
        {
          MultiplyConstantInternal<int16_t, true>(image, factor);
        }
        else
        {
          MultiplyConstantInternal<int16_t, false>(image, factor);
        }
        return;

      default:
        throw OrthancException(ErrorCode_NotImplemented);
    }
  }


  static bool IsIdentityRescaling(float offset,
                                  float scaling)
  {
    return (std::abs(offset) <= 10.0f * std::numeric_limits<float>::epsilon() &&
            std::abs(scaling - 1.0f) <= 10.0f * std::numeric_limits<float>::epsilon());
  }
  

  void ImageProcessing::ShiftScale2(ImageAccessor& image,
                                    float offset,
                                    float scaling,
                                    bool useRound)
  {
    // We compute "a * x + b"
    const float a = scaling;
    const float b = offset;

    if (IsIdentityRescaling(offset, scaling))
    {
      return;
    }
    
    switch (image.GetFormat())
    {
      case PixelFormat_Grayscale8:
        if (useRound)
        {
          ShiftScaleIntegerInternal<uint8_t, uint8_t, true, false>(image, image, a, b);
        }
        else
        {
          ShiftScaleIntegerInternal<uint8_t, uint8_t, false, false>(image, image, a, b);
        }
        return;

      case PixelFormat_Grayscale16:
        if (useRound)
        {
          ShiftScaleIntegerInternal<uint16_t, uint16_t, true, false>(image, image, a, b);
        }
        else
        {
          ShiftScaleIntegerInternal<uint16_t, uint16_t, false, false>(image, image, a, b);
        }
        return;

      case PixelFormat_SignedGrayscale16:
        if (useRound)
        {
          ShiftScaleIntegerInternal<int16_t, int16_t, true, false>(image, image, a, b);
        }
        else
        {
          ShiftScaleIntegerInternal<int16_t, int16_t, false, false>(image, image, a, b);
        }
        return;

      case PixelFormat_Float32:
        if (useRound)
        {
          ShiftScaleFloatInternal<float>(image, image, a, b);
        }
        else
        {
          ShiftScaleFloatInternal<float>(image, image, a, b);
        }
        return;

      default:
        throw OrthancException(ErrorCode_NotImplemented);
    }
  }


  void ImageProcessing::ShiftScale2(ImageAccessor& target,
                                    const ImageAccessor& source,
                                    float offset,
                                    float scaling,
                                    bool useRound)
  {
    // We compute "a * x + b"
    const float a = scaling;
    const float b = offset;
    
    if (target.GetFormat() == source.GetFormat() &&
        IsIdentityRescaling(offset, scaling))
    {
      Copy(target, source);
      return;
    }
    
    switch (target.GetFormat())
    {
      case PixelFormat_Grayscale8:

        switch (source.GetFormat())
        {
          case PixelFormat_Float32:
            if (useRound)
            {
              ShiftScaleIntegerInternal<uint8_t, float, true, false>(target, source, a, b);
            }
            else
            {
              ShiftScaleIntegerInternal<uint8_t, float, false, false>(target, source, a, b);
            }
            return;

          default:
            throw OrthancException(ErrorCode_NotImplemented);
        }
        
      default:
        throw OrthancException(ErrorCode_NotImplemented);
    }
  }


  void ImageProcessing::ShiftScale(ImageAccessor& image,
                                   float offset,
                                   float scaling,
                                   bool useRound)
  {
    // Rewrite "(x + offset) * scaling" as "a * x + b"

    const float a = scaling;
    const float b = offset * scaling;
    ShiftScale2(image, b, a, useRound);
  }


  void ImageProcessing::ShiftScale(ImageAccessor& target,
                                   const ImageAccessor& source,
                                   float offset,
                                   float scaling,
                                   bool useRound)
  {
    // Rewrite "(x + offset) * scaling" as "a * x + b"

    const float a = scaling;
    const float b = offset * scaling;
    ShiftScale2(target, source, b, a, useRound);
  }



  void ImageProcessing::Invert(ImageAccessor& image, int64_t maxValue)
  {
    const unsigned int width = image.GetWidth();
    const unsigned int height = image.GetHeight();

    switch (image.GetFormat())
    {
      case PixelFormat_Grayscale16:
      {
        uint16_t maxValueUint16 = (uint16_t)(std::min(maxValue, static_cast<int64_t>(std::numeric_limits<uint16_t>::max())));

        for (unsigned int y = 0; y < height; y++)
        {
          uint16_t* p = reinterpret_cast<uint16_t*>(image.GetRow(y));

          for (unsigned int x = 0; x < width; x++, p++)
          {
            *p = maxValueUint16 - (*p);
          }
        }

        return;
      }
      case PixelFormat_Grayscale8:
      {
        uint8_t maxValueUint8 = (uint8_t)(std::min(maxValue, static_cast<int64_t>(std::numeric_limits<uint8_t>::max())));

        for (unsigned int y = 0; y < height; y++)
        {
          uint8_t* p = reinterpret_cast<uint8_t*>(image.GetRow(y));

          for (unsigned int x = 0; x < width; x++, p++)
          {
            *p = maxValueUint8 - (*p);
          }
        }

        return;
      }

      default:
        throw OrthancException(ErrorCode_NotImplemented);
    }

  }

  void ImageProcessing::Invert(ImageAccessor& image)
  {
    switch (image.GetFormat())
    {
      case PixelFormat_Grayscale8:
        return Invert(image, 255);
      default:
        throw OrthancException(ErrorCode_NotImplemented); // you should use the Invert(image, maxValue) overload
    }
  }



  namespace
  {
    template <PixelFormat Format>
    class BresenhamPixelWriter
    {
    private:
      typedef typename PixelTraits<Format>::PixelType  PixelType;

      ImageAccessor&  image_;
      PixelType                value_;

      void PlotLineLow(int x0,
                       int y0,
                       int x1,
                       int y1)
      {
        int dx = x1 - x0;
        int dy = y1 - y0;
        int yi = 1;

        if (dy < 0)
        {
          yi = -1;
          dy = -dy;
        }

        int d = 2 * dy - dx;
        int y = y0;

        for (int x = x0; x <= x1; x++)
        {
          Write(x, y);
          
          if (d > 0)
          {
            y = y + yi;
            d = d - 2 * dx;
          }

          d = d + 2*dy;
        }
      }
      
      void PlotLineHigh(int x0,
                        int y0,
                        int x1,
                        int y1)
      {
        int dx = x1 - x0;
        int dy = y1 - y0;
        int xi = 1;

        if (dx < 0)
        {
          xi = -1;
          dx = -dx;
        }

        int d = 2 * dx - dy;
        int x = x0;

        for (int y = y0; y <= y1; y++)
        {
          Write(x, y);
          
          if (d > 0)
          {
            x = x + xi;
            d = d - 2 * dy;
          }

          d = d + 2 * dx;
        }
      }

    public:
      BresenhamPixelWriter(ImageAccessor& image,
                           int64_t value) :
        image_(image),
        value_(PixelTraits<Format>::IntegerToPixel(value))
      {
      }

      BresenhamPixelWriter(ImageAccessor& image,
                           const PixelType& value) :
        image_(image),
        value_(value)
      {
      }

      void Write(int x,
                 int y)
      {
        if (x >= 0 &&
            y >= 0 &&
            static_cast<unsigned int>(x) < image_.GetWidth() &&
            static_cast<unsigned int>(y) < image_.GetHeight())
        {
          PixelType* p = reinterpret_cast<PixelType*>(image_.GetRow(y));
          p[x] = value_;
        }
      }

      void DrawSegment(int x0,
                       int y0,
                       int x1,
                       int y1)
      {
        // This is an implementation of Bresenham's line algorithm
        // https://en.wikipedia.org/wiki/Bresenham%27s_line_algorithm#All_cases

        if (abs(y1 - y0) < abs(x1 - x0))
        {
          if (x0 > x1)
          {
            PlotLineLow(x1, y1, x0, y0);
          }
          else
          {
            PlotLineLow(x0, y0, x1, y1);
          }
        }
        else
        {
          if (y0 > y1)
          {
            PlotLineHigh(x1, y1, x0, y0);
          }
          else
          {
            PlotLineHigh(x0, y0, x1, y1);
          }
        }
      }
    };
  }

  
  void ImageProcessing::DrawLineSegment(ImageAccessor& image,
                                        int x0,
                                        int y0,
                                        int x1,
                                        int y1,
                                        int64_t value)
  {
    switch (image.GetFormat())
    {
      case PixelFormat_Grayscale8:
      {
        BresenhamPixelWriter<PixelFormat_Grayscale8> writer(image, value);
        writer.DrawSegment(x0, y0, x1, y1);
        break;
      }

      case PixelFormat_Grayscale16:
      {
        BresenhamPixelWriter<PixelFormat_Grayscale16> writer(image, value);
        writer.DrawSegment(x0, y0, x1, y1);
        break;
      }

      case PixelFormat_SignedGrayscale16:
      {
        BresenhamPixelWriter<PixelFormat_SignedGrayscale16> writer(image, value);
        writer.DrawSegment(x0, y0, x1, y1);
        break;
      }

      default:
        throw OrthancException(ErrorCode_NotImplemented);
    }
  }

  
  void ImageProcessing::DrawLineSegment(ImageAccessor& image,
                                        int x0,
                                        int y0,
                                        int x1,
                                        int y1,
                                        uint8_t red,
                                        uint8_t green,
                                        uint8_t blue,
                                        uint8_t alpha)
  {
    switch (image.GetFormat())
    {
      case PixelFormat_BGRA32:
      {
        PixelTraits<PixelFormat_BGRA32>::PixelType pixel;
        pixel.red_ = red;
        pixel.green_ = green;
        pixel.blue_ = blue;
        pixel.alpha_ = alpha;

        BresenhamPixelWriter<PixelFormat_BGRA32> writer(image, pixel);
        writer.DrawSegment(x0, y0, x1, y1);
        break;
      }

      case PixelFormat_RGBA32:
      {
        PixelTraits<PixelFormat_RGBA32>::PixelType pixel;
        pixel.red_ = red;
        pixel.green_ = green;
        pixel.blue_ = blue;
        pixel.alpha_ = alpha;

        BresenhamPixelWriter<PixelFormat_RGBA32> writer(image, pixel);
        writer.DrawSegment(x0, y0, x1, y1);
        break;
      }

      case PixelFormat_RGB24:
      {
        PixelTraits<PixelFormat_RGB24>::PixelType pixel;
        pixel.red_ = red;
        pixel.green_ = green;
        pixel.blue_ = blue;

        BresenhamPixelWriter<PixelFormat_RGB24> writer(image, pixel);
        writer.DrawSegment(x0, y0, x1, y1);
        break;
      }

      default:
        throw OrthancException(ErrorCode_NotImplemented);
    }
  }

  void ComputePolygonExtent(int32_t& left, int32_t& right, int32_t& top, int32_t& bottom, const std::vector<ImageProcessing::ImagePoint>& points)
  {
    left = std::numeric_limits<int32_t>::max();
    right = std::numeric_limits<int32_t>::min();
    top = std::numeric_limits<int32_t>::max();
    bottom = std::numeric_limits<int32_t>::min();

    for (size_t i = 0; i < points.size(); i++)
    {
      const ImageProcessing::ImagePoint& p = points[i];
      left = std::min(p.GetX(), left);
      right = std::max(p.GetX(), right);
      bottom = std::max(p.GetY(), bottom);
      top = std::min(p.GetY(), top);
    }
  }


  namespace
  {
#define USE_POLYGON_FRACTIONS 1

    class PolygonEdge
    {
    private:
      int   yUpper;

#if USE_POLYGON_FRACTIONS == 1
      int  x;
      int  xOffset;
      int  dxPerScanNumerator;
      int  dxPerScanDenominator;
#else
      float xIntersect;
      float dxPerScan;
#endif

    public:
      PolygonEdge(const ImageProcessing::ImagePoint& lower,
                  const ImageProcessing::ImagePoint& upper,
                  int yComp)
      {
        // cf. "makeEdgeRec()" in textbook

        assert(upper.GetY() != lower.GetY());

#if USE_POLYGON_FRACTIONS == 1
        x = lower.GetX();
        xOffset = 0;
        dxPerScanNumerator = upper.GetX() - lower.GetX();
        dxPerScanDenominator = upper.GetY() - lower.GetY();
#else
        dxPerScan = (static_cast<float>(upper.GetX() - lower.GetX()) /
                     static_cast<float>(upper.GetY() - lower.GetY()));
        xIntersect = lower.GetX();
#endif
    
        if (upper.GetY() < yComp)
        {
          yUpper = upper.GetY() - 1;
        }
        else
        {
          yUpper = upper.GetY();
        }
      }

      void NextScanLine()
      {
#if USE_POLYGON_FRACTIONS == 1
        xOffset += dxPerScanNumerator;

        while (xOffset >= dxPerScanDenominator)
        {
          x++;
          xOffset -= dxPerScanDenominator;
        }

        while (xOffset < 0)
        {
          x--;
          xOffset += dxPerScanDenominator;
        }

#else
        xIntersect += dxPerScan;
#endif
      }


      int GetEnterX() const
      {
#if USE_POLYGON_FRACTIONS == 1
        assert(xOffset >= 0 && xOffset < dxPerScanDenominator);
        if (xOffset == 0)
        {
          return x;
        }
        else
        {
          return x + 1;
        }
#else
        return static_cast<int>(std::ceil(xIntersect));
#endif
      }
  
      int GetExitX() const
      {
#if USE_POLYGON_FRACTIONS == 1
        assert(xOffset >= 0 && xOffset < dxPerScanDenominator);
        return x;
#else
        return static_cast<int>(std::floor(xIntersect));
#endif
      }
  
      int GetUpperY() const
      {
        return yUpper;
      }

      bool operator< (const PolygonEdge& other) const
      {
#if USE_POLYGON_FRACTIONS == 1
        assert(xOffset >= 0 && xOffset < dxPerScanDenominator);
        assert(other.xOffset >= 0 && other.xOffset < other.dxPerScanDenominator);
        return x < other.x;
#else
        // cf. "insertEdge()" in textbook
        return (xIntersect < other.xIntersect);
#endif
      }
    };
  }
  

  // For an index, return y-coordinate of next nonhorizontal line
  static int GetPolygonNextY(const std::vector<ImageProcessing::ImagePoint>& points,
                             size_t k)
  {
    // cf. "yNext()" in textbook
    size_t j = k;
  
    for (;;)
    {
      j++;
      if (j == points.size())
      {
        j = 0;
      }

      if (points[k].GetY() != points[j].GetY())
      {
        return points[j].GetY();
      }
    }
  }


  static int GetPolygonPreviousY(const std::vector<ImageProcessing::ImagePoint>& points,
                                 size_t k)
  {
    size_t j = k;

    for (;;)
    {
      if (j > 0)
      {
        j --;
      }
      else
      {
        j = points.size() - 1;
      }

      if (points[k].GetY() != points[j].GetY())
      {
        return points[j].GetY();
      }
    }
  }



  void ImageProcessing::FillPolygon(IPolygonFiller& filler,
                                    const std::vector<ImagePoint>& points)
  {
    /**
     * This implementation is a C++ adaption of Section 3.11 (pages
     * 117-124) of textbook "Computer Graphics - C Version (2nd
     * Edition)" by Hearn and Baker, 1997.
     **/
  
    typedef std::map<int, std::list<PolygonEdge> > EdgeTable;

    if (points.size() < 2)
    {
      return;
    }

    bool onlyHorizontalSegments = true;
    for (size_t i = 1; i < points.size(); i++)
    {
      if (points[0].GetY() != points[i].GetY())
      {
        onlyHorizontalSegments = false;
        break;
      }
    }

    if (onlyHorizontalSegments)
    {
      // Degenerate case: There are only horizontal lines. If this is
      // the case, "GetPolygonPreviousY()" would be an infinite loop
      int x1 = points[0].GetX();
      int x2 = x1;
      for (size_t i = 1; i < points.size(); i++)
      {
        assert(points[i].GetY() == points[0].GetY());

        const int x = points[i].GetX();
        x1 = std::min(x1, x);
        x2 = std::max(x2, x);
      }
      filler.Fill(points[0].GetY(), x1, x2);
      return;
    }

    EdgeTable globalEdgeTable;

    // cf. "buildEdgeList()" in textbook

    // Error in the textbook: we use "GetPolygonPreviousY()" instead of "points.size() - 2"
    int yPrev = GetPolygonPreviousY(points, points.size() - 1);
    ImagePoint v1(points[points.size() - 1]);

    for (size_t i = 0; i < points.size(); i++)
    {
      ImagePoint v2(points[i]);

      if (v1.GetY() != v2.GetY())
      {
        // Non-horizontal line
        if (v1.GetY() < v2.GetY())
        {
          // Up-going edge
          PolygonEdge edge(v1, v2, GetPolygonNextY(points, i));
          globalEdgeTable[v1.GetY()].push_back(edge);
        }
        else if (v1.GetY() > v2.GetY())
        {
          // Down-going edge
          PolygonEdge edge(v2, v1, yPrev);
          globalEdgeTable[v2.GetY()].push_back(edge);
        }

        // Error in the textbook: "yPrev" must NOT be updated on horizontal lines
        yPrev = v1.GetY();
      }
      
      v1 = v2;
    }
    
    assert(!globalEdgeTable.empty());

    std::vector<PolygonEdge> activeEdges;

    for (EdgeTable::const_iterator it = globalEdgeTable.begin(); it != globalEdgeTable.end(); ++it)
    {
      // cf. "buildActiveList()" in textbook
      activeEdges.reserve(activeEdges.size() + it->second.size());
      for (std::list<PolygonEdge>::const_iterator it2 = it->second.begin(); it2 != it->second.end(); ++it2)
      {
        activeEdges.push_back(*it2);
      }

      assert(!activeEdges.empty());

      EdgeTable::const_iterator next = it;
      ++next;

      int rampEnd;
      if (next == globalEdgeTable.end())
      {
        rampEnd = activeEdges[0].GetUpperY() + 1;

        for (size_t i = 1; i < activeEdges.size(); i++)
        {
          rampEnd = std::max(rampEnd, activeEdges[i].GetUpperY() + 1);
        }
      }
      else
      {
        rampEnd = next->first;
      }

      for (int y = it->first; y < rampEnd; y++)
      {
        // cf. "updateActiveList()" in textbook
        std::vector<PolygonEdge> stillActive;
        stillActive.reserve(activeEdges.size());

        for (size_t i = 0; i < activeEdges.size(); i++)
        {
          if (y <= activeEdges[i].GetUpperY())
          {
            stillActive.push_back(activeEdges[i]);
          }
        }

        activeEdges.swap(stillActive);

        assert(activeEdges.size() % 2 == 0);
        std::sort(activeEdges.begin(), activeEdges.end());

        // cf. "fillScan()" in textbook
        for (size_t k = 0; k + 1 < activeEdges.size(); )
        {
          int a = activeEdges[k].GetExitX();
          int b = activeEdges[k + 1].GetEnterX();

          // Fix wrt. the textbook: merge overlapping segments
          k += 2;
          while (k + 1 < activeEdges.size() &&
                 activeEdges[k].GetExitX() == b)
          {
            assert(a <= b);
            b = activeEdges[k + 1].GetEnterX();
            k += 2;
          }

          assert(a <= b);
          filler.Fill(y, a, b);
        }

        // cf. "updateActiveList()" in textbook
        for (size_t k = 0; k < activeEdges.size(); k++)
        {
          activeEdges[k].NextScanLine();
        }
      }
    }
  }

  
  void ImageProcessing::FillPolygon(ImageAccessor& image,
                                    const std::vector<ImagePoint>& points,
                                    int64_t value)
  {
    class Filler : public IPolygonFiller
    {
    private:
      ImageAccessor& image_;
      int64_t        value_;

    public:
      Filler(ImageAccessor& image,
             int64_t value) :
        image_(image),
        value_(value)
      {
      }
      
      virtual void Fill(int y,
                        int x1,
                        int x2) ORTHANC_OVERRIDE
      {
        assert(x1 <= x2);
    
        if (x1 < static_cast<int>(image_.GetWidth()) &&
            x2 >= 0 &&
            y >= 0 &&
            y < static_cast<int>(image_.GetHeight()))
        {
          unsigned int yy = static_cast<unsigned int>(y);
          unsigned int a = static_cast<unsigned int>(std::max(0, x1));
          unsigned int b = static_cast<unsigned int>(std::min(x2, static_cast<int>(image_.GetWidth()) - 1));

          assert(a <= b);

          ImageAccessor region;
          image_.GetRegion(region, a, yy, b - a + 1, 1);
          Set(region, value_);
        }
      }
    };

    
    if (image.GetFormat() == PixelFormat_Grayscale8 ||
        image.GetFormat() == PixelFormat_Grayscale16 ||
        image.GetFormat() == PixelFormat_SignedGrayscale16)
    {
      Filler filler(image, value);
      FillPolygon(filler, points);
    }
    else
    {
      throw OrthancException(ErrorCode_NotImplemented);
    }
  }


  template <PixelFormat Format>
  static void ResizeInternal(ImageAccessor& target,
                             const ImageAccessor& source)
  {
    assert(target.GetFormat() == source.GetFormat() &&
           target.GetFormat() == Format);
      
    const unsigned int sourceWidth = source.GetWidth();
    const unsigned int sourceHeight = source.GetHeight();
    const unsigned int targetWidth = target.GetWidth();
    const unsigned int targetHeight = target.GetHeight();

    if (targetWidth == 0 || targetHeight == 0)
    {
      return;
    }

    if (sourceWidth == 0 || sourceHeight == 0)
    {
      // Avoids division by zero below
      ImageProcessing::Set(target, 0);
      return;
    }
      
    const float scaleX = static_cast<float>(sourceWidth) / static_cast<float>(targetWidth);
    const float scaleY = static_cast<float>(sourceHeight) / static_cast<float>(targetHeight);


    /**
     * Create two lookup tables to quickly know the (x,y) position
     * in the source image, given the (x,y) position in the target
     * image.
     **/
      
    std::vector<unsigned int>  lookupX(targetWidth);
      
    for (unsigned int x = 0; x < targetWidth; x++)
    {
      int sourceX = static_cast<int>(std::floor((static_cast<float>(x) + 0.5f) * scaleX));
      if (sourceX < 0)
      {
        sourceX = 0;  // Should never happen
      }
      else if (sourceX >= static_cast<int>(sourceWidth))
      {
        sourceX = sourceWidth - 1;
      }

      lookupX[x] = static_cast<unsigned int>(sourceX);
    }
      
    std::vector<unsigned int>  lookupY(targetHeight);
      
    for (unsigned int y = 0; y < targetHeight; y++)
    {
      int sourceY = static_cast<int>(std::floor((static_cast<float>(y) + 0.5f) * scaleY));
      if (sourceY < 0)
      {
        sourceY = 0;  // Should never happen
      }
      else if (sourceY >= static_cast<int>(sourceHeight))
      {
        sourceY = sourceHeight - 1;
      }

      lookupY[y] = static_cast<unsigned int>(sourceY);
    }


    /**
     * Actual resizing
     **/
      
    for (unsigned int targetY = 0; targetY < targetHeight; targetY++)
    {
      unsigned int sourceY = lookupY[targetY];

      for (unsigned int targetX = 0; targetX < targetWidth; targetX++)
      {
        unsigned int sourceX = lookupX[targetX];

        typename ImageTraits<Format>::PixelType pixel;
        ImageTraits<Format>::GetPixel(pixel, source, sourceX, sourceY);
        ImageTraits<Format>::SetPixel(target, pixel, targetX, targetY);
      }
    }            
  }



  void ImageProcessing::Resize(ImageAccessor& target,
                               const ImageAccessor& source)
  {
    if (source.GetFormat() != target.GetFormat())
    {
      throw OrthancException(ErrorCode_IncompatibleImageFormat);
    }

    if (source.GetWidth() == target.GetWidth() &&
        source.GetHeight() == target.GetHeight())
    {
      Copy(target, source);
      return;
    }

    switch (source.GetFormat())
    {
      case PixelFormat_Grayscale8:
        ResizeInternal<PixelFormat_Grayscale8>(target, source);
        break;

      case PixelFormat_Float32:
        ResizeInternal<PixelFormat_Float32>(target, source);
        break;

      case PixelFormat_RGB24:
        ResizeInternal<PixelFormat_RGB24>(target, source);
        break;

      default:
        throw OrthancException(ErrorCode_NotImplemented);
    }
  }


  ImageAccessor* ImageProcessing::Halve(const ImageAccessor& source,
                                        bool forceMinimalPitch)
  {
    std::unique_ptr<Image> target(new Image(source.GetFormat(), source.GetWidth() / 2,
                                            source.GetHeight() / 2, forceMinimalPitch));
    Resize(*target, source);
    return target.release();
  }

    
  template <PixelFormat Format>
  static void FlipXInternal(ImageAccessor& image)
  {     
    const unsigned int height = image.GetHeight();
    const unsigned int width = image.GetWidth();

    for (unsigned int y = 0; y < height; y++)
    {
      for (unsigned int x1 = 0; x1 < width / 2; x1++)
      {
        unsigned int x2 = width - 1 - x1;
          
        typename ImageTraits<Format>::PixelType a, b;
        ImageTraits<Format>::GetPixel(a, image, x1, y);
        ImageTraits<Format>::GetPixel(b, image, x2, y);
        ImageTraits<Format>::SetPixel(image, a, x2, y);
        ImageTraits<Format>::SetPixel(image, b, x1, y);
      }
    }        
  }

    
  void ImageProcessing::FlipX(ImageAccessor& image)
  {
    switch (image.GetFormat())
    {
      case PixelFormat_Grayscale8:
        FlipXInternal<PixelFormat_Grayscale8>(image);
        break;

      case PixelFormat_RGB24:
        FlipXInternal<PixelFormat_RGB24>(image);
        break;

      default:
        throw OrthancException(ErrorCode_NotImplemented);
    }
  }

    
  template <PixelFormat Format>
  static void FlipYInternal(ImageAccessor& image)
  {     
    const unsigned int height = image.GetHeight();
    const unsigned int width = image.GetWidth();

    for (unsigned int y1 = 0; y1 < height / 2; y1++)
    {
      unsigned int y2 = height - 1 - y1;
        
      for (unsigned int x = 0; x < width; x++)
      {
        typename ImageTraits<Format>::PixelType a, b;
        ImageTraits<Format>::GetPixel(a, image, x, y1);
        ImageTraits<Format>::GetPixel(b, image, x, y2);
        ImageTraits<Format>::SetPixel(image, a, x, y2);
        ImageTraits<Format>::SetPixel(image, b, x, y1);
      }
    }        
  }

    
  void ImageProcessing::FlipY(ImageAccessor& image)
  {
    switch (image.GetFormat())
    {
      case PixelFormat_Grayscale8:
        FlipYInternal<PixelFormat_Grayscale8>(image);
        break;

      case PixelFormat_RGB24:
        FlipYInternal<PixelFormat_RGB24>(image);
        break;

      default:
        throw OrthancException(ErrorCode_NotImplemented);
    }
  }


  // This is a slow implementation of horizontal convolution on one
  // individual channel, that checks for out-of-image values
  template <typename RawPixel, unsigned int ChannelsCount>
  static float GetHorizontalConvolutionFloatSecure(const ImageAccessor& source,
                                                   const std::vector<float>& horizontal,
                                                   size_t horizontalAnchor,
                                                   unsigned int x,
                                                   unsigned int y,
                                                   float leftBorder,
                                                   float rightBorder,
                                                   unsigned int channel)
  {
    const RawPixel* row = reinterpret_cast<const RawPixel*>(source.GetConstRow(y)) + channel;

    float p = 0;

    for (unsigned int k = 0; k < horizontal.size(); k++)
    {
      float value;

      if (x + k < horizontalAnchor)   // Negation of "x - horizontalAnchor + k >= 0"
      {
        value = leftBorder;
      }
      else if (x + k >= source.GetWidth() + horizontalAnchor)   // Negation of "x - horizontalAnchor + k < width"
      {
        value = rightBorder;
      }
      else
      {
        // The value lies within the image
        value = row[(x - horizontalAnchor + k) * ChannelsCount];
      }

      p += value * horizontal[k];
    }

    return p;
  }
  

  
  // This is an implementation of separable convolution that uses
  // floating-point arithmetics, and an intermediate Float32
  // image. The out-of-image values are taken as the border
  // value. Further optimization is possible.
  template <typename RawPixel, unsigned int ChannelsCount, bool UseRound>
  static void SeparableConvolutionFloat(ImageAccessor& image /* inplace */,
                                        const std::vector<float>& horizontal,
                                        size_t horizontalAnchor,
                                        const std::vector<float>& vertical,
                                        size_t verticalAnchor,
                                        float normalization)
  {
    // WARNING - "::min()" should be replaced by "::lowest()" if
    // dealing with float or double (which is not the case so far)
    assert(sizeof(RawPixel) <= 2);  // Safeguard to remember about "float/double"

    const unsigned int width = image.GetWidth();
    const unsigned int height = image.GetHeight();
    

    /**
     * Horizontal convolution
     **/

    Image tmp(PixelFormat_Float32, ChannelsCount * width, height, false);

    for (unsigned int y = 0; y < height; y++)
    {
      const RawPixel* row = reinterpret_cast<const RawPixel*>(image.GetConstRow(y));

      float leftBorder[ChannelsCount], rightBorder[ChannelsCount];
      
      for (unsigned int c = 0; c < ChannelsCount; c++)
      {
        leftBorder[c] = row[c];
        rightBorder[c] = row[ChannelsCount * (width - 1) + c];
      }

      float* p = static_cast<float*>(tmp.GetRow(y));

      if (width < horizontal.size())
      {
        // It is not possible to have the full kernel within the image, use the direct implementation
        for (unsigned int x = 0; x < width; x++)
        {
          for (unsigned int c = 0; c < ChannelsCount; c++, p++)
          {
            *p = GetHorizontalConvolutionFloatSecure<RawPixel, ChannelsCount>
              (image, horizontal, horizontalAnchor, x, y, leftBorder[c], rightBorder[c], c);
          }
        }
      }
      else
      {
        // Deal with the left border
        for (unsigned int x = 0; x < horizontalAnchor; x++)
        {
          for (unsigned int c = 0; c < ChannelsCount; c++, p++)
          {
            *p = GetHorizontalConvolutionFloatSecure<RawPixel, ChannelsCount>
              (image, horizontal, horizontalAnchor, x, y, leftBorder[c], rightBorder[c], c);
          }
        }

        // Deal with the central portion of the image (all pixel values
        // scanned by the kernel lie inside the image)

        for (unsigned int x = 0; x < width - horizontal.size() + 1; x++)
        {
          for (unsigned int c = 0; c < ChannelsCount; c++, p++)
          {
            *p = 0;
            for (unsigned int k = 0; k < horizontal.size(); k++)
            {
              *p += static_cast<float>(row[(x + k) * ChannelsCount + c]) * horizontal[k];
            }
          }
        }

        // Deal with the right border
        for (unsigned int x = static_cast<unsigned int>(
               horizontalAnchor + width - horizontal.size() + 1); x < width; x++)
        {
          for (unsigned int c = 0; c < ChannelsCount; c++, p++)
          {
            *p = GetHorizontalConvolutionFloatSecure<RawPixel, ChannelsCount>
              (image, horizontal, horizontalAnchor, x, y, leftBorder[c], rightBorder[c], c);
          }
        }
      }
    }


    /**
     * Vertical convolution
     **/

    std::vector<const float*> rows(vertical.size());

    for (unsigned int y = 0; y < height; y++)
    {
      for (unsigned int k = 0; k < vertical.size(); k++)
      {
        if (y + k < verticalAnchor)
        {
          rows[k] = reinterpret_cast<const float*>(tmp.GetConstRow(0));   // Use top border
        }
        else if (y + k >= height + verticalAnchor)
        {
          rows[k] = reinterpret_cast<const float*>(tmp.GetConstRow(height - 1));  // Use bottom border
        }
        else
        {
          rows[k] = reinterpret_cast<const float*>(tmp.GetConstRow(static_cast<unsigned int>(y + k - verticalAnchor)));
        }
      }

      RawPixel* p = reinterpret_cast<RawPixel*>(image.GetRow(y));
        
      for (unsigned int x = 0; x < width; x++)
      {
        for (unsigned int c = 0; c < ChannelsCount; c++, p++)
        {
          float accumulator = 0;
        
          for (unsigned int k = 0; k < vertical.size(); k++)
          {
            accumulator += rows[k][ChannelsCount * x + c] * vertical[k];
          }

          accumulator *= normalization;

          if (accumulator <= static_cast<float>(std::numeric_limits<RawPixel>::min()))
          {
            *p = std::numeric_limits<RawPixel>::min();
          }
          else if (accumulator >= static_cast<float>(std::numeric_limits<RawPixel>::max()))
          {
            *p = std::numeric_limits<RawPixel>::max();
          }
          else
          {
            if (UseRound)
            {
              assert(sizeof(RawPixel) < sizeof(int));
              *p = static_cast<RawPixel>(boost::math::iround(accumulator));
            }
            else
            {
              *p = static_cast<RawPixel>(accumulator);
            }
          }
        }
      }
    }
  }


  void ImageProcessing::SeparableConvolution(ImageAccessor& image /* inplace */,
                                             const std::vector<float>& horizontal,
                                             size_t horizontalAnchor,
                                             const std::vector<float>& vertical,
                                             size_t verticalAnchor,
                                             bool useRound)
  {
    if (horizontal.size() == 0 ||
        vertical.size() == 0 ||
        horizontalAnchor >= horizontal.size() ||
        verticalAnchor >= vertical.size())
    {
      throw OrthancException(ErrorCode_ParameterOutOfRange);
    }
    
    if (image.GetWidth() == 0 ||
        image.GetHeight() == 0)
    {
      return;
    }
    
    /**
     * Compute normalization
     **/
    
    float sumHorizontal = 0;
    for (size_t i = 0; i < horizontal.size(); i++)
    {
      sumHorizontal += horizontal[i];
    }
    
    float sumVertical = 0;
    for (size_t i = 0; i < vertical.size(); i++)
    {
      sumVertical += vertical[i];
    }

    if (fabsf(sumHorizontal) <= std::numeric_limits<float>::epsilon() ||
        fabsf(sumVertical) <= std::numeric_limits<float>::epsilon())
    {
      throw OrthancException(ErrorCode_ParameterOutOfRange, "Singular convolution kernel");
    }      

    const float normalization = 1.0f / (sumHorizontal * sumVertical);

    switch (image.GetFormat())
    {
      case PixelFormat_Grayscale8:
        if (useRound)
        {
          SeparableConvolutionFloat<uint8_t, 1u, true>
            (image, horizontal, horizontalAnchor, vertical, verticalAnchor, normalization);
        }
        else
        {
          SeparableConvolutionFloat<uint8_t, 1u, false>
            (image, horizontal, horizontalAnchor, vertical, verticalAnchor, normalization);
        }
        break;

      case PixelFormat_RGB24:
        if (useRound)
        {
          SeparableConvolutionFloat<uint8_t, 3u, true>
            (image, horizontal, horizontalAnchor, vertical, verticalAnchor, normalization);
        }
        else
        {
          SeparableConvolutionFloat<uint8_t, 3u, false>
            (image, horizontal, horizontalAnchor, vertical, verticalAnchor, normalization);
        }
        break;

      default:
        throw OrthancException(ErrorCode_NotImplemented);
    }
  }


  void ImageProcessing::SmoothGaussian5x5(ImageAccessor& image,
                                          bool useRound)
  {
    std::vector<float> kernel(5);
    kernel[0] = 1;
    kernel[1] = 4;
    kernel[2] = 6;
    kernel[3] = 4;
    kernel[4] = 1;

    SeparableConvolution(image, kernel, 2, kernel, 2, useRound);
  }


  void ImageProcessing::FitSize(ImageAccessor& target,
                                const ImageAccessor& source)
  {
    if (target.GetWidth() == 0 ||
        target.GetHeight() == 0)
    {
      return;
    }

    if (source.GetWidth() == target.GetWidth() &&
        source.GetHeight() == target.GetHeight())
    {
      Copy(target, source);
      return;
    }

    Set(target, 0);

    // Preserve the aspect ratio
    float cw = static_cast<float>(source.GetWidth());
    float ch = static_cast<float>(source.GetHeight());
    float r = std::min(
      static_cast<float>(target.GetWidth()) / cw,
      static_cast<float>(target.GetHeight()) / ch);

    unsigned int sw = std::min(static_cast<unsigned int>(boost::math::iround(cw * r)), target.GetWidth());  
    unsigned int sh = std::min(static_cast<unsigned int>(boost::math::iround(ch * r)), target.GetHeight());

    Image resized(target.GetFormat(), sw, sh, false);
  
    //ImageProcessing::SmoothGaussian5x5(source);
    ImageProcessing::Resize(resized, source);

    assert(target.GetWidth() >= resized.GetWidth() &&
           target.GetHeight() >= resized.GetHeight());
    unsigned int offsetX = (target.GetWidth() - resized.GetWidth()) / 2;
    unsigned int offsetY = (target.GetHeight() - resized.GetHeight()) / 2;

    ImageAccessor region;
    target.GetRegion(region, offsetX, offsetY, resized.GetWidth(), resized.GetHeight());
    ImageProcessing::Copy(region, resized);
  }


  ImageAccessor* ImageProcessing::FitSize(const ImageAccessor& source,
                                          unsigned int width,
                                          unsigned int height)
  {
    std::unique_ptr<ImageAccessor> target(new Image(source.GetFormat(), width, height, false));
    FitSize(*target, source);
    return target.release();
  }

    
  ImageAccessor* ImageProcessing::FitSizeKeepAspectRatio(const ImageAccessor& source,
                                                         unsigned int width,
                                                         unsigned int height)
  {
    std::unique_ptr<ImageAccessor> target(new Image(source.GetFormat(), width, height, false));
    Set(*target, 0);

    if (width != 0 &&
        height != 0 &&
        source.GetWidth() != 0 &&
        source.GetHeight() != 0)
    {
      float ratio = std::min(static_cast<float>(width) / static_cast<float>(source.GetWidth()),
                             static_cast<float>(height) / static_cast<float>(source.GetHeight()));

      unsigned int resizedWidth = static_cast<unsigned int>(
        boost::math::iround(ratio * static_cast<float>(source.GetWidth())));

      unsigned int resizedHeight = static_cast<unsigned int>(
        boost::math::iround(ratio * static_cast<float>(source.GetHeight())));

      std::unique_ptr<ImageAccessor> resized(FitSize(source, resizedWidth, resizedHeight));

      ImageAccessor region;
      target->GetRegion(region, (width - resizedWidth) / 2,
                        (height - resizedHeight) / 2, resizedWidth, resizedHeight);
      Copy(region, *resized);
    }

    return target.release();
  }


  void ImageProcessing::ConvertJpegYCbCrToRgb(ImageAccessor& image)
  {
    // http://dicom.nema.org/medical/dicom/current/output/chtml/part03/sect_C.7.6.3.html#sect_C.7.6.3.1.2
    // https://en.wikipedia.org/wiki/YCbCr#JPEG_conversion
    
    // TODO - Check out the outcome of Mathieu's discussion about
    // truncation of YCbCr-to-RGB conversion:
    // https://groups.google.com/forum/#!msg/comp.protocols.dicom/JHuGeyWbTz8/ARoTWrJzAQAJ

    const unsigned int width = image.GetWidth();
    const unsigned int height = image.GetHeight();
    const unsigned int pitch = image.GetPitch();
    uint8_t* buffer = reinterpret_cast<uint8_t*>(image.GetBuffer());
        
    if (image.GetFormat() != PixelFormat_RGB24 ||
        pitch < 3 * width)
    {
      throw OrthancException(ErrorCode_IncompatibleImageFormat);
    }

    for (unsigned int y = 0; y < height; y++)
    {
      uint8_t* p = buffer + y * pitch;
          
      for (unsigned int x = 0; x < width; x++, p += 3)
      {
        const float Y  = p[0];
        const float Cb = p[1];
        const float Cr = p[2];

        const float result[3] = {
          Y                             + 1.402f    * (Cr - 128.0f),
          Y - 0.344136f * (Cb - 128.0f) - 0.714136f * (Cr - 128.0f),
          Y + 1.772f    * (Cb - 128.0f)
        };

        for (uint8_t i = 0; i < 3 ; i++)
        {
          if (result[i] < 0)
          {
            p[i] = 0;
          }
          else if (result[i] > 255)
          {
            p[i] = 255;
          }
          else
          {
            p[i] = static_cast<uint8_t>(result[i]);
          }
        }    
      }
    }
  }


  void ImageProcessing::SwapEndianness(ImageAccessor& image /* inplace */)
  {
    const unsigned int width = image.GetWidth();
    const unsigned int height = image.GetHeight();
    
    switch (image.GetFormat())
    {
      case PixelFormat_Grayscale8:
      case PixelFormat_RGB24:
      case PixelFormat_RGBA32:
      case PixelFormat_BGRA32:
        // No swapping required
        break;

      case PixelFormat_Grayscale16:
      case PixelFormat_SignedGrayscale16:
        for (unsigned int y = 0; y < height; y++)
        {
          uint8_t* t = reinterpret_cast<uint8_t*>(image.GetRow(y));
          for (unsigned int x = 0; x < width; x++)
          {
            uint8_t a = t[0];
            t[0] = t[1];
            t[1] = a;
            t += 2;
          }
        }    
        break;

      case PixelFormat_Grayscale32:
      case PixelFormat_Float32:
        for (unsigned int y = 0; y < height; y++)
        {
          uint8_t* t = reinterpret_cast<uint8_t*>(image.GetRow(y));
          for (unsigned int x = 0; x < width; x++)
          {
            uint8_t a = t[0];
            uint8_t b = t[1];
            t[0] = t[3];
            t[1] = t[2];
            t[2] = b;
            t[3] = a;
            t += 4;
          }
        }    
        break;

      case PixelFormat_RGB48:  // uint16_t per channel
        for (unsigned int y = 0; y < height; y++)
        {
          uint8_t* t = reinterpret_cast<uint8_t*>(image.GetRow(y));
          for (unsigned int x = 0; x < 3 * width; x++)
          {
            uint8_t a = t[0];
            t[0] = t[1];
            t[1] = a;
            t += 2;
          }
        }    
        break;

      default:
        throw OrthancException(ErrorCode_NotImplemented);
    }          
  }


  template <typename PixelType,
            typename Functor>
  static void ApplyImageOntoImage(Functor f,
                                  ImageAccessor& image /* inout */,
                                  const ImageAccessor& other)
  {
    const unsigned int width = image.GetWidth();
    const unsigned int height = image.GetHeight();
    
    if (width != other.GetWidth() ||
        height != other.GetHeight())
    {
      throw OrthancException(ErrorCode_IncompatibleImageSize);
    }
    else if (image.GetFormat() != other.GetFormat() ||
             GetBytesPerPixel(image.GetFormat()) != sizeof(PixelType))
    {
      throw OrthancException(ErrorCode_IncompatibleImageFormat);
    }
    else
    {
      for (unsigned int y = 0; y < height; y++)
      {
        PixelType* p = reinterpret_cast<PixelType*>(image.GetRow(y));
        const PixelType* q = reinterpret_cast<const PixelType*>(other.GetConstRow(y));
        
        for (unsigned int x = 0; x < width; x++, p++, q++)
        {
          f(*p, *q);
        }
      }
    }
  }


  namespace
  {
    // For older version of gcc, templated functors cannot be defined
    // as types internal to functions, hence the anonymous namespace
    
    struct MaximumFunctor
    {
      void operator() (uint8_t& a, const uint8_t& b)
      {
        a = std::max(a, b);
      }

      void operator() (uint16_t& a, const uint16_t& b)
      {
        a = std::max(a, b);
      }
    };
  }
  

  void ImageProcessing::Maximum(ImageAccessor& image,
                                const ImageAccessor& other)
  {
    switch (image.GetFormat())
    {
      case PixelFormat_Grayscale8:
        ApplyImageOntoImage<uint8_t, MaximumFunctor>(MaximumFunctor(), image, other);
        return;

      case PixelFormat_Grayscale16:
        ApplyImageOntoImage<uint16_t, MaximumFunctor>(MaximumFunctor(), image, other);
        return;

      default:
        throw OrthancException(ErrorCode_NotImplemented);
    }
  }
}