#include "trop.h"

#include "tpixel.h"

#include "tpixelutils.h"

// calls to _mm_* functions disabled in code for now (marked as comment)

// so disable include <emmintrin.h></emmintrin.h>

/*

#if defined(_WIN32) && defined(x64)

  #define ENABLE_SSE2

#endif

#ifdef ENABLE_SSE2

#include <emmintrin.h>  // per SSE2</emmintrin.h>

#endif

 */

namespace {

inline double luminance(TPixel32 *pix) {

  return 0.2126 * pix->r + 0.7152 * pix->g + 0.0722 * pix->b;

}

inline double luminance(TPixel64 *pix) {

  return 0.2126 * pix->r + 0.7152 * pix->g + 0.0722 * pix->b;

}

}  // namespace

//-----------------------------------------------------------------------------

#define FOR_EACH_PIXEL_BEGIN_LOOP(UpType, up, DownType, down, OutType, out)    \

  {                                                                            \

    int upWrap   = up->getWrap();                                              \

    int downWrap = down->getWrap();                                            \

    int outWrap  = out->getWrap();                                             \

                                                                               \

    up->lock();                                                                \

    down->lock();                                                              \

    out->lock();                                                               \

    UpType *upPix = 0, *upRow   = up->pixels();                                \

    DownType *downPix, *downRow = down->pixels();                              \

    OutType *outPix, *outRow    = out->pixels();                               \

    UpType *endPix;                                                            \

    int upLx        = up->getLx();                                             \

    UpType *lastPix = upRow + upWrap * (up->getLy() - 1) + upLx;               \

    while (upPix < lastPix) {                                                  \

      upPix   = upRow;                                                         \

      downPix = downRow;                                                       \

      outPix  = outRow;                                                        \

      endPix  = upPix + upLx;                                                  \

      while (upPix < endPix) {

//-----------------------------------------------------------------------------

#define FOR_EACH_PIXEL_END_LOOP(up, down, out)                                 \

  ++upPix;                                                                     \

  ++downPix;                                                                   \

  ++outPix;                                                                    \

  }                                                                            \

  upRow += upWrap;                                                             \

  downRow += downWrap;                                                         \

  outRow += outWrap;                                                           \

  }                                                                            \

  up->unlock();                                                                \

  down->unlock();                                                              \

  out->unlock();                                                               \

  }

//-----------------------------------------------------------------------------

#define FOR_EACH_PIXEL_32_BEGIN_LOOP                                           \

  assert(up32 &&down32 &&out32);                                               \

  FOR_EACH_PIXEL_BEGIN_LOOP(TPixelRGBM32, up32, TPixelRGBM32, down32,          \

                            TPixelRGBM32, out32)

//-----------------------------------------------------------------------------

#define FOR_EACH_PIXEL_32_END_LOOP                                             \

  assert(up32 &&down32 &&out32);                                               \

  FOR_EACH_PIXEL_END_LOOP(up32, down32, out32)

//-----------------------------------------------------------------------------

#define FOR_EACH_PIXEL_64_BEGIN_LOOP                                           \

  assert(up64 &&down64 &&out64);                                               \

  FOR_EACH_PIXEL_BEGIN_LOOP(TPixelRGBM64, up64, TPixelRGBM64, down64,          \

                            TPixelRGBM64, out64)

//-----------------------------------------------------------------------------

#define FOR_EACH_PIXEL_64_END_LOOP                                             \

  assert(up64 &&down64 &&out64);                                               \

  FOR_EACH_PIXEL_END_LOOP(up64, down64, out64)

//-----------------------------------------------------------------------------

#define FOR_EACH_PIXEL_F_BEGIN_LOOP                                            \

  assert(upF &&downF &&outF);                                                  \

  FOR_EACH_PIXEL_BEGIN_LOOP(TPixelF, upF, TPixelF, downF, TPixelF, outF)

//-----------------------------------------------------------------------------

#define FOR_EACH_PIXEL_F_END_LOOP                                              \

  assert(upF &&downF &&outF);                                                  \

  FOR_EACH_PIXEL_END_LOOP(upF, downF, outF)

//-----------------------------------------------------------------------------

#define FOR_EACH_PIXEL_8_BEGIN_LOOP                                            \

  assert(up8 &&down8 &&out8);                                                  \

  FOR_EACH_PIXEL_BEGIN_LOOP(TPixelGR8, up8, TPixelGR8, down8, TPixelGR8, out8)

