#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) {

		FOR_EACH_PIXEL_BEGIN_LOOP(TPixelRGBM32, source32, TPixelRGBM32, matte32, TPixelRGBM32, out32)

		if (downPix->m == 255)

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

		else if (downPix->m == 0)

			*outPix = *upPix;

		else {

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

			UINT fac = MAGICFAC * (255 - 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 == 65535)

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

		else if (downPix->m == 0)

			*outPix = *upPix;

		else {

			double fac = (65535 - 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);