#include "blend.h"
 
//TPoint structure
#include "tgeometry.h"
 
//Palette - pixel functions
#include "tpalette.h"
#include "tpixelutils.h"
 
#include <vector>
 
//=================================================================================
 
//===========================
//    Blur pattern class
//---------------------------
 
//! The BlurPattern class delineates the idea of a 'blur'
//! pattern from a number of random sample points taken
//! in a neighbourhood of the blurred pixel. The pattern
//! develops in a radial manner if specified, so that possible
//! 'obstacles' in the blur can be identified.
 
class BlurPattern
{
public:
	typedef std::vector<TPoint> SamplePath;
 
	std::vector<TPoint> m_samples;
	std::vector<SamplePath> m_samplePaths;
 
	BlurPattern(double distance, unsigned int samplesCount, bool radial);
	~BlurPattern() {}
};
 
//---------------------------------------------------------------------------------
 
//Builds the specified number of samples count, inside the specified distance
//from the origin. If the pattern is radial, paths to the samples points are
//calculated.
BlurPattern::BlurPattern(double distance, unsigned int samplesCount, bool radial)
{
	const double randFactor = 2.0 * distance / RAND_MAX;
 
	m_samples.resize(samplesCount);
 
	//Build the samples
	unsigned int i;
	for (i = 0; i < samplesCount; ++i) {
		//NOTE: The following method ensures a perfectly flat probability distribution.
 
		TPoint candidatePoint(tround(rand() * randFactor - distance), tround(rand() * randFactor - distance));
		double distanceSq = sq(distance);
		while (sq(candidatePoint.x) + sq(candidatePoint.y) > distanceSq)
			candidatePoint = TPoint(tround(rand() * randFactor - distance), tround(rand() * randFactor - distance));
 
		m_samples[i] = candidatePoint;
	}
 
	m_samplePaths.resize(samplesCount);
 
	//If necessary, build the paths
	if (radial) {
		for (i = 0; i < samplesCount; ++i) {
			TPoint &sample = m_samples[i];
 
			int l = tmax(abs(sample.x), abs(sample.y));
 
			m_samplePaths[i].reserve(l);
 
			double dx = sample.x / (double)l;
			double dy = sample.y / (double)l;
 
			double x, y;
			int j;
			for (j = 0, x = dx, y = dy; j < l; x += dx, y += dy, ++j)
				m_samplePaths[i].push_back(TPoint(tround(x), tround(y)));
		}
	}
}
 
//=================================================================================
 
//=================================
//    Raster Selection classes
//---------------------------------
 
struct SelectionData {
	UCHAR m_selectedInk : 1;
	UCHAR m_selectedPaint : 1;
	UCHAR m_pureInk : 1;
	UCHAR m_purePaint : 1;
};
 
//=================================================================================
 
// Implements an array of selection infos using bitfields. It seems that bitfields are more optimized than
// using raw bits and bitwise operators, and use just the double of the space required with bit arrays.
class SelectionArrayPtr
{
	SelectionData *m_buffer;
 
public:
	SelectionArrayPtr() : m_buffer(0) {}
 
	inline void allocate(unsigned int count)
	{
		m_buffer = new SelectionData[count];
		memset(m_buffer, 0, count * sizeof(SelectionData));
	}
 
	inline void destroy()
	{
		delete[] m_buffer;
		m_buffer = 0;
	}
 
	inline SelectionData *data() const
	{
		return m_buffer;
	}
 
	inline SelectionData *data()
	{
		return m_buffer;
	}
};
 
//=================================================================================
 
// Bitmap used to store blend color selections and pure color informations.
class SelectionRaster
{
	SelectionArrayPtr m_selection;
 
	int m_wrap;
 
public:
	SelectionRaster(TRasterCM32P cm);
 
	void updateSelection(TRasterCM32P cm, const BlendParam &param);
 
	SelectionData *data() const
	{
		return m_selection.data();
	}
 
