weight+=mask[xk+radius];
			}
			Mediate(x,y)=temp/weight;
		}
		Progress(y/2);
	}
//column process
	double* Column = new double[nHeight];
	for(x=0;x<nWidth;x++)
	{
		for(y=radius;y<nHeight-radius;y++)
		{
			temp=0;
			for(yk=-radius;yk<radius+1;yk++)
				temp+=Mediate(x,y+yk)*mask[yk+radius];
			Column[y]=temp;
		}
		for(y=0;y<radius;y++)
		{
			temp=0;
			weight=0;
			for(yk=-y;yk<radius+1;yk++)
			{
				temp+=Mediate(x,y+yk)*mask[yk+radius];
				weight+=mask[yk+radius];
			}
			Column[y]=temp/weight;
		}
        for(y=nHeight-radius;y<nHeight;y++)
		{
			temp=0;
			weight=0;
			for(yk=-radius;yk<nHeight-y;yk++)
			{
				temp+=Mediate(x,y+yk)*mask[yk+radius];
				weight+=mask[yk+radius];
			}
			Column[y]=temp/weight;
		}

		for(y=0;y<nHeight;y++)
		{
			Out(x,y)=(int)Column[y];
		}
		Progress(nHeight/2+x*nHeight/2/nWidth);
	}
	delete []Column;
	delete []mask;
//    Laplacian(nWidth,nHeight,lpOutput,lpMediate);
	delete []lpMediate;
//general method(hide)
	PutPoints(nWidth,nHeight,lpOut,lpOutput);
	delete lpInput;
	delete lpOutput;
}

int Conv2(int nWidth,int nHeight,double *lpOutput,BYTE *lpInput,double mask[],int windowX,int windowY)
{
	int x,y,xk,yk;
	int radiusX=windowX/2;
	int radiusY=windowY/2;
	double temp;
    for(y=radiusY;y<nHeight-radiusY;y++)
	{
		for(x=radiusX;x<nWidth-radiusX;x++)
		{
			temp=0;
			for(xk=-radiusX;xk<radiusX+1;xk++)
				for(yk=-radiusY;yk<radiusY+1;yk++)
					temp+=In(x+xk,y+yk)*mask[xk+radiusX+(yk+radiusY)*windowX];
			Out(x,y)=temp;
		}
	}    
	return 0;
}

int Zerocross(int nWidth,int nHeight,BYTE *lpOutput,double *lpInput,double thresh)
{
	int x,y;
//zerocross cal
	for(y=1;y<nHeight-1;y++)
	{
		for(x=1;x<nWidth-1;x++)
		{
			if(In(x,y)<0 && In(x+1,y)>0 &&
				fabs(In(x,y)-In(x+1,y))>thresh)
			{
				Out(x,y)=255;
			}
			else if(In(x-1,y)>0 && In(x,y)<0 &&
				fabs(In(x-1,y)-In(x,y))>thresh)
			{
				Out(x,y)=255;
			}
			else if(In(x,y)<0 && In(x,y+1)>0 &&
				fabs(In(x,y)-In(x,y+1))>thresh)
			{
				Out(x,y)=255;
			}
			else if(In(x,y-1)>0 && In(x,y)<0 &&
				fabs(In(x,y)-In(x,y-1))>thresh)
			{
				Out(x,y)=255;
			}
		/*	if((In(x,y)<0 && In(x+1,y)>0)||
			   (In(x,y)<0 && In(x,y+1)>0)||
			   (In(x,y)>0 && In(x+1,y)<0)||
			   (In(x,y)>0 && In(x,y+1)<0))
			{
				Out(x,y)=255;
			}
		*/	else if(In(x,y)==0)
			{
				if(In(x-1,y)<0 && In(x+1,y)>0 &&
				   fabs(In(x-1,y)-In(x+1,y))>2*thresh)
				{
					Out(x,y)=255;
				}
			    else if(In(x-1,y)>0 && In(x+1,y)<0 &&
				    fabs(In(x-1,y)-In(x+1,y))>2*thresh)
				{
				    Out(x,y)=255;
				}
			    else if(In(x,y-1)>0 && In(x,y+1)<0 &&
			    	fabs(In(x,y-1)-In(x,y+1))>2*thresh)
				{
				    Out(x,y)=255;
				}
			    else if(In(x,y+1)>0 && In(x,y-1)<0 &&
				    fabs(In(x,y+1)-In(x,y-1))>2*thresh)
				{
				    Out(x,y)=255;
				}
				else 
				{
				    Out(x,y)=0;
				}
			}
			else 
			{
				Out(x,y)=0;
			}
		}
	}
	return 0;
}

