process.cpp

这是VC++ 2003.net图像处理的光盘源程序!!!非常好的

void KRP_op(BYTE **list0,BYTE **list1,int Dx,int Dy,int *pp,int bb)
{          //  带方向的边缘检测共用部分 (Kirsch，Robinson，Prewitt)
  BYTE t[8];
  int  i,j,k,p,g,m,tab[8][3][3];
   
  for (k=0;k<8;k++)
    for (i=0;i<3;i++)
	  for (j=0;j<3;j++)
		tab[k][i][j]=*pp++;

  for (i=1;i<Dy-1;i++) {
    for (j=1;j<Dx-1;j++) {
		for (p=0;p<8;p++) {	               
	    g = (tab[p][0][0]*list1[i-1][j-1]+tab[p][0][1]*list1[i-1][j]
		    +tab[p][0][2]*list1[i-1][j+1]+tab[p][1][0]*list1[i][j-1]
		    +tab[p][1][1]*list1[i][j]    +tab[p][1][2]*list1[i][j+1]
		    +tab[p][2][0]*list1[i+1][j-1]+tab[p][2][1]*list1[i+1][j]
		    +tab[p][2][2]*list1[i+1][j+1])/Scale;     
		if (g<0)  g=0;	                       
		else if (g>0xff)  g=0xff;
		t[p]= (BYTE) g;	       
	  }
	  m=0;
	  for (p=0;p<8;p++) {	   
	    if (t[p]>m)  m=t[p];
	  }
      if ((bb==0)||(m==t[bb-1]))
	    list0[i][j] = (BYTE) m;         
 	  else  list0[i][j] = 0;
	}
  }
}

void WINAPI Kirsch(BYTE **list0,BYTE **list1,int Dx,int Dy,int direction)
{                                      //  Kirsch 带方向的边缘检测
   Scale=8;                                        
   KRP_op(list0,list1,Dx,Dy,Tkirsch,direction);    
}

void WINAPI Robinson(BYTE **list0,BYTE **list1,int Dx,int Dy,int direction)
{                                    //  Robinson 带方向的边缘检测
   Scale=2;                                        
   KRP_op(list0,list1,Dx,Dy,Trobinson,direction);  
}

void WINAPI Prewitt(BYTE **list0,BYTE **list1,int Dx,int Dy,int direction)
{                                     //  Prewitt 带方向的边缘检测
   Scale=2;                                        
   KRP_op(list0,list1,Dx,Dy,Tprewitt,direction);   
}

void WINAPI TemplateB(BYTE **list0,BYTE **list1,int Dx,int Dy,int t)
{                                          	//  t x t 邻域处理共用程序 (t = 3 or 5)
	int   i,j,g,m,n,*p;
 
	for (i=t/2;i<Dy-t/2;i++) {
		for (j=t/2;j<Dx-t/2;j++) {
			g=0;
			p=&Mask5[0][0];
			for (m=i-t/2;m<=i+t/2;m++) {	
				for (n=j-t/2;n<=j+t/2;n++) 
					g+=(*p++)*list1[m][n];
 
            }
			g=g/Scale+Offset;	
			if (g>0xff)   g=0xff;	
			else if (g<0) g=0;
			list0[i][j] = g;           

		}
	}
}

void WINAPI LoG(BYTE **list0,BYTE **list1,int Dx,int Dy)
{                                                //  高斯型拉普拉斯算法
	int  mask[25]={ 0, 0,-1, 0, 0,
		            0,-1,-2,-1, 0,
		           -1,-2,16,-2,-1,
		            0,-1,-2,-1, 0,
				    0, 0,-1, 0, 0};

  Scale=1;		Offset=254;
  memcpy(Mask5,mask,25*sizeof(int));
  TemplateB(list0,list1,Dx,Dy,5);
  Offset=0;
}

