fig2.cpp

图像处理源代码VC

 **********************************************************************/
float Gauss( float mean, float deviation)
{
    int        u1,u2;
    float      v1,v2,S;
    double     tmpv;
    static int FirstFlag = 1;

    if ( FirstFlag)
    {
//	randomize();    /*  or use srand() function  */
		srand( (unsigned)time( NULL ) );
	FirstFlag = 0;
    }

    do
    {
      u1 = rand();
      u2 = rand();
      v1 = 2 * (float)u1 / RAND_MAX - 1;
      v2 = 2 * (float)u2 / RAND_MAX - 1;
      S  = v1 * v1 + v2 * v2;
    } while (S>1);

    tmpv = v1 * sqrt(-2*log(S)/S);
    /*  now tmpv is a random variable of standard narmal N(0,1)  */

    tmpv = tmpv * deviation + mean;

    return(tmpv);
}
void addnoise(unsigned char **Image, unsigned char **Imageout,
              int /*row*/col, int /*col*/row, int var)
{
	   int   i, j;
	   float temp;
		for(i=0;i<row;i++)
      for(j=0;j<col;j++)
	   {
		   temp = Gauss(0, sqrt(var));
		   if(temp+Image[i][j]<0)
		       Imageout[i][j] = 0;
		   if(temp+Image[i][j]>255)
		       Imageout[i][j] = 255;
		   else
		       Imageout[i][j] = temp + Image[i][j];
		}
}


/****************************************************************
 *     Function  :  Statistic
 *      Parameter :  Image, row, col    : image to be statisticed
 *                   x,y                : startpoint of the window to be statisticed
 *                   width              : width of the window
 *                   mean               : mean of the window
 *                   var                : variation of the window
 ****************************************************************/
void Statistic(unsigned char **Image, int row, int col,
		       int x, int y, int width, int *mean, int *var)

{
	int    i, j;
	float  sum;

	sum = 0.;
	for (i=x-width/2; i<=x+width/2; i++){
	     if(i<=0 || i>=row)
		 continue;
	     for(j=y-width/2; j<=y+width/2; j++){
		if (j<=0 || j>=col)
		    continue;
		sum += Image[i][j];
	     }
	}

	*mean = (int)(sum/(width*width));

	sum = 0.;
	for (i=x-width/2; i<=x+width/2; i++){
	     if(i<=0 || i>=row)
		 continue;
	     for(j=y-width/2; j<=y+width/2; j++){
		if(j<=0 || j>=col)
		   continue;
		sum += ((int)Image[i][j]-*mean)*((int)Image[i][j]-*mean);
	     }
	}

	*var = (int)(sum/(width*width-1));
}
/****************************************************************
 *      Function  :     Average
 *      Parameter :     Source, row, col : image to be smoothed
 *                      x, y             : startpoint of the window to be smoothed
 *                      width            : width of the window
 *                      Dest             : destination image
 *      Return    :     None
 ****************************************************************/
unsigned char Average(unsigned char **Source, int row, int col,
		    int x, int y, int width)

{
	 int      i, j;
	 float    sum;
	 unsigned char adest;
	 sum = 0.;

	 for (i=x-width/2; i<=x+width/2; i++)
	 {
	     if (i<=0 || i>=row)
		continue;
	     for(j=y-width/2; j<=y+width/2; j++)
	     {
		  if (j<=0 || j>=col)
		     continue;
		  sum+=Source[i][j];
	     }
	 }

	 adest=(unsigned char)(sum/width/width);
	 return(adest);
}
/****************************************************************
 *        Function   :   Smoothing                              *
 *                                                              *
 *        Parameters :   char **Source                          *
 *                       int   row                              *
 *                       int   col                              *
 *                       char **Dest                            *
 *        Return     :   None                                   *
 *                                                              *
 ****************************************************************/