void LoG(int nWidth,int nHeight,BYTE *lpIn,BYTE *lpOut,void (* Progress)(int Pos))
{

	BYTE *lpInput=new BYTE[nWidth*nHeight];
	BYTE *lpOutput=new BYTE[nWidth*nHeight];
	GetPoints(nWidth,nHeight,lpIn,lpInput);
	double thresh=0,sigma=1;

    int x,y,xx,yy;
	int radius,window; //window of the converlutin kernel(window=2*radius+1)
	double dev_inv=0.5/(sigma*sigma);  //1/(2*sigma^2)
	radius = (int)ceil(3*sigma);
	window = radius*2+1;
    //generate mask of convolution (matlab code)
	//std2 = std*std;
    //h1 = exp(-(x.*x + y.*y)/(2*std2));
    
	double* mask = new double [window*window];
	double sum=0;
	double std2=2*sigma*sigma;
    for(y=0;y<radius;y++)
	{
		for(x=0;x<radius;x++)
		{
			yy=(y-radius)*(y-radius);
			xx=(x-radius)*(x-radius);
			mask[x+y*window]=exp(-(xx+yy)*dev_inv);
//			mask[window-1-x+y*window]=mask[x+y*window];
//			mask[window-1-x+(window-1-y)*window]=mask[x+y*window];
//			mask[x+(window-1-y)*window]=mask[x+y*window];
			sum=sum+4*mask[x+y*window];
		}
	}
	x=radius;
	for(y=0;y<radius;y++)
	{
		yy=(y-radius)*(y-radius);
		mask[x+y*window]=exp(-yy*dev_inv);
//		mask[x+(window-1-y)*window]=mask[x+y*window];
		sum=sum+2*mask[x+y*window];
	}
	y=radius;
	for(x=0;x<radius;x++)
	{
		xx=(x-radius)*(x-radius);
		mask[x+y*window]=exp(-xx*dev_inv);
//		mask[window-1-x+y*window]=mask[x+y*window];
		sum=sum+2*mask[x+y*window];
	}
	mask[radius+radius*window]=1;
	sum=sum+mask[radius+radius*window];
//h = h1.*(x.*x + y.*y - 2*std2)/(2*pi*(std^6));	
	double denominator=2*3.14159265*pow(sigma,6)*sum;
	sum=0;
	for(x=0;x<radius;x++)
	{
		for(y=0;y<radius;y++)
		{
			yy=(y-radius)*(y-radius);
			xx=(x-radius)*(x-radius);
			mask[x+y*window]=mask[x+y*window]*(xx+yy-2*std2)/denominator;
//			mask[window-1-x+y*window]=mask[x+y*window];
//			mask[window-1-x+(window-1-y)*window]=mask[x+y*window];
//			mask[x+(window-1-y)*window]=mask[x+y*window];
			sum=sum+4*mask[x+y*window];
		}
	}
	x=radius;
	for(y=0;y<radius;y++)
	{
		yy=(y-radius)*(y-radius);
		mask[x+y*window]=mask[x+y*window]*(yy-2*std2)/denominator;
//		mask[x+(window-1-y)*window]=mask[x+y*window];
		sum=sum+2*mask[x+y*window];
	}
	y=radius;
	for(x=0;x<radius;x++)
	{
		xx=(x-radius)*(x-radius);
		mask[x+y*window]=mask[x+y*window]*(xx-2*std2)/denominator;
//		mask[window-1-x+y*window]=mask[x+y*window];
		sum=sum+2*mask[x+y*window];
	}
	mask[radius+radius*window]=mask[radius+radius*window]*(-2*std2)/denominator;
	sum=sum+mask[radius+radius*window];
//h = h - sum(h(:))/prod(size(h)); % make the filter sum to zero
	double average=sum/(window*window);
	for(x=0;x<radius;x++)
	{
		for(y=0;y<radius;y++)
		{
			mask[x+y*window]=mask[x+y*window]-average;
			mask[window-1-x+y*window]=mask[x+y*window];
			mask[window-1-x+(window-1-y)*window]=mask[x+y*window];
			mask[x+(window-1-y)*window]=mask[x+y*window];
		}
	}
	x=radius;
	for(y=0;y<radius;y++)
	{
		mask[x+y*window]=mask[x+y*window]-average;
		mask[x+(window-1-y)*window]=mask[x+y*window];
	}
	y=radius;
	for(x=0;x<radius;x++)
	{
		mask[x+y*window]=mask[x+y*window]-average;
		mask[window-1-x+y*window]=mask[x+y*window];
	}
	mask[radius+radius*window]=mask[radius+radius*window]-average;

    double* lpMediate = new double[nWidth*nHeight];
	//2D convolution
	if(Conv2(nWidth,nHeight,lpMediate,lpInput,mask,window,window))
	{
		TRACE("error in Conv2");
		return;
	};
//calcaulate the thresh for default
//  thresh = .75*mean2(abs(b(rr,cc)));
	sum = 0;
    for(y=radius;y<nHeight-radius;y++)
	{
		for(x=radius;x<nWidth-radius;x++)
		{
			sum = sum + fabs(Mediate(x,y));
		}
	}
	thresh=0.5*sum/((nWidth-2*radius)*(nHeight-2*radius));
//zerocross cal
    if(Zerocross(nWidth,nHeight,lpOutput,lpMediate,thresh))
	{
		TRACE("error in Zerocross");
		return;
	}
	delete mask;
	PutPoints(nWidth,nHeight,lpOut,lpOutput);
	delete lpInput;
	delete lpOutput;
}