	SelectionData *data()
	{
		return m_selection.data();
	}
 
	void destroy()
	{
		m_selection.destroy();
	}
 
	bool isSelectedInk(int xy) const
	{
		return (m_selection.data() + xy)->m_selectedInk;
	}
 
	bool isSelectedInk(int x, int y) const
	{
		return isSelectedInk(x + y * m_wrap);
	}
 
	bool isSelectedPaint(int xy) const
	{
		return (m_selection.data() + xy)->m_selectedPaint;
	}
 
	bool isSelectedPaint(int x, int y) const
	{
		return isSelectedPaint(x + y * m_wrap);
	}
 
	bool isPureInk(int xy) const
	{
		return (m_selection.data() + xy)->m_pureInk;
	}
 
	bool isPureInk(int x, int y) const
	{
		return isPureInk(x + y * m_wrap);
	}
 
	bool isPurePaint(int xy) const
	{
		return (m_selection.data() + xy)->m_purePaint;
	}
 
	bool isPurePaint(int x, int y) const
	{
		return isPurePaint(x + y * m_wrap);
	}
 
	bool isToneColor(int xy) const
	{
		return !(isPureInk(xy) || isPurePaint(xy));
	}
 
	bool isToneColor(int x, int y) const
	{
		return isToneColor(x + y * m_wrap);
	}
};
 
//---------------------------------------------------------------------------------
 
inline UCHAR linearSearch(const int *v, unsigned int vSize, int k)
{
	const int *vEnd = v + vSize;
	for (; v < vEnd; ++v)
		if (*v == k)
			return 1;
	return 0;
}
 
//---------------------------------------------------------------------------------
 
// I've seen the std::binary_search go particularly slow... perhaps it was the debug mode,
// but I guess this is the fastest version possible.
inline UCHAR binarySearch(const int *v, unsigned int vSize, int k)
{
	//NOTE: v.size() > 0 due to external restrictions. See SelectionRaster's constructor.
 
	int a = -1, b, c = vSize;
 
	for (b = c >> 1; b != a; b = (a + c) >> 1) {
		if (v[b] == k)
			return 1;
		else if (k < v[b])
			c = b;
		else
			a = b;
	}
 
	return 0;
}
 
//---------------------------------------------------------------------------------
 
SelectionRaster::SelectionRaster(TRasterCM32P cm)
{
	unsigned int lx = cm->getLx(), ly = cm->getLy(), wrap = cm->getWrap();
	unsigned int size = lx * ly;
 
	m_wrap = lx;
 
	m_selection.allocate(size);
 
	cm->lock();
	TPixelCM32 *pix, *pixBegin = (TPixelCM32 *)cm->getRawData();
 
	SelectionData *selData = data();
 
	unsigned int i, j;
	for (i = 0; i < ly; ++i) {
		pix = pixBegin + i * wrap;
		for (j = 0; j < lx; ++j, ++pix, ++selData) {
			selData->m_pureInk = pix->getTone() == 0;
			selData->m_purePaint = pix->getTone() == 255;
		}
	}
 
	cm->unlock();
}
 
//---------------------------------------------------------------------------------
 
