fig2.cpp

图像处理源代码VC

   }

/*
   cout<<endl<<"Do you want filted Image without edge tracing?(Y/N)";
   cin>>YesNo;
   if(YesNo=='Y'||YesNo=='y')
	 {
     cout<<endl<<"Filted image name:"; cin>>outfile;
     for(i=0;i<row;i++) for(j=0;j<col;j++)
       img[i][j]=(unsigned char)(((edgedimg[i][j]-min)*255)/(max-min));
     OutputImageWithName(img,row,col,outfile);
   }

  cout<<endl<<"Do you want binary image?(Y/N)";
  cin>>YesNo;
  if(YesNo=='Y'||YesNo=='y')
  {
    cout<<endl<<"Binary image name:"; cin>>outfile;
    for(i=0;i<row;i++)  for(j=0;j<col;j++)
       img[i][j]=edgedimg[i][j]<0?0:255;
    OutputImageWithName(img,row,col,outfile);
  }
*/
//  cout<<endl<<"Do you want edged image?(Y/N)"; cin>>YesNo;
//  if(YesNo=='Y'||YesNo=='y')
  {
//   cout<<endl<<"Edge Image  name:"; cin>>outfile;
     for(i=1;i<row;i++)    for(j=1;j<col;j++)
		{
	    img[i][j]=0;
		 if(edgedimg[i][j]>0&&(edgedimg[i-1][j]<0||edgedimg[i][j-1]<0))
		   img[i][j]=255;
		 if(edgedimg[i][j]<0&&(edgedimg[i-1][j]>0||edgedimg[i][j-1]>0))
		   img[i][j]=255;
      }
//    OutputImageWithName(img,row,col,outfile);
  }// end edged image output.

  ffree_2d((void**)edgedimg,row);

  return 1;
}



struct compl{
   float x,y;
};


/*****************************  FILTER.C  ******************************/

/***************************  comp_mult  *******************************/
/*  Complex multiplition  */
void comp_mult(struct compl a,struct compl b,struct compl *pc)
{
  pc->x = a.x*b.x-a.y*b.y;
  pc->y = a.x*b.y+a.y*b.x;
}

/***************************  fft  **********************************/
/*  Do Fourier transformation of one line  */
void fft(struct compl *ca,int n,char inv)
{
  int i,j,k,l,le,le1,ip,sign;
  struct compl t,u,w;

  j = 1;
  for(i=1;i<=n-1;i++) {
    if (i<j) {
      t = ca[j-1];
      ca[j-1] = ca[i-1];
      ca[i-1] = t;
    }
    k = n/2;
    while(k<j) {
      j -= k;
      k = k/2;
    }
    j += k;
  }
  if (inv=='t') sign = -1; else sign = 1;
  for(l=1;l<=(log((double)n)/log(2.0)+0.5);l++) {
    le = exp(l*log(2.0))+0.5;
    le1 = le/2;
    u.x = 1;
    u.y = 0;
    w.x = cos(3.1415927/le1);
    w.y = sign*sin(3.1415927/le1);
    for(j=1;j<=le1;j++) {
      i = j;
      while(i<=n) {
	ip = i+le1;
	comp_mult(ca[ip-1],u,&t);
	ca[ip-1].x = ca[i-1].x-t.x;
	ca[ip-1].y = ca[i-1].y-t.y;
	ca[i-1].x = ca[i-1].x+t.x;
	ca[i-1].y = ca[i-1].y+t.y;
	i += le;
      }
      comp_mult(u,w,&t);
      u = t;
    }
  }
}

