/////////////////////////////////////////////////////////////////////////////
//
// Image class for filtering functions
//
// Copyright © 1999
//
// IMM, Department of Mathematical Modelling
// Technical University of Denmark, Building 321
// DK-2800 Lyngby, Denmark
// http://www.imm.dtu.dk/~diva
//
// author: Rune Fisker & Lars Gunder Knudsen 
//
// 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"


//  Generates a Sobel filter kernel
//  Input:
//		seqFilters:  list to add new filter in
//  Output:
//		Filter is added to seqFilters
//
//	author:	Lars Gunder Knudsen, 990126
//	modified: Rune Fisker 3/9-99
template <class TPixel>
void CDImageConv<TPixel>::GenFilterSobel(CVisSequence<float> &seqFilters){

	CDImageConv<float> mFilter1(3,3);
	mFilter1.SetName("3x3 Sobel");

	mFilter1.Pixel(0,0)=-1.0;mFilter1.Pixel(0,1)= 0.0;mFilter1.Pixel(0,2)= 1.0;
	mFilter1.Pixel(1,0)=-2.0;mFilter1.Pixel(1,1)= 0.0;mFilter1.Pixel(1,2)= 2.0;
	mFilter1.Pixel(2,0)=-1.0;mFilter1.Pixel(2,1)= 0.0;mFilter1.Pixel(2,2)= 1.0;

	seqFilters.Insert(0,mFilter1);

	mFilter1.Pixel(0,0)=-1.0;mFilter1.Pixel(0,1)=-2.0;mFilter1.Pixel(0,2)=-1.0;
	mFilter1.Pixel(1,0)= 0.0;mFilter1.Pixel(1,1)= 0.0;mFilter1.Pixel(1,2)= 0.0;
	mFilter1.Pixel(2,0)= 1.0;mFilter1.Pixel(2,1)= 2.0;mFilter1.Pixel(2,2)= 1.0;

	seqFilters.Insert(1,mFilter1);
}

//  Generates a Prewitt filter kernel
//  Input:
//		seqFilters:  list to add new filter in
//  Output:
//		Filter is added to seqFilters
//
//	author:	Lars Gunder Knudsen, 990126
template <class TPixel>
void CDImageConv<TPixel>::GenFilterPrewitt(CVisSequence<float> &seqFilters){

	CDVisImage<float> mFilter1(3,3);
	mFilter1.Pixel(0,0)=-1.0;mFilter1.Pixel(0,1)= 0.0;mFilter1.Pixel(0,2)= 1.0;
	mFilter1.Pixel(1,0)=-1.0;mFilter1.Pixel(1,1)= 0.0;mFilter1.Pixel(1,2)= 1.0;
	mFilter1.Pixel(2,0)=-1.0;mFilter1.Pixel(2,1)= 0.0;mFilter1.Pixel(2,2)= 1.0;

	seqFilters.Insert(0,mFilter1);

	mFilter1.Pixel(0,0)=-1.0;mFilter1.Pixel(0,1)=-1.0;mFilter1.Pixel(0,2)=-1.0;
	mFilter1.Pixel(1,0)= 0.0;mFilter1.Pixel(1,1)= 0.0;mFilter1.Pixel(1,2)= 0.0;
	mFilter1.Pixel(2,0)= 1.0;mFilter1.Pixel(2,1)= 1.0;mFilter1.Pixel(2,2)= 1.0;

	seqFilters.Insert(1,mFilter1);
}

//  Generates a Laplace (4n) filter kernel
//  Input:
//		seqFilters:  list to add new filter in
//  Output:
//		Filter is added to seqFilters
//
//	author:	Lars Gunder Knudsen, 990126
//	modified: Rune Fisker 3/9-99
template <class TPixel>
void CDImageConv<TPixel>::GenFilterLaplace4n(CVisSequence<float> &seqFilters)
{
	CDImage<float> mFilter1(3,3);

	mFilter1.Pixel(0,0)= 0.0;mFilter1.Pixel(0,1)=-1.0;mFilter1.Pixel(0,2)= 0.0;
	mFilter1.Pixel(1,0)=-1.0;mFilter1.Pixel(1,1)= 4.0;mFilter1.Pixel(1,2)=-1.0;
	mFilter1.Pixel(2,0)= 0.0;mFilter1.Pixel(2,1)=-1.0;mFilter1.Pixel(2,2)= 0.0;

	seqFilters.Insert(0,mFilter1);
}