void smoothing(unsigned char **Source, unsigned char **Dest,
               int /*row*/col, int row/*col*/, int NoiseVar)
{
	int     i, j;
	int     mean, var;

	for (i=0; i<row; i++)
	    for (j=0; j<col; j++)
	    {
		 Statistic((unsigned char **)Source, row, col, i, j, 9, &mean, &var);
		  /* WINWID=9    the big-smooth-window width. */
		 //printf("\npoint(%d, %d), mean: %d, var: %d", i, j, mean, var);
		 if (var>1.2*NoiseVar)
			 Dest[i][j]=Average((unsigned char **)Source, row, col, i, j, 3);

		     /*  AVGWID=3 the small-smooth-window width. */
		 else
			 Dest[i][j]=Average((unsigned char **)Source, row, col, i, j, 9);
		   /*WINWID=9  the big-smooth-window width. */
	    }
}

/******************************FUNCTION:DrawHough()************************************/
/*					FUNCTION:draw detected lines after detect;                        */
/*					ARGUMENT:unsigned char **hough,**buf;                         */
/*                      int row,col,rohmax,DigitalMax;                           */
/*					RETURN:  none;                                                    */
/**************************************************************************************/
#define EPS  0.1
void DrawHough(unsigned char **buf,	unsigned char **hough,int row,int col,
            	int DigitalMax,int rohmax,int LineTh,int maxm,int minm)
{
	register int i,j,k,m;
	double theta,pi;

	pi=atan(1.0)*4.0;

	//statistic the min and max points detected
	minm=buf[0][0];
	maxm=minm;
	for(i=0;i<rohmax;i++)
		for(j=0;j<DigitalMax;j++)
		{
			if(buf[i][j]>maxm) maxm=buf[i][j];
			if(buf[i][j]<minm) minm=buf[i][j];
		}
      //printf("\nthe detected minm and maxm points is: %d---%d",minm,maxm);
      //printf("\ninput your LineThresold (40 to my image):");
	//input your choice line threshold
       //scanf("%d",&LineTh);
	for(i=0;i<rohmax;i++)
		for(j=0;j<DigitalMax;j++)
		{
			if(buf[i][j]<LineTh)continue;
		     //	printf("\nthe detected points of your threshold: %d",buf[i][j]);
			//	if detected then save lines
			theta=(double)((float)j/(float)DigitalMax)*1.5*pi;
			//decide the vertical line
			if(fabs(theta-pi/2)<EPS){
				for(k=0;k<row;k++)
				hough[k][i]=255;}
			//decide the horizontical line
			else if(fabs(theta-pi)<EPS){
				for(k=0;k<col;k++)
				hough[i][k]=255; }
			//decide others lines
			else{
				for(k=0;k<col;k++)
				{
					m=k*tan(theta)-i/cos(theta);
					if((m>0)&&(m<row))
						hough[m][k]=255;
				}
				}
		 }
}
/****************************FUNCTION:HoughDetectLine()**********************************/
/*	             FUNCTION:hough detect line;                                           */
/*              ARGUMENT:unsigned char **image,**hough;                               */
/*						       int row,col;                                                 */
/*						       int DigitalMax;                                              */
/*			   	 RETURN  :none;                                                        */
/**************************************************************************************/
#define EPS  0.1
int HoughDetectLine(unsigned char **image, unsigned char **hough,
                     int row,int col,int DigitalMax,int LineTh,int maxm,int minm)
{
	register int x,y,i;
	int k,rohmax;
	unsigned char **buf;
	double theta,pi;

	//polar equation of line: roh=x*cos(theta)+y*sin(theta)
	rohmax=(int)(sqrt((double)row*row+(double)col*col));
	buf=(unsigned char **)fspace_2d(rohmax,DigitalMax,sizeof(unsigned char));
	if(!buf){
	  ffree_2d((void **)buf,rohmax);
	  return 0;
   }
	pi=atan(1.0)*4.0;

     //	printf("please wait...");
	for(x=0;x<rohmax;x++)
		for(y=0;y<col;y++)
			buf[x][y]=0;
	//statics the points of the origin image
	for(y=0;y<row;y++)
		for(x=0;x<col;x++)
		{
			if(!image[y][x]) continue;
			{
				for(i=0;i<DigitalMax;i++)
				{
					theta=(double)((float)i/(float)DigitalMax)*1.5*pi;
					//	k,i is on the same line	:
					k=(int)(-y*cos(theta)+x*sin(theta)+0.5);
					//if roh and theta fit the poly equation :buf[][]++
					if((k>0)&&(k<rohmax))
						buf[k][i]++;
				}

			}
		 }
	//	creat and save hough detect image
	DrawHough(buf,(unsigned char **)hough,row,col,DigitalMax,rohmax,LineTh,maxm,minm);
	ffree_2d((void **)buf,rohmax);

	return 1;
}