//-----------------------------------------------------------------------------

#define FOR_EACH_PIXEL_8_END_LOOP                                              \

  assert(up8 &&down8 &&out8);                                                  \

  FOR_EACH_PIXEL_END_LOOP(up32, down32, out32)

//-----------------------------------------------------------------------------

void TRop::add(const TRasterP &rup, const TRasterP &rdown, const TRasterP &rout,

               double v) {

  TRaster32P up32   = rup;

  TRaster32P down32 = rdown;

  TRaster32P out32  = rout;

  if (up32 && down32 && out32) {

    FOR_EACH_PIXEL_32_BEGIN_LOOP

    TINT32 r, g, b, m;

    if (upPix->m == 0)

      *outPix = *downPix;

    else {

      r         = downPix->r + tround(upPix->r * v);

      g         = downPix->g + tround(upPix->g * v);

      b         = downPix->b + tround(upPix->b * v);

      m         = downPix->m + tround(upPix->m * v);

      outPix->r = (UCHAR)tcrop<tint32>(r, (TINT32)0, (TINT32)255);</tint32>

      outPix->g = (UCHAR)tcrop<tint32>(g, (TINT32)0, (TINT32)255);</tint32>

      outPix->b = (UCHAR)tcrop<tint32>(b, (TINT32)0, (TINT32)255);</tint32>

      outPix->m = (UCHAR)tcrop<tint32>(m, (TINT32)0, (TINT32)255);</tint32>

    }

    FOR_EACH_PIXEL_32_END_LOOP

    return;

  }

  TRaster64P up64   = rup;

  TRaster64P down64 = rdown;

  TRaster64P out64  = rout;

  if (up64 && down64 && out64) {

    FOR_EACH_PIXEL_64_BEGIN_LOOP

    TINT32 r, g, b, m;

    r = downPix->r + tround(upPix->r * v);

    g = downPix->g + tround(upPix->g * v);

    b = downPix->b + tround(upPix->b * v);

    m = downPix->m + tround(upPix->m * v);

    outPix->r = (USHORT)tcrop<tint32>(r, 0, 0xffff);</tint32>

    outPix->g = (USHORT)tcrop<tint32>(g, 0, 0xffff);</tint32>

    outPix->b = (USHORT)tcrop<tint32>(b, 0, 0xffff);</tint32>

    outPix->m = (USHORT)tcrop<tint32>(m, 0, 0xffff);</tint32>

    FOR_EACH_PIXEL_64_END_LOOP

    return;

  }

  TRasterFP upF   = rup;

  TRasterFP downF = rdown;

  TRasterFP outF  = rout;

  if (upF && downF && outF) {

    FOR_EACH_PIXEL_F_BEGIN_LOOP

    outPix->r = downPix->r + upPix->r * v;

    outPix->g = downPix->g + upPix->g * v;

    outPix->b = downPix->b + upPix->b * v;

    outPix->m = tcrop<float>(downPix->m + upPix->m * v, 0.f, 1.f);</float>

    FOR_EACH_PIXEL_F_END_LOOP

    return;

  }

  TRasterGR8P up8   = rup;

  TRasterGR8P down8 = rdown;

  TRasterGR8P out8  = rout;

  if (up8 && down8 && out8) {

    FOR_EACH_PIXEL_8_BEGIN_LOOP

    USHORT value = troundp(upPix->value * v) + downPix->value;

    outPix->value = (UCHAR)tcrop<ushort>(value, 0, 255);</ushort>

    FOR_EACH_PIXEL_8_END_LOOP

    return;

  }

  throw TRopException("TRop::add invalid raster combination");

}

//-----------------------------------------------------------------------------

