fig2.cpp

图像处理源代码VC

  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);
    }
  }

//-----------微分f(x,y)-f(x-1,y)<-->F(u,v)*(1-exp(-j*6.28u/N))---------

  for (i=0;i<TrimRow;i++) 
  for (j=0;j<TrimCol;j++){
	 pu[0].x = 1-cos(6.28315307*i/TrimCol);
	 pu[0].y = 1+sin(6.28315307*i/TrimRow);
	 comp_mult(pv[i][j],pu[0],&pv[i][j]);
  }

/********************************IFFT*********************************/
	//逆变换后即得 f(x,y)-f(x-1,y)
   b = 'f';
//   centerflag = 'y';
   centerflag = 'n';

   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;
	 }

   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+256)/2;	//整体提高亮度,以显示小于0的点
//	pv2[i][j] = pv[starty+i][startx+j].x;
/*
	float max=-100000,min=+100000;
	for(i=0;i<TrimRow;i++)
	for(j=0;j<TrimCol;j++){
		if(pv[i][j].x>max)
			max=pv[i][j].x;
		if(pv[i][j].x<min)
			min=pv[i][j].x;
	}

	for(i=0;i<row;i++)
	for(j=0;j<col;j++)
		pv2[i][j] = ((pv[starty+i][startx+j].x-min)/(max-min))*255.0;
*/

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


void JFFT( 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';

  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);
    }
  }

//-----------积分 f(x,y)+f(x+1,y)<-->F(u,v)*(1+exp(-j*6.28u/N))---------
  for (i=0;i<TrimRow;i++) 
  for (j=0;j<TrimCol;j++){
	 pu[0].x = 1+cos(6.28315307*i/TrimCol);
	 pu[0].y = 1-sin(6.28315307*i/TrimRow);
	 comp_mult(pv[i][j],pu[0],&pv[i][j]);
  }


/********************************IFFT*********************************/
	//逆变换后即得 f(x,y)+f(x-1,y)
   b = 'f';
//   centerflag = 'y';
   centerflag = 'n';

   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/2;

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

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

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


/*-------------------Otsu segment------------------------------------------*/
unsigned long binarize(unsigned char **ImgSrc, unsigned char **ImgDes,
											 unsigned ImgSizeX, unsigned ImgSizeY, unsigned thresh)
{
	int i,j;
  unsigned long onum;

  onum=0l;
	for (i=0;i<ImgSizeY;i++)
  {
		for (j=0;j<ImgSizeX;j++)
    {
			if (ImgSrc[i][j]>thresh)
      {
  	    onum++;
				ImgDes[i][j]=GL1;
      }
      else
			  ImgDes[i][j]=0;
    }
  }
  return(onum);
}
//       test_grey---->some grey level.
float BetweenClassVariance(float histogram[256],int test_grey)
{
	 int i;
	 float w0=0.0,w1=0.0;
	 float u0=0.0,u1,ut=0.0;
	 float sigmaB,sigmaT=0.0;

	 /*give class possibility*/
	 for(i=0;i<=test_grey;i++)
		 w0=w0+histogram[i];             //  C0类概率w0
	 w1=1.0-w0;                        //  C1类概率w1

	 if(w0==0.0)
	   return(0.0);
	 if(w0==1.0)
		 return(0.0);
	 
	 /*give means of two class*/
	 for(i=0;i<GL;i++)
		 ut=ut+i*histogram[i];         // 总体均值ut
	 for(i=0;i<=test_grey;i++)
		 u0=u0+i*histogram[i]/w0;		   // C0类均值u0
	 u1=(ut-u0)/(1.0-w0);            // C1类均值u1
	 sigmaB=w0*w1*(u1-u0)*(u1-u0);   // 类间方差
	 for(i=0;i<GL;i++)               // 总体方差
		 sigmaT=sigmaT+(i-ut)*(i-ut)*histogram[i];

	 return(sigmaB/sigmaT);        	 /*return ostu threshold*/
}
// 直方图
long *histo(unsigned char **ImgSrc,
						unsigned ImgSizeX, unsigned ImgSizeY )
{
  register int i,j;
  static long hist[GL];

  for (i=0;i<GL;i++)
    hist[i]=0L;
	for (i=0;i<ImgSizeY;i++)
  {
		for (j=0;j<ImgSizeX;j++)
      hist[ImgSrc[i][j]]++;
	}
  return(hist);
}
void OtsuSegment(unsigned char **ImgSrc, unsigned char **ImgDes,
								 unsigned ImgSizeX, unsigned ImgSizeY)
{
	long *histo(unsigned char **ImgSrc,unsigned ImgSizeX, unsigned ImgSizeY);
  long *hist;
	float HistPTR[GL1],b_variance[GL1];
	int thresh,i;
	float temp;
	float max_variance=0.0;

	hist= histo(ImgSrc,ImgSizeX,ImgSizeY);
	/* histogram normalize */
	temp=(float)ImgSizeX*ImgSizeY;
	for(i=0;i<GL;i++)
    HistPTR[i]=hist[i]/temp;

  /* call subprogram of BetweenClassVariance calculating */
	for(i=0;i<GL;i++)
	  b_variance[i]=BetweenClassVariance(HistPTR,i);

	/* give grey level of segment when max between_class variance */
	for(i=0;i<GL;i++)
	  if(b_variance[i]>=max_variance)
	  {
		  max_variance=b_variance[i];
			thresh=i;
	  }

	binarize(ImgSrc, ImgDes, ImgSizeX, ImgSizeY, thresh);
}
/*------------------------------------------------------------------------*/


void KJFFT( 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);
    }
  }
  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);
    }
  }

//---------空间域平移____F(u,v)*(-1)^(u+v)-->f(x-N/2,y-N/2)-----------
  for(i=0;i<TrimRow;i++)
  for(j=0;j<TrimCol;j++){
     pv[i][j].x *= pow(-1,j+i);
     pv[i][j].y *= pow(-1,j+i);

  }


/********************************IFFT*********************************/
	//逆变换后即得 f(x,y)-f(x-1,y)
   b = 'f';
   centerflag = 'y';
//   centerflag = 'n';

   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;
	 }

   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;
*/
    for(i=0;i<=row/2;i++){
    	for(j=0;j<=col/2;j++)
			pv2[i][j] = pv[i][j].x;
		for(j=col-1;j>=col/2;j--)
			pv2[i][j] = pv[i][FFTDataSize-col+j].x;
  	}
	for(i=row-1;i>=row/2;i--){
		for(j=0;j<=col/2;j++)
			pv2[i][j] = pv[FFTDataSize-row+i][j].x;
		for(j=col-1;j>=col/2;j--)
			pv2[i][j] = pv[FFTDataSize-row+i][FFTDataSize-col+j].x;
	}



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



void fatal_error(char* msg)
{
	AfxMessageBox(msg);
}