void LineDetectA(BYTE **list0,BYTE **list1,int Dx,int Dy,int *pp,int bb)
{          //  带方向的边缘检测共用部分 (Kirsch，Robinson，Prewitt)
  BYTE t[8];
  int  i,j,k,p,g,m,tab[8][3][3];
   
  for (k=0;k<4;k++)
    for (i=0;i<3;i++)
	  for (j=0;j<3;j++)
		tab[k][i][j]=*pp++;

  for (i=1;i<Dy-1;i++) {
    for (j=1;j<Dx-1;j++) {
		for (p=0;p<4;p++) {	               
	    g = (tab[p][0][0]*list1[i-1][j-1]+tab[p][0][1]*list1[i-1][j]
		    +tab[p][0][2]*list1[i-1][j+1]+tab[p][1][0]*list1[i][j-1]
		    +tab[p][1][1]*list1[i][j]    +tab[p][1][2]*list1[i][j+1]
		    +tab[p][2][0]*list1[i+1][j-1]+tab[p][2][1]*list1[i+1][j]
		    +tab[p][2][2]*list1[i+1][j+1])/Scale;     
		if (g<0)  g=0;	                       
		else if (g>0xff)  g=0xff;
		t[p]= (BYTE) g;	       
	  }
	  m=0;
	  for (p=0;p<4;p++) {	   
	    if (t[p]>m)  m=t[p];
	  }
      if (bb==0)
	    list0[i][j] = (BYTE) m;         
	  else if (m==t[bb-1]) 
		list0[i][j] = (BYTE) m;         
	  else  list0[i][j] = 0;
	}
  }
}

void WINAPI LineDetect(BYTE **list0,BYTE **list1,int Dx,int Dy)
{                                      
	Scale=2;                                        
	LineDetectA(list0,list1,Dx,Dy,TLineDetect,0);   
}

void WINAPI SobelThinning(BYTE **list0,BYTE **list1,int Dx,int Dy) 	//  Sobel边缘细化
{
  int  i,j,g;

  Sobel(list0,list1,Dx,Dy); 	 
  for (i=1;i<Dy-1;i++) {
    for (j=1;j<Dx-1;j++) {
      g = list1[i][j]-list0[i][j]; 	 
      if ((list1[i][j-1]+list1[i][j+1]==0)|| 	 
          (list1[i-1][j]+list1[i+1][j]==0)) {	 
        list0[i][j]=list1[i][j];	 
      }
      else {
        if (g<0) list0[i][j]=0; 	 
        else     list0[i][j]=g; 	 
	  }
	}
  }
}

void WINAPI SobelThin(BYTE **list0,BYTE **list1,int Dx,int Dy) 	//  Sobel边缘细化
{
  BYTE *bufu,*bufv;
  BYTE **listu,**listv;
  int  i,j,g;
 
  bufu =(BYTE*)  malloc(Dx*Dy*sizeof(BYTE));
  for (i=0;i<Dy;i++) 
	listu=(BYTE**) malloc(Dy*sizeof(BYTE*));
  for(i=0;i<Dy;i++) listu[i]=bufu+i*Dx;

  bufv =(BYTE*)  malloc(Dx*Dy*sizeof(BYTE));
  for (i=0;i<Dy;i++) 
	listv=(BYTE**) malloc(Dy*sizeof(BYTE*));
  for(i=0;i<Dy;i++) listv[i]=bufv+i*Dx;

  SobelThinning(list0,list1,Dx,Dy);
  SobelThinning(listu,list0,Dx,Dy);
  SobelThinning(listv,listu,Dx,Dy);

  for (i=1;i<Dy-1;i++) {
	for (j=0;j<Dx-1;j++) {
	  g=list0[i][j];
	  g=(g+listu[i][j]+listv[i][j])/3;
	  listu[i][j]=g;
	}
  }
  OnePWidthW(list0,listu,Dx,Dy);

  free(bufu);  
  free(listu);  
  free(bufv);  
  free(listv);  
}

//--------------------------------------------------------------------------------