//  Generates a Laplace (8n) filter kernel
//  Input:
//		seqFilters:  list to add new filter in
//  Output:
//		Filter is added to seqFilters
//
//	author:	Lars Gunder Knudsen, 990126
//	modified: Rune Fisker 3/9-99
template <class TPixel>
void CDImageConv<TPixel>::GenFilterLaplace8n(CVisSequence<float> &seqFilters)
{
	CDImage<float> mFilter1(3,3);

	mFilter1.Pixel(0,0)=-1.0;mFilter1.Pixel(0,1)=-1.0;mFilter1.Pixel(0,2)=-1.0;
	mFilter1.Pixel(1,0)=-1.0;mFilter1.Pixel(1,1)= 8.0;mFilter1.Pixel(1,2)=-1.0;
	mFilter1.Pixel(2,0)=-1.0;mFilter1.Pixel(2,1)=-1.0;mFilter1.Pixel(2,2)=-1.0;

	seqFilters.Insert(0,mFilter1);
}

//  Generates a Mean filter kernel
//  Input:
//		seqFilters:  list to add new filter in
//  Output:
//		Filter is added to seqFilters
//
//	author:	Lars Gunder Knudsen, 990126
//	modified: Rune Fisker 3/9-99
template <class TPixel>
void CDImageConv<TPixel>::GenFilterMean(CVisSequence<float> &seqFilters,int nX,int nY){

	CDImage<float> mFilter1(nX,nY);
	mFilter1FillPixels(1.0/(nX*nY));
	seqFilters.Insert(0,mFilter1);
}

//  Generates a Gauss filter kernel
//  Input:
//		seqFilters:  list to add new filter in
//  Output:
//		Filter is added to seqFilters
//
//	author:	Lars Gunder Knudsen, 990126
//	modified: Rune Fisker 3/9-99
template <class TPixel>
void CDImageConv<TPixel>::GenFilterGauss(CVisSequence<float> &seqFilters,int nX,int nY, double dStd){

	CDImage<float> mFilter(nX,nY);

	if (dStd == -1)
	{
		dStd = (min(nX,nY)-1)/2;
	}

	double dbl_y,dbl_px,dbl_py,dMean=0;
	double d2SqrStd = 2*dStd*dStd;
	double dbl_cx = (double)(nX-1)/2.0;
	double dbl_cy = (double)(nY-1)/2.0;

	// create filter
	for(int i=0;i<nX;i++){
		for(int j=0;j<nY;j++){
			// calculating gauss:
			dbl_px = dbl_cx-(double)i;
			dbl_py = dbl_cy-(double)j;
			dbl_y = exp( - (dbl_px*dbl_px + dbl_py*dbl_py) / d2SqrStd );
			mFilter.Pixel(i,j) = dbl_y;
			dMean += dbl_y;
		}
	}

	// allways make the filter sum to 1
	mFilter /= dMean;

	seqFilters.Insert(0,mFilter);
}

//  Generates a LoG filter kernel
//  Input:
//		seqFilters:  list to add new filter in
//  Output:
//		Filter is added to seqFilters
//
//	author:	Lars Gunder Knudsen, 990126
//	modified: Rune Fisker 3/9-99
template <class TPixel>
void CDImageConv<TPixel>::GenFilterLoG(CVisSequence<float> &seqFilters,int nX,int nY, double dParam1){

	CDImage<float> mFilter1(nX,nY);

	int i,j;
	double dbl_y,dbl_cx,dbl_cy,dbl_px,dbl_py,dbl_Param1Sq;
	dbl_Param1Sq = dParam1*dParam1;
	dbl_cx = (double)(nX-1)/2.0;
	dbl_cy = (double)(nY-1)/2.0;
	for(i=0;i<nX;i++){
		for(j=0;j<nY;j++){
			// calculating gauss:
			dbl_px = dbl_cx-(double)i;
			dbl_py = dbl_cy-(double)j;
			dbl_y = (-1.0/(3.1415926535*dbl_Param1Sq*dbl_Param1Sq))*
				(1-(dbl_px*dbl_px+dbl_py*dbl_py)/(2*dbl_Param1Sq))*
				exp(-(dbl_px*dbl_px+dbl_py*dbl_py)/(2*dbl_Param1Sq));
			mFilter1.Pixel(i,j) = dbl_y;
		}
	}

	seqFilters.Insert(0,mFilter1);
}

		
/////////////////////////////////////////////////////////////////////////////