void SelectionRaster::updateSelection(TRasterCM32P cm, const BlendParam &param)
{
	//Make a hard copy of color indexes. We do so since we absolutely prefer
	//having them SORTED!
	std::vector<int> cIndexes = param.colorsIndexes;
	std::sort(cIndexes.begin(), cIndexes.end());
 
	unsigned int lx = cm->getLx(), ly = cm->getLy(), wrap = cm->getWrap();
 
	//Scan each cm pixel, looking if its ink or paint is in param's colorIndexes.
	cm->lock();
	TPixelCM32 *pix, *pixBegin = (TPixelCM32 *)cm->getRawData();
 
	SelectionData *selData = data();
 
	const int *v = &cIndexes[0]; //NOTE: cIndexes.size() > 0 due to external check.
	unsigned int vSize = cIndexes.size();
	unsigned int i, j;
 
	//NOTE: It seems that linear searches are definitely best for small color indexes.
	if (vSize > 50) {
		for (i = 0; i < ly; ++i) {
			pix = pixBegin + i * wrap;
			for (j = 0; j < lx; ++j, ++pix, ++selData) {
				selData->m_selectedInk = binarySearch(v, vSize, pix->getInk());
				selData->m_selectedPaint = binarySearch(v, vSize, pix->getPaint());
			}
		}
	} else {
		for (i = 0; i < ly; ++i) {
			pix = pixBegin + i * wrap;
			for (j = 0; j < lx; ++j, ++pix, ++selData) {
				selData->m_selectedInk = linearSearch(v, vSize, pix->getInk());
				selData->m_selectedPaint = linearSearch(v, vSize, pix->getPaint());
			}
		}
	}
 
	cm->unlock();
}
 
//=================================================================================
 
//========================
//    Blend functions
//------------------------
 
// Pixel whose channels are doubles. Used to store intermediate values for pixel blending.
struct DoubleRGBMPixel {
	double r;
	double g;
	double b;
	double m;
 
	DoubleRGBMPixel() : r(0.0), g(0.0), b(0.0), m(0.0) {}
};
 
//---------------------------------------------------------------------------------
 
const double maxTone = TPixelCM32::getMaxTone();
 
// Returns the ink & paint convex factors associated with passed tone.
inline void getFactors(int tone, double &inkFactor, double &paintFactor)
{
	paintFactor = tone / maxTone;
	inkFactor = (1.0 - paintFactor);
}
 
//---------------------------------------------------------------------------------
 
// Copies the cmIn paint and ink colors to the output rasters.
void buildLayers(
	const TRasterCM32P &cmIn, const std::vector<TPixel32> &palColors,
	TRaster32P &inkRaster, TRaster32P &paintRaster)
{
	//Separate cmIn by copying the ink & paint colors directly to the layer rasters.
	TPixelCM32 *cmPix, *cmBegin = (TPixelCM32 *)cmIn->getRawData();
	TPixel32 *inkPix = (TPixel32 *)inkRaster->getRawData();
	TPixel32 *paintPix = (TPixel32 *)paintRaster->getRawData();
 
	unsigned int i, j, lx = cmIn->getLx(), ly = cmIn->getLy(), wrap = cmIn->getWrap();
	for (i = 0; i < ly; ++i) {
		cmPix = cmBegin + i * wrap;
 
		for (j = 0; j < lx; ++j, ++cmPix, ++inkPix, ++paintPix) {
			*inkPix = palColors[cmPix->getInk()];
			*paintPix = palColors[cmPix->getPaint()];
 
			//Should pure colors be checked...?
		}
	}
}
 
//---------------------------------------------------------------------------------
 
// Returns true or false whether the selectedColor is the only selectable color
// in the neighbourhood. If so, the blend copies it to the output layer pixel directly.
inline bool isFlatNeighbourhood(
	int selectedColor,
	const TRasterCM32P &cmIn, const TPoint &pos,
	const SelectionRaster &selRas,
	const BlurPattern &blurPattern)
{
	TPixelCM32 &pix = cmIn->pixels(pos.y)[pos.x];
	int lx = cmIn->getLx(), ly = cmIn->getLy();
	unsigned int xy;
 
	TPoint samplePix;
 
	const TPoint *samplePoint = blurPattern.m_samples.empty() ? 0 : &blurPattern.m_samples[0];
 
	//Read the samples to determine if they only have posSelectedColor
	unsigned int i, samplesCount = blurPattern.m_samples.size();
	for (i = 0; i < samplesCount; ++i, ++samplePoint) {
		//Make sure the sample is inside the image
		samplePix.x = pos.x + samplePoint->x;
		samplePix.y = pos.y + samplePoint->y;
 
		xy = samplePix.x + lx * samplePix.y;
 
		if (samplePix.x < 0 || samplePix.y < 0 || samplePix.x >= lx || samplePix.y >= ly)
			continue;
 
		if (!selRas.isPurePaint(xy) && selRas.isSelectedInk(xy))
			if (cmIn->pixels(samplePix.y)[samplePix.x].getInk() != selectedColor)
				return false;
 
		if (!selRas.isPureInk(xy) && selRas.isSelectedPaint(xy))
			if (cmIn->pixels(samplePix.y)[samplePix.x].getPaint() != selectedColor)
				return false;
	}
 
	return true;
}
 