void filter(unsigned char **pv1,unsigned char **pv2,
				int row,int col, int choice, int n, double r0)
{
  char centerflag,saveflag,b;
  int i,j,n1,n2,nn1,nn2,ii,jj;
  float max;
  struct compl temcompl;
  struct compl **pv,*pu;
  double  dis,dis0;
  int TrimRow, TrimCol;
  void fatal_error(char *msg);

  int MaxImgSize=(row>col)?row:col;
  int FFTDataSize=2;
  int PowPara=0;
	while(FFTDataSize<MaxImgSize){
	 FFTDataSize=pow(2.0,PowPara);
	 PowPara++;
  }
  int startx,starty;
  startx=(FFTDataSize-col)/2;
  starty=(FFTDataSize-row)/2;

/*
  char msg[100];
  sprintf(msg,"Size=%d",FFTDataSize);
  fatal_error(msg);
*/

  b = 't';           /* the flag of FFT */
  centerflag = 'y';  /* move to center point */

  TrimRow=FFTDataSize;//row;
  TrimCol=FFTDataSize;//col;

  pv = (struct compl **)fspace_2d(TrimRow,TrimCol,sizeof(struct compl));
  pu = (struct compl *)calloc(TrimRow,sizeof(struct compl));

  for(i=0;i<TrimRow;i++)
  for(j=0;j<TrimCol;j++){
	  pv[i][j].x = 0;
	  pv[i][j].y = 0;
  }

  for(i=0;i<row;i++)
  for(j=0;j<col;j++){
	  if (((centerflag=='y')||(centerflag=='Y'))&&((i+j)%2!=0))
		 pv[starty+i][startx+j].x = -pv1[i][j]; /* move to center point */
		else
		 pv[starty+i][startx+j].x = pv1[i][j];
	}

/******************************* FFT***************************/
  for(i=0;i<TrimRow;i++)
  {
    fft(pv[i],TrimCol,b);
    for(j=0;j<TrimCol;j++)
    {
	    pv[i][j].x = pv[i][j].x/sqrt((double)TrimCol);
		  pv[i][j].y = pv[i][j].y/sqrt((double)TrimCol);
    }
  }
  for(j=0;j<TrimCol;j++)
  {
    for(i=0;i<TrimRow;i++)
      pu[i]=pv[i][j];
    fft(pu,TrimRow,b);
    for(i=0;i<TrimRow;i++)
    {
	    pv[i][j].x=pu[i].x/sqrt((double)TrimRow);
	    pv[i][j].y=pu[i].y/sqrt((double)TrimRow);
    }
  }

/*****************filtering through diffirent filter********************/

  /*  printf("Please input your choice:");
    printf("\n 1).BHPF\n 2).BLPF\n 3).EHPF\n 4).ELPF\n 5).IHPF\n 6).ILPF\n");
    scanf("%d",&choice);   */

  n1=TrimRow;n2=TrimCol;
  switch ( choice )
  {
    case 1:      /*   Butterwrothy high pass filter */
	  {
			for(i=0;i<n1;i++)
	    for(j=0;j<n2;j++)
		    if((i!=n1/2)&&(j!=n2/2))
		    {
					ii=(double)(i);jj=(double)(j);nn1=(double)(n1);nn2=(double)(n2);
		      dis=sqrt(((ii-nn1/2.0)*(ii-nn1/2.0)+(jj-nn2/2.0)*(jj-nn2/2.0)));
		      dis0=r0/dis;
		      pv[i][j].x=pv[i][j].x/(1+0.414*pow(dis0,2*n));
		      pv[i][j].y=pv[i][j].y/(1+0.414*pow(dis0,2*n));
		    }
		    else
		    { 
		      pv[i][j].x=0.0;
		      pv[i][j].y=0.0;
		    }
	  }
	  break;

	  case 2:    /* Butterworthy low pass filter */
	  {
	    for(i=0;i<n1;i++)
		  for(j=0;j<n2;j++)
	    {
		    ii=(double)(i);jj=(double)(j);nn1=(double)(n1);nn2=(double)(n2);
		    dis=sqrt(((ii-nn1/2.0)*(ii-nn1/2.0)+(jj-nn2/2.0)*(jj-nn2/2.0)));
		    dis0=dis/r0;
		    pv[i][j].x=pv[i][j].x/(1+0.414*pow(dis0,2*n));
		    pv[i][j].y=pv[i][j].y/(1+0.414*pow(dis0,2*n));
		  }
	  }
	  break;

	  case 3:   /* Exp high pass filter */
	  {
		  for(i=0;i<n1;i++)
		  for(j=0;j<n2;j++)
		    if((i!=n1/2)&&(j!=n2/2))
		    {
		      ii=(double)(i);jj=(double)(j);nn1=(double)(n1);nn2=(double)(n2);
		      dis=sqrt(((ii-nn1/2.0)*(ii-nn1/2.0)+(jj-nn2/2.0)*(jj-nn2/2.0)));
		      dis0=r0/dis;
					pv[i][j].x=pv[i][j].x/exp(0.347*pow(dis0,n));
					pv[i][j].y=pv[i][j].y/exp(0.347*pow(dis0,n));
		    }
		    else
		    {
		      pv[i][j].x=0.0;
		      pv[i][j].y=0.0;
		    }
	  }
	  break;

	  case 4:    /* Exp low pass filter */
	  {
		  for(i=0;i<n1;i++)
	    for(j=0;j<n2;j++)
	    {
	      ii=(double)(i);jj=(double)(j);nn1=(double)(n1);nn2=(double)(n2);
	      dis=sqrt(((ii-nn1/2.0)*(ii-nn1/2.0)+(jj-nn2/2.0)*(jj-nn2/2.0)));
	      dis0=dis/r0;
				pv[i][j].x=pv[i][j].x/exp(0.347*pow(dis0,n));
				pv[i][j].y=pv[i][j].y/exp(0.347*pow(dis0,n));
      }
	  }
	  break;

	  case 5:   /* Ideal high pass filter */
	  {
		  for(i=0;i<n1;i++)
	    for(j=0;j<n2;j++)
	    {
		    ii=(double)(i);jj=(double)(j);nn1=(double)(n1);nn2=(double)(n2);
		    dis=sqrt(((ii-nn1/2.0)*(ii-nn1/2.0)+(jj-nn2/2.0)*(jj-nn2/2.0)));
			  if( dis<=r0 )
			  {
				  pv[i][j].x=0;pv[i][j].y=0;
			  }
		  }
		}
	  break;

	  case 6:   /* Ideal low pass filter */
	  {
	 	  for(i=0;i<n1;i++)
		  for(j=0;j<n2;j++)
		  {
		    ii=(double)(i);jj=(double)(j);nn1=(double)(n1);nn2=(double)(n2);
		    dis=sqrt(((ii-nn1/2.0)*(ii-nn1/2.0)+(jj-nn2/2.0)*(jj-nn2/2.0)));
				if( dis>r0 )
			  {
				  pv[i][j].x=0;pv[i][j].y=0;
			  }
		  }
	  }
	  break;

		default:
	    ffree_2d((void **)pv,TrimRow);
	    free((void *)pu);
      return;
  }

/********************************IFFT*********************************/
   b = 'f';
   centerflag = 'y';

   n1=TrimRow;
   n2=TrimCol;

   if ((centerflag=='y')||(centerflag=='Y')) /*move to center */
	 for(i=0;i<n1/2;i++)
	 for(j=0;j<n2/2;j++)
	 {
	   temcompl           = pv[i][j];
	   pv[i][j]           = pv[i+n1/2][j+n2/2];
	   pv[i+n1/2][j+n2/2] = temcompl;
	   temcompl           = pv[i][j+n2/2];
	   pv[i][j+n2/2]      = pv[i+n1/2][j];
	   pv[i+n1/2][j]      = temcompl;
	 }

/***************************IFFT************************/
   for(i=0;i<n1;i++)
   {
      fft(pv[i],n2,b);
      for(j=0;j<n2;j++)
      {
			pv[i][j].x = pv[i][j].x/sqrt((double)n2);
			pv[i][j].y = pv[i][j].y/sqrt((double)n2);
      }
   }

   for(j=0;j<n2;j++)
	 {
      for(i=0;i<n1;i++)
			  pu[i]=pv[i][j];
      fft(pu,n1,b);
      for(i=0;i<n1;i++)
      {
			  pv[i][j].x=pu[i].x/sqrt((double)n1);
			  pv[i][j].y=pu[i].y/sqrt((double)n1);
      }
   }

   for(i=0;i<row;i++)
   for(j=0;j<col;j++)
	   pv2[i][j] = pv[starty+i][startx+j].x+0.5;

	 ffree_2d((void **)pv,TrimRow);
	 free((void *)pu);
}