void BiValue1(int nWidth,int nHeight,BYTE *lpInput,BYTE *lpOutput)
{
	BYTE *lpPoints=new BYTE[nWidth*nHeight];
	GetPoints(nWidth,nHeight,lpInput,lpPoints);
	double t=GetNumber("請(qǐng)輸入閾值",128,255);
	BYTE threshold=(BYTE)t;
	int x,y;
	for(y=0;y<nHeight;y++)
	{
		for(x=0;x<nWidth;x++)
		{
			if (Point(x,y)>threshold) Point(x,y)=255;
			else Point(x,y)=0;
		}
	}

	PutPoints(nWidth,nHeight,lpOutput,lpPoints);
	delete lpPoints;
}



void Sobel(int nWidth,int nHeight,BYTE *lpInput,BYTE *lpOutput,void (* Progress)(int Pos))
{
	int x,y,x1,y1,i;
	BYTE *lpPoints=new BYTE[nWidth*nHeight];
	BYTE *lpPoints1=new BYTE[nWidth*nHeight];
	memset(lpPoints1,0,nWidth*nHeight);
	GetPoints(nWidth,nHeight,lpInput,lpPoints);

	int d,max;
	static s[8][9]={
		{-1,-2,-1,0,0,0,1,2,1},
		{0,-1,-2,1,0,-1,2,1,0},
		{1,0,-1,2,0,-2,1,0,-1},
		{2,1,0,1,0,-1,0,-1,-2},
		{1,2,1,0,0,0,-1,-2,-1},
		{0,1,2,-1,0,1,-2,-1,0},
		{-1,0,1,-2,0,2,-1,0,1},
		{-2,-1,0,-1,0,1,0,1,2}
	};
	for(y=1;y<nHeight-1;y++)
	{
		for(x=1;x<nWidth-1;x++)
		{
			max=0;
			for(i=0;i<8;i++)
			{
				d=0;
				for(y1=0;y1<3;y1++)
				for(x1=0;x1<3;x1++)
				{
					d+=s[i][x1+y1*3]*Point(x+x1-1,y+y1-1);
				}
				if (d>max) max=d;
			}
			if (max>255) max=255;
			Point1(x,y)=(BYTE)max;
		}
		Progress(y);
	}
	PutPoints(nWidth,nHeight,lpOutput,lpPoints1);
	delete lpPoints;
	delete lpPoints1;
}

