#define EXTRA_NAME "@gausssmooth."
#include "loadbmp.h"

#define In(x,y) lpInput[(x)+(y)*nWidth]
#define Out(x,y) lpOutput[(x)+(y)*nWidth]
#define Mediate(x,y) lpMediate[(x)+(y)*nWidth]
#define Point(x,y) lpPoints[(x)+(y)*nWidth]
#define sigma 1

void Gauss()
{
	int x,y,xx,xk,yk;
	BYTE *lpInput=new BYTE[nWidth*nHeight];
	BYTE *lpOutput=new BYTE[nWidth*nHeight];
	GetPoints(lpInput);
	int radius,window; 
	double dev_inv=0.5/(sigma*sigma);
	radius = (int)ceil(3*sigma);
	window = radius*2+1;
	memcpy(lpOutput,lpInput,nWidth*nHeight);
	double* lpMediate=new double[nWidth*nHeight];

	double* mask=new double[window];
	double sum=0;
	double temp;
	double weight;
	for(x=0;x<radius;x++)
	{
		xx=(x-radius)*(x-radius);
		mask[x]=exp(-xx*dev_inv);
		mask[window-1-x]=mask[x];
		sum+=2*mask[x];
	}
	mask[radius]=1;
	sum+=1;
	for(x=0;x<window;x++)
	{
		mask[x]/=sum;
	}
	for(y=0;y<nHeight;y++)
	{
		for(x=radius;x<nWidth-radius;x++)
		{
			temp=0;
			for(xk=-radius;xk<radius+1;xk++)
				temp+=In(x+xk,y)*mask[xk+radius];
			Mediate(x,y)=temp;
		}
		for(x=0;x<radius;x++)
		{
			temp=0;
			weight=0;
			for(xk=-x;xk<radius+1;xk++)
			{
				temp+=In(x+xk,y)*mask[xk+radius];
				weight+=mask[xk+radius];
			}
			Mediate(x,y)=temp/weight;
		}
		for(x=nWidth-radius;x<nWidth;x++)
		{
			temp=0;
			weight=0;
			for(xk=-radius;xk<nWidth-x;xk++)
			{
				temp+=In(x+xk,y)*mask[xk+radius];
				weight+=mask[xk+radius];
			}
			Mediate(x,y)=temp/weight;
		}
	}
	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];
		}
	}
	delete []Column;
	delete []mask;
	delete []lpMediate;
	PutPoints(lpOutput);
	delete lpInput;
	delete lpOutput;
}

void main(int argc, char *argv[])
{
	if(argc==2)
		FileName=argv[1];
	else
		return;	
	OpenFile();
	Gauss();
	SaveAs();
}