void TRop::add(const TRasterP &rup, const TRasterP &rdown,

               const TRasterP &rout) {

  TRaster32P up32   = rup;

  TRaster32P down32 = rdown;

  TRaster32P out32  = rout;

  if (up32 && down32 && out32) {

    FOR_EACH_PIXEL_32_BEGIN_LOOP

    USHORT r, g, b, m;

    r = downPix->r + upPix->r;

    g = downPix->g + upPix->g;

    b = downPix->b + upPix->b;

    m = downPix->m + upPix->m;

    outPix->r = (UCHAR)tcrop<ushort>(r, 0, 255);</ushort>

    outPix->g = (UCHAR)tcrop<ushort>(g, 0, 255);</ushort>

    outPix->b = (UCHAR)tcrop<ushort>(b, 0, 255);</ushort>

    outPix->m = (UCHAR)tcrop<ushort>(m, 0, 255);</ushort>

    FOR_EACH_PIXEL_32_END_LOOP

    return;

  }

  TRaster64P up64   = rup;

  TRaster64P down64 = rdown;

  TRaster64P out64  = rout;

  if (up64 && down64 && out64) {

    FOR_EACH_PIXEL_64_BEGIN_LOOP

    TINT32 r, g, b, m;

    r = downPix->r + upPix->r;

    g = downPix->g + upPix->g;

    b = downPix->b + upPix->b;

    m = downPix->m + upPix->m;

    outPix->r = (USHORT)tcrop<tint32>(r, 0, 0xffff);</tint32>

    outPix->g = (USHORT)tcrop<tint32>(g, 0, 0xffff);</tint32>

    outPix->b = (USHORT)tcrop<tint32>(b, 0, 0xffff);</tint32>

    outPix->m = (USHORT)tcrop<tint32>(m, 0, 0xffff);</tint32>

    FOR_EACH_PIXEL_64_END_LOOP

    return;

  }

  TRasterFP upF   = rup;

  TRasterFP downF = rdown;

  TRasterFP outF  = rout;

  if (upF && downF && outF) {

    FOR_EACH_PIXEL_F_BEGIN_LOOP

    outPix->r = downPix->r + upPix->r;

    outPix->g = downPix->g + upPix->g;

    outPix->b = downPix->b + upPix->b;

    outPix->m = tcrop<float>(downPix->m + upPix->m, 0.f, 1.f);</float>

    FOR_EACH_PIXEL_F_END_LOOP

    return;

  }

  TRasterGR8P up8   = rup;

  TRasterGR8P down8 = rdown;

  TRasterGR8P out8  = rout;

  if (up8 && down8 && out8) {

    FOR_EACH_PIXEL_8_BEGIN_LOOP

    USHORT value = upPix->value + downPix->value;

    outPix->value = (UCHAR)tcrop<ushort>(value, 0, 255);</ushort>

    FOR_EACH_PIXEL_8_END_LOOP

    return;

  }

  throw TRopException("TRop::add invalid raster combination");

}

//-----------------------------------------------------------------------------

void TRop::colordodge(const TRasterP &rup, const TRasterP &rdown,

                      const TRasterP &rout) {

  TRaster32P up32   = rup;

  TRaster32P down32 = rdown;

  TRaster32P out32  = rout;

  if (up32 && down32 && out32) {

    FOR_EACH_PIXEL_32_BEGIN_LOOP

    USHORT r, g, b, m;

    r = (USHORT)((downPix->r << 8) / (256.0 - upPix->r));

    g = (USHORT)((downPix->g << 8) / (256.0 - upPix->g));

    b = (USHORT)((downPix->b << 8) / (256.0 - upPix->b));

    m = downPix->m + upPix->m;

    outPix->r = (UCHAR)tcrop<ushort>(r, 0, 255);</ushort>

    outPix->g = (UCHAR)tcrop<ushort>(g, 0, 255);</ushort>

    outPix->b = (UCHAR)tcrop<ushort>(b, 0, 255);</ushort>

    outPix->m = (UCHAR)tcrop<ushort>(m, 0, 255);</ushort>

    FOR_EACH_PIXEL_32_END_LOOP

  } else {

    TRaster64P up64   = rup;

    TRaster64P down64 = rdown;

    TRaster64P out64  = rout;

    if (up64 && down64 && out64) {

      FOR_EACH_PIXEL_64_BEGIN_LOOP

      TINT32 r, g, b, m;

      r = (TINT32)(65536.0 * (downPix->r / (65536.0 - upPix->r)));

      g = (TINT32)(65536.0 * (downPix->g / (65536.0 - upPix->g)));

      b = (TINT32)(65536.0 * (downPix->b / (65536.0 - upPix->b)));

      m = downPix->m + upPix->m;

      outPix->r = (USHORT)tcrop<tint32>(r, 0, 0xffff);</tint32>

      outPix->g = (USHORT)tcrop<tint32>(g, 0, 0xffff);</tint32>

      outPix->b = (USHORT)tcrop<tint32>(b, 0, 0xffff);</tint32>

      outPix->m = (USHORT)tcrop<tint32>(m, 0, 0xffff);</tint32>

      FOR_EACH_PIXEL_64_END_LOOP

    } else {

      TRasterFP upF   = rup;

      TRasterFP downF = rdown;

      TRasterFP outF  = rout;

      if (upF && downF && outF) {

        FOR_EACH_PIXEL_F_BEGIN_LOOP

        outPix->r = downPix->r / (1.f - upPix->r);

        outPix->g = downPix->g / (1.f - upPix->g);

        outPix->b = downPix->b / (1.f - upPix->b);

        outPix->m = tcrop<float>(downPix->m + upPix->m, 0.f, 1.f);</float>

        FOR_EACH_PIXEL_F_END_LOOP

      } else {

        TRasterGR8P up8   = rup;

        TRasterGR8P down8 = rdown;

        TRasterGR8P out8  = rout;

        if (up8 && down8 && out8) {

          FOR_EACH_PIXEL_8_BEGIN_LOOP

          USHORT value;

          if (downPix->value)

            value = (USHORT)((downPix->value << 8) / (255.0 - upPix->value));

          outPix->value = (UCHAR)tcrop<ushort>(value, 0, 255);</ushort>

          FOR_EACH_PIXEL_8_END_LOOP

        } else

          throw TRopException("TRop::color dodge invalid raster combination");

      }

    }

  }

}

