#include "trop.h"

#include "tpixel.h"

#include "tpixelutils.h"

#ifdef _WIN32

#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_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

  } 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 = 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

    } else {

      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

      } else

        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

  } 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 = 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

    } else {

      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

      } else

        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 {

      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

      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

  } else {

    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

    } 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::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

    } else {

      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

      } else {

        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

        } else

          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

    } else {

      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

      } else {

        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

        } else

          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,

                          0.0f, outMf);

        outPix->m = outMf;

      } else

        *outPix = TPixel32::Transparent;

      FOR_EACH_PIXEL_32_END_LOOP

    }

    return;

  }

  // 64-bit images case

  TRaster64P up64 = rup, down64 = rdown, out64 = rout;

  if (up64 && down64 && out64) {

    static const double maxChannelF = double(TPixel64::maxChannelValue);

    static const USHORT maxChannelC = USHORT(TPixel64::maxChannelValue);

    double vf =

        v * (TPixel64::maxChannelValue / double(TPixel32::maxChannelValue));

    if (matte) {

      double dnMf, upMf_norm, outMf;

      FOR_EACH_PIXEL_64_BEGIN_LOOP

      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.0, outMf);

      outPix->g =

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

      outPix->b =

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

      outPix->m = outMf;

      FOR_EACH_PIXEL_64_END_LOOP

    } else {

      double umf_norm, dmf_norm, umdmf_norm, outMf;

      double mSumf, uf, df, ufdf, normalizer;

      FOR_EACH_PIXEL_64_BEGIN_LOOP

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

      if (mSumf > 0.0) {

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

        outMf = upPix->m + (1.0 - umf_norm) * downPix->m;

        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.0, 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.0, 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,

                          0.0, outMf);

        outPix->m = outMf;

      } else

        *outPix = TPixel64::Transparent;

      FOR_EACH_PIXEL_64_END_LOOP

    }

    return;

  }

  // According to the specifics, throw an exception. I think it's not

  // appropriate, though.

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

}

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

void TRop::ropin(const TRasterP &source, const TRasterP &matte,

                 const TRasterP &rout) {

  TRaster32P source32 = source;

  TRaster32P matte32  = matte;

  TRaster32P out32    = rout;

  TRaster64P source64 = source;

  TRaster64P matte64  = matte;

  TRaster64P out64    = rout;

  if (source32 && matte32 && out32) {

    FOR_EACH_PIXEL_BEGIN_LOOP(TPixelRGBM32, source32, TPixelRGBM32, matte32,

                              TPixelRGBM32, out32)

    if (downPix->m == 0)

      outPix->r = outPix->g = outPix->b = outPix->m = 0;

    else if (downPix->m == 255)

      *outPix = *upPix;

    else {

      /*

__m128i zeros = _mm_setzero_si128();

__m128i upPix_packed_i= _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(DWORD*)upPix),

zeros);

__m128  upPix_packed  = _mm_cvtepi32_ps(_mm_unpacklo_epi16(upPix_packed_i,

zeros));

float fac = downPix->m / 255.0;

__m128 fac_packed = _mm_load1_ps(&fac);

upPix_packed = _mm_mul_ps(upPix_packed, fac_packed);

__m128i outPix_packed_i = _mm_cvtps_epi32(upPix_packed);

outPix_packed_i = _mm_packs_epi32(outPix_packed_i, zeros);

outPix_packed_i = _mm_packus_epi16(outPix_packed_i, zeros);

*(DWORD*)(outPix) = _mm_cvtsi128_si32(outPix_packed_i);

*/

      const int MAGICFAC = (257U * 256U + 1U);

      UINT fac           = MAGICFAC * downPix->m;

      outPix->r = (UINT)(upPix->r * fac + (1U << 23)) >> 24;

      outPix->g = (UINT)(upPix->g * fac + (1U << 23)) >> 24;

      outPix->b = (UINT)(upPix->b * fac + (1U << 23)) >> 24;

      outPix->m = (UINT)(upPix->m * fac + (1U << 23)) >> 24;

    }

    FOR_EACH_PIXEL_END_LOOP(source32, matte32, out32)

  } else if (source64 && matte64 && out64)

  {

    FOR_EACH_PIXEL_BEGIN_LOOP(TPixelRGBM64, source64, TPixelRGBM64, matte64,

                              TPixelRGBM64, out64)

    if (downPix->m == 0)

      outPix->r = outPix->g = outPix->b = outPix->m = 0;

    else if (downPix->m == 65535)

      *outPix = *upPix;

    else {

      /*

__m128i zeros = _mm_setzero_si128();

__m128i upPix_packed_i= _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(DWORD*)upPix),

zeros);

__m128  upPix_packed  = _mm_cvtepi32_ps(_mm_unpacklo_epi16(upPix_packed_i,

zeros));

float fac = downPix->m / 255.0;

__m128 fac_packed = _mm_load1_ps(&fac);

upPix_packed = _mm_mul_ps(upPix_packed, fac_packed);

__m128i outPix_packed_i = _mm_cvtps_epi32(upPix_packed);

outPix_packed_i = _mm_packs_epi32(outPix_packed_i, zeros);

outPix_packed_i = _mm_packus_epi16(outPix_packed_i, zeros);

*(DWORD*)(outPix) = _mm_cvtsi128_si32(outPix_packed_i);

*/

      double fac = downPix->m / 65535.0;

      outPix->r = (USHORT)(upPix->r * fac);

      outPix->g = (USHORT)(upPix->g * fac);

      outPix->b = (USHORT)(upPix->b * fac);

      outPix->m = (USHORT)(upPix->m * fac);

    }

    FOR_EACH_PIXEL_END_LOOP(source64, matte64, out64)

  } else

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

}

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

void TRop::ropout(const TRasterP &source, const TRasterP &matte,

                  const TRasterP &rout) {

  TRaster32P source32 = source;

  TRaster32P matte32  = matte;

  TRaster32P out32    = rout;

  TRaster64P source64 = source;

  TRaster64P matte64  = matte;

  TRaster64P out64    = rout;

  if (source32 && matte32 && out32) {