/*************************FUNCTION:HoughDetectCircle()***********************/
/*			FUNCTION:detect circle;                            */
/*			ARGUMENT:unsigned char **image,**hough;            */
/* 		 	int row,col,DigitalMax;                           */
/*			RETURN:  none;                                     */
/****************************************************************************/
int HoughDetectCircle(unsigned char **image,unsigned char **hough,
        int row,int col,int DigitalMax,int CircleTh,int minm,int maxm)
{
	register int i,j,k,m;
	int x,y,rmax;
	double pi,theta,*PtrSin,*PtrCos;
	unsigned char **buf;
	//arguments of circle
	struct circle{
		int CenterX;
		int CenterY;
		int r;
		//to every radius there is a max and min value of detected points
		int max;
		int min;
		} *CirPtr;

       //	rmax=min(row,col)/2;  ***********
	if (row<col) rmax=row/2;
        else rmax=col/2;
	pi=atan(1.0)*4.0;
	CirPtr=(struct circle*)calloc(rmax,sizeof(struct circle));
	buf=(unsigned char **)fspace_2d(row,col,sizeof(unsigned char));
	PtrSin=(double *)calloc(DigitalMax,sizeof(double));
	PtrCos=(double *)calloc(DigitalMax,sizeof(double));
	if(!CirPtr || !buf || !PtrSin || !PtrCos){
	  free(CirPtr);
	  free(PtrCos);
	  free(PtrSin);
	  ffree_2d((void **)buf,row);
	  return 0;
   }

     //@@@@@@**->*
    //	printf("Please wait...");
	//calculate sin and cos values from 0 to 2PI
	for(i=0;i<DigitalMax;i++)
		{
			theta=i*pi/90;
			PtrSin[i]=sin(theta);
			PtrCos[i]=cos(theta);
			}
	// statics points of the origin image
	//k is radius,i and j is the location of the midcircle points
	int test = 0;//test
	for(k=0;k<rmax;k++)
	{
		for(i=0;i<row;i++)
			for(j=0;j<col;j++)
				buf[i][j]=0;
		for(i=0;i<row;i++)
			for(j=0;j<col;j++)
			{
				if(!image[i][j])continue;
				{
					//the poly equation of circle is
					//x=x0+r*cos(theta),y=y0+r*sin(theta)
					for(m=0;m<DigitalMax;m++)
					{
						x=(int)(j+k*PtrCos[m]+0.5);
						y=(int)(i+k*PtrSin[m]+0.5);
						//if x,y fit the equation of circle then buf[][]++
						if((x>=0)&&(x<col)&&(y>=0)&&(y<row))
						buf[y][x]++;
					}
				}
			}
		CirPtr[k].CenterX=0;
		CirPtr[k].CenterY=0;
		CirPtr[k].r=0;
		//statics the max and min values of every raidus
		CirPtr[k].max=buf[0][0];
		CirPtr[k].min=buf[0][0];
		for(i=0;i<row;i++)
			for(j=0;j<col;j++)
			{
				if(buf[i][j]<CirPtr[k].min)
					CirPtr[k].min=buf[i][j];
				if(buf[i][j]>CirPtr[k].max){
					CirPtr[k].CenterX=j;
					CirPtr[k].CenterY=i;
					CirPtr[k].r=k;
					CirPtr[k].max=buf[i][j];}
			}
	}
	//to select a threshold of circles we set max and min of all radius
	//then you can select CircleTh between the min and max values
	minm=CirPtr[0].min;
	maxm=CirPtr[0].max;
	for(k=0;k<rmax;k++)
	{
		 if(CirPtr[k].max>maxm)maxm=CirPtr[k].max;
		 if(CirPtr[k].min<minm)minm=CirPtr[k].min;
	}
      //printf("\nmin and max detected points=%d---%d",minm,maxm);
       //printf("\ninput your circle threshold(110 to my image):");
	//input your choice threshold
       //scanf("%d",&CircleTh);
	//compare to the threshold then decide circle
	for(k=0;k<rmax;k++)
	{
		if(CirPtr[k].max>CircleTh)
		{
		       //printf("\ndetected radius=%d and points=%d",k,CirPtr[k].max);
			for(m=0;m<DigitalMax;m++)
			{
				//return origin location and create detected circles
				x=(int)(CirPtr[k].CenterX+CirPtr[k].r*PtrCos[m]);
				y=(int)(CirPtr[k].CenterY+CirPtr[k].r*PtrSin[m]);
				if((x>=0)&&(x<col)&&(y>=0)&&(y<row))
				  hough[y][x]=255;
				else   // czy add
				{
					if(y < 0) y = 0;
					if(x < 0) x = 0;
					if(x > col) x = col;
					if(y > row) y = row;
					
					hough[y][x]=0;      // czy add
				}
			}
		}
	}
	 //free the created space
	 free(CirPtr);
	 free(PtrCos);
	 free(PtrSin);
	 ffree_2d((void **)buf,row);
	 return 1;
}