//-----------------------------------------------------------------------------

void TRop::colorburn(const TRasterP &rup, const TRasterP &rdown,

                     const TRasterP &rout) {

  TRaster32P up32   = rup;

  TRaster32P down32 = rdown;

  TRaster32P out32  = rout;

  if (up32 && down32 && out32) {

    FOR_EACH_PIXEL_32_BEGIN_LOOP

    //  downPix->r=0;

    //  downPix->g=255;

    //  downPix->b=0;

    double r, g, b;

    if (upPix->m) {

      if (downPix->r == 0 || downPix->r == 255)

        r = downPix->r;

      else if (upPix->r)

        r = 255 - (((255 - downPix->r) << 8) / (double)upPix->r);

      else

        r = 0;

      if (downPix->g == 0 || downPix->g == 255)

        g = downPix->g;

      else if (upPix->g)

        g = 255 - (((255 - downPix->g) << 8) / (double)upPix->g);

      else

        g = 0;

      if (downPix->b == 0 || downPix->b == 255)

        b = downPix->b;

      else if (upPix->b)

        b = 255 - (((255 - downPix->b) << 8) / (double)upPix->b);

      else

        b = 0;

      if (upPix->m != 255) {

        TPixel32 tmpPix;

        tmpPix.r = (UCHAR)tcrop<double>(r, .0, 255.0);</double>

        tmpPix.g = (UCHAR)tcrop<double>(g, .0, 255.0);</double>

        tmpPix.b = (UCHAR)tcrop<double>(b, .0, 255.0);</double>

        tmpPix.m = upPix->m;

        ;

        overPix<tpixel32, uchar="">(*outPix, *downPix, tmpPix);</tpixel32,>

      } else {

        outPix->r = (UCHAR)tcrop<double>(r, .0, 255.0);</double>

        outPix->g = (UCHAR)tcrop<double>(g, .0, 255.0);</double>

        outPix->b = (UCHAR)tcrop<double>(b, .0, 255.0);</double>

        outPix->m = downPix->m;

      }

    } else {

      outPix = downPix;

    }

    FOR_EACH_PIXEL_32_END_LOOP

  } else {

    TRaster64P up64   = rup;

    TRaster64P down64 = rdown;

    TRaster64P out64  = rout;

    if (up64 && down64 && out64) {

      FOR_EACH_PIXEL_64_BEGIN_LOOP

      double r, g, b;

      if (upPix->m) {

        if (downPix->r == 0 || downPix->r == 65535)

          r = downPix->r;

        else if (upPix->r)

          r = 65535 - 65536 * ((65535.0 - downPix->r) / (double)upPix->r);

        else

          r = 0;

        if (downPix->g == 0 || downPix->g == 65535)

          g = downPix->g;

        else if (upPix->g)

          g = 65535 - 65536 * ((65535.0 - downPix->g) / (double)upPix->g);

        else

          g = 0;

        if (downPix->b == 0 || downPix->b == 65535)

          b = downPix->b;

        else if (upPix->b)

          b = 65535 - 65536 * ((65535.0 - downPix->b) / (double)upPix->b);

        else

          b = 0;

        if (upPix->m != 65535) {

          TPixel64 tmpPix;

          tmpPix.r = (USHORT)tcrop<double>(r, .0, 65535.0);</double>

          tmpPix.g = (USHORT)tcrop<double>(g, .0, 65535.0);</double>

          tmpPix.b = (USHORT)tcrop<double>(b, .0, 65535.0);</double>

          tmpPix.m = upPix->m;

          overPix<tpixel64, ushort="">(*outPix, *downPix, tmpPix);</tpixel64,>

        } else {

          outPix->r = (USHORT)tcrop<double>(r, .0, 65535.0);</double>

          outPix->g = (USHORT)tcrop<double>(g, .0, 65535.0);</double>

          outPix->b = (USHORT)tcrop<double>(b, .0, 65535.0);</double>

          outPix->m = downPix->m;

        }

      } else {

        outPix = downPix;

      }

      FOR_EACH_PIXEL_64_END_LOOP

    } else {

      TRasterFP upF   = rup;

      TRasterFP downF = rdown;

      TRasterFP outF  = rout;

      if (upF && downF && outF) {

        FOR_EACH_PIXEL_F_BEGIN_LOOP

        float r, g, b;

        if (upPix->m > 0.f) {

          if (downPix->r <= 0.f || downPix->r >= 1.f)

            r = downPix->r;

          else if (upPix->r)

            r = 1.f - (1.f - downPix->r) / upPix->r;

          else

            r = 0.f;

          if (downPix->g <= 0.f || downPix->g >= 1.f)

            g = downPix->g;

          else if (upPix->g)

            g = 1.f - (1.f - downPix->g) / upPix->g;

          else

            g = 0.f;

          if (downPix->b <= 0.f || downPix->b >= 1.f)

            b = downPix->b;

          else if (upPix->b)

            b = 1.f - (1.f - downPix->b) / upPix->b;

          else

            b = 0.f;

          if (upPix->m < 1.f) {

            overPix<tpixelf, float="">(*outPix, *downPix,</tpixelf,>

                                    TPixelF(r, g, b, upPix->m));

          } else {

            outPix->r = r;

            outPix->g = g;

            outPix->b = b;

            outPix->m = downPix->m;

          }

        } else {

          outPix = downPix;

        }

        FOR_EACH_PIXEL_F_END_LOOP

      } else

        throw TRopException("TRop::color burn invalid raster combination");

    }

  }

}