/////////////////////////////////////////////////////////////////////////////
//	Convolving image with filter kernel 
//
//	author:		Lars Gunder Knudsen, 990126
//	modified:	Rune Fisker, 3/5-1999
template <class TPixel>
void CDImageConv<TPixel>::Conv2DReal(const CVisImage<float>& mFilter)
{
	// height and width of tmp image
	int nWidth=mFilter.Width()-1+Right()-Left();
	int nHeight=mFilter.Height()-1+Bottom()-Top();

	// distance to border of image...
	int n_dTop=0,n_dBottom=0,n_dLeft=0,n_dRight=0;
	// 'copy rect' - coordinates
	int n_cTop,n_cBottom,n_cLeft,n_cRight;
	// backup of RECT for orig. image
	int n_bTop,n_bBottom,n_bLeft,n_bRight;


	n_bLeft = Left();
	n_bTop = Top();
	n_bRight = Right();
	n_bBottom = Bottom();

	n_cLeft = mFilter.Width()/2;
	n_cTop = mFilter.Height()/2;
	n_cRight = nWidth-mFilter.Width()/2;
	n_cBottom = nHeight-mFilter.Height()/2;


	if(Left()>0){
		if(Left()>mFilter.Width()/2){
			n_dLeft = mFilter.Width()/2;
			n_cLeft = 0;
		} else {
			n_dLeft = /*n_cLeft - */Left();
			n_cLeft -=  Left();
		}
	}
	if(Top()>0){
		if(Top()>mFilter.Height()/2){
			n_dTop = mFilter.Height()/2;
			n_cTop = 0;
		} else {
			n_dTop = /*n_cTop - */Top();
			n_cTop -=  Top();
		}
	}
	if(Right()<=MemoryShape().right){
		if( (MemoryShape().right-Right()) >= (mFilter.Width()/2)){
			n_dRight = mFilter.Width()/2;
			n_cRight = nWidth;
		} else {
			n_dRight =  (mFilter.Width()/2) - (MemoryShape().right-Right());
			n_cRight +=  MemoryShape().right-Right();
		}
	}
	if(Bottom()<=MemoryShape().bottom){
		if( (MemoryShape().bottom-Bottom()) >= (mFilter.Height()/2)){
			n_dBottom = mFilter.Height()/2;
			n_cBottom = nHeight;
		} else {
			n_dBottom =  (mFilter.Height()/2) - (MemoryShape().bottom-Bottom());
			n_cBottom +=  MemoryShape().bottom-Bottom();
		}
	}

	// get source an destination ROI's for pixel copying
/*	CRect rSrc,rDest,rROI,rImage;

	rROI   = Rect();
	rImage = MemoryRect();

	// source ROI, i.e. copy all relavant pixels inside the image
	rSrc.left  = max(0, rROI.left - mFilter.Width()/2);
	rSrc.right = min(rImage.right, rROI.right + mFilter.Width()/2);
	rSrc.top   = max(0, rROI.top - mFilter.Height()/2);
	rSrc.bottom= min(rImage.bottom, rROI.bottom + mFilter.Height()/2);

	// destination ROI
	rDest = rSrc; // always same size
	rDest.OffsetRect(-rDest.TopLeft()); // topleft is zero
	rDest.OffsetRect(-rROI.left + mFilter.Width()/2 + rSrc.left, - rROI.top + mFilter.Height()/2 + rSrc.top);
*/
	SetRect(n_bLeft-n_dLeft,n_bTop-n_dTop,n_bRight+n_dRight,n_bBottom+n_dBottom);

	// add 1 in height and width to avoid problems with Pad
	CDImageConv<TPixel> imageConvSrc(nWidth+1,nHeight+1);
		
	imageConvSrc.SetRect(n_cLeft,n_cTop,n_cRight,n_cBottom);
	CopyPixelsTo(imageConvSrc,evisnormalizeCopyBytesSameType);

	// pad if nesesary
	if (((n_cRight-n_cLeft) != nWidth) || ((n_cBottom-n_cTop) != nHeight))
	{
		imageConvSrc.SetRect(imageConvSrc.MemoryRect());
		imageConvSrc.Pad(BorderMode(), n_cLeft, n_cTop, n_cRight, n_cBottom);
	}

	imageConvSrc.SetRect(0,0,nWidth,nHeight);
	SetRect(n_bLeft,n_bTop,n_bRight,n_bBottom);

	double dConv;
	const float *pdFilter;
	int u,v,x,y;
	TPixel *pPixel,*pPixelDest;

	CRect rConvROI(mFilter.Width()/2,mFilter.Height()/2,nWidth-mFilter.Width()/2,nHeight-mFilter.Height()/2);
	int iYDest = n_bTop; 

	// for every pixel in the ROI area
	for(y=rConvROI.top;y<rConvROI.bottom;y++,iYDest++)
	{
		pPixelDest = &(Pixel(n_bLeft,iYDest)); // pointer to first destination pixel in row iYDest
		for(x=rConvROI.left;x<rConvROI.right;x++)
		{
			dConv = 0;
			for(v=0;v<mFilter.Height();v++)
			{
				pPixel = &(imageConvSrc.Pixel(x-rConvROI.left,y+v-rConvROI.top)); // pointer to first source pixel in row y+v-rConvROI.top
				pdFilter = mFilter.RowPointer(v); // pointer the first column in row v in the filter
				for(u=0;u<mFilter.Width();u++)
				{
					dConv += *pPixel++ * (*pdFilter++); // do convolution
				}
			}
			*pPixelDest++ = (TPixel)dConv; // insert res. in dest.
		}
	}

	// add to history
	CString strName(mFilter.Name());
	AddToHistory("DIVA: Filtering by " + strName + " (Conv2DReal)");
}