//---------------------------------------------------------------------------------
 
// Calculates the estimate of blend selection in the neighbourhood specified by
// blurPattern.
inline void addSamples(
	const TRasterCM32P &cmIn, const TPoint &pos,
	const TRaster32P &inkRas, const TRaster32P &paintRas,
	const SelectionRaster &selRas,
	const BlurPattern &blurPattern,
	DoubleRGBMPixel &pixSum, double &factorsSum)
{
	double inkFactor, paintFactor;
	unsigned int xy, j, l;
	int lx = cmIn->getLx(), ly = cmIn->getLy();
	TPixel32 *color;
	TPoint samplePos, pathPos;
 
	const TPoint *samplePoint = blurPattern.m_samples.empty() ? 0 : &blurPattern.m_samples[0];
	const TPoint *pathPoint;
 
	unsigned int i, blurSamplesCount = blurPattern.m_samples.size();
	for (i = 0; i < blurSamplesCount; ++i, ++samplePoint) {
		//Add each samples contribute to the sum
		samplePos.x = pos.x + samplePoint->x;
		samplePos.y = pos.y + samplePoint->y;
		if (samplePos.x < 0 || samplePos.y < 0 || samplePos.x >= lx || samplePos.y >= ly)
			continue;
 
		//Ensure that each pixel on the sample's path (if any) is selected
		l = blurPattern.m_samplePaths[i].size();
		pathPoint = blurPattern.m_samplePaths[i].empty() ? 0 : &blurPattern.m_samplePaths[i][0];
		for (j = 0; j < l; ++j, ++pathPoint) {
			pathPos.x = pos.x + pathPoint->x;
			pathPos.y = pos.y + pathPoint->y;
			xy = pathPos.x + lx * pathPos.y;
 
			if (!(selRas.isPurePaint(xy) || selRas.isSelectedInk(xy)))
				break;
 
			if (!(selRas.isPureInk(xy) || selRas.isSelectedPaint(xy)))
				break;
		}
 
		if (j < l)
			continue;
 
		xy = samplePos.x + lx * samplePos.y;
 
		if (selRas.isSelectedInk(xy) && !selRas.isPurePaint(xy)) {
			getFactors(cmIn->pixels(samplePos.y)[samplePos.x].getTone(), inkFactor, paintFactor);
 
			color = &inkRas->pixels(samplePos.y)[samplePos.x];
			pixSum.r += inkFactor * color->r;
			pixSum.g += inkFactor * color->g;
			pixSum.b += inkFactor * color->b;
			pixSum.m += inkFactor * color->m;
			factorsSum += inkFactor;
		}
 
		if (selRas.isSelectedPaint(xy) && !selRas.isPureInk(xy)) {
			getFactors(cmIn->pixels(samplePos.y)[samplePos.x].getTone(), inkFactor, paintFactor);
 
			color = &paintRas->pixels(samplePos.y)[samplePos.x];
			pixSum.r += paintFactor * color->r;
			pixSum.g += paintFactor * color->g;
			pixSum.b += paintFactor * color->b;
			pixSum.m += paintFactor * color->m;
			factorsSum += paintFactor;
		}
	}
}
 
//---------------------------------------------------------------------------------
 
typedef std::pair<TRaster32P, TRaster32P> RGBMRasterPair;
 