//-----------------------------------------------------------------------------

void TRop::screen(const TRasterP &rup, const TRasterP &rdown,

                  const TRasterP &rout) {

  TRaster32P up32   = rup;

  TRaster32P down32 = rdown;

  TRaster32P out32  = rout;

  if (up32 && down32 && out32) {

    FOR_EACH_PIXEL_32_BEGIN_LOOP

    double r, g, b;

    r = 256 - ((256 - upPix->r) * (256 - downPix->r) >> 8);

    g = 256 - ((256 - upPix->g) * (256 - downPix->g) >> 8);

    b = 256 - ((256 - upPix->b) * (256 - downPix->b) >> 8);

    if (upPix->m != 255) {

      double m;

      m = 256 - ((256 - upPix->m) * (256 - downPix->m) >> 8);

      TPixel32 tmpPix;

      tmpPix.r = (UCHAR)tcrop<double>(r, 0, 255);</double>

      tmpPix.g = (UCHAR)tcrop<double>(g, 0, 255);</double>

      tmpPix.b = (UCHAR)tcrop<double>(b, 0, 255);</double>

      tmpPix.m = (UCHAR)tcrop<double>(m, 0, 255);</double>

      overPix<tpixel32, uchar="">(*outPix, *downPix, tmpPix);</tpixel32,>

    } else {

      outPix->r = (UCHAR)tcrop<double>(r, 0, 255);</double>

      outPix->g = (UCHAR)tcrop<double>(g, 0, 255);</double>

      outPix->b = (UCHAR)tcrop<double>(b, 0, 255);</double>

      outPix->m = upPix->m;

    }

    FOR_EACH_PIXEL_32_END_LOOP

    return;

  }

  TRaster64P up64   = rup;

  TRaster64P down64 = rdown;

  TRaster64P out64  = rout;

  if (up64 && down64 && out64) {

    FOR_EACH_PIXEL_64_BEGIN_LOOP

    double r, g, b;

    r = 65536 - (65536 - upPix->r) * ((65536 - downPix->r) / 65536.0);

    g = 65536 - (65536 - upPix->g) * ((65536 - downPix->g) / 65536.0);

    b = 65536 - (65536 - upPix->b) * ((65536 - downPix->b) / 65536.0);

    if (upPix->m != 65535) {

      double m;

      m = 65536 - (65536 - upPix->m) * ((65536 - downPix->m) / 65536.0);

      TPixel64 tmpPix;

      tmpPix.r = (USHORT)tcrop<double>(r, 0, 65535);</double>

      tmpPix.g = (USHORT)tcrop<double>(g, 0, 65535);</double>

      tmpPix.b = (USHORT)tcrop<double>(b, 0, 65535);</double>

      tmpPix.m = (USHORT)tcrop<double>(m, 0, 65535);</double>

      overPix<tpixel64, ushort="">(*outPix, *downPix, tmpPix);</tpixel64,>

    } else {

      outPix->r = (USHORT)tcrop<double>(r, 0, 65535);</double>

      outPix->g = (USHORT)tcrop<double>(g, 0, 65535);</double>

      outPix->b = (USHORT)tcrop<double>(b, 0, 65535);</double>

      outPix->m = upPix->m;

    }

    FOR_EACH_PIXEL_64_END_LOOP

    return;

  }

  TRasterFP upF   = rup;

  TRasterFP downF = rdown;

  TRasterFP outF  = rout;

  if (upF && downF && outF) {

    FOR_EACH_PIXEL_F_BEGIN_LOOP

    float r, g, b;

    r = 1.f - (1.f - upPix->r) * (1.f - downPix->r);

    g = 1.f - (1.f - upPix->g) * (1.f - downPix->g);

    b = 1.f - (1.f - upPix->b) * (1.f - downPix->b);

    if (upPix->m <= 1.f) {

      float m;

      m = 1.f - (1.f - upPix->m) * (1.f - downPix->m);

      TPixelF tmpPix(r, g, b, tcrop<float>(m, 0.f, 1.f));</float>

      overPix<tpixelf, float="">(*outPix, *downPix, tmpPix);</tpixelf,>

    } else {

      outPix->r = r;

      outPix->g = g;

      outPix->b = b;

      outPix->m = upPix->m;

    }

    FOR_EACH_PIXEL_F_END_LOOP

    return;

  }

  TRasterGR8P up8   = rup;

  TRasterGR8P down8 = rdown;

  TRasterGR8P out8  = rout;

  if (up8 && down8 && out8) {

    FOR_EACH_PIXEL_8_BEGIN_LOOP

    USHORT value;

    if (downPix->value)

      value = (USHORT)((downPix->value << 8) / (255.0 - upPix->value));

    outPix->value = (UCHAR)tcrop<ushort>(value, 0, 255);</ushort>

    FOR_EACH_PIXEL_8_END_LOOP

    return;

  }

  throw TRopException("TRop::color dodge invalid raster combination");

}

