最大方差阈值分割（vc实现）

void CISLSView::OnThresholdOtsu()
 {
//程序编制：李立宗  lilizong@gmail.com
//2012-8-14
  if(myImage1.IsNull())
   OnOpenResourceFile();
  if(!myImage2.IsNull())
   myImage2.Destroy();
  if(myImage2.IsNull()){
   myImage2.Create(myImage1.GetWidth(),myImage1.GetHeight(),24,0);
  }
  //COLORREF pixel;
  int maxY = myImage1.GetHeight();
  int maxX=myImage1.GetWidth();
  byte* pRealData;
  byte* pRealData2;
  pRealData=(byte*)myImage1.GetBits();
  pRealData2=(byte*)myImage2.GetBits();
  int pit=myImage1.GetPitch();
  int pit2=myImage2.GetPitch();
  //需要注意,pit和pit2的值并不一样,所以如果使用一个值,会导致不同的结果出现
  //CString str;
  //str.Format(TEXT("%d"),pit);
  //MessageBox(str);
  //str.Format(TEXT("%d"),pit2);
  //MessageBox(str);
  int bitCount=myImage1.GetBPP()/8;
  int bitCount2=myImage2.GetBPP()/8;
  int tempR,tempG,tempB;
  //float temp,tempX,tempY;
  int temp;
  int T=128;
  //int pixel[4];
  float u0,u1;   //均值
  float w0,w1;   //概率
  float sum0,sum1;  //像素和
  int optIndex,optT;   //最优阈值，及其所在像素的值
  float fVaria,fMaxVaria=0;   //临时方差,最大方差
  int i;        //循环变量
  //int pixelR[256],pixelG[256],pixelB[256];
  int pixel[256]={0};   //不要忘记初始化
  for (int y=0; y<maxY-1; y++) {
   for (int x=0; x<maxX-1; x++) {
    temp=*(pRealData+pit*(y)+(x)*bitCount);
    if(bitCount==3)
    {
     tempR=*(pRealData+pit*(y)+(x)*bitCount);
     tempG=*(pRealData+pit*(y)+(x)*bitCount+1);
     tempB=*(pRealData+pit*(y)+(x)*bitCount+2);
     temp=(int)(tempR*0.49+tempG*0.31+tempB*0.2);
     //temp=(int)((tempR+tempG+tempB)/3);
    }
    
     *(pRealData2+pit2*(y)+(x)*bitCount2)=temp;
     *(pRealData2+pit2*(y)+(x)*bitCount2+1)=temp;
     *(pRealData2+pit2*(y)+(x)*bitCount2+2)=temp;
   }
  }
  for (int y=0; y<maxY; y++) {
   for (int x=0; x<maxX; x++) {
    temp=*(pRealData2+pit2*(y)+(x)*bitCount2);
    pixel[temp]++;

   }
  }
  //CString str;
  //str.Format(TEXT("%d"),pixel[3]);
  //MessageBox(str);

  for(optIndex=0;optIndex<256;optIndex++)
  {
   u0=0;
   sum0=0;
   for(i=0;i<=optIndex;i++)
   {
    u0+=i*pixel[i];
    sum0+=pixel[i];
   }
   u0/=sum0;
   w0=(float)(sum0)/(maxX*maxY);
   u1=0;
   sum1=0;
   for(i=optIndex;i<=255;i++)
   {
    u1+=i*pixel[i];
    sum1+=pixel[i];
   }
   u1/=sum1;
   w1=1-w0;
   fVaria=w0*w1*(u0-u1)*(u0-u1);
   if(fVaria>fMaxVaria)
   {
    fMaxVaria=fVaria;
    optT=optIndex;
   }
  }

  for (int y=0; y<maxY; y++) {
   for (int x=0; x<maxX; x++) {
    temp=*(pRealData2+pit2*(y)+(x)*bitCount2);
    if(temp>optT)
     tempR=255;
    else
     tempR=0;
     tempG=tempR;
     tempB=tempR;    
    *(pRealData2+pit2*y+x*bitCount2)=tempR;
    *(pRealData2+pit2*y+x*bitCount2+1)=tempG;
    *(pRealData2+pit2*y+x*bitCount2+2)=tempB;
   }
  }
  Invalidate();
 }