//---------------------------------------------------------------------------------
 
// Performs a single color blending. This function can be repeatedly invoked to
// perform multiple color blending.
inline void doBlend(
	const TRasterCM32P &cmIn,
	RGBMRasterPair &inkLayer, RGBMRasterPair &paintLayer,
	const SelectionRaster &selRas,
	const std::vector<BlurPattern> &blurPatterns)
{
	//Declare some vars
	unsigned int blurPatternsCount = blurPatterns.size();
	int lx = cmIn->getLx(), ly = cmIn->getLy();
	double totalFactor;
 
	TPixelCM32 *cmPix, *cmBegin = (TPixelCM32 *)cmIn->getRawData();
 
	TPixel32
		*inkIn = (TPixel32 *)inkLayer.first->getRawData(),
		*inkOut = (TPixel32 *)inkLayer.second->getRawData(),
		*paintIn = (TPixel32 *)paintLayer.first->getRawData(),
		*paintOut = (TPixel32 *)paintLayer.second->getRawData();
 
	const BlurPattern *blurPattern, *blurPatternsBegin = &blurPatterns[0];
	bool builtSamples = false;
 
	DoubleRGBMPixel samplesSum;
 
	//For every cmIn pixel
	TPoint pos;
	SelectionData *selData = selRas.data();
	cmPix = cmBegin;
	for (pos.y = 0; pos.y < ly; ++pos.y, cmPix = cmBegin + pos.y * cmIn->getWrap())
		for (pos.x = 0; pos.x < lx; ++pos.x, ++inkIn, ++inkOut, ++paintIn, ++paintOut, ++selData, ++cmPix) {
			blurPattern = blurPatternsBegin + (rand() % blurPatternsCount);
 
			//Build the ink blend color
			if (!selData->m_purePaint && selData->m_selectedInk) {
				if (!builtSamples) {
					//Build samples contributes
					totalFactor = 1.0;
					samplesSum.r = samplesSum.g = samplesSum.b = samplesSum.m = 0.0;
 
					if (!isFlatNeighbourhood(cmPix->getInk(), cmIn, pos, selRas, *blurPattern))
						addSamples(cmIn, pos, inkLayer.first, paintLayer.first, selRas, *blurPattern,
								   samplesSum, totalFactor);
 
					builtSamples = true;
				}
 
				//Output the blended pixel
				inkOut->r = (samplesSum.r + inkIn->r) / totalFactor;
				inkOut->g = (samplesSum.g + inkIn->g) / totalFactor;
				inkOut->b = (samplesSum.b + inkIn->b) / totalFactor;
				inkOut->m = (samplesSum.m + inkIn->m) / totalFactor;
			} else {
				//If the color is not blended, then just copy the old layer pixel
				*inkOut = *inkIn;
			}
 
			//Build the paint blend color
			if (!selData->m_pureInk && selData->m_selectedPaint) {
				if (!builtSamples) {
					//Build samples contributes
					totalFactor = 1.0;
					samplesSum.r = samplesSum.g = samplesSum.b = samplesSum.m = 0.0;
 
					if (!isFlatNeighbourhood(cmPix->getPaint(), cmIn, pos, selRas, *blurPattern))
						addSamples(cmIn, pos, inkLayer.first, paintLayer.first, selRas, *blurPattern,
								   samplesSum, totalFactor);
 
					builtSamples = true;
				}
 
				//Output the blended pixel
				paintOut->r = (samplesSum.r + paintIn->r) / totalFactor;
				paintOut->g = (samplesSum.g + paintIn->g) / totalFactor;
				paintOut->b = (samplesSum.b + paintIn->b) / totalFactor;
				paintOut->m = (samplesSum.m + paintIn->m) / totalFactor;
			} else {
				//If the color is not blended, then just copy the old layer pixel
				*paintOut = *paintIn;
			}
 
			builtSamples = false;
		}
}
 