// flag=1: move to center
void FFT( unsigned char **pv1,unsigned char **pv2,
					int col,int row,char flag)
{
  int i,j;
  struct compl temcompl;
  struct compl **pv,*pu;
  int TrimRow, TrimCol;
  void fatal_error(char *msg);

  int MaxImgSize=(row>col)?row:col;
  int FFTDataSize=2;
  int PowPara=0;
	while(FFTDataSize<MaxImgSize){
		FFTDataSize=pow(2.0,PowPara);
	  PowPara++;
  }
  int startx,starty;
	startx=(FFTDataSize-col)/2;
  starty=(FFTDataSize-row)/2;

	char b = 't';           /* the flag of FFT */

  TrimRow=FFTDataSize;//row;
  TrimCol=FFTDataSize;//col;

  pv = (struct compl **)fspace_2d(TrimRow,TrimCol,sizeof(struct compl));
  pu = (struct compl *)calloc(TrimRow,sizeof(struct compl));

	for(i=0;i<row;i++)
  for(j=0;j<col;j++){
		if (flag==1 && (i+j)%2!=0)
		 pv[starty+i][startx+j].x = -pv1[i][j]; /* move to center point */
		else
		 pv[starty+i][startx+j].x = pv1[i][j];
	}

  for(i=0;i<TrimRow;i++)
  {
    fft(pv[i],TrimCol,b);
    for(j=0;j<TrimCol;j++)
    {
	    pv[i][j].x = pv[i][j].x/sqrt((double)TrimCol);
		  pv[i][j].y = pv[i][j].y/sqrt((double)TrimCol);
    }
  }
	for(j=0;j<TrimCol;j++)
  {
    for(i=0;i<TrimRow;i++)
      pu[i]=pv[i][j];
    fft(pu,TrimRow,b);
    for(i=0;i<TrimRow;i++)
    {
	    pv[i][j].x=pu[i].x/sqrt((double)TrimRow);
	    pv[i][j].y=pu[i].y/sqrt((double)TrimRow);
    }
  }

  float max=0;
	for(i=0;i<TrimRow;i++)
	for(j=0;j<TrimCol;j++){
		pv[i][j].x = log(10*sqrt(pv[i][j].x*pv[i][j].x+pv[i][j].y*pv[i][j].y)+1);
		if(pv[i][j].x>max)
			max=pv[i][j].x;
	}

	if(flag==1){
		for(i=0;i<row;i++)
		for(j=0;j<col;j++)
			pv2[i][j] = (pv[starty+i][startx+j].x/max)*255.0;
	}
	else{
	  for(i=0;i<=row/2;i++){
    for(j=0;j<=col/2;j++)
			pv2[i][j] = (pv[i][j].x/max)*255.0;
		for(j=col-1;j>col/2;j--)
			pv2[i][j] = (pv[i][FFTDataSize-col+j].x/max)*255.0;
  	}
	  for(i=row-1;i>row/2;i--){
		for(j=0;j<=col/2;j++)
			pv2[i][j] = (pv[FFTDataSize-row+i][j].x/max)*255.0;
		for(j=col-1;j>col/2;j--)
			pv2[i][j] = (pv[FFTDataSize-row+i][FFTDataSize-col+j].x/max)*255.0;
		}
  }

	ffree_2d((void **)pv,TrimRow);
	free((void *)pu);
}