template <class TPixel>
void CDImageConv<TPixel>::Conv2D(const CVisImage<float>& imgFilter)
{
	try // IPL
	{
		// only float images are supported for float kernels
		Convolve2DFPIpl(imgFilter);	
	}
	catch(CVisError error)
	{
		Conv2DReal(imgFilter);
	}
}


template <class TPixel>
void CDImageConv<TPixel>::Conv2D(const CVisSequence<float>& seqFilters, EDivaFilterCombine edivafiltercombine)
{
	switch (edivafiltercombine)
	{
		case edivafiltercombineSUM		:	Convolve2DFPIpl(seqFilters, IPL_SUM); break;
		case edivafiltercombineSUMSQ	:	Convolve2DFPIpl(seqFilters, IPL_SUMSQ); break;
		case edivafiltercombineSUMSQROOT:	Convolve2DFPIpl(seqFilters, IPL_SUMSQROOT); break;
		case edivafiltercombineMAX		:	Convolve2DFPIpl(seqFilters, IPL_MAX); break;
		case edivafiltercombineMIN		:	Convolve2DFPIpl(seqFilters, IPL_MIN); break;

		default: throw CVisError("This function does not support this combine method",eviserrorPixFmt,"Conv2D",__FILE__, __LINE__);
	}

	// no alternative to IPL yet
}


template <class TPixel>
void CDImageConv<TPixel>::MeanFilter(int nCols, int nRows)
{
	try
	{
		// try IPL blur
		BlurIpl(nCols,nRows);
	}
	catch(CVisError error)
	{
		// use our conv
//		CVisDMatrix mFilter(nRows,nCols);
		CDImageConv<float> mFilter(nRows,nCols);
		mFilter.FillPixels(1.0/(nRows*nCols));
		Conv2DReal(mFilter);
	}		
}

