/////////////////////////////////////////////////////////////////////////////
//
// Image class containning basic image analysis functions
//
// Copyright © 1999
//
// DTU Image Viever and Analyser (DIVA)
// Department of Mathematical Modelling
// Technical University of Denmark (DTU), Building 321
// DK-2800 Lyngby, Denmark
// http://www.imm.dtu.dk/~diva
//
// author: Rune Fisker 
//
// Disclaimer:
//
// No guarantees of performance accompany this software,
// nor is any responsibility assumed on the part of the author(s).
// The software has been tested extensively and every effort has been
// made to insure its reliability.
//
// This software is provided by IMM and the contributor(s) ``as is'' and
// any express or implied warranties, including, but not limited to, the
// implied warranties of merchantability and fitness for a particular purpose
// are disclaimed.  In no event shall IMM or the contributor(s) be liable
// for any direct, indirect, incidental, special, exemplary, or consequential
// damages (including, but not limited to, procurement of substitute goods
// or services; loss of use, data, or profits; or business interruption)
// however caused and on any theory of liability, whether in contract, strict
// liability, or tort (including negligence or otherwise) arising in any way
// out of the use of this software, even if advised of the possibility of
// such damage.
//
// This software is partly based on the Microsoft Vision Software Developers Kit VisSDK
//
/////////////////////////////////////////////////////////////////////////////

//#include "stdafx.h"

/////////////////////////////////////////////////////////////////////////////
// Performs a simple threshold, i.e. all pixel values with a value above or equal to the
// threshold is set to 255 and the rest is set to zero
//
// input:
//	tValue: threshold value
//
// author:	Rune Fisker 8/9-1998
template <class TPixel>
void CDImageBasic<TPixel>::Threshold(TPixel tValue) 
{
	assert(NBands() == 1);

	for (int r = Top();r < Bottom(); r++)
	{
		TPixel *ptpixel = RowPointer(r);
		for (int c = Left();c < Right();c++)
		{
			if (ptpixel[c] >= tValue)
				ptpixel[c] = (TPixel) 1;
			else
				ptpixel[c] = 0;
		}
	}

	// update history
	CString sHistory;
	sHistory.Format("DIVA: Threshold by t = %f (Threshold)",(double)tValue);
	AddToHistory(sHistory);
}


/////////////////////////////////////////////////////////////////////////////
// Performs simple inverted threshold, i.e. corresponds to a simple threshold 
// where the result is inverted.
//
// input:
//	tValue: threshold value
//
// author:	Rune Fisker 26/1-1999
template <class TPixel>
void CDImageBasic<TPixel>::ThresholdInv(TPixel tValue) 
{
	assert(NBands() == 1);

	for (int r = Top();r < Bottom(); r++)
	{
		TPixel *ptpixel = RowPointer(r);
		for (int c = Left();c < Right();c++)
		{
			if (ptpixel[c] < tValue)
				ptpixel[c] = (TPixel) 1;
			else
				ptpixel[c] = 0;
		}
	}

	// update history
	CString sHistory;
	sHistory.Format("DIVA: Inverse threshold by t = %f (ThresholdInv)",(double)tValue);
	AddToHistory(sHistory);

}

/////////////////////////////////////////////////////////////////////////////
// Performs soft threshold, i.e. all pixels which have a value, which is below
// the threshold is set to the threshold value
//
// input:
//	tValue: threshold value
//
// author:	Rune Fisker 8/4-1999
template <class TPixel>
void CDImageBasic<TPixel>::ThresholdSoft(TPixel tValue) 
{
	assert(NBands() == 1);

	for (int r = Top();r < Bottom(); r++)
	{
		TPixel *ptpixel = RowPointer(r);
		for (int c = Left();c < Right();c++)
		{
			if (ptpixel[c] < tValue)
				ptpixel[c] = tValue;
		}
	}

	// update history
	CString sHistory;
	sHistory.Format("DIVA: Soft threshold by t = %f (ThresholdSoft)",(double)tValue);
	AddToHistory(sHistory);

}

/////////////////////////////////////////////////////////////////////////////
// Performs inverse soft threshold, i.e. all pixels which have a value, which is above 
// to the threshold is set to the threshold value
//
// input:
//	tValue: threshold value
//
// author:	Rune Fisker 26/1-1999
template <class TPixel>
void CDImageBasic<TPixel>::ThresholdSoftInv(TPixel tValue) 
{
	assert(NBands() == 1);

	for (int r = Top();r < Bottom(); r++)
	{
		TPixel *ptpixel = RowPointer(r);
		for (int c = Left();c < Right();c++)
		{
			if (ptpixel[c] > tValue)
				ptpixel[c] = tValue;
		}
	}

	// update history
	CString sHistory;
	sHistory.Format("DIVA: Inverse soft threshold by t = %f (ThresholdSoftInv)",(double)tValue);
	AddToHistory(sHistory);

}


/////////////////////////////////////////////////////////////////////////////
// takes square of each pixel. Supports only gray images
//
// author:	Rune Fisker 8/4-1999
template <class TPixel>
void CDImageBasic<TPixel>::Sqr() 
{
	assert(NBands() == 1);

	TPixel *ppixel;
	for (int y = Top();y < Bottom(); y++)
	{
		ppixel = &(RowPointer(y)[Left()]);
		for (int x = Left();x < Right();x++)
		{
			*ppixel++ *= *ppixel;
		}
	}

	// update history
	AddToHistory("DIVA: Sqaure of each pixel (Sqr)");
}

/////////////////////////////////////////////////////////////////////////////
//	takes the square root of each pixel. Supports only gray images
//
//	author:	Rune Fisker, 29/4-1998
template <class TPixel>
void CDImageBasic<TPixel>::Sqrt()
{
	TPixel t = 0;
	TPixel tMinRange = VisRangeMin(t,false);

	TPixel *ppixel;
	for (int y = Top();y < Bottom(); y++)
	{
		ppixel = &(RowPointer(y)[Left()]);
		for (int x = Left();x < Right(); x++)
		{
			*ppixel++ = sqrt(*ppixel);
		}
	}

	// update history
	AddToHistory("DIVA: Sqaure root of each pixel (Sqrt)");

}

/////////////////////////////////////////////////////////////////////////////
//	takes log of each pixel. Supports only gray images
//
//	author:	Rune Fisker, 29/4-1998
template <class TPixel>
void CDImageBasic<TPixel>::Log()
{
	TPixel *ppixel;
	for (int y = Top();y < Bottom(); y++)
	{
		ppixel = &(RowPointer(y)[Left()]);
		for (int x = Left();x < Right(); x++)
		{
			*ppixel++ = log(*ppixel);
		}
	}

	// update history
	AddToHistory("DIVA: Log of each pixel (Log)");
}



/////////////////////////////////////////////////////////////////////////////
//	Vertical flip 
//
//	author:	Rune Fisker, 29/9-1998
template <class TPixel>
void CDImageBasic<TPixel>::FlipVertical() 
{
	assert(NBands() == 1);

	int nRowSize = Width()*sizeof(TPixel);
	
	BYTE *ptDown,*ptUp,*ptTmp;
	ptTmp = new BYTE [nRowSize];
	
	// swap the rows
	for (int r = Top();r < Top()+Height()/2; r++)
	{
		ptDown = (BYTE *) RowPointer(r) + Left()*sizeof(TPixel);
		ptUp = (BYTE *) RowPointer(Bottom()-r+Top()-1) + Left()*sizeof(TPixel);
		memcpy(ptTmp,ptDown,nRowSize);
		memcpy(ptDown,ptUp,nRowSize);
		memcpy(ptUp,ptTmp,nRowSize);
	}

	delete ptTmp;

	// update history
	AddToHistory("DIVA: Vertical flip (FlipVertical)");

}

/////////////////////////////////////////////////////////////////////////////
//	Horizontal flip 
//
//	author:	Rune Fisker, 29/9-1998
template <class TPixel>
void CDImageBasic<TPixel>::FlipHorizontal() 
{
	TPixel tpixelTmp;
	int nColLeft;

	for (int r = Top();r < Bottom(); r++)
	{
		TPixel *ptpixel = RowPointer(r);
		for (int c = Left();c < Left()+Width()/2; c++)
		{
			nColLeft = Right()-c+Left()-1;
			
			tpixelTmp = ptpixel[c];
			ptpixel[c] = ptpixel[nColLeft];
			ptpixel[nColLeft] = tpixelTmp;
		}
	}

	// update history
	AddToHistory("DIVA: Horizontal flip (FlipHorizontal)");
}

/////////////////////////////////////////////////////////////////////////////
//	Horizontal flip 
//
//	author:	Rune Fisker, 29/9-1998
template <class TPixel>
LPVOID CDImageBasic<TPixel>::Image2DIB(int &nSizeDIB)
{
	if ((PixFmt() == evispixfmtGrayByte) ||
		(PixFmt() == evispixfmtRGBAByte))
	{
		// create bitmapinfo
		struct
		{
			BITMAPINFOHEADER header;
			RGBQUAD mpbul[256];
		}			bmi;

		// file bitmapinforheader
		FillBitmapInfoHeader(&bmi.header,0);
		
		// width have to be corrected
		bmi.header.biWidth = Width();

		// May need to fill in color map values.
		if (bmi.header.biBitCount <= 8)
		{			
			if ((UseColorMap()) && (ColorMap().CbData() >= 1024))
			{
				memcpy(&(bmi.mpbul),ColorMap().PbData(),sizeof(bmi.mpbul));
				bmi.header.biClrUsed = ColorMap().CbData()/sizeof(RGBQUAD);
			}
			else if (bmi.header.biBitCount > 4) 
			{
				memcpy(&(bmi.mpbul),s_mpbulByteColorMapGrayBasic,sizeof(bmi.mpbul));
				bmi.header.biClrUsed = 256;	
			}
			else if (UseColorMap())
			{
				memcpy(&(bmi.mpbul),s_mpbulNibbleColorMapDefaultBasic,sizeof(bmi.mpbul));
				bmi.header.biClrUsed = 16;
			}
			else
			{
				memcpy(&(bmi.mpbul),s_mpbulNibbleColorMapGrayBasic,sizeof(bmi.mpbul));
				bmi.header.biClrUsed = 16;
			}
		}

		// unfortunately a large number of programs can handle negative height,
		bmi.header.biHeight = -bmi.header.biHeight;
		// so the image need to be vertical flipped
		FlipVertical();

		// set WORD alignment
		int nWordAlign;
		if ((bmi.header.biBitCount/8*bmi.header.biWidth) % sizeof(DWORD))
			nWordAlign = sizeof(DWORD)-(bmi.header.biBitCount/8*bmi.header.biWidth) % sizeof(DWORD);
		else
			nWordAlign = 0;

		// calc. sizes
		int nSizeImage = bmi.header.biBitCount/8 * bmi.header.biHeight * (bmi.header.biWidth+nWordAlign);
		int nSizeBMI = bmi.header.biSize + bmi.header.biClrUsed*sizeof(RGBQUAD);
		nSizeDIB = nSizeBMI + nSizeImage;

		LPVOID *pDIB;		
		pDIB = new LPVOID [nSizeDIB];

		if (pDIB!=NULL) {
			// copy bmi
			memcpy(pDIB, &bmi.header, nSizeBMI);				

			BYTE *pbSrc,*pbDest;
			int nSizeRow = bmi.header.biBitCount/8 * bmi.header.biWidth;
			// copy ROI image
			for (int r = Top();r < Bottom(); r++)
			{
				pbSrc = (BYTE*) RowPointer(r) + Left()*bmi.header.biBitCount/8;
				pbDest = (BYTE*)pDIB + (nSizeRow+nWordAlign)*(r - Top()) + nSizeBMI;
				
				memcpy(pbDest,pbSrc,nSizeRow);
			}
		}
		else
			return NULL;

		// flip back
		FlipVertical();

		return pDIB;
	}
	return NULL;
}