void Prewitte(int nWidth,int nHeight,BYTE *lpInput,BYTE *lpOutput,void (* Progress)(int Pos))
{
	int x,y,x1,y1,i;
	BYTE *lpPoints=new BYTE[nWidth*nHeight];
	BYTE *lpPoints1=new BYTE[nWidth*nHeight];
	memset(lpPoints1,0,nWidth*nHeight);
	GetPoints(nWidth,nHeight,lpInput,lpPoints);

	int d,max;
	static s[8][9]={
		{-1,-1,-1,0,0,0,1,1,1},
		{0,-1,-1,1,0,-1,1,1,0},
		{1,0,-1,1,0,-1,1,0,-1},
		{1,1,0,1,0,-1,0,-1,-1},
		{1,1,1,0,0,0,-1,-1,-1},
		{0,1,1,-1,0,1,-1,-1,0},
		{-1,0,1,-1,0,1,-1,0,1},
		{-1,-1,0,-1,0,1,0,1,1}
	};
	for(y=1;y<nHeight-1;y++)
	{
		for(x=1;x<nWidth-1;x++)
		{
			max=0;
			for(i=0;i<8;i++)
			{
				d=0;
				for(y1=0;y1<3;y1++)
				for(x1=0;x1<3;x1++)
				{
					d+=s[i][x1+y1*3]*Point(x+x1-1,y+y1-1);
				}
				if (d>max) max=d;
			}
			if (max>255) max=255;
			Point1(x,y)=(BYTE)max;
		}
		Progress(y);
	}
	PutPoints(nWidth,nHeight,lpOutput,lpPoints1);
	delete lpPoints;
	delete lpPoints1;
}

void Roberts(int nWidth,int nHeight,BYTE *lpInput,BYTE *lpOutput,void (* Progress)(int Pos))
{
	int x,y;
	BYTE *lpPoints=new BYTE[nWidth*nHeight];
	BYTE *lpPoints1=new BYTE[nWidth*nHeight];
	memset(lpPoints1,0,nWidth*nHeight);
	GetPoints(nWidth,nHeight,lpInput,lpPoints);

	int d;
	for(y=0;y<nHeight-1;y++)
	{
		for(x=0;x<nWidth-1;x++)
		{
			d=2*(abs(Point(x,y)-Point(x+1,y+1))+abs(Point(x+1,y)-Point(x,y+1)));
			if (d>255) d=255;
			Point1(x,y)=(BYTE)d;
		}
		Progress(y);
	}
	PutPoints(nWidth,nHeight,lpOutput,lpPoints1);
	delete lpPoints;
	delete lpPoints1;
}