template <class TPixel>
void CDImageConv<TPixel>::FixedFilter(EDivaFixedFilter edivafixedfilter)
{
	CVisSequence<float> seqFilters;

	switch (edivafixedfilter)
	{
		case edivafixedfilterLAPLACIAN3x3: FixedFilterIpl(IPL_LAPLACIAN_3x3); AddToHistory("DIVA: Filtering by 3x3 Laplacain (FixedFilter)"); break;
		case edivafixedfilterLAPLACIAN5x5: FixedFilterIpl(IPL_LAPLACIAN_5x5); AddToHistory("DIVA: Filtering by 5x5 Laplacain (FixedFilter)"); break;
		case edivafixedfilterGAUSSIAN3x3: FixedFilterIpl(IPL_GAUSSIAN_3x3); AddToHistory("DIVA: Filtering by 3x3 Gaussian (FixedFilter)"); break;
		case edivafixedfilterGAUSSIAN5x5: FixedFilterIpl(IPL_GAUSSIAN_5x5); AddToHistory("DIVA: Filtering by 5x5 Gaussian (FixedFilter)"); break;
		case edivafixedfilterHIGHPASS3x3: FixedFilterIpl(IPL_HIPASS_3x3); AddToHistory("DIVA: Filtering by 3x3 Highpass (FixedFilter)"); break;
		case edivafixedfilterHIGHPASS5x5: FixedFilterIpl(IPL_HIPASS_5x5); AddToHistory("DIVA: Filtering by 5x5 Highpass (FixedFilter)"); break;
		case edivafixedfilterSHARPEN3x3: FixedFilterIpl(IPL_SHARPEN_3x3); AddToHistory("DIVA: Filtering by 3x3 Sharpen (FixedFilter)"); break;
		case edivafixedfilterSOBEL: GenFilterSobel(seqFilters); Conv2D(seqFilters,edivafiltercombineSUMSQROOT); break;

		default: throw CVisError("Can not handle this filter",eviserrorUnknown,"FixedFilter","DVisConvClass.h", __LINE__); 
	}
}

/////////////////////////////////////////////////////////////////////////////
//	computes the cross correlation between image and image mask
//	See documentation for further explanation
//
//	author:	Rune Fisker, 18/5-1999
template <class TPixel>
void CDImageConv<TPixel>::Corr2D(CDImageConv<float>& imageDest, const CVisImage<float>& imageMask, CDImageConv<float>* pimageSqr)
{
	CDImageConv<float> imageCares(imageMask.Rect());
	imageCares.FillPixels(1);

	Corr2D(imageDest, imageMask, imageCares, pimageSqr);
}

/////////////////////////////////////////////////////////////////////////////
//	computes the cross correlation between image and image mask
//	See documentation for further explanation
//
//	author:	Rune Fisker, 18/5-1999
template <class TPixel>
void CDImageConv<TPixel>::Corr2D(CDImageConv<float>& imageDest, const CVisImage<float>& imageMask, const float flDontCare, CDImageConv<float>* pimageSqr)
{
	// create matrix indicating cares (pixels that are not don't cares)
	// 1 = cares and 0 = dont cares
	CDImageConv<float> imageCares(imageMask.Rect());
	imageCares.FillPixels(1);

	for (int y=imageMask.Top();y<imageMask.Bottom();y++)
	{
		for (int x=imageMask.Left();x<imageMask.Right();x++)
		{
			if (imageMask.Pixel(x,y) == flDontCare)
			{
				imageCares.Pixel(x,y) = 0;
			}
		}
	}

	Corr2D(imageDest, imageMask, imageCares, pimageSqr);
}