template <class TPixel>
bool CDImageBasic<TPixel>::DIB2Image(LPVOID pDIB)
{
	if (pDIB == NULL)
		return false;

	bool bUseBitmapCore;

	BITMAPINFOHEADER *bmih = (BITMAPINFOHEADER *)pDIB;
	BITMAPCOREHEADER *bmch = (BITMAPCOREHEADER *)pDIB;

	if (bmih->biSize == sizeof(BITMAPINFOHEADER)) 
	{
		bUseBitmapCore=false;

		if (bmih->biCompression != BI_RGB) 
			return false;
	}
	else if (bmch->bcSize == sizeof(BITMAPCOREHEADER)) 
	{
		bUseBitmapCore=true;
	}
				
	switch (bmih->biBitCount) {
	case 8 :
		if (PixFmt() != evispixfmtGrayByte)
			return false;
		break;
	case 24 :
	case 32 :
		if (PixFmt() != evispixfmtRGBAByte)
			return false;
		break;
	default :
		return false;
	}

	if (bmih->biPlanes != 1) 
		return false;
	
	Deallocate();

	if (bmih->biHeight > 0)
		Allocate(bmih->biWidth,bmih->biHeight);
	else
		Allocate(bmih->biWidth,-bmih->biHeight);

	// set WORD alignment
	int nWordAlign;
	if ((bmih->biBitCount/8*bmih->biWidth) % sizeof(DWORD))
		nWordAlign = sizeof(DWORD)-(bmih->biBitCount/8*bmih->biWidth) % sizeof(DWORD);
	else
		nWordAlign = 0;

	int nSizeBMI = bmih->biSize + bmih->biClrUsed*sizeof(RGBQUAD);

	BYTE *pbSrc,*pbDest;
	int nSizeRow = bmih->biBitCount/8 * Width();
	
	// copy image
	for (int r = Top();r < Bottom(); r++)
	{
		pbDest = (BYTE*) RowPointer(r) + Left()*bmih->biBitCount/8;
		pbSrc = (BYTE*)pDIB + (nSizeRow+nWordAlign)*(r - Top()) + nSizeBMI;

		if (bmih->biBitCount == 24)
		{
			for (int c = 0; c < Width(); c++ )
			{
				memcpy(pbDest,pbSrc,nSizeRow);
				pbSrc += 3;
				pbDest += 3;
				*pbDest++ = 255;
			}
		}
		else 
		{
			memcpy(pbDest,pbSrc,nSizeRow);
		}
	}

	if (bmih->biClrUsed)
	{
		if (bUseBitmapCore==TRUE) {
			RGBTRIPLE* pcolormap = (RGBTRIPLE*) ((BYTE*)pDIB + bmih->biSize);
			bool bCorrect = false;

			for (int i=0;i<bmih->biClrUsed;i++) {
				if (pcolormap[i].rgbtBlue != pcolormap[i].rgbtGreen) bCorrect=true; 
				if (pcolormap[i].rgbtRed != pcolormap[i].rgbtGreen) bCorrect=true; 
			}		
		} else {
			RGBQUAD* pcolormap = (RGBQUAD*) ((BYTE*)pDIB + bmih->biSize);
			bool bCorrect = false;

			for (int i=0;i<bmih->biClrUsed;i++) {
				if (pcolormap[i].rgbBlue != pcolormap[i].rgbGreen) bCorrect=true; 
				if (pcolormap[i].rgbRed != pcolormap[i].rgbGreen) bCorrect=true; 
			}

		}
	}

	if (bmih->biHeight > 0)
		FlipVertical();

	// update history
	ClearHistory();
	AddToHistory("DIVA: Paste from clipboard (DIB2Image)");

	return true;
}


/////////////////////////////////////////////////////////////////////////////
//	Draw a rectangle into the image
//
//	input:
//		rFrame: rectangle
//		tVal: pixel value
//
//	author:	Rune Fisker, 8/4-1999
//	modified:	Rune Fisker, 24/11-1999	correction following error in DrawLine
//	DIVAUPDATE
template <class TPixel>
void CDImageBasic<TPixel>::DrawRect(CRect rFrame, TPixel tVal)
{
	DrawLine(rFrame.left, rFrame.top, rFrame.right, rFrame.top  , tVal);
	DrawLine(rFrame.right, rFrame.top ,rFrame.right ,rFrame.bottom  , tVal);
	DrawLine(rFrame.right, rFrame.bottom, rFrame.left, rFrame.bottom , tVal);
	DrawLine(rFrame.left, rFrame.bottom, rFrame.left, rFrame.top  ,  tVal);

	// update history
	CString str;
	str.Format("DIVA: draw rectangle (l=%i,t=%i,r=%i,b=%i) (DrawRect)",rFrame.left ,rFrame.top ,rFrame.right ,rFrame.bottom);
	AddToHistory(str);
}

/////////////////////////////////////////////////////////////////////////////
//	Draw a line into the image
//
//	input:
//		x1,y1,x2,y2: draw line from (x1,y1) to (x2,y2)
//		tVal: pixel value
//
//	author:		Rune Fisker, 8/4-1999
//	modified:	Rune Fisker, 17/9-1999	coordinates swap, ie. drawed points from (y,x) instead of (x,y)
//	DIVAUPDATE
template <class TPixel>
void CDImageBasic<TPixel>::DrawLine(int x1 ,int y1 ,int x2 ,int y2, TPixel tVal, bool fInsideCheck)
{
//	int     d,dx,dy;
//	int     aincr,bincr,yincr;
//	int     xx,yy;
	
	bool     flip=true;
	
	if (abs(x2-x1) < abs(y2-y1)) 
	{
		Swap(x1,y1);
		Swap(x2,y2);
		flip = false;
	}
	
	if (x1>x2) 
	{
		Swap(x1,x2);
		Swap(y1,y2);
	}
	
	int yincr;

	if (y2>y1) 
		yincr =  1;
	else
		yincr = -1;
	
	double dx=x2-x1;
	double dy=abs(y2-y1);
	double d=2*dy-dx;
	int aincr=2*(dy-dx);
	int bincr=2*dy;
	int xx=x1;
	int yy=y1;
	int nPoints = 0;


	if (fInsideCheck)
	{
		if (!flip) 
		{
			if (ContainsPoint(yy,xx))
				Pixel(yy,xx)=tVal;
	//		Pixel(xx,yy)=tVal;
		}
		else 
		{
			if (ContainsPoint(xx,yy))
				Pixel(xx,yy)=tVal;
	//		Pixel(yy,xx)=tVal;
		}

		for (xx=x1+1;xx<=x2;xx++) 
		{
			if (d>=0) 
			{
				yy+=yincr;
				d+=aincr;
			}
			else 
				d+=bincr;

			if (!flip) 
			{
				if (ContainsPoint(yy,xx))
					Pixel(yy,xx)=tVal;
	//			Pixel(xx,yy)=tVal;
			}
			else 
			{
				if (ContainsPoint(xx,yy))
					Pixel(xx,yy)=tVal;
	//			Pixel(yy,xx)=tVal;
			}
		}       
	}
	else
	{
		if (!flip) 
		{
			Pixel(yy,xx)=tVal;
	//		Pixel(xx,yy)=tVal;
		}
		else 
		{
			Pixel(xx,yy)=tVal;
	//		Pixel(yy,xx)=tVal;
		}

		for (xx=x1+1;xx<=x2;xx++) 
		{
			if (d>=0) 
			{
				yy+=yincr;
				d+=aincr;
			}
			else 
				d+=bincr;
			
			if (!flip) 
			{
				Pixel(yy,xx)=tVal;
				//			Pixel(xx,yy)=tVal;
			}
			else 
			{
				Pixel(xx,yy)=tVal;
				//			Pixel(yy,xx)=tVal;
			}
		}       
		
	}

	// update history
	CString str;
	str.Format("DIVA: draw line from (l=%i,t=%i) to (r=%i,b=%i) (DrawLine)",x1,y1,x2,y2);
	AddToHistory(str);
}

/////////////////////////////////////////////////////////////////////////////
// operator overload of +
//
// author:	Rune Fisker 9/4-1999
template <class TPixel>
inline CDImageBasic<TPixel> CDImageBasic<TPixel>::operator+(const double dbl)
{
	CDImageBasic<TPixel> image(*this);

	return (image += dbl);
}

/////////////////////////////////////////////////////////////////////////////
// operator overload of +
//
// author:	Rune Fisker 9/4-1999
template <class TPixel>
inline CDImageBasic<TPixel> CDImageBasic<TPixel>::operator+(const CDImageBasic<TPixel>& imageSrc) const
{
	assert((Width()-image.Width()) <= 0);
	assert((Height()-image.Height()) <= 0);

	CDImageBasic<TPixel> imageDest(*this);

	int yScr = imageDest.Top();

	for (int y = Top();y < Bottom(); y++)
	{
		TPixel *ppixelDest = &(imageDest.Pixel(Left(),y));
		TPixel *ppixelSrc = &(image.Pixel(Left(),yScr++));
		
		for (int x = Left();x < Right();x++)
		{
			*ppixelDest++ += *ppixelSrc++;
		}
	}

	// update history
	CString strName(image.Name(true));
	imageDest.AddToHistory("DIVA: add image (operator+ image)" + strName);

	return imageDest;
}


/////////////////////////////////////////////////////////////////////////////
// operator overload of -
//
// author:	Rune Fisker 9/4-1999
template <class TPixel>
inline CDImageBasic<TPixel> CDImageBasic<TPixel>::operator-(const double dbl)
{
	CDImageBasic<TPixel> image(*this);

	// update history
	AddToHistory("DIVA: Change sign (operator-)" + image.Name(true));

	return (image -= dbl);
}

/////////////////////////////////////////////////////////////////////////////
// operator overload of *
//
// author:	Rune Fisker 9/4-1999
template <class TPixel>
inline CDImageBasic<TPixel> CDImageBasic<TPixel>::operator*(const double dbl)
{
	CDImageBasic<TPixel> image(*this);

	return (image *= dbl);
}

/////////////////////////////////////////////////////////////////////////////
// operator overload of /
//
// author:	Rune Fisker 9/4-1999
template <class TPixel>
inline CDImageBasic<TPixel> CDImageBasic<TPixel>::operator/(const double dbl)
{
	assert(dbl != 0);

	CDImageBasic<TPixel> image(*this);

	return (image /= dbl);
}

/////////////////////////////////////////////////////////////////////////////
// operator overload of /=
//
// author:	Rune Fisker 21/2-1999
template <class TPixel>
inline CDImageBasic<TPixel>& CDImageBasic<TPixel>::operator/=(const double dbl)
{
	assert(dbl != 0);

	for (int y = Top();y < Bottom(); y++)
	{
		TPixel *ppixel = &(Pixel(Left(),y));
		for (int x = Left();x < Right();x++)
		{
			*ppixel++ /= dbl;
		}
	}

	// update history
	CString str;
	str.Format("DIVA: Divide by %f (operator/=)",dbl);
	AddToHistory(str);

	return *this;
}
	