void PseudoColor1(int nWidth,int nHeight,BYTE *lpInput,BYTE *lpOutput,void (* Progress)(int Pos))
{
	int x,y;
	BYTE *lpPoints=new BYTE[nWidth*nHeight];
	int nByteWidth=nWidth*3;
	if (nByteWidth%4) nByteWidth+=4-(nByteWidth%4);
	GetPoints(nWidth,nHeight,lpInput,lpPoints);
	BYTE b,R,G,B;
	for(y=0;y<nHeight;y++)
	{
		for(x=0;x<nWidth;x++)
		{
			b=Point(x,y);
			if (b<64)
			{
				B=255;
				G=4*b;
				R=0;
			}
			else if (b<128)
			{
				B=(127-b)*4;
				G=255;
				R=0;
			}
			else if (b<192)
			{
				B=0;
				G=255;
				R=(b-128)*4;
			}
			else
			{
				B=0;
				G=(255-b)*4;
				R=255;
			}

			lpOutput[x*3+nByteWidth*y]=B;
			lpOutput[x*3+1+nByteWidth*y]=G;
			lpOutput[x*3+2+nByteWidth*y]=R;
		}
		Progress(y);
	}
	delete lpPoints;
}

BYTE psf2(BYTE x)
{
	double x1=x;
	return (BYTE)((5-x1/42.5)*x1);
}

void PseudoColor2(int nWidth,int nHeight,BYTE *lpInput,BYTE *lpOutput,void (* Progress)(int Pos))
{
	int x,y;
	BYTE *lpPoints=new BYTE[nWidth*nHeight];
	int nByteWidth=nWidth*3;
	if (nByteWidth%4) nByteWidth+=4-(nByteWidth%4);
	GetPoints(nWidth,nHeight,lpInput,lpPoints);
	BYTE b,R,G,B;
	for(y=0;y<nHeight;y++)
	{
		for(x=0;x<nWidth;x++)
		{
			b=Point(x,y);
			if (b<85)
			{
				B=psf2(84-b);
				G=psf2(b);
				R=0;
			}
			else if (b<170)
			{
				B=0;
				G=psf2(169-b);
				R=psf2(b-85);
			}
			else
			{
				B=psf2(b-170);
				G=0;
				R=psf2(255-b);
			}

			lpOutput[x*3+nByteWidth*y]=B;
			lpOutput[x*3+1+nByteWidth*y]=G;
			lpOutput[x*3+2+nByteWidth*y]=R;
		}
		Progress(y);
	}
	delete lpPoints;
}



typedef struct MaxValue{
	int Value;
	int Dist;
	int AngleNumber;
}MaxValue;

