log.cpp

关于月球车的边缘算子处理程序

#define EXTRA_NAME "@logedge."
#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]

int Conv2(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(BYTE *lpOutput,double *lpInput,double thresh)
{
	int x,y;
	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;
			}
			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()
{

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

    int x,y,xx,yy;
	int radius,window;
	double dev_inv=0.5/(sigma*sigma);
	radius = (int)ceil(3*sigma);
	window = radius*2+1;
    
	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);
			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);
		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);
		sum=sum+2*mask[x+y*window];
	}
	mask[radius+radius*window]=1;
	sum=sum+mask[radius+radius*window];	
	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;
			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;
		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;
		sum=sum+2*mask[x+y*window];
	}
	mask[radius+radius*window]=mask[radius+radius*window]*(-2*std2)/denominator;
	sum=sum+mask[radius+radius*window];
	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];
	if(Conv2(lpMediate,lpInput,mask,window,window))
	{
		return;
	};
	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));
    if(Zerocross(lpOutput,lpMediate,thresh))
	{
		return;
	}
	delete mask;
	PutPoints(lpOutput);
	delete lpInput;
	delete lpOutput;
}

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