/////////////////////////////////////////////////////////////////////////////
// operator overload of /=
//
// author:	Rune Fisker 21/2-1999
template <class TPixel>
inline CDImageBasic<TPixel>& CDImageBasic<TPixel>::operator/=(const CDImageBasic<TPixel>& image)
{
	assert((Width()-image.Width()) <= 0);
	assert((Height()-image.Height()) <= 0);

	int yScr = image.Top();;
	for (int y = Top();y < Bottom(); y++)
	{
		TPixel *ppixelDest = &(Pixel(Left(),y));
		const TPixel *ppixelSrc = &(image.Pixel(Left(),yScr++));
		
		for (int x = Left();x < Right();x++)
		{
			*ppixelDest++ /= *ppixelSrc++;
		}
	}

	// update history
	CString strName(image.Name(true));
	AddToHistory("DIVA: Divde by image (operator/=)" + strName);

	return *this;
}

/////////////////////////////////////////////////////////////////////////////
// operator overload of +=
//
// author:	Rune Fisker 21/2-1999
template <class TPixel>
inline CDImageBasic<TPixel>& CDImageBasic<TPixel>::operator+=(const double dbl)
{

	for (int y = Top();y < Bottom(); y++)
	{
		TPixel *ppixel = &(Pixel(Left(),y));
		for (int x = Left();x < Right();x++)
		{
			*ppixel++ += dbl;
		}
	}

	// add to history
	CString str;
	str.Format("DIVA: Add %f (operator+=)",dbl);
	AddToHistory(str);

	return *this;
}

/////////////////////////////////////////////////////////////////////////////
// operator overload of +=
//
// author:	Rune Fisker 21/2-1999
template <class TPixel>
inline CDImageBasic<TPixel>& CDImageBasic<TPixel>::operator+=(const CDImageBasic<TPixel>& image)
{
	assert((Width()-image.Width()) <= 0);
	assert((Height()-image.Height()) <= 0);

	int yScr = image.Top();;
	for (int y = Top();y < Bottom(); y++)
	{
		TPixel *ppixelDest = &(Pixel(Left(),y));
		const TPixel *ppixelSrc = &(image.Pixel(Left(),yScr++));
		
		for (int x = Left();x < Right();x++)
		{
			*ppixelDest++ += *ppixelSrc++;
		}
	}

	// update history
	CString strName(image.Name(true));
	AddToHistory("DIVA: Add image (operator+=)" + strName);

	return *this;
}

/////////////////////////////////////////////////////////////////////////////
// operator overload of -=
//
// author:	Rune Fisker 21/2-1999
template <class TPixel>
inline CDImageBasic<TPixel>& CDImageBasic<TPixel>::operator-=(const double dbl)
{
	for (int y = Top();y < Bottom(); y++)
	{
		TPixel *ppixel = &(Pixel(Left(),y));
		for (int x = Left();x < Right();x++)
		{
			*ppixel++ -= dbl;
		}
	}

	// update history
	CString str;
	str.Format("DIVA: Subtract %f (operator-=)",dbl);
	AddToHistory(str);

	return *this;
}
	
/////////////////////////////////////////////////////////////////////////////
// operator overload of -=
//
// author:	Rune Fisker 21/2-1999
template <class TPixel>
inline CDImageBasic<TPixel>& CDImageBasic<TPixel>::operator-=(const CDImageBasic<TPixel>& image)
{
	assert((Width()-image.Width()) <= 0);
	assert((Height()-image.Height()) <= 0);

	int yScr = image.Top();;
	for (int y = Top();y < Bottom(); y++)
	{
		TPixel *ppixelDest = &(Pixel(Left(),y));
		const TPixel *ppixelSrc = &(image.Pixel(Left(),yScr++));
		
		for (int x = Left();x < Right();x++)
		{
			*ppixelDest++ -= *ppixelSrc++;
		}
	}

	// update history
	CString strName(image.Name(true));
	AddToHistory("DIVA: Subtract image (operator*=)" + strName);

	return *this;
}
	
/////////////////////////////////////////////////////////////////////////////
// operator overload of *=
//
// author:	Rune Fisker 22/2-1999
template <class TPixel>
inline CDImageBasic<TPixel>& CDImageBasic<TPixel>::operator*=(const double dbl)
{
	for (int y = Top();y < Bottom(); y++)
	{
		TPixel *ppixel = &(Pixel(Left(),y));
		for (int x = Left();x < Right();x++)
		{
			*ppixel++ *= dbl;
		}
	}

	// update history
	CString str;
	str.Format("DIVA: Multiply by %f (operator*=)",dbl);
	AddToHistory(str);

	return *this;
}

/////////////////////////////////////////////////////////////////////////////
// operator overload of *=
//
// author:	Rune Fisker 21/2-1999
template <class TPixel>
inline CDImageBasic<TPixel>& CDImageBasic<TPixel>::operator*=(const CDImageBasic<TPixel>& image)
{
	assert((Width()-image.Width()) <= 0);
	assert((Height()-image.Height()) <= 0);

	int yScr = image.Top();;
	for (int y = Top();y < Bottom(); y++)
	{
		TPixel *ppixelDest = &(Pixel(Left(),y));
		TPixel *ppixelSrc = &(image.Pixel(Left(),yScr++));
		
		for (int x = Left();x < Right();x++)
		{
			*ppixelDest++ *= *ppixelSrc++;
		}
	}

	// update history
	AddToHistory("DIVA: Multiply by image (operator*=)" + image.Name(true));

	return *this;
}

/////////////////////////////////////////////////////////////////////////////
//	Label a binary image using 4 neighbourhood
//
//	author:	Rune Fisker 11/1-1999
//	modified: Rune Fisker 13/8-1999
//		check for "more labels that pixel range" inserted
template <class TPixel>
void CDImageBasic<TPixel>::Label() 
{
	// label counter
	int iLabel = 1;

	// initial max number of labels for equavilance tabel. 
	int nMaxLabels = Right()+Bottom();
	
	// Label equavilance tabel. 
	// Holds the lowest number of the labels, which are equavilent to this label
	int* piLabelEqual = new int[nMaxLabels];
	
	if (piLabelEqual == NULL)
		throw CVisError("Out of memory",eviserrorMemory,"Label","DVisBasic.h", __LINE__);

	piLabelEqual[0] = 0;

	// handle top left border case
	if (Pixel(Left(),Top()) > 0)
	{
		Pixel(Left(),Top()) = iLabel;
		piLabelEqual[iLabel] = iLabel++;
	}

	// handle top border case
	for (int x = Left()+1; x < Right(); x++)
	{
		if (Pixel(x,Top()) > 0)
		{
			if (Pixel(x-1,Top()) == 0)
			{
				Pixel(x,Top()) = iLabel;
				piLabelEqual[iLabel] = iLabel++;
			}
			else
			{
				Pixel(x,Top()) = Pixel(x-1,Top());
			}
		}
	}

	// handle left border case
	for (int y = Top()+1; y < Bottom(); y++)
	{
		if (Pixel(Left(),y) > 0)
		{
			if (Pixel(Left(),y-1) == 0)
			{
				Pixel(Left(),y) = iLabel;
				piLabelEqual[iLabel] = iLabel++;
			}
			else
			{
				Pixel(Left(),y) = Pixel(Left(),y-1);
			}
		}
	}
	
	// handle rest of the image
	for (y = Top()+1; y < Bottom(); y++)
	{ 
		for (x = Left()+1; x < Right(); x++)
		{			
			if (Pixel(x,y) > 0)
			{
				// assign label
				if ((Pixel(x-1,y) == 0) & (Pixel(x,y-1) == 0))
				{
					Pixel(x,y) = iLabel;
					piLabelEqual[iLabel] = iLabel++;

					// throw error if more labels that pixel range
					if (iLabel > VisRangeMax(Pixel(x,y),false))
					{
						delete piLabelEqual;

						Threshold(1);

						throw CVisError("More labels than the pixel range. Convert image to higher pixel range.",eviserrorUnknown,"Label","DVisBasic.h", __LINE__);
					}

					// handle if array to small
					if(iLabel == nMaxLabels)
					{
						nMaxLabels *= 2;
		
						int* piLabelEqualOld = piLabelEqual;

						piLabelEqual = new int[nMaxLabels];
						
						if (piLabelEqual == NULL)
							throw CVisError("Out of memory",eviserrorMemory,"Label","DVisBasic.h", __LINE__);
					
						memcpy(piLabelEqual,piLabelEqualOld,sizeof(int)*nMaxLabels/2);

						delete piLabelEqualOld;
					}
				}
				else if ((Pixel(x-1,y) != 0) & (Pixel(x,y-1) == 0))
				{
					Pixel(x,y) = Pixel(x-1,y);  
				}
				else if ((Pixel(x-1,y) == 0) & (Pixel(x,y-1) != 0))
				{
					Pixel(x,y) = Pixel(x,y-1);
				}
				else if (Pixel(x-1,y) == Pixel(x,y-1))
				{
					Pixel(x,y) = Pixel(x-1,y);
				}
				else
				{
					Pixel(x,y) = max(Pixel(x-1,y),Pixel(x,y-1));
					
					// update correspondance tabel
					int nChangeLabel = max(piLabelEqual[(int) Pixel(x-1,y)],piLabelEqual[(int) Pixel(x,y-1)]);				
					piLabelEqual[(int)Pixel(x,y)] = min(piLabelEqual[(int) Pixel(x-1,y)],piLabelEqual[(int) Pixel(x,y-1)]);

					// update all other elements in the correspondance tabel, 
					// which are equal to nChangeLabel
					if (piLabelEqual[(int)Pixel(x,y)] != nChangeLabel)
					{
						for (int i=1;i<iLabel;i++)
						{
							if (nChangeLabel == piLabelEqual[i])
							{
								piLabelEqual[i] = piLabelEqual[(int) Pixel(x,y)];
							}
						}
					}
				}
			}
		}
	}
	
	// order the labels 
	int nextlabel = 1;

	for (int i=1; i < iLabel; i++)
	{
		if (piLabelEqual[i] == i)
		{
			piLabelEqual[i] = nextlabel++;
		}
		else
		{
			piLabelEqual[i] = piLabelEqual[piLabelEqual[i]];
		}		
	}

	// relabel the image
	for (y = Top(); y < Bottom(); y++) 
	{
		for (x = Left(); x < Right(); x++)
		{ 	
			Pixel(x,y) = piLabelEqual[(int) Pixel(x,y)];               
		}
	}

	delete piLabelEqual;

	// update history
	AddToHistory("DIVA: Label (label)");
}


/////////////////////////////////////////////////////////////////////////////
// Enlarge the image 
//
// input:
//	nFactor: the factor to enlarge the width and height of the image
//
// author:	Rune Fisker 22/2-1999
template <class TPixel>
void CDImageBasic<TPixel>::EnlargePyr(int nFactor)
{
	assert(nFactor > 1);
	
	CDImageBasic<TPixel> imageDest(Width()*nFactor, Height()*nFactor);
	int yDest = 0;

	for (int ySrc=Top();ySrc<Bottom();ySrc++,yDest += nFactor)
	{
		TPixel *ppixelSrc = &( Pixel(Left(),ySrc) );

		for (int xSrc=Left();xSrc<Right();xSrc++)
		{				
			// insert pixel value in nFactor x nFactor window
			for (int i=0;i<nFactor;i++)
			{
				TPixel *ppixelDest = &( imageDest.Pixel(nFactor*(xSrc-Left()),yDest+i) );

				for (int j=0;j<nFactor;j++)
				{
					*ppixelDest++ = *ppixelSrc;
				}
			}

			ppixelSrc++;
		}
	}


	Deallocate();
	Allocate(imageDest.Width(), imageDest.Height());

	imageDest.CopyPixelsTo(*this,evisnormalizeCopyBytesSameType);

	// update history
	CString str;
	str.Format("DIVA: enlarge by a factor %i (EnlargePyr)",nFactor);
	AddToHistory(str);

}