//////////////////////////////////////////////////////////////
bool HoughTF(long lWidth,long lHeight,BYTE *lpInput,BYTE *lpOutput)
{
	BYTE *lpDIBBits;

	BYTE *lpSrc; //point to the source bitmap
	BYTE *lpDst;  //point to the template bitmap
	BYTE *lpTrans; //point to the the transforming bitmap
	BYTE *lpNewDIBBits; //point to the memory bitmap

	int *lpTransArea; //point to the transforming area
	int iMaxDist; //the maximize distance of the transforming area
	int iMaxAngleNumber;// the maximize angle of the transforming area

	int iDist; //the transforming coordinate distance
	int iAngleNumber;// the transforming coordinate angle

	BYTE pixel;

	MaxValue MaxValue1; //store the first two biggest (r,a)
	MaxValue MaxValue2;
 

	lpDIBBits=(BYTE *)malloc(lWidth*lHeight);
	BiValue1(lWidth,lHeight,lpInput,lpOutput);
	GetPoints(lWidth,lHeight,lpOutput,lpDIBBits);
	

	lpNewDIBBits=(BYTE *)malloc(lWidth*lHeight);
	if(lpNewDIBBits==NULL)
		return false;

	//initialize the allocated memory with value 255
	lpDst=(BYTE *)lpNewDIBBits;
	//memset(lpDst,(BYTE)255,lWidth*lHeight);
	memcpy(lpDst,lpDIBBits,lWidth*lHeight);

	//calculate the transforming area
	iMaxDist=(int)(sqrt((double)lWidth*lWidth+lHeight*lHeight)+0.5);
    //angle from 0 to 180 degree,2 degree per grid
	iMaxAngleNumber=90;

	//allocate memory for the transforming area and initialized with 0
	lpTransArea=(int *)malloc(iMaxAngleNumber*iMaxDist*sizeof(int));
	if(lpTransArea==NULL)
		return false;
	lpTrans=(BYTE *)lpTransArea;
	memset(lpTrans,0,iMaxAngleNumber*iMaxDist*sizeof(int));

	for(int j=0;j<lHeight;j++)
	{
		for(int i=0;i<lWidth;i++)
		{
			lpSrc=(BYTE *)lpDIBBits+lWidth*j+i;
			pixel=(BYTE)*lpSrc;

			if(pixel!=255 && *lpSrc!=0)//the grey value is supposed  to be 0 or 255
				return false;
			if(pixel==255)//the white point
			{
				for(iAngleNumber=0;iAngleNumber<iMaxAngleNumber;iAngleNumber++)
				{
					iDist=(int)fabs(i*cos(iAngleNumber*2*pi/180.0)+j*sin(iAngleNumber*2*pi/180.0));
					if(iDist>0)
					  *(lpTransArea+iDist*iMaxAngleNumber+iAngleNumber)=\
					   *(lpTransArea+iDist*iMaxAngleNumber+iAngleNumber)+1;
				}
			}
		}
	}

	//find out the first two biggest points
	MaxValue1.Value=0;
	MaxValue2.Value=0;

	for(iDist=0;iDist<iMaxDist;iDist++)
	{
		for(iAngleNumber=0;iAngleNumber<iMaxAngleNumber;iAngleNumber++)
		{
			if((int)*(lpTransArea+iDist*iMaxAngleNumber+iAngleNumber)>MaxValue1.Value)
			{
				MaxValue1.Value=(int)*(lpTransArea+iDist*iMaxAngleNumber+iAngleNumber);
				MaxValue1.Dist=iDist;
				MaxValue1.AngleNumber=iAngleNumber;
			}


		}
	}

	//clear the nearby points
	for(iDist=-9;iDist<10;iDist++)
	{
		for(iAngleNumber=-1;iAngleNumber<2;iAngleNumber++)
			if(iDist+MaxValue1.Dist>0 && iDist+MaxValue1.Dist<iMaxDist \
				&& iAngleNumber+MaxValue1.AngleNumber>=0\
				&&iAngleNumber+MaxValue1.AngleNumber<iMaxAngleNumber)
			{
				*(lpTransArea+(iDist+MaxValue1.Dist)*iMaxAngleNumber+\
					(iAngleNumber + MaxValue1.AngleNumber))=0;
			}
	}
	

	for(iDist=0;iDist<iMaxDist;iDist++)
	{
		for(iAngleNumber=0;iAngleNumber<iMaxAngleNumber;iAngleNumber++)
		{
			if((int)*(lpTransArea+iDist*iMaxAngleNumber+iAngleNumber)>MaxValue2.Value\
				&&iDist!=MaxValue1.Dist )
			{
				MaxValue2.Value=(int)*(lpTransArea+iDist*iMaxAngleNumber+iAngleNumber);
				MaxValue2.Dist=iDist;
				MaxValue2.AngleNumber=iAngleNumber;
			}


		}
	}

	if(abs(MaxValue1.AngleNumber-MaxValue2.AngleNumber)<=2)
	{
		for(long j=0;j<lHeight;j++)
		{
			for(long i=0;i<lWidth;i++)
			{
				lpDst=(BYTE *)lpNewDIBBits+lWidth*j+i;
				//set the color value white if the pixel in the first line
				iDist=(int)fabs(i*cos(MaxValue1.AngleNumber*2*pi/180.0)+\
					j*sin(MaxValue1.AngleNumber*2*pi/180.0));
				if(iDist==MaxValue1.Dist)
					*lpDst=(BYTE)255;
				//set the color value white if the pixel in the second line
				iDist=(int)fabs(i*cos(MaxValue2.AngleNumber*2*pi/180.0)+\
					j*sin(MaxValue1.AngleNumber*2*pi/180.0));
				if(iDist==MaxValue2.Dist)
					*lpDst=(BYTE)255;