//-----------------------------------------------------------------------------

void TRop::sub(const TRasterP &rup, const TRasterP &rdown, const TRasterP &rout,

               bool matte) {

  TRaster32P up32   = rup;

  TRaster32P down32 = rdown;

  TRaster32P out32  = rout;

  if (matte) {

    if (up32 && down32 && out32) {

      FOR_EACH_PIXEL_32_BEGIN_LOOP

      SHORT r = downPix->r - upPix->r;

      SHORT g = downPix->g - upPix->g;

      SHORT b = downPix->b - upPix->b;

      SHORT m = downPix->m - upPix->m;

      outPix->r = (UCHAR)tcrop<short>(r, 0, 255);</short>

      outPix->g = (UCHAR)tcrop<short>(g, 0, 255);</short>

      outPix->b = (UCHAR)tcrop<short>(b, 0, 255);</short>

      outPix->m = (UCHAR)tcrop<short>(m, 0, 255);</short>

      FOR_EACH_PIXEL_32_END_LOOP

      return;

    }

    TRaster64P up64   = rup;

    TRaster64P down64 = rdown;

    TRaster64P out64  = rout;

    if (up64 && down64 && out64) {

      FOR_EACH_PIXEL_64_BEGIN_LOOP

      TINT32 r  = downPix->r - upPix->r;

      TINT32 g  = downPix->g - upPix->g;

      TINT32 b  = downPix->b - upPix->b;

      TINT32 m  = downPix->m - upPix->m;

      outPix->r = (USHORT)tcrop<tint32>(r, 0, 0xffff);</tint32>

      outPix->g = (USHORT)tcrop<tint32>(g, 0, 0xffff);</tint32>

      outPix->b = (USHORT)tcrop<tint32>(b, 0, 0xffff);</tint32>

      outPix->m = (USHORT)tcrop<tint32>(m, 0, 0xffff);</tint32>

      FOR_EACH_PIXEL_64_END_LOOP

      return;

    }

    TRasterFP upF   = rup;

    TRasterFP downF = rdown;

    TRasterFP outF  = rout;

    if (upF && downF && outF) {

      FOR_EACH_PIXEL_F_BEGIN_LOOP

      outPix->r = downPix->r - upPix->r;

      outPix->g = downPix->g - upPix->g;

      outPix->b = downPix->b - upPix->b;

      outPix->m = tcrop<float>(downPix->m - upPix->m, 0.f, 1.f);</float>

      FOR_EACH_PIXEL_F_END_LOOP

      return;

    }

    TRasterGR8P up8   = rup;

    TRasterGR8P down8 = rdown;

    TRasterGR8P out8  = rout;

    if (up8 && down8 && out8) {

      FOR_EACH_PIXEL_8_BEGIN_LOOP

      SHORT value   = upPix->value - downPix->value;

      outPix->value = (UCHAR)tcrop<short>(value, 0, 255);</short>

      FOR_EACH_PIXEL_8_END_LOOP

      return;

    }

    throw TRopException("TRop::sub invalid raster combination");

  } else {

    if (up32 && down32 && out32) {

      FOR_EACH_PIXEL_32_BEGIN_LOOP

      SHORT r = downPix->r - upPix->r;

      SHORT g = downPix->g - upPix->g;

      SHORT b = downPix->b - upPix->b;

      SHORT m = downPix->m;  // - upPix->m;

      outPix->r = (UCHAR)tcrop<short>(r, 0, 255);</short>

      outPix->g = (UCHAR)tcrop<short>(g, 0, 255);</short>

      outPix->b = (UCHAR)tcrop<short>(b, 0, 255);</short>

      outPix->m = (UCHAR)tcrop<short>(m, 0, 255);</short>

      FOR_EACH_PIXEL_32_END_LOOP

      return;

    }

    TRaster64P up64   = rup;

    TRaster64P down64 = rdown;

    TRaster64P out64  = rout;

    if (up64 && down64 && out64) {

      FOR_EACH_PIXEL_64_BEGIN_LOOP

      TINT32 r  = downPix->r - upPix->r;

      TINT32 g  = downPix->g - upPix->g;

      TINT32 b  = downPix->b - upPix->b;

      TINT32 m  = downPix->m;  // - upPix->m;

      outPix->r = (USHORT)tcrop<tint32>(r, 0, 0xffff);</tint32>

      outPix->g = (USHORT)tcrop<tint32>(g, 0, 0xffff);</tint32>

      outPix->b = (USHORT)tcrop<tint32>(b, 0, 0xffff);</tint32>

      outPix->m = (USHORT)tcrop<tint32>(m, 0, 0xffff);</tint32>

      FOR_EACH_PIXEL_64_END_LOOP

      return;

    }

    TRasterFP upF   = rup;

    TRasterFP downF = rdown;

    TRasterFP outF  = rout;

    if (upF && downF && outF) {

      FOR_EACH_PIXEL_F_BEGIN_LOOP

      outPix->r = downPix->r - upPix->r;

      outPix->g = downPix->g - upPix->g;

      outPix->b = downPix->b - upPix->b;

      outPix->m = tcrop<float>(downPix->m, 0.f, 1.f);</float>

      FOR_EACH_PIXEL_F_END_LOOP

      return;

    }

    TRasterGR8P up8   = rup;

    TRasterGR8P down8 = rdown;

    TRasterGR8P out8  = rout;

    if (up8 && down8 && out8) {

      FOR_EACH_PIXEL_8_BEGIN_LOOP

      SHORT value   = upPix->value - downPix->value;

      outPix->value = (UCHAR)tcrop<short>(value, 0, 255);</short>

      FOR_EACH_PIXEL_8_END_LOOP

      return;

    }

    throw TRopException("TRop::sub invalid raster combination");

  }

}