/////////////////////////////////////////////////////////////////////////////
// Reduce the image 
//
// input:
//	nFactor: the factor to reduce the width and height of the image
//
// author:	Rune Fisker 22/2-1999
template <class TPixel>
void CDImageBasic<TPixel>::ReducePyr(int nFactor)
{
	assert(nFactor > 1);
	
	CDImageBasic<TPixel> imageDest(Width()/nFactor, Height()/nFactor);

	double dSum;
	int nDiv = nFactor*nFactor;
	int yDest = 0;
	int nNewRight = nFactor*imageDest.Width()+Left();
	int nNewBottom = nFactor*imageDest.Height()+Top();


	for (int ySrc=Top();ySrc<nNewBottom;yDest++,ySrc+= nFactor)
	{
		TPixel *ppixelDest = &( imageDest.Pixel(0,yDest) );

		for (int xSrc=Left();xSrc<nNewRight;xSrc += nFactor)
		{				
			dSum = 0;

			// clac. mean in nFactor x nFactor window
			for (int i=0;i<nFactor;i++)
			{
				TPixel *ppixelSrc = &( Pixel(xSrc,ySrc+i) );

				for (int j=0;j<nFactor;j++)
				{
					dSum += *ppixelSrc++;
				}
			}

			*ppixelDest++ = dSum/nDiv;
		}
	}


	Deallocate();
	Allocate(imageDest.Width(), imageDest.Height());

	imageDest.CopyPixelsTo(*this,evisnormalizeCopyBytesSameType);

	// update history
	CString str;
	str.Format("DIVA:  reduce by a factor %i (ReducePyr)",nFactor);
	AddToHistory(str);

}



/////////////////////////////////////////////////////////////////////////////
// Remove labeled blobs closer to the border that a given distance
//
// input:
//	nDist: distance to the border
//
// author:	Rune Fisker 11/1-1999
template <class TPixel>
void CDImageBasic<TPixel>::BlobRemoveBorder(int nDist) 
{
	// find max label
	int	nMaxLabel = Max();

	// create array for new labels
	TPixel* ptpNewLabel = new TPixel[nMaxLabel+1];

	if (ptpNewLabel == NULL)
		throw CVisError("Out of memory",eviserrorMemory,"BlobRemoveBorder","DVisBasic.h", __LINE__);

	for (int i=0;i<=nMaxLabel;i++)
		ptpNewLabel[i] = 1;

	ptpNewLabel[0] = 0;

	// check left border
	for (int x=Left();x<Left()+nDist;x++)
	{
		for (int y=Top(); y<Bottom();y++)
		{
			if (Pixel(x,y) > 0)
				ptpNewLabel[(int) Pixel(x,y)] = 0;
		}
	}

	// check Right border
	for (x=Right()-nDist;x<Right();x++)
	{
		for (int y=Top(); y<Bottom();y++)
		{
			if (Pixel(x,y) > 0)
				ptpNewLabel[(int) Pixel(x,y)] = 0;
		}
	}

	
	// check top border
	for (x=Left();x<Right();x++)
	{
		for (int y=Top(); y<Top()+nDist;y++)
		{
			if (Pixel(x,y) > 0)
				ptpNewLabel[(int)Pixel(x,y)] = 0;
		}
	}

	// check bottom border
	for (x=Left();x<Right();x++)
	{
		for (int y=Bottom()-nDist; y<Bottom();y++)
		{
			if (Pixel(x,y) > 0)
				ptpNewLabel[(int)Pixel(x,y)] = 0;
		}
	}

	// opdated new label tabel, so the rest of the structures are label correct
	int nextlabel = 1;
	
	for (i=1; i<=nMaxLabel;i++)
	{
		if (ptpNewLabel[i])
		{
			ptpNewLabel[i] = nextlabel++;
		}
	}

	// set all the labeled structures to their new label
	for (x=Left();x<Right();x++)
	{
		for (int y=Top(); y<Bottom();y++)
		{      
			Pixel(x,y) = ptpNewLabel[(int) Pixel(x,y)];
		}
	}

	delete ptpNewLabel;

	// update history
	CString str;
	str.Format("DIVA: remove blobs close to the border than %i (BlobRemoveBorder)",nDist);
	AddToHistory(str);
}

/////////////////////////////////////////////////////////////////////////////
// Calc. the first order moment (center of mass) for labeled blobs
//
// input:
//	vXC: x-coordinate for first order moment
//	vYC: y-coordinate for first order moment
//
// author:	Rune Fisker 11/1-1999
template <class TPixel>
void CDImageBasic<TPixel>::BlobFirstMoment(CVisDVector& vXC, CVisDVector& vYC) 
{
	// get max label
	int nMaxLabel = Max();
	
	// resize vectors
	vXC.Resize(nMaxLabel+1);
	vYC.Resize(nMaxLabel+1);

	// create counter array
	CVisDVector vNPixels(nMaxLabel+1);

	// init arays
	vXC = 0;
	vYC = 0;
	vNPixels = 0;

	// sum over all pixels
	for (int x = Left(); x < Right(); x++)
	{
		for (int y = Top(); y < Bottom(); y++)
		{ 
			if (Pixel(x,y) != 0)
			{				
				vXC[(int)Pixel(x,y)] += x;
				vYC[(int)Pixel(x,y)] += y;
				
				vNPixels[(int)Pixel(x,y)] += 1;
			}
		}
	}

	// find the mean
	for (int i=1;i<=nMaxLabel;i++)
	{
		if (vNPixels[i] > 0)
		{
			vXC[i] /= vNPixels[i];
			vYC[i] /= vNPixels[i];
			
		}
	}
}

// sfj code begins

/////////////////////////////////////////////////////////////////////////////
// Interpolation function, allows READ ONLY access to the image between 
// the grid points
//
// author:	Søren Falch Jørgensen 14/7-1999


template <class TPixel>
inline double CDImageBasic<TPixel>::PixelI(float x, float y) const 
{
  // if image is of type gray make Byte display image
 //  if ((PixFmt() & evispixfmtGray))

  return (double) Pixel3(x, y) ;	// Bicubic interpolation: best but slowest
}

/////////////////////////////////////////////////////////////////////////////
// Interpolation function, allows READ ONLY access to the image between 
// the grid points
//
// author:	Søren Falch Jørgensen 14/7-1999
template <class TPixel>
inline TPixel CDImageBasic<TPixel>::Pixel0(float x, float y) const
{
	int Y = (int)(y+0.5);
	int X = (int)(x+0.5);

// MB
	//assert (X>0 && X<Width() && Y>0 && Y<Height());	
// mb

	return Pixel(X,Y);

}

/////////////////////////////////////////////////////////////////////////////
// Interpolation function, allows READ ONLY access to the image between 
// the grid points
//
// author:	Søren Falch Jørgensen 14/7-1999
template <class TPixel>
inline double CDImageBasic<TPixel>::Pixel1(float x, float y) const
{
	
	register int Y = (int)y;
	register int X = (int)x;
		
	register double c=y-Y;
	register double d=x-X;
	
	register double t1 = CB(Pixel(X,  Y), Pixel(X,  Y+1), c);
	register double t2 = CB(Pixel(X+1,Y), Pixel(X+1,Y+1), c);

	return CB(t1, t2, d);

}
//	
//   Linear interpolation
//
template <class TPixel>
inline double CDImageBasic<TPixel>::CB(TPixel y1, TPixel y2, double d) const
{
	return d*(y2-y1)+y1;		// = (1-d)*y1+d*y2;
}

// MB
// HUGE HACK!
typedef double TSum;
// mb

/////////////////////////////////////////////////////////////////////////////
// Interpolation function, allows READ ONLY access to the image between 
// the grid points
//
// author:	Søren Falch Jørgensen 14/7-1999
template <class TPixel>
inline double CDImageBasic<TPixel>::Pixel2(float x, float y) const
{
	if (y<0.5 || x<0.5 || y>= Height()-1.5 || x>= Width()-1.5)
		return 0;

	int Y= (int)(y+0.5); // Nearest Neighbour
	int X= (int)(x+0.5);

	double c= y-Y; // Distance form NN to point
	double d= x-X;

	TSum t1=CQ(Pixel(X-1,Y-1), Pixel(X,Y-1), Pixel(X+1,Y-1), d);
	TSum t2=CQ(Pixel(X-1,Y  ), Pixel(X,Y  ), Pixel(X+1,Y  ), d);
	TSum t3=CQ(Pixel(X-1,Y+1), Pixel(X,Y+1), Pixel(X+1,Y+1), d);

	return CQ(t1,t2,t3,c);
}

template <class TPixel>
inline double CDImageBasic<TPixel>::CQ(TPixel y1, TPixel y2, TPixel y3, double d) const
{

	return y1*((1+d)*(1+d) - 3.0*(1+d) + 2.25)/2.0 + (0.75 - d*d)*y2 + y3*((1-d)*(1-d) - 3.0*(1-d) + 2.25)/2.0;

//	return y1-d*((3.0*y1-4.0*y2+y3)-d*(y1-2.0*y2+y3))/2.0;
}

/////////////////////////////////////////////////////////////////////////////
// Interpolation function, allows READ ONLY access to the image between 
// the grid points
//
// author:	Søren Falch Jørgensen 14/7-1999
template <class TPixel>
inline double CDImageBasic<TPixel>::Pixel3(float x, float y) const
{

	int X= (int)(x);
	int Y= (int)(y);

	if (Y<1 || X<1 || Y+2>=Height() || X+2>=Width())	// Handle surrounding area as 0
		return 0;


	double c= y-Y;
	double d= x-X;

	double t1=CC2(Pixel(X-1,Y-1), Pixel(X,Y-1), Pixel(X+1,Y-1), Pixel(X+2,Y-1),d);
	double t2=CC2(Pixel(X-1,Y  ), Pixel(X,Y  ), Pixel(X+1,Y  ), Pixel(X+2,Y  ),d);
	double t3=CC2(Pixel(X-1,Y+1), Pixel(X,Y+1), Pixel(X+1,Y+1), Pixel(X+2,Y+1),d);
	double t4=CC2(Pixel(X-1,Y+2), Pixel(X,Y+2), Pixel(X+1,Y+2), Pixel(X+2,Y+2),d);

	return CC2(t1,t2,t3,t4,c);

}

//
//	Cubic interpolation
//
// A "hard" bicubic interpolation

template <class TPixel>
inline double CDImageBasic<TPixel>::CC(TPixel y1, TPixel y2, TPixel y3, TPixel y4, double d) const
{
	return (TSum)y2+d*((y3-y1)+d*((2.0*(y1-y2) + y3 - y4)+d*(y2-y1-y3+y4)));
}

//--------------------------------------------------------------------------
// A "soft" bicubic interpolation
template <class TPixel>
inline double CDImageBasic<TPixel>::CC2(TPixel y1, TPixel y2, TPixel y3, TPixel y4, double d) const {
	return (double)y2+d*((-0.5*y1+0.5*y3)+d*((y1-2.5*y2+2.0*y3-0.5*y4)+d*(-0.5*y1-1.5*y3+1.5*y2+0.5*y4)));
}

// SFJ code ends



////////////////////////////////////////////////////////////////////////////
//
//	colormaps copied from CVisImageBase.cpp
//
////////////////////////////////////////////////////////////////////////////