void HistogSmooth(long *pg,int n)
{
   int  i,j,k;

   for (j=0;j<n;j++) {
      k=pg[0];
	  for (i=1;i<255;i++) 
		 pg[i-1]=(pg[i-1]+pg[i]+pg[i+1])/3;
      for (i=254;i>0;i--) pg[i]=pg[i-1];
	  pg[0]=k;
   }
}

int  MaxMin(double *tab,int flag)
{
   double  max,min;
   int     i,p,q;

   max=min=tab[128];		
   p=q=128;
   for (i=0;i<256;i++) {
      if (tab[i]>max) {
	     max=tab[i];	
		 p=i;
      }
	  if (tab[i]<min) {
		 min=tab[i];	
		 q=i;
	  }
   }
   if (flag==1) {
      for (i=0;i<256;i++) 
         tab[i]=120*(max-tab[i])/(max-min);
   }
   else {         
      for (i=0;i<256;i++) 
         tab[i]=120*(tab[i]-min)/(max-min);
	  p=q;
   } 
   return(p);
}

int  ValleyAnalyse(long *pg)
{
   int  i,j,k,t,tab[10][2],n,min,g;

   HistogSmooth(pg,2);
   g=2;
   min=pg[0];      k=0;
   for (i=1;i<256;i++) {
	  if (pg[i]<min) {
	     k=i;
		 min=pg[i];
	  }
   }
   t=i=j=k;
   while((pg[i]<=g)&&(i<255)) i++;
   while((pg[j]<=g)&&(j>0))   j--;
   if ((i<250)||(j>5)) k=(i+j)/2;
   tab[1][0]=j;		tab[1][1]=i;

   tab[2][0]=0;     tab[2][1]=255;
   n=i+10;
   if (n>220) n=220;
   min=pg[n];      k=n;
   for (i=n+1;i<256;i++) {
	  if (pg[i]<min) {
	     k=i;
		 min=pg[i];
	  }
   }
   i=j=k;
   while((pg[i]-min<=g)&&(i<255)) i++;
   while((pg[j]-min<=g)&&(j>0))   j--;
   if ((i<250)||(j>5)) k=(i+j)/2;
    if ((i<224)||(j>5)) {
	  k=(i+j)/2;
      tab[2][0]=j;		tab[2][1]=i;
   }

   tab[0][0]=0;         tab[0][1]=255;
   n=t-10;
   if (n<20) n=20;
   min=pg[n];     k=n;
   for (i=n-1;i>0;i--) {
	  if (pg[i]<min) {
	     k=i;
		 min=pg[i];
	  }
   }
   i=j=k;
   while((pg[i]-min<=g)&&(i<255)) i++;
   while((pg[j]-min<=g)&&(j>0))   j--;
   if ((i<250)||(j>5)) k=(i+j)/2;
    if ((i<224)||(j>5)) {
	  k=(i+j)/2;
      tab[0][0]=j;		tab[0][1]=i;
   }

   if (tab[0][1]-tab[0][0]>80) tab[0][1]=255;
   if (tab[2][1]-tab[2][0]>80) tab[2][0]=0;
   if (tab[0][0]==0)   tab[0][1]=255;
   if (tab[2][1]==255) tab[2][0]=0;
   n=1;
   if (tab[0][1]<tab[1][0]) 
      if (tab[0][1]-tab[0][0]>tab[1][1]-tab[1][0]) n=0; 
   if (tab[2][0]>tab[1][1]) 
      if (tab[2][1]-tab[2][0]>tab[1][1]-tab[1][0]) n=2; 
   k=(tab[n][1]+tab[n][0])/2;
   return(k);
}

int  WINAPI Otsu(long *pg)      //  大津法取阈值
{
   int  i,j,p;
   double m0,m1,M0,M1,u,v,w[256];

   M0=M1=0;
   for (i=0;i<256;i++) {
      M0+=pg[i];     M1+=pg[i]*i;
   }
   for (j=0;j<256;j++) {
      m0=m1=0;
      for (i=0;i<=j;i++) {
	 m0+=pg[i];	 m1+=pg[i]*i;
      }
      if (m0) u=m1/m0;
      else   u=0;
      if (M0-m0) v=(M1-m1)/(M0-m0);
      else      v=0;
      w[j]=m0*(M0-m0)*(u-v)*(u-v);
   }

   p=MaxMin(w,1);
   return(p);
}