//---------------------------------------------------------------------------------
 
typedef std::vector<BlurPattern> BlurPatternContainer;
 
//---------------------------------------------------------------------------------
 
/*! This function performs a group of <a> spatial color blending <\a> operations on Toonz Images.
    The BlendParam structure stores the blend options recognized by this function; it includes
    a list of the palette indexes involved in the blend operation, plus:
    \li \b Intensity represents the \a radius of the blur operation between blend colors.
    \li \b Smoothness is the number of samples per pixel used to approximate the blur.
    <li> <b> Stop at Contour <\b> specifies if lines from pixels to neighbouring samples
         should not trespass color indexes not included in the blend operation <\li>
    The succession of input blend parameters are applied in the order.
*/
 
template <typename PIXEL>
void blend(TToonzImageP ti, TRasterPT<PIXEL> rasOut, const std::vector<BlendParam> &params)
{
	assert(ti->getRaster()->getSize() == rasOut->getSize());
 
	//Extract the interesting raster. It should be the savebox of passed cmap, plus - if
	//some param has the 0 index as blending color - the intensity of that blend param.
	unsigned int i, j;
	TRect saveBox(ti->getSavebox());
 
	int enlargement = 0;
	for (i = 0; i < params.size(); ++i)
		for (j = 0; j < params[i].colorsIndexes.size(); ++j)
			if (params[i].colorsIndexes[j] == 0)
				enlargement = tmax(enlargement, tceil(params[i].intensity));
	saveBox = saveBox.enlarge(enlargement);
 
	TRasterCM32P cmIn(ti->getRaster()->extract(saveBox));
	TRasterPT<PIXEL> rasOutExtract = rasOut->extract(saveBox);
 
	//Ensure that cmIn and rasOut have the same size
	unsigned int lx = cmIn->getLx(), ly = cmIn->getLy();
 
	//Build the pure colors infos
	SelectionRaster selectionRaster(cmIn);
 
	//Now, build a little group of BlurPatterns - and for each, one for passed param.
	//A small number of patterns per param is needed to make the pattern look not ever the same.
	const int blurPatternsPerParam = 10;
	std::vector<BlurPatternContainer> blurGroup(params.size());
 
	for (i = 0; i < params.size(); ++i) {
		BlurPatternContainer &blurContainer = blurGroup[i];
		blurContainer.reserve(blurPatternsPerParam);
 
		for (j = 0; j < blurPatternsPerParam; ++j)
			blurContainer.push_back(BlurPattern(params[i].intensity, params[i].smoothness, params[i].stopAtCountour));
	}
 
	//Build the palette
	TPalette *palette = ti->getPalette();
	std::vector<TPixel32> paletteColors;
	paletteColors.resize(palette->getStyleCount());
	for (i = 0; i < paletteColors.size(); ++i)
		paletteColors[i] = premultiply(palette->getStyle(i)->getAverageColor());
 
	//Build the 4 auxiliary rasters for the blending procedure: they are ink / paint versus input / output in the blend.
	//The output raster is reused to spare some memory - it should be, say, the inkLayer's second at the end of the overall
	//blending procedure. It could be the first, without the necessity of clearing it before blending the layers, but things
	//get more complicated when PIXEL is TPixel64...
	RGBMRasterPair inkLayer, paintLayer;
 
	TRaster32P rasOut32P_1(lx, ly, lx, (TPixel32 *)rasOut->getRawData(), false);
	inkLayer.first = (params.size() % 2) ? rasOut32P_1 : TRaster32P(lx, ly);
	inkLayer.second = (params.size() % 2) ? TRaster32P(lx, ly) : rasOut32P_1;
 
	if (PIXEL::maxChannelValue >= TPixel64::maxChannelValue) {
		TRaster32P rasOut32P_2(lx, ly, lx, ((TPixel32 *)rasOut->getRawData()) + lx * ly, false);
		paintLayer.first = (params.size() % 2) ? rasOut32P_2 : TRaster32P(lx, ly);
		paintLayer.second = (params.size() % 2) ? TRaster32P(lx, ly) : rasOut32P_2;
	} else {
		paintLayer.first = TRaster32P(lx, ly);
		paintLayer.second = TRaster32P(lx, ly);
	}
 
	inkLayer.first->clear();
	inkLayer.second->clear();
	paintLayer.first->clear();
	paintLayer.second->clear();
 
	//Now, we have to perform the blur of each of the cm's pixels.
	cmIn->lock();
	rasOut->lock();
 
	inkLayer.first->lock();
	inkLayer.second->lock();
	paintLayer.first->lock();
	paintLayer.second->lock();
 
	//Convert the initial cmIn to fullcolor ink - paint layers
	buildLayers(cmIn, paletteColors, inkLayer.first, paintLayer.first);
 
	//Perform the blend on separated ink - paint layers
	for (i = 0; i < params.size(); ++i) {
		if (params[i].colorsIndexes.size() == 0)
			continue;
 
		selectionRaster.updateSelection(cmIn, params[i]);
		doBlend(cmIn, inkLayer, paintLayer, selectionRaster, blurGroup[i]);
 
		tswap(inkLayer.first, inkLayer.second);
		tswap(paintLayer.first, paintLayer.second);
	}
 
	//Release the unnecessary rasters
	inkLayer.second->unlock();
	paintLayer.second->unlock();
	inkLayer.second = TRaster32P();
	paintLayer.second = TRaster32P();
 
	//Clear rasOut - since it was reused to spare space...
	rasOut->clear();
 
	//Add the ink & paint layers on the output raster
	double PIXELmaxChannelValue = PIXEL::maxChannelValue;
	double toPIXELFactor = PIXELmaxChannelValue / (double)TPixel32::maxChannelValue;
	double inkFactor, paintFactor;
	TPoint pos;
 
	PIXEL *outPix, *outBegin = (PIXEL *)rasOutExtract->getRawData();
	TPixelCM32 *cmPix, *cmBegin = (TPixelCM32 *)cmIn->getRawData();
	int wrap = rasOutExtract->getWrap();
 
	TPixel32 *inkPix = (TPixel32 *)inkLayer.first->getRawData();
	TPixel32 *paintPix = (TPixel32 *)paintLayer.first->getRawData();
 
	for (i = 0; i < ly; ++i) {
		outPix = outBegin + wrap * i;
		cmPix = cmBegin + wrap * i;
		for (j = 0; j < lx; ++j, ++outPix, ++cmPix, ++inkPix, ++paintPix) {
			getFactors(cmPix->getTone(), inkFactor, paintFactor);
 
			outPix->r = tcrop(toPIXELFactor * (inkFactor * inkPix->r + paintFactor * paintPix->r), 0.0, PIXELmaxChannelValue);
			outPix->g = tcrop(toPIXELFactor * (inkFactor * inkPix->g + paintFactor * paintPix->g), 0.0, PIXELmaxChannelValue);
			outPix->b = tcrop(toPIXELFactor * (inkFactor * inkPix->b + paintFactor * paintPix->b), 0.0, PIXELmaxChannelValue);
			outPix->m = tcrop(toPIXELFactor * (inkFactor * inkPix->m + paintFactor * paintPix->m), 0.0, PIXELmaxChannelValue);
		}
	}
 
	inkLayer.first->unlock();
	paintLayer.first->unlock();
 
	cmIn->unlock();
	rasOut->unlock();
 
	//Destroy the auxiliary bitmaps
	selectionRaster.destroy();
}
 
template void blend<TPixel32>(TToonzImageP cmIn, TRasterPT<TPixel32> rasOut, const std::vector<BlendParam> &params);
template void blend<TPixel64>(TToonzImageP cmIn, TRasterPT<TPixel64> rasOut, const std::vector<BlendParam> &params);