// Default alpha value used in color maps.
//enum { knVisAlphaDefault = 0 };
enum { knVisAlphaDefault = 0xff000000 };


// @mdata:(IMPL) const unsigned long | CVisImageBase | s_mpbulNibbleColorMapGray |
//
// The default colormap for 4-bit grayscale images.
// @comm This may change in the final release of the SDK.
template <class TPixel>
const unsigned long CDImageBasic<TPixel>::s_mpbulNibbleColorMapGrayBasic[16] =
{
	(const unsigned long) (0x000000 | knVisAlphaDefault), // 0000
	(const unsigned long) (0x111111 | knVisAlphaDefault), // 0x01
	(const unsigned long) (0x222222 | knVisAlphaDefault), // 0x02
	(const unsigned long) (0x333333 | knVisAlphaDefault), // 0x03
	(const unsigned long) (0x444444 | knVisAlphaDefault), // 0x04
	(const unsigned long) (0x555555 | knVisAlphaDefault), // 0x05
	(const unsigned long) (0x666666 | knVisAlphaDefault), // 0x06
	(const unsigned long) (0x777777 | knVisAlphaDefault), // 0x07
	(const unsigned long) (0x888888 | knVisAlphaDefault), // 0x08
	(const unsigned long) (0x999999 | knVisAlphaDefault), // 0x09
	(const unsigned long) (0xaaaaaa | knVisAlphaDefault), // 0x0a
	(const unsigned long) (0xbbbbbb | knVisAlphaDefault), // 0x0b
	(const unsigned long) (0xcccccc | knVisAlphaDefault), // 0x0c
	(const unsigned long) (0xdddddd | knVisAlphaDefault), // 0x0d
	(const unsigned long) (0xeeeeee | knVisAlphaDefault), // 0x0e
	(const unsigned long) (0xffffff | knVisAlphaDefault), // 0x0f
};


// @mdata:(IMPL) const unsigned long | CVisImageBase | s_mpbulNibbleColorMapDefault |
//
// The default colormap for 4-bit color images.
// @comm This may change in the final release of the SDK.
// LATER:  Find the standard VGA values for these colors
template <class TPixel>
const unsigned long CDImageBasic<TPixel>::s_mpbulNibbleColorMapDefaultBasic[16] =
{
	(const unsigned long) (0x000000 | knVisAlphaDefault), // 0000
	(const unsigned long) (0x0000ff | knVisAlphaDefault), // 0x01
	(const unsigned long) (0x00ff00 | knVisAlphaDefault), // 0x02
	(const unsigned long) (0xff0000 | knVisAlphaDefault), // 0x03
	(const unsigned long) (0x00ffff | knVisAlphaDefault), // 0x04
	(const unsigned long) (0xff00ff | knVisAlphaDefault), // 0x05
	(const unsigned long) (0xffff00 | knVisAlphaDefault), // 0x06
	(const unsigned long) (0xffffff | knVisAlphaDefault), // 0x07
	(const unsigned long) (0x808080 | knVisAlphaDefault), // 0x08
	(const unsigned long) (0x000080 | knVisAlphaDefault), // 0x09
	(const unsigned long) (0x008000 | knVisAlphaDefault), // 0x0a
	(const unsigned long) (0x800000 | knVisAlphaDefault), // 0x0b
	(const unsigned long) (0x008080 | knVisAlphaDefault), // 0x0c
	(const unsigned long) (0x800080 | knVisAlphaDefault), // 0x0d
	(const unsigned long) (0x808000 | knVisAlphaDefault), // 0x0e
	(const unsigned long) (0xff0080 | knVisAlphaDefault), // 0x0f
};