int  WINAPI Minimum(long *pg) //  最小误差法取阈值
{
   double u1,u2,tab[256];
   double s1,s2,P1,P2,PP1,PP2,SS1,SS2;
   int    i,g,p;

   for (i=0;i<256;i++) {
	   for (g=0,  P1=0;g<=i; g++) 
		   P1 += (double) pg[g];
	   for (g=i+1,P2=0;g<256;g++) 
		   P2 += (double) pg[g];
	   for (g=0,  u1=0;g<=i; g++) 
		   u1 += (double) pg[g]*g;
	   for (g=i+1,u2=0;g<256;g++) 
		   u2 += (double) pg[g]*g;
       u1 /= P1;		
	   u2 /= P2;
	   for (g=0,  s1=0;g<=i; g++) 
		   s1 += (double) pg[g]*(g-u1)*(g-u1);
	   for (g=i+1,s2=0;g<256;g++)
		   s2 += (double) pg[g]*(g-u2)*(g-u2);
	   s1 = sqrt(s1/P1);	
	   s2 = sqrt(s2/P2);
	   SS1=SS2=PP1=PP2=0;
       if (s1) SS1=P1*log(s1);
       if (s2) SS2=P2*log(s2);
       if (P1) PP1=P1*log(P1);
       if (P2) PP2=P2*log(P2);
       tab[i]=1+2*(SS1+SS2)-2*(PP1+PP2);
   }

   p=MaxMin(tab,0);
   return(p);
}

int  WINAPI Ptile(long *pg,double nn)       //  分位数法
{
   int  i;
   double mm,kk;

   for (i=0,mm=0;i<256;i++) mm+=pg[i];
   kk=0;
   for (i=0;i<256;i++) {
	  kk+=(double) pg[i]/mm;
      if (kk>=nn)  break;
   }
   return i;
}

int  WINAPI KSW_Entropic(long *pg)            //  熵法取阈值
{
   long   i,t,s;
   double p[256],Pt[256],Ht[256],HT,H[256];
   double A,B,C;
  
   for (i=0,s=0;i<256;i++) s+=pg[i];
   for (i=0;i<256;i++) p[i]=((double) pg[i])/s;

   Pt[0]=p[0];
   for (i=1;i<256;i++) Pt[i]=Pt[i-1]+p[i];

   for (i=0;i<256;i++) {
      if (p[i]==0) Ht[i]=0;
	  else  Ht[i]=-p[i]*log(p[i]);
   }
   for (i=1;i<256;i++) Ht[i]+=Ht[i-1];
   HT=Ht[255];

   for (i=0;i<256;i++) {
	  A=Pt[i]*(1-Pt[i]);
      if (A>0) A=log(A);
	  B=Pt[i];
	  if (B>0) B=Ht[i]/B;
	  C=1-Pt[i];
	  if (C>0) C=(HT-Ht[i])/C;
	  H[i]=A+B+C;
   }

   t=MaxMin(H,1);
   return(t);
}

int  WINAPI Moment(long *pg)                //  矩法取阈值
{
   long   i,t;
   double m0,m1,m2,m3,p0,C0,C1;

   m0=m1=m2=m3=0;
   for (i=0;i<256;i++) {
      m0+=(double) pg[i];		
	  m1+=(double) pg[i]*i;
      m2+=(double) pg[i]*i*i;
	  m3+=(double) pg[i]*i*i*i;
   }
   C0=(m1*m3-m2*m2)/(m0*m2-m1*m1);
   C1=(m1*m2-m0*m3)/(m0*m2-m1*m1);
   p0=0.5-(m1/m0+C1/2)/sqrt(C1*C1-4*C0);
   t=Ptile(pg,p0);
   return(t);
}