/////////////////////////////////////////////////////////////////////////////
//	computes the cross correlation between image and image mask
//	See documentation for further explanation
//
//	author:	Rune Fisker, 18/5-1999
template <class TPixel>
void CDImageConv<TPixel>::Corr2D(CDImageConv<float>& imageDest, const CVisImage<float>& imageMask, const CDImageConv<float>& imageCares, CDImageConv<float>* pimageSqr)
{
	// assume same size
	assert(Width() == imageDest.Width());
	assert(Height() == imageDest.Height());

	// save ROI
	CRect rROI = Rect();
	SetRect(MemoryRect());

	// number of dont cares in the filter
	int nCoeff = 0;

	// matrix used to actual conv. with the image (all dont cares are set to zero)
	CDImageConv<float> imageMaskReal(imageMask.Rect());
	imageMaskReal.FillPixels(0);

	/////////////////////////////////////////////////////////////////////////////
	// calc. Sum_ij M and N * Sum_ij M^2
	// and check for dont cares
	float dSumM = 0, dSumMSqr = 0;

	for (int y=imageMask.Top();y<imageMask.Bottom();y++)
	{
		for (int x=imageMask.Left();x<imageMask.Right();x++)
		{
			if (imageCares.Pixel(x,y) == 1)
			{
				dSumM += imageMask.Pixel(x,y);
				dSumMSqr += imageMask.Pixel(x,y)*imageMask.Pixel(x,y);
				imageMaskReal.Pixel(x,y) = imageMask.Pixel(x,y);
				nCoeff++;
			}
		}
	}

	// is there any dont cares
	bool fDontCares = false;
	if (nCoeff != (imageMask.NPoints()))
	{
		fDontCares = true;
	}

	float dMDen = nCoeff * dSumMSqr - dSumM*dSumM;

	if (dMDen == 0)
		throw CVisError("The standard deviation of the filter coefficients is zero, so the croos correlation is undefined",edviserrorDivisionByZero,"Corr2D()",__FILE__, __LINE__);

	/////////////////////////////////////////////////////////////////////////////
	// generate N * Sum_ij I^2
	CDImageConv<float> imageNSumISqr(MemoryRect());
	imageNSumISqr.SetBorderMode(BorderMode());

	if (pimageSqr != NULL)
	{
		assert(Rect() == pimageSqr->Rect());
		
		// set rect to full image
		CRect rROISqr = pimageSqr->Rect();
		pimageSqr->SetRect(pimageSqr->MemoryRect());

		// copy squared values
		pimageSqr->CopyPixelsTo(imageNSumISqr,evisnormalizeCopyBytesSameType);

		// set ROI back
		pimageSqr->SetRect(rROISqr);
	}
	else
	{
		// copy pixels
		CopyPixelsTo(imageNSumISqr,evisnormalizeCopyBytesSameType);

		// square pixels
		imageNSumISqr.Sqr();
	}

	imageNSumISqr.SetRect(rROI);
	
	// clac. local sum
	if (fDontCares == false)
	{
		imageNSumISqr.MeanFilter(imageMask.Width(),imageMask.Height());
		imageNSumISqr *= nCoeff;
	}
	else
	{
		imageNSumISqr.Conv2D(imageCares);
	}

	imageNSumISqr *= nCoeff;

	/////////////////////////////////////////////////////////////////////////////
	// generate Sum_ij I
	CDImageConv<float> imageSumI(MemoryRect());
	imageSumI.SetBorderMode(BorderMode());

	CopyPixelsTo(imageSumI,evisnormalizeCopyBytesSameType);

	imageSumI.SetRect(rROI);
	
	// clac. local sum
	if (fDontCares == false)
	{
		imageSumI.MeanFilter(imageMask.Width(),imageMask.Height());
		imageSumI *= nCoeff;
	}
	else
	{
		imageSumI.Conv2D(imageCares);	
	}

	
	/////////////////////////////////////////////////////////////////////////////
	// generate N Sum_ij I*M (I*M = convolution)
	CDImageConv<float> imageNSumIM(MemoryRect());
	imageNSumISqr.SetBorderMode(BorderMode());

	CopyPixelsTo(imageNSumIM,evisnormalizeCopyBytesSameType);

	imageNSumIM.SetRect(rROI);
	imageNSumIM.Conv2D(imageMaskReal);	

	imageNSumIM *= nCoeff;

	/////////////////////////////////////////////////////////////////////////////
	// calc. correlation

	float *ppixelSumI,*ppixelNSumISqr,*ppixelNSumIM,flIDen,*ppixelDest;

	// restore ROI's
	SetRect(rROI);

	for (y=imageDest.Top();y<imageDest.Bottom();y++)
	{
		ppixelDest = &(imageDest.RowPointer(y)[imageDest.Left()]);
		ppixelSumI = &(imageSumI.RowPointer(y)[imageSumI.Left()]);
		ppixelNSumISqr = &(imageNSumISqr.RowPointer(y)[imageNSumISqr.Left()]);
		ppixelNSumIM = &(imageNSumIM.RowPointer(y)[imageNSumIM.Left()]);

		for (int x=imageDest.Left();x<imageDest.Right();x++)
		{
			flIDen = *ppixelNSumISqr++ - *ppixelSumI*(*ppixelSumI);

			if (flIDen > 0)
			{
				*ppixelDest++ = (*ppixelNSumIM++ - *ppixelSumI * dSumM)/sqrt(flIDen*dMDen);
			}
			else
			{
				*ppixelDest++ = FLT_EPSILON;
				*ppixelNSumIM++;
			}
			
			*ppixelSumI++;
		}
	}

}