// @mdata:(IMPL) const unsigned long | CVisImageBase | s_mpbulByteColorMapGray |
//
// The default colormap for 8-bit grayscale images.
// @comm This may change in the final release of the SDK.
template <class TPixel>
const unsigned long CDImageBasic<TPixel>::s_mpbulByteColorMapGrayBasic[256] =
{
	(const unsigned long) (0x000000 | knVisAlphaDefault), // 0000
	(const unsigned long) (0x010101 | knVisAlphaDefault), // 0x01
	(const unsigned long) (0x020202 | knVisAlphaDefault), // 0x02
	(const unsigned long) (0x030303 | knVisAlphaDefault), // 0x03
	(const unsigned long) (0x040404 | knVisAlphaDefault), // 0x04
	(const unsigned long) (0x050505 | knVisAlphaDefault), // 0x05
	(const unsigned long) (0x060606 | knVisAlphaDefault), // 0x06
	(const unsigned long) (0x070707 | knVisAlphaDefault), // 0x07
	(const unsigned long) (0x080808 | knVisAlphaDefault), // 0x08
	(const unsigned long) (0x090909 | knVisAlphaDefault), // 0x09
	(const unsigned long) (0x0a0a0a | knVisAlphaDefault), // 0x0a
	(const unsigned long) (0x0b0b0b | knVisAlphaDefault), // 0x0b
	(const unsigned long) (0x0c0c0c | knVisAlphaDefault), // 0x0c
	(const unsigned long) (0x0d0d0d | knVisAlphaDefault), // 0x0d
	(const unsigned long) (0x0e0e0e | knVisAlphaDefault), // 0x0e
	(const unsigned long) (0x0f0f0f | knVisAlphaDefault), // 0x0f
	(const unsigned long) (0x101010 | knVisAlphaDefault), // 0x10
	(const unsigned long) (0x111111 | knVisAlphaDefault), // 0x11
	(const unsigned long) (0x121212 | knVisAlphaDefault), // 0x12
	(const unsigned long) (0x131313 | knVisAlphaDefault), // 0x13
	(const unsigned long) (0x141414 | knVisAlphaDefault), // 0x14
	(const unsigned long) (0x151515 | knVisAlphaDefault), // 0x15
	(const unsigned long) (0x161616 | knVisAlphaDefault), // 0x16
	(const unsigned long) (0x171717 | knVisAlphaDefault), // 0x17
	(const unsigned long) (0x181818 | knVisAlphaDefault), // 0x18
	(const unsigned long) (0x191919 | knVisAlphaDefault), // 0x19
	(const unsigned long) (0x1a1a1a | knVisAlphaDefault), // 0x1a
	(const unsigned long) (0x1b1b1b | knVisAlphaDefault), // 0x1b
	(const unsigned long) (0x1c1c1c | knVisAlphaDefault), // 0x1c
	(const unsigned long) (0x1d1d1d | knVisAlphaDefault), // 0x1d
	(const unsigned long) (0x1e1e1e | knVisAlphaDefault), // 0x1e
	(const unsigned long) (0x1f1f1f | knVisAlphaDefault), // 0x1f
	(const unsigned long) (0x202020 | knVisAlphaDefault), // 0x20
	(const unsigned long) (0x212121 | knVisAlphaDefault), // 0x21
	(const unsigned long) (0x222222 | knVisAlphaDefault), // 0x22
	(const unsigned long) (0x232323 | knVisAlphaDefault), // 0x23
	(const unsigned long) (0x242424 | knVisAlphaDefault), // 0x24
	(const unsigned long) (0x252525 | knVisAlphaDefault), // 0x25
	(const unsigned long) (0x262626 | knVisAlphaDefault), // 0x26
	(const unsigned long) (0x272727 | knVisAlphaDefault), // 0x27
	(const unsigned long) (0x282828 | knVisAlphaDefault), // 0x28
	(const unsigned long) (0x292929 | knVisAlphaDefault), // 0x29
	(const unsigned long) (0x2a2a2a | knVisAlphaDefault), // 0x2a
	(const unsigned long) (0x2b2b2b | knVisAlphaDefault), // 0x2b
	(const unsigned long) (0x2c2c2c | knVisAlphaDefault), // 0x2c
	(const unsigned long) (0x2d2d2d | knVisAlphaDefault), // 0x2d
	(const unsigned long) (0x2e2e2e | knVisAlphaDefault), // 0x2e
	(const unsigned long) (0x2f2f2f | knVisAlphaDefault), // 0x2f
	(const unsigned long) (0x303030 | knVisAlphaDefault), // 0x30
	(const unsigned long) (0x313131 | knVisAlphaDefault), // 0x31
	(const unsigned long) (0x323232 | knVisAlphaDefault), // 0x32
	(const unsigned long) (0x333333 | knVisAlphaDefault), // 0x33
	(const unsigned long) (0x343434 | knVisAlphaDefault), // 0x34
	(const unsigned long) (0x353535 | knVisAlphaDefault), // 0x35
	(const unsigned long) (0x363636 | knVisAlphaDefault), // 0x36
	(const unsigned long) (0x373737 | knVisAlphaDefault), // 0x37
	(const unsigned long) (0x383838 | knVisAlphaDefault), // 0x38
	(const unsigned long) (0x393939 | knVisAlphaDefault), // 0x39
	(const unsigned long) (0x3a3a3a | knVisAlphaDefault), // 0x3a
	(const unsigned long) (0x3b3b3b | knVisAlphaDefault), // 0x3b
	(const unsigned long) (0x3c3c3c | knVisAlphaDefault), // 0x3c
	(const unsigned long) (0x3d3d3d | knVisAlphaDefault), // 0x3d
	(const unsigned long) (0x3e3e3e | knVisAlphaDefault), // 0x3e
	(const unsigned long) (0x3f3f3f | knVisAlphaDefault), // 0x3f
	(const unsigned long) (0x404040 | knVisAlphaDefault), // 0x40
	(const unsigned long) (0x414141 | knVisAlphaDefault), // 0x41
	(const unsigned long) (0x424242 | knVisAlphaDefault), // 0x42
	(const unsigned long) (0x434343 | knVisAlphaDefault), // 0x43
	(const unsigned long) (0x444444 | knVisAlphaDefault), // 0x44
	(const unsigned long) (0x454545 | knVisAlphaDefault), // 0x45
	(const unsigned long) (0x464646 | knVisAlphaDefault), // 0x46
	(const unsigned long) (0x474747 | knVisAlphaDefault), // 0x47
	(const unsigned long) (0x484848 | knVisAlphaDefault), // 0x48
	(const unsigned long) (0x494949 | knVisAlphaDefault), // 0x49
	(const unsigned long) (0x4a4a4a | knVisAlphaDefault), // 0x4a
	(const unsigned long) (0x4b4b4b | knVisAlphaDefault), // 0x4b
	(const unsigned long) (0x4c4c4c | knVisAlphaDefault), // 0x4c
	(const unsigned long) (0x4d4d4d | knVisAlphaDefault), // 0x4d
	(const unsigned long) (0x4e4e4e | knVisAlphaDefault), // 0x4e
	(const unsigned long) (0x4f4f4f | knVisAlphaDefault), // 0x4f
	(const unsigned long) (0x505050 | knVisAlphaDefault), // 0x50
	(const unsigned long) (0x515151 | knVisAlphaDefault), // 0x51
	(const unsigned long) (0x525252 | knVisAlphaDefault), // 0x52
	(const unsigned long) (0x535353 | knVisAlphaDefault), // 0x53
	(const unsigned long) (0x545454 | knVisAlphaDefault), // 0x54
	(const unsigned long) (0x555555 | knVisAlphaDefault), // 0x55
	(const unsigned long) (0x565656 | knVisAlphaDefault), // 0x56
	(const unsigned long) (0x575757 | knVisAlphaDefault), // 0x57
	(const unsigned long) (0x585858 | knVisAlphaDefault), // 0x58
	(const unsigned long) (0x595959 | knVisAlphaDefault), // 0x59
	(const unsigned long) (0x5a5a5a | knVisAlphaDefault), // 0x5a
	(const unsigned long) (0x5b5b5b | knVisAlphaDefault), // 0x5b
	(const unsigned long) (0x5c5c5c | knVisAlphaDefault), // 0x5c
	(const unsigned long) (0x5d5d5d | knVisAlphaDefault), // 0x5d
	(const unsigned long) (0x5e5e5e | knVisAlphaDefault), // 0x5e
	(const unsigned long) (0x5f5f5f | knVisAlphaDefault), // 0x5f
	(const unsigned long) (0x606060 | knVisAlphaDefault), // 0x60
	(const unsigned long) (0x616161 | knVisAlphaDefault), // 0x61
	(const unsigned long) (0x626262 | knVisAlphaDefault), // 0x62
	(const unsigned long) (0x636363 | knVisAlphaDefault), // 0x63
	(const unsigned long) (0x646464 | knVisAlphaDefault), // 0x64
	(const unsigned long) (0x656565 | knVisAlphaDefault), // 0x65
	(const unsigned long) (0x666666 | knVisAlphaDefault), // 0x66
	(const unsigned long) (0x676767 | knVisAlphaDefault), // 0x67
	(const unsigned long) (0x686868 | knVisAlphaDefault), // 0x68
	(const unsigned long) (0x696969 | knVisAlphaDefault), // 0x69
	(const unsigned long) (0x6a6a6a | knVisAlphaDefault), // 0x6a
	(const unsigned long) (0x6b6b6b | knVisAlphaDefault), // 0x6b
	(const unsigned long) (0x6c6c6c | knVisAlphaDefault), // 0x6c
	(const unsigned long) (0x6d6d6d | knVisAlphaDefault), // 0x6d
	(const unsigned long) (0x6e6e6e | knVisAlphaDefault), // 0x6e
	(const unsigned long) (0x6f6f6f | knVisAlphaDefault), // 0x6f
	(const unsigned long) (0x707070 | knVisAlphaDefault), // 0x70
	(const unsigned long) (0x717171 | knVisAlphaDefault), // 0x71
	(const unsigned long) (0x727272 | knVisAlphaDefault), // 0x72
	(const unsigned long) (0x737373 | knVisAlphaDefault), // 0x73
	(const unsigned long) (0x747474 | knVisAlphaDefault), // 0x74
	(const unsigned long) (0x757575 | knVisAlphaDefault), // 0x75
	(const unsigned long) (0x767676 | knVisAlphaDefault), // 0x76
	(const unsigned long) (0x777777 | knVisAlphaDefault), // 0x77
	(const unsigned long) (0x787878 | knVisAlphaDefault), // 0x78
	(const unsigned long) (0x797979 | knVisAlphaDefault), // 0x79
	(const unsigned long) (0x7a7a7a | knVisAlphaDefault), // 0x7a
	(const unsigned long) (0x7b7b7b | knVisAlphaDefault), // 0x7b
	(const unsigned long) (0x7c7c7c | knVisAlphaDefault), // 0x7c
	(const unsigned long) (0x7d7d7d | knVisAlphaDefault), // 0x7d
	(const unsigned long) (0x7e7e7e | knVisAlphaDefault), // 0x7e
	(const unsigned long) (0x7f7f7f | knVisAlphaDefault), // 0x7f
	(const unsigned long) (0x808080 | knVisAlphaDefault), // 0x80
	(const unsigned long) (0x818181 | knVisAlphaDefault), // 0x81
	(const unsigned long) (0x828282 | knVisAlphaDefault), // 0x82
	(const unsigned long) (0x838383 | knVisAlphaDefault), // 0x83
	(const unsigned long) (0x848484 | knVisAlphaDefault), // 0x84
	(const unsigned long) (0x858585 | knVisAlphaDefault), // 0x85
	(const unsigned long) (0x868686 | knVisAlphaDefault), // 0x86
	(const unsigned long) (0x878787 | knVisAlphaDefault), // 0x87
	(const unsigned long) (0x888888 | knVisAlphaDefault), // 0x88
	(const unsigned long) (0x898989 | knVisAlphaDefault), // 0x89
	(const unsigned long) (0x8a8a8a | knVisAlphaDefault), // 0x8a
	(const unsigned long) (0x8b8b8b | knVisAlphaDefault), // 0x8b
	(const unsigned long) (0x8c8c8c | knVisAlphaDefault), // 0x8c
	(const unsigned long) (0x8d8d8d | knVisAlphaDefault), // 0x8d
	(const unsigned long) (0x8e8e8e | knVisAlphaDefault), // 0x8e
	(const unsigned long) (0x8f8f8f | knVisAlphaDefault), // 0x8f
	(const unsigned long) (0x909090 | knVisAlphaDefault), // 0x90
	(const unsigned long) (0x919191 | knVisAlphaDefault), // 0x91
	(const unsigned long) (0x929292 | knVisAlphaDefault), // 0x92
	(const unsigned long) (0x939393 | knVisAlphaDefault), // 0x93
	(const unsigned long) (0x949494 | knVisAlphaDefault), // 0x94
	(const unsigned long) (0x959595 | knVisAlphaDefault), // 0x95
	(const unsigned long) (0x969696 | knVisAlphaDefault), // 0x96
	(const unsigned long) (0x979797 | knVisAlphaDefault), // 0x97
	(const unsigned long) (0x989898 | knVisAlphaDefault), // 0x98
	(const unsigned long) (0x999999 | knVisAlphaDefault), // 0x99
	(const unsigned long) (0x9a9a9a | knVisAlphaDefault), // 0x9a
	(const unsigned long) (0x9b9b9b | knVisAlphaDefault), // 0x9b
	(const unsigned long) (0x9c9c9c | knVisAlphaDefault), // 0x9c
	(const unsigned long) (0x9d9d9d | knVisAlphaDefault), // 0x9d
	(const unsigned long) (0x9e9e9e | knVisAlphaDefault), // 0x9e
	(const unsigned long) (0x9f9f9f | knVisAlphaDefault), // 0x9f
	(const unsigned long) (0xa0a0a0 | knVisAlphaDefault), // 0xa0
	(const unsigned long) (0xa1a1a1 | knVisAlphaDefault), // 0xa1
	(const unsigned long) (0xa2a2a2 | knVisAlphaDefault), // 0xa2
	(const unsigned long) (0xa3a3a3 | knVisAlphaDefault), // 0xa3
	(const unsigned long) (0xa4a4a4 | knVisAlphaDefault), // 0xa4
	(const unsigned long) (0xa5a5a5 | knVisAlphaDefault), // 0xa5
	(const unsigned long) (0xa6a6a6 | knVisAlphaDefault), // 0xa6
	(const unsigned long) (0xa7a7a7 | knVisAlphaDefault), // 0xa7
	(const unsigned long) (0xa8a8a8 | knVisAlphaDefault), // 0xa8
	(const unsigned long) (0xa9a9a9 | knVisAlphaDefault), // 0xa9
	(const unsigned long) (0xaaaaaa | knVisAlphaDefault), // 0xaa
	(const unsigned long) (0xababab | knVisAlphaDefault), // 0xab
	(const unsigned long) (0xacacac | knVisAlphaDefault), // 0xac
	(const unsigned long) (0xadadad | knVisAlphaDefault), // 0xad
	(const unsigned long) (0xaeaeae | knVisAlphaDefault), // 0xae
	(const unsigned long) (0xafafaf | knVisAlphaDefault), // 0xaf
	(const unsigned long) (0xb0b0b0 | knVisAlphaDefault), // 0xb0
	(const unsigned long) (0xb1b1b1 | knVisAlphaDefault), // 0xb1
	(const unsigned long) (0xb2b2b2 | knVisAlphaDefault), // 0xb2
	(const unsigned long) (0xb3b3b3 | knVisAlphaDefault), // 0xb3
	(const unsigned long) (0xb4b4b4 | knVisAlphaDefault), // 0xb4
	(const unsigned long) (0xb5b5b5 | knVisAlphaDefault), // 0xb5
	(const unsigned long) (0xb6b6b6 | knVisAlphaDefault), // 0xb6
	(const unsigned long) (0xb7b7b7 | knVisAlphaDefault), // 0xb7
	(const unsigned long) (0xb8b8b8 | knVisAlphaDefault), // 0xb8
	(const unsigned long) (0xb9b9b9 | knVisAlphaDefault), // 0xb9
	(const unsigned long) (0xbababa | knVisAlphaDefault), // 0xba
	(const unsigned long) (0xbbbbbb | knVisAlphaDefault), // 0xbb
	(const unsigned long) (0xbcbcbc | knVisAlphaDefault), // 0xbc
	(const unsigned long) (0xbdbdbd | knVisAlphaDefault), // 0xbd
	(const unsigned long) (0xbebebe | knVisAlphaDefault), // 0xbe
	(const unsigned long) (0xbfbfbf | knVisAlphaDefault), // 0xbf
	(const unsigned long) (0xc0c0c0 | knVisAlphaDefault), // 0xc0
	(const unsigned long) (0xc1c1c1 | knVisAlphaDefault), // 0xc1
	(const unsigned long) (0xc2c2c2 | knVisAlphaDefault), // 0xc2
	(const unsigned long) (0xc3c3c3 | knVisAlphaDefault), // 0xc3
	(const unsigned long) (0xc4c4c4 | knVisAlphaDefault), // 0xc4
	(const unsigned long) (0xc5c5c5 | knVisAlphaDefault), // 0xc5
	(const unsigned long) (0xc6c6c6 | knVisAlphaDefault), // 0xc6
	(const unsigned long) (0xc7c7c7 | knVisAlphaDefault), // 0xc7
	(const unsigned long) (0xc8c8c8 | knVisAlphaDefault), // 0xc8
	(const unsigned long) (0xc9c9c9 | knVisAlphaDefault), // 0xc9
	(const unsigned long) (0xcacaca | knVisAlphaDefault), // 0xca
	(const unsigned long) (0xcbcbcb | knVisAlphaDefault), // 0xcb
	(const unsigned long) (0xcccccc | knVisAlphaDefault), // 0xcc
	(const unsigned long) (0xcdcdcd | knVisAlphaDefault), // 0xcd
	(const unsigned long) (0xcecece | knVisAlphaDefault), // 0xce
	(const unsigned long) (0xcfcfcf | knVisAlphaDefault), // 0xcf
	(const unsigned long) (0xd0d0d0 | knVisAlphaDefault), // 0xd0
	(const unsigned long) (0xd1d1d1 | knVisAlphaDefault), // 0xd1
	(const unsigned long) (0xd2d2d2 | knVisAlphaDefault), // 0xd2
	(const unsigned long) (0xd3d3d3 | knVisAlphaDefault), // 0xd3
	(const unsigned long) (0xd4d4d4 | knVisAlphaDefault), // 0xd4
	(const unsigned long) (0xd5d5d5 | knVisAlphaDefault), // 0xd5
	(const unsigned long) (0xd6d6d6 | knVisAlphaDefault), // 0xd6
	(const unsigned long) (0xd7d7d7 | knVisAlphaDefault), // 0xd7
	(const unsigned long) (0xd8d8d8 | knVisAlphaDefault), // 0xd8
	(const unsigned long) (0xd9d9d9 | knVisAlphaDefault), // 0xd9
	(const unsigned long) (0xdadada | knVisAlphaDefault), // 0xda
	(const unsigned long) (0xdbdbdb | knVisAlphaDefault), // 0xdb
	(const unsigned long) (0xdcdcdc | knVisAlphaDefault), // 0xdc
	(const unsigned long) (0xdddddd | knVisAlphaDefault), // 0xdd
	(const unsigned long) (0xdedede | knVisAlphaDefault), // 0xde
	(const unsigned long) (0xdfdfdf | knVisAlphaDefault), // 0xdf
	(const unsigned long) (0xe0e0e0 | knVisAlphaDefault), // 0xe0
	(const unsigned long) (0xe1e1e1 | knVisAlphaDefault), // 0xe1
	(const unsigned long) (0xe2e2e2 | knVisAlphaDefault), // 0xe2
	(const unsigned long) (0xe3e3e3 | knVisAlphaDefault), // 0xe3
	(const unsigned long) (0xe4e4e4 | knVisAlphaDefault), // 0xe4
	(const unsigned long) (0xe5e5e5 | knVisAlphaDefault), // 0xe5
	(const unsigned long) (0xe6e6e6 | knVisAlphaDefault), // 0xe6
	(const unsigned long) (0xe7e7e7 | knVisAlphaDefault), // 0xe7
	(const unsigned long) (0xe8e8e8 | knVisAlphaDefault), // 0xe8
	(const unsigned long) (0xe9e9e9 | knVisAlphaDefault), // 0xe9
	(const unsigned long) (0xeaeaea | knVisAlphaDefault), // 0xea
	(const unsigned long) (0xebebeb | knVisAlphaDefault), // 0xeb
	(const unsigned long) (0xececec | knVisAlphaDefault), // 0xec
	(const unsigned long) (0xededed | knVisAlphaDefault), // 0xed
	(const unsigned long) (0xeeeeee | knVisAlphaDefault), // 0xee
	(const unsigned long) (0xefefef | knVisAlphaDefault), // 0xef
	(const unsigned long) (0xf0f0f0 | knVisAlphaDefault), // 0xf0
	(const unsigned long) (0xf1f1f1 | knVisAlphaDefault), // 0xf1
	(const unsigned long) (0xf2f2f2 | knVisAlphaDefault), // 0xf2
	(const unsigned long) (0xf3f3f3 | knVisAlphaDefault), // 0xf3
	(const unsigned long) (0xf4f4f4 | knVisAlphaDefault), // 0xf4
	(const unsigned long) (0xf5f5f5 | knVisAlphaDefault), // 0xf5
	(const unsigned long) (0xf6f6f6 | knVisAlphaDefault), // 0xf6
	(const unsigned long) (0xf7f7f7 | knVisAlphaDefault), // 0xf7
	(const unsigned long) (0xf8f8f8 | knVisAlphaDefault), // 0xf8
	(const unsigned long) (0xf9f9f9 | knVisAlphaDefault), // 0xf9
	(const unsigned long) (0xfafafa | knVisAlphaDefault), // 0xfa
	(const unsigned long) (0xfbfbfb | knVisAlphaDefault), // 0xfb
	(const unsigned long) (0xfcfcfc | knVisAlphaDefault), // 0xfc
	(const unsigned long) (0xfdfdfd | knVisAlphaDefault), // 0xfd
	(const unsigned long) (0xfefefe | knVisAlphaDefault), // 0xfe
	(const unsigned long) (0xffffff | knVisAlphaDefault)  // 0xff
};



