最大方差阈值分割（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();

 }