/*
#include <stdio.h>
#include <alloc.h>
#include <dos.h>
#include <stdlib.h>
#include <conio.h>
#include <math.h>
#include "comlib.h"
*/
/***********************************************************
Function:		convolving image stored in array img.
Parameter:
					<img>:the array of the source image;
					<row>:the row of the array img;
					<col>:the column of the array img;
					<sigma>:the space constant of Gauss function.
Call by:			main().
************************************************************/
int **Convolve2D(unsigned char **img,int row,int col,float sigma)
{
	int M,M2,K=4232,i,j,k,l,ki,lj;
	long sum;
	int **Ilog;
	int **Mask;
	M=(int)ceil(3*1.414*sigma)*2-1;
	M2=(M-1)/2;
	Ilog=(int **)fspace_2d(row,col,sizeof(int));
	if(Ilog==NULL)
	{
	      //	printf("Alloc error");
		return NULL;
	}
	Mask=(int **)fspace_2d(M,M,sizeof(unsigned int));
	if(Mask==NULL)
	{
		//printf("Alloc error");
		return NULL;
	}
	for(i=-M2;i<=M2;i++)
		for(j=-M2;j<=M2;j++)
			Mask[i+M2][j+M2]=(2-(i*i+j*j)/(sigma*sigma))*(K*exp(-(i*i+j*j)/(2*sigma*sigma)));
	for(k=0;k<row;k++)
		for(l=0;l<col;l++)
		{
			sum=0;
			for(i=-M2;i<=M2;i++)
				for(j=-M2;j<=M2;j++)
				{
					if(k+i>=0&&l+j>=0&&k+i<row&&l+j<col)
						sum+=((long)img[k+i][l+j]*(long)Mask[i+M2][j+M2])/((long)M);
					else
					{
						if((ki=k+i)<0) ki=0;
						if((lj=l+j)<0) lj=0;
						if(ki>=row) ki=row-1;
						if(lj>=col) lj=col-1;
						sum+=((long)img[ki][lj]*(long)Mask[i+M2][j+M2])/((long)M);
					}
				}
			Ilog[k][l]=sum/((long)M);
		     //	printf("\r%d  %d  %d",k,l,Ilog[k][l]);
		}
	ffree_2d((void **)Mask,M);
	return(Ilog);
}
/***********************************************************
Function:		convolving image stored in array img.
Parameter:
					<img>:the array of the source image;
					<row>:the row of the array img;
					<col>:the column of the array img;
					<sigma>:the space constant of Gauss function.
Call by :		main().
************************************************************/
int **Convolve1D(unsigned char **img,int row,int col,float sigma)
{
	int M,M2,SK,i,j,k,l,ki,lj;
	int **Ilog;
	long sum,sum1;
	int *RMask,*CMask;
	M=(int)ceil(3*1.414*sigma)*2-1;
	M2=(M-1)/2;
	SK=4232;
	Ilog=(int **)fspace_2d(row,col,sizeof(int));
	if(Ilog==NULL)
	{
	     //	printf("Alloc error");
		return NULL;
	}
	RMask=(int *)calloc(M,sizeof(int));