/////////////////////////////////////////////////////////////////////////////
//	Takes exp of each pixel. Supports only gray images
//
//	author:	Rune Fisker, 29/4-1998
template <class TPixel>
void CDImageBasic<TPixel>::Exp()
{
	for (int y = Top();y < Bottom(); y++)
	{
		for (int x = Left();x < Right(); x++)
		{
			Pixel(x,y) = exp(Pixel(x,y));
		}
	}

	// update history
	AddToHistory("DIVA: Exp of each pixel (Exp)");

}



/////////////////////////////////////////////////////////////////////////////
//	Inverts the image. Supports only gray images
//
//	author:	Klaus B. Hilger, 12/7-1999
template <class TPixel>
void CDImageBasic<TPixel>::Invert()
{
	//code for invert
	TPixel tmaxpixel = Max();
	TPixel tminpixel = Min();

	for (int y = Top();y < Bottom(); y++)
	{
		TPixel *pT = &Pixel(Left(),y);

		for (int x = Left();x < Right(); x++)
		{
			*pT++ = -(*pT) + (tmaxpixel+tminpixel);
		}
	}

	// update history
	AddToHistory("DIVA: Inverting the image (Invert)");
}

/////////////////////////////////////////////////////////////////////////////
//	Regulates the Brightness of an image. Supports only gray images
//
//	author:	Klaus B. Hilger, 12/7-1999
template <class TPixel>
void CDImageBasic<TPixel>::Brightness()
{
	//code for invert
	TPixel maxpixel = Max();

	for (int y = Top();y < Bottom(); y++)
	{
		TPixel *pT = &Pixel(Left(),y);

		for (int x = Left();x < Right(); x++)
		{
			*pT++ = -(*pT) + maxpixel;
		}
	}

	// update history
	AddToHistory("DIVA: Regulating the brightness of the image (Brightness)");
}



/////////////////////////////////////////////////////////////////////////////
//	Linear Mappings. Supports only gray images


/////////////////////////////////////////////////////////////////////////////
//	Genneral Linear Mapping. Supports only gray images
//
//	author:	Klaus B. Hilger, 12/7-1999
template <class TPixel>
void CDImageBasic<TPixel>::GeneralLinMap(const double& dG, const double& dB)
{
 
	for (int y = Top();y < Bottom(); y++)
	{
		TPixel *pT = &Pixel(Left(),y);

		for (int x = Left();x < Right(); x++)
		{
			*pT++=dG*(*pT)+dB;
		}
	}

	// update history
	
	CString sHistory;
	sHistory.Format("DIVA: General linear mapping using gain = %f, bias = %f (GeneralLinMap)",(double)dG, (double) dB);
	AddToHistory(sHistory);  //RUNE: how to move to Dlg OnDestroy

}

/////////////////////////////////////////////////////////////////////////////
//	Linear Mapping to desired Mean and STD. Supports only gray images
//
//	author:	Klaus B. Hilger, 12/7-1999
template <class TPixel>
void CDImageBasic<TPixel>::MeanSTDLinMap(const double& dDesiredMean, const double& dDesiredSTD)
{


	double dMean;
	double dSTD;
	Mean(dMean);
	Std(dSTD);

	if (dSTD<0) throw CVisError("Current STD is zero!!!",eviserrorPixFmt,"MeanSTDLinMap","DImageMyClass.inl", __LINE__);


	for (int y = Top();y < Bottom(); y++)
	{
		TPixel *pT = &Pixel(Left(),y);

		for (int x = Left();x < Right(); x++)
		{
			
			*pT++=(dDesiredSTD/dSTD)*((*pT)-dMean)+dDesiredMean;
	
		}
	}

	// update history
	
	CString sHistory;
	sHistory.Format("DIVA: MeanStd mapping with desired mean = %f, std = %f (MeanSTDLinMap)",(double)dDesiredMean, (double) dDesiredSTD);
	AddToHistory(sHistory);

}


/////////////////////////////////////////////////////////////////////////////
//	Linear Mapping to desired Max and Min. Supports only gray images
//
//	author:	Klaus B. Hilger, 12/7-1999
template <class TPixel>
void CDImageBasic<TPixel>::MaxMinLinMap(const double& dDesiredMax, const double& dDesiredMin)
{


	double dMaxPixel = Max();  
	double dMinPixel = Min();
	double dSpan = dMaxPixel-dMinPixel;
    if (dSpan<0) throw CVisError("Pixel span equal zero!!!",eviserrorPixFmt,"MaxMinLinMap","DImageMyClass.inl", __LINE__);
	double dDesiredSpan = dDesiredMax - dDesiredMin;


	for (int y = Top();y < Bottom(); y++)
	{
		TPixel *pT = &Pixel(Left(),y);

		for (int x = Left();x < Right(); x++)
		{
			
			*pT++=(dDesiredSpan/dSpan)*((*pT)-dMinPixel)+dDesiredMin;
		}
	}

	// update history
	
	CString sHistory;
	sHistory.Format("DIVA: MinMax mapping with desired min = %f, max = %f (MaxMinLinMap)",(double)dDesiredMin, (double) dDesiredMax);
	AddToHistory(sHistory);

}



/////////////////////////////////////////////////////////////////////////////
//	Non-Linear Mappings. Supports only gray images



/////////////////////////////////////////////////////////////////////////////
//	Gamma mapping. Supports only gray images
//
//	author:	Klaus B. Hilger, 12/7-1999
//
//          modified to use both RGBA and gray scale images by MBS 021401
//

template <class TPixel>
void CDImageBasic<TPixel>::GammaMapping(const double& dG)
{

    if (dG<0) throw CVisError("Parameter must be larger than zero!!!",eviserrorPixFmt,"GammaMapping","DImageMyClass.inl", __LINE__);

    TPixel pixVal;

    if ( PixFmt() & evispixfmtRGBA == false ) {

        // gray scale version

        // find the actual ranges an stretch between these
        TPixel pixMax = this->Max();  
        TPixel pixMin = this->Min();  
	    TPixel pixSpan = pixMax-pixMin;        
	    for (int y = Top();y < Bottom(); y++)
	    {
		    TPixel *ptPixel = RowPointer(y);
		    for (int x = Left();x < Right(); x++)
		    {
                pixVal = (ptPixel[x]-pixMin)/pixSpan;			
			    pixVal = pow(pixVal, dG);
			    ptPixel[x] = pixSpan*pixVal+pixMin;	
		    }
	    }

    } else {

        // RGBA version

        // stretch between [0;255] 
        // (implicitly assumes a cast conversion from a byte RGBA image)
        for (int y = Top();y < Bottom(); y++)
	    {		
            TPixel *ptPixel = RowPointer(y);
		    for (int x = Left();x < Right(); x++)
		    {
			    pixVal = ptPixel[x]/255.0;            
			    pixVal = pow(pixVal, dG);
			    Pixel(x,y) = 255.0*pixVal;	
		    }
	    }
    }

	// update history	
	CString sHistory;
	sHistory.Format("DIVA: Gamma mapping using g = %f (GammaMapping)",(double)dG);
	AddToHistory(sHistory);
}

/////////////////////////////////////////////////////////////////////////////
//	Hyperbolic mapping. Supports only gray images
//
//	author:	Klaus B. Hilger, 13/7-1999
template <class TPixel>
void CDImageBasic<TPixel>::HyperbolicMapping(const double& dK)
{

	if (dK<0) throw CVisError("Parameter must be larger than zero!!!",eviserrorPixFmt,"HyperbolicMapping","DImageMyClass.inl", __LINE__);

	double dMaxPixel = Max();  
	double dMinPixel = Min();
	double dSpan = dMaxPixel-dMinPixel;
 
    double dVal;
	
	for (int y = Top();y < Bottom(); y++)
	{
		TPixel *pT = &Pixel(Left(),y);
		
		for (int x = Left();x < Right(); x++)
		{
			dVal = (*pT - dMinPixel)/dSpan;
			if (dVal==0)
				*pT++=dMinPixel;
			else 
				if (dVal==1)
					*pT++=dMaxPixel;
				else
					*pT++= (dVal/((1-dK)*dVal+dK))*dSpan+dMinPixel;
		}
	}

	// update history
	
	CString sHistory;
	sHistory.Format("DIVA: Hyperbolic mapping using k = %f (HyperbolicMapping)",(double)dK);
	AddToHistory(sHistory);

}