//-----------------------------------------------------------------------------

void TRop::mult(const TRasterP &rup, const TRasterP &rdown,

                const TRasterP &rout, int v, bool matte) {

  /*

Let U be 'up', D be 'down' and M be 'multiplied' (the result), and suppose for

the moment

that pixels are both NOT premultiplied and normalized to [0, 1].

The additional value v is used to add to RGB components of the U pixel:

U'_rgb = (U_rgb + v / 255)

The matte component is either U_m D_m in case (matte == true), or D_m in case

(matte == false).

Please, observe that in the case (matte == false) that choice makes the product

NOT COMMUTATIVE,

but I think that's justified - it's simply the most obvious use case.

In case (matte == true), each channel is multiplied independently, and that's

it.

The matter is more complicated when (matte == false). The problem in this case

is dealing with rgb

components when U and D both have some transparency.  When U is fully

transparent, we expect the

result to be D, and vice-versa, which is non-trivial.

We REQUIRE that (let's denote r only here):

M_r = M_r_u = U_r (1 - D_m) + U_r D_r D_m,   when  U_m == 1    (When U is fully

opaque, M_r is a D_m-linear

                                                                combination of

U_r and U_r D_r)

M_r = M_r_d = D_r (1 - U_m) + U_r D_r U_m,   when  D_m == 1    (Vice-versa, when

it's D that is fully opaque)

Finally, we're building a weighted sum, by U_m and D_m of the two above:

M_r = (M_r_u * U_m + M_r_d * D_m) / (U_m + D_m) =

    = (...) =

    = [ U_r U_m (1 - D_m) + D_r D_m (1 - U_m) + 2 U_r U_m  D_r D_m ] / (U_m +

D_m)

*/

  // 32-bit images case

  TRaster32P up32 = rup, down32 = rdown, out32 = rout;

  if (up32 && down32 && out32) {

    static const float maxChannelF = float(TPixel32::maxChannelValue);

    static const UCHAR maxChannelC = UCHAR(TPixel32::maxChannelValue);

    float vf = v;

    if (matte) {

      float dnMf, upMf_norm, outMf;

      FOR_EACH_PIXEL_32_BEGIN_LOOP  // Awful... should be explicit...

          dnMf  = downPix->m;

      upMf_norm = upPix->m / maxChannelF;

      outMf = downPix->m * upMf_norm;

      outPix->r =

          tcrop((upPix->r / upMf_norm + vf) * (downPix->r / dnMf), 0.0f, outMf);

      outPix->g =

          tcrop((upPix->g / upMf_norm + vf) * (downPix->g / dnMf), 0.0f, outMf);

      outPix->b =

          tcrop((upPix->b / upMf_norm + vf) * (downPix->b / dnMf), 0.0f, outMf);

      outPix->m = outMf;

      // NOTE:  + 0.5f in the crop arguments could take care of rounding...

      FOR_EACH_PIXEL_32_END_LOOP

    } else {

      float umf_norm, dmf_norm, umdmf_norm, outMf;

      float mSumf, uf, df, ufdf, normalizer;

      FOR_EACH_PIXEL_32_BEGIN_LOOP

      mSumf = upPix->m + float(downPix->m);

      if (mSumf > 0.0f) {

        umf_norm = upPix->m / maxChannelF, dmf_norm = downPix->m / maxChannelF;

        outMf = upPix->m +

                (1.0f - umf_norm) *

                    downPix->m;  // umf_norm should be ensured in [0.0, 1.0].

        // Convex combination should be in the conversion range.

        normalizer = outMf / (maxChannelF * mSumf);

        umdmf_norm = umf_norm * dmf_norm;

        uf = upPix->r + vf * umdmf_norm, df = downPix->r, ufdf = uf * df;

        outPix->r = tcrop((uf * (maxChannelC - downPix->m) +

                           df * (maxChannelC - upPix->m) + ufdf + ufdf) *

                              normalizer,

                          0.0f, outMf);

        uf = upPix->g + vf * umdmf_norm, df = downPix->g, ufdf = uf * df;

        outPix->g = tcrop((uf * (maxChannelC - downPix->m) +

                           df * (maxChannelC - upPix->m) + ufdf + ufdf) *

                              normalizer,

                          0.0f, outMf);

        uf = upPix->b + vf * umdmf_norm, df = downPix->b, ufdf = uf * df;

        outPix->b = tcrop((uf * (maxChannelC - downPix->m) +

                           df * (maxChannelC - upPix->m) + ufdf + ufdf) *

                              normalizer,