void WFFT( unsigned char **pv1,unsigned char **pv2,
				  int col,int row)
{
  char centerflag,saveflag,b;
  int i,j,n1,n2,nn1,nn2,ii,jj;
  struct compl temcompl;
  struct compl **pv,*pu;
  double  dis,dis0;
  int TrimRow, TrimCol;
  void fatal_error(char *msg);

  int MaxImgSize=(row>col)?row:col;
  int FFTDataSize=2;
  int PowPara=0;
  while(FFTDataSize<MaxImgSize){
	 FFTDataSize=pow(2.0,PowPara);
	 PowPara++;
  }
  int startx,starty;
  startx=(FFTDataSize-col)/2;
  starty=(FFTDataSize-row)/2;

  b = 't';           /* the flag of FFT */
//  centerflag = 'y';  /* move to center point */
  centerflag = 'n';  /* move to center point */

  TrimRow=FFTDataSize;//row;
  TrimCol=FFTDataSize;//col;

  pv = (struct compl **)fspace_2d(TrimRow,TrimCol,sizeof(struct compl));
  pu = (struct compl *)calloc(TrimRow,sizeof(struct compl));

  for(i=0;i<TrimRow;i++)
  for(j=0;j<TrimCol;j++){
	  pv[i][j].x = 0;
	  pv[i][j].y = 0;
  }

  for(i=0;i<row;i++)
  for(j=0;j<col;j++){
	  if (((centerflag=='y')||(centerflag=='Y'))&&((i+j)%2!=0))
		 pv[starty+i][startx+j].x = -pv1[i][j]; /* move to center point */
	  else
		 pv[starty+i][startx+j].x = pv1[i][j];
  }

/******************************* FFT***************************/
  for(i=0;i<TrimRow;i++)
  {
    fft(pv[i],TrimCol,b);
	 for(j=0;j<TrimCol;j++)
    {
	    pv[i][j].x = pv[i][j].x/sqrt((double)TrimCol);
		  pv[i][j].y = pv[i][j].y/sqrt((double)TrimCol);
    }
  }