/////////////////////////////////////////////////////////////////////////////
//	Logarithmic mapping. Supports only gray images
//
//	author:	Klaus B. Hilger, 13/7-1999
template <class TPixel>
void CDImageBasic<TPixel>::LogarithmicMapping(const double& dK)
{
	if (dK<0) throw CVisError("Parameter must be larger than zero!!!",eviserrorPixFmt,"LogarithmicMapping","DImageMyClass.inl", __LINE__);

	double dMaxPixel = Max();  
	double dMinPixel = Min();
	double dSpan = dMaxPixel-dMinPixel;
 
    double dVal;
	
	for (int y = Top();y < Bottom(); y++)
	{
		TPixel *pT = &Pixel(Left(),y);
		
		for (int x = Left();x < Right(); x++)
		{
			dVal = (*pT - dMinPixel)/dSpan;
			if (dVal==0)
				*pT++=dMinPixel;
			else 
				if (dVal==1)
					*pT++=dMaxPixel;
				else
					*pT++= (log(1+(exp(dK)-1)*dVal)/dK)*dSpan+dMinPixel;
		}
	}

	// update history
	
	CString sHistory;
	sHistory.Format("DIVA: Logarithmic mapping using k = %f (LogarithmicMapping)",(double)dK);
	AddToHistory(sHistory);

}

/////////////////////////////////////////////////////////////////////////////
//	Exponential mapping. Supports only gray images
//
//	author:	Klaus B. Hilger, 13/7-1999
template <class TPixel>
void CDImageBasic<TPixel>::ExponentialMapping(const double& dK)
{
	if (dK<0) throw CVisError("Parameter must be larger than zero!!!",eviserrorPixFmt,"ExponentialMapping","DImageMyClass.inl", __LINE__);

	double dMaxPixel = Max();  
	double dMinPixel = Min();
	double dSpan = dMaxPixel-dMinPixel;
 
    double dVal;
	
	for (int y = Top();y < Bottom(); y++)
	{
		TPixel *pT = &Pixel(Left(),y);
		
		for (int x = Left();x < Right(); x++)
		{
			dVal = (*pT - dMinPixel)/dSpan;
			if (dVal==0)
				*pT++=dMinPixel;
			else 
				if (dVal==1)
					*pT++=dMaxPixel;
				else
					*pT++= ((exp(dK*dVal)-1)/(exp(dK)-1))*dSpan+dMinPixel;
		}
	}

	// update history
	
	CString sHistory;
	sHistory.Format("DIVA: Exponential mapping using k = %f (ExponentialMapping)",(double)dK);
	AddToHistory(sHistory);

}

/////////////////////////////////////////////////////////////////////////////
//	Histogram. Supports only gray images
//
//	author:	Klaus B. Hilger, 16/7-1999
template <class TPixel>
void CDImageBasic<TPixel>::Histogram(const int& iBins, const double& dMin, const double& dMax, CDVector& vHistogram)
{
	if (iBins<0) throw CVisError("Number of bins must be larger than zero!!!",eviserrorPixFmt,"Histogram","DImageMyClass.inl", __LINE__);

	vHistogram.Resize(iBins);
	vHistogram = 0;

	double dBinSpan=(dMax-dMin)/iBins;
	int index;
	
	for (int y = Top();y < Bottom(); y++)
	{
		TPixel *pT = &Pixel(Left(),y);
		for (int x = Left();x < Right(); x++)
		{
			if (*pT >= dMin && *pT <= dMax) 
			{
				index = floor((*pT-dMin)/dBinSpan);
				if (*pT==dMax) index -= 1;
				vHistogram[index] += 1;
			}
			pT++;
		}
	}	
}


/////////////////////////////////////////////////////////////////////////////
//	BetaHistogramMatching. Supports only gray images
//
//	author:	Klaus B. Hilger, 16/7-1999
template <class TPixel>
void CDImageBasic<TPixel>::BetaHistogramMatching(const double& dMin, const double& dMax,const CDVector& vOrgHistogram, const double& dAlpha, const double& dBeta)
{
	// vHistogram is the current histogram of the image between dMin and dMax
	// alpha, beta describes the desired Beta destribution between dMin and dMax

	CDVector vHistogram = vOrgHistogram;
	
	int iLength = vHistogram.Length();
	if (iLength<=0) throw CVisError("Number of bins must be larger than zero!!!",eviserrorPixFmt,"BetaHistogramMatching","DImageMyClass.inl", __LINE__);


	int index;
	double dSpan =dMax-dMin;
	double dBinSpan = dSpan/iLength;

	// Calculate the current normalized cummulative histogram

	double dSum=0;
	for (int iH=0; iH<iLength; iH++)
	{
		dSum += vHistogram[iH];
		vHistogram[iH] = dSum; ///iNPix;
	}
	int iNPix=dSum;

	vHistogram.ToMatlab("DIVA.m","AH"," ",false);

	// Calculate desired normalized cummmulative Beta destribution
	CDVector vBetaHistogram;
	vBetaHistogram.Resize(iLength);
	vBetaHistogram = 0;

	if (dAlpha==1 && dBeta==1)
	{
		for (iH=0,dSum=0; iH<iLength; iH++)
			vBetaHistogram[iH] += (iH+1)*((double) iNPix)/iLength;
	}
	else
	{
		for (iH=0, dSum=0; iH<iLength; iH++)
		{
			double dX= ((double) iH + 0.5)/iLength;
			if (dX!=0 && dX!=1) dSum += exp((dAlpha-1)*log(dX))*exp((dBeta-1)*log(1-dX));
			vBetaHistogram[iH] = dSum;
		}
		vBetaHistogram /= (dSum/iNPix);
	}
	vBetaHistogram.ToMatlab("DIVA.m","AB"," ",true);
	
	
	// Creating LookUpTable
	CDVector vLookUpTable;
	vLookUpTable.Resize(iLength);
	vLookUpTable = 0;
	int iTargetIndex = 0;
	double dIndexBase; 
	int iI, iJ;
	iI = iJ = 0;
	
	for (iH=0; iH<iLength; iH++)
	{
		double dTarget = vHistogram[iH];  // Determine Target to find in the vBetaHistogram
		double test;
		while (vBetaHistogram[iJ]<dTarget)	{test=vBetaHistogram[iJ]; iJ++; }

		double testI=vBetaHistogram[iI];
		double testJ=vBetaHistogram[iJ];

		if (fabs(vBetaHistogram[iI]-dTarget)<=fabs(vBetaHistogram[iJ]-dTarget))
		{
			iTargetIndex=iI;
			iJ=iI;
		}
		else
		{
			if (fabs(vBetaHistogram[iJ]-dTarget)>fabs(vBetaHistogram[iJ-1]-dTarget))
			{
				iJ--;
				while ((vBetaHistogram[iJ]==vBetaHistogram[iJ-1]) && (iJ!=0)) iJ--;
			}
			iTargetIndex=iJ;
			iI=iJ;
			
		}
		test=vBetaHistogram[iJ-1];
		testI=vBetaHistogram[iJ];
		testJ=vBetaHistogram[iJ+1];
		
		if (iH==0) dIndexBase = iTargetIndex;
		vLookUpTable[iH]= (double) (dSpan/((iLength-1)-dIndexBase)*(iTargetIndex-dIndexBase)+dMin);  
	//	vLookUpTable[iH]= (double) (dSpan/(iLength-1)*iTargetIndex+dMin);  
	}

	
	vLookUpTable.ToMatlab("DIVA.m","L"," ",true);
	
	// Updating the pixel values	

	for (int y = Top();y < Bottom(); y++)
	{
		TPixel *pT = &Pixel(Left(),y);
		for (int x = Left();x < Right(); x++)
		{
			if (*pT>=dMin && *pT <= dMax) 
			{
				index = floor((*pT-dMin)/dBinSpan);
				if (*pT==dMax) index -= 1;
				*pT= (double) vLookUpTable[index];
			}
			pT++;
		}
	}	
}


/////////////////////////////////////////////////////////////////////////////
//	BetaHistogram. Supports only gray images
//
//	author:	Klaus B. Hilger, 16/7-1999
template <class TPixel>
void CDImageBasic<TPixel>::BetaHistogram(CDVector& vHistogram, CDVector& vCumHistogram, const int iBins, const double& dAlpha, const double& dBeta)
{
	
	vHistogram.Resize(iBins);	
	vHistogram = 0;
	vCumHistogram.Resize(iBins);	
	vCumHistogram = 0;
	double dSum=0;
	
	int iNPix=NPixels();
	if (dAlpha==1 && dBeta==1)
	{
		vHistogram = ((double) iNPix)/iBins;
		for (int iH=0; iH<iBins; iH++)
			vCumHistogram[iH] = (iH+1)*((double) iNPix)/iBins;;
	}
	else
	{
		for (int iH=0; iH<iBins; iH++)
		{
			double dX= ((double) iH + 0.5)/iBins;
			vHistogram = exp((dAlpha-1)*log(dX))*exp((dBeta-1)*log(1-dX));
			dSum += vHistogram[iH];
			vCumHistogram[iH] = dSum;
		}
		vHistogram /= (dSum/iNPix);
		vCumHistogram /= (dSum/iNPix);
	}




}


/////////////////////////////////////////////////////////////////////////////
//	HistogramMatching. Supports only gray images
//
//	author:	Klaus B. Hilger, 16/7-1999
template <class TPixel>
void CDImageBasic<TPixel>::HistogramMatching(CDVector& vRefHistogram, CDVector& vCumRefHistogram, CDVector& vImageHistogram, const double& dAlpha, const double& dBeta)
{
	
	double dERSum = 0;
	double dEISum = 0;
	double dERPercent;
	double dEIPercent;
	double dEMin;
	double dEDiff;
	
	int iBins = vRefHistogram.Length();

	double dRPop = vCumRefHistogram[iBins-1];
	int iIPop = NPixels();
	int iRMin;
	
	CDVector vERCum;
	vERCum.Resize(iBins);
	CDVector vEICum;
	vEICum.Resize(iBins);
	
	CDVector vTable;
	vTable.Resize(iBins);
	
	for (int iH = 0; iH < 256 ; iH++)
	{
		dERPercent = ((double) vRefHistogram[iH])/dRPop;
		dERSum = dERSum+dERPercent;
		vERCum[iH]=dERSum;
		
		dEIPercent = ((double) vImageHistogram[iH])/iIPop;
		dEISum = dEISum+dEIPercent;
		vEICum[iH]=dEISum;
	}
	
	int iRLast = 0;
	for (int iX = 0; iX < 256 ; iX++)
	{
		dEMin=2000;  // set larger than largest possible!!!!
		for (int iK = iRLast;iK<256; iK++) 
		{
			dEDiff=fabs(vEICum[iX]-vERCum[iK]);
			if (dEDiff<dEMin)
			{
				iRMin = iK;
				dEMin = dEDiff;
			}
		}
		iRLast = iRMin;
		vTable[iX]=iRMin;
	}
/*
	vEICum.ToMatlab("DIVAH2.m","CI"," ",true);
	vERCum.ToMatlab("DIVAH2.m","CR"," ",true);
	vTable.ToMatlab("\DIVAH2.m","L"," ",true);
*/
	double dSpan =255-0;
	double dBinSpan = dSpan/iBins;

	for (int y = Top();y < Bottom(); y++)
	{
		TPixel *pT = &Pixel(Left(),y);
		for (int x = Left();x < Right(); x++)
		{
			int index = floor((*pT-0)/dBinSpan);
			if (*pT==255) index -= 1;
			*pT++= (double) vTable[index];
		}
	}		
	
}