int  WINAPI BiasNormal(long *pg)
{
   int  i,j,p;
   double u,u1,u2,t,t1,t2,tab[256];
   double k1,k2;

   t=u=0;
   for (i=0;i<256;i++) {
      t+=pg[i];		u+=i*pg[i];
   }
   for (j=0;j<256;j++) {
      t1=u1=0;
      for (i=0;i<=j;i++) {
	     t1+=pg[i];		u1+=i*pg[i];
      }
      t2=t-t1;
      u2=u-u1;
      if (t1)	u1=u1/t1;
      else		u1=0;
      if (t2)	u2=u2/t2;
      else		u2=0;

	  for (i=0,k1=0;i<j;i++)   
		  k1+=(i-u1)*(i-u1)*(i-u1)*pg[i];
      for (i=j,k2=0;i<256;i++) 
		  k2+=(i-u2)*(i-u2)*(i-u2)*pg[i];
	  k1=k1*t1;
	  k2=k2*t2;

	  tab[j]=fabs(k1)+fabs(k2);
   }

   p=MaxMin(tab,0);
   return(p);
}
int  WINAPI ColligationThreshold(long *pg)
{
	int		thre,list[256],buf[20],bn;
	int		i,j,k,m,n,tab[10];

	thre=FormAnalysis(pg,list,buf,&bn);
	if (bn<4) return(thre);

	tab[0]=Otsu(pg);
	tab[1]=KSW_Entropic(pg);
	tab[2]=Minimum(pg);
	tab[3]=BiasNormal(pg);
	tab[4]=Moment(pg);

	m=5;
	for (i=0;i<m-1;i++) {                  
       k=i;      
	   n=tab[i];
       for (j=i+1;j<m;j++) {
	      if (tab[j]<n) {
	          k=j;	
			  n=tab[j];
		  }
       }
       tab[k]=tab[i];
       tab[i]=n;      
	}
	k=tab[m/2];
	if ((bn>3)&&(bn<6)&&(buf[3]-buf[1]>100)) {
		for (i=2;i<bn;i+=2) {
			if (abs(buf[i]-k)<40) k=buf[i];
		}
	}
	return(k);
}

int  WINAPI SimpleThreshold(long *pg)
{
	if (pg[254]>pg[1]) return(255);
	else return(1);
}

void HistogramSmooth(long *pg,int n)
{
	TableSmooth(pg,256,n);
}

int  GrayImageType(long *pg)
{
	int i,k,t;
	long pga[256],max,a,b;

	memcpy(pga,pg,256*sizeof(long));
	pga[0]=pga[255]=0;
	HistogSmooth(pga,15);
	max=pga[0];
	k=0;
	for (i=0;i<256;i++) {
		if (pga[i]>max) {
			max=pga[i];
			k=i;
		}
	}
	a=b=0;
	for (i=0;i<k;i++) a+=pg[i];
	for (i=k+1;i<256;i++) b+=pg[i];
	if (a>b) t=0;
	else t=1;
	return(t);
}

int  WINAPI FormAnalysis(long *pg,int *tab,int *buf,int *bn)
{
	int		num=20,h=120,a=10,b=5,g=20; 
 	long	pga[256];
	int		i,j,m,n,k,f,v,s,t;
	int		thre,buf2[20];
 
 	memcpy(pga,pg,256*sizeof(long));
	pga[0]=pga[255]=0;
	HistogSmooth(pga,num);

	m=SearchScale(pga,5);
	k=m/h+1;
	for(i=0;i<256;i++) 
	{
 		n=pga[i]/k;
		if (n>h) n=h;
		tab[i]=n;
 	}

	m=n=0;
	buf[0]=0;
	k=tab[0];
	f=j=1;
	for (i=0;i<255;i++) {
		if (tab[i]==k) n++;
		else {
    		if (f==1) {
				if (tab[i]<k) {
					buf[j++]=m+n/2;
					f=0;
				}
			}
			else {
				if (tab[i]>k) {
					buf[j++]=m+n/2;
					f=1;
				}
			}
			m=i;
			n=1;
			k=tab[i];