精华内容
下载资源
问答
  • OTSU 算法

    2020-03-26 08:49:47
    处理双峰图像,Otsu算法试图寻找阈值(t)使类内方差最小化 其中: 变量解释:t表示阈值,P(i)表示未归一化的直方图,灰度值i出现的概率。 实例说明可参考:...

    OTSU:

    处理双峰图像,Otsu算法试图寻找阈值(t)使类内方差最小化

    其中:

    变量解释:t表示阈值,P(i)表示未归一化的直方图,灰度值i出现的概率。

    实例说明可参考:http://www.labbookpages.co.uk/software/imgProc/otsuThreshold.html

     

     

    展开全文
  • Otsu算法

    千次阅读 2018-11-06 15:42:28
    版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/kksc1099054857/article/details/78321776 1....nbsp;   &...一维Otsu算法有计算简...
    版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/kksc1099054857/article/details/78321776

    1.简介:

           一维Otsu算法有计算简洁、稳定、自适应强等优点,被广泛用于图像分割中。但一维Otsu算法没有考虑图像像素点之间的关系,当图像中有噪声时,会导致分割的效果不理想。因此,刘健庄等人在1993年提出了二维的Otsu算法,提升了算法的抗噪声能力。

    2.算法思想:

           同时考虑像素的灰度值分布和它们邻域像素的平均灰度值分布,因此形成的阈值是一个二维矢量,最佳的阈值在一个二维的测度准则下确定最大值时得到。

    3.算法过程:

    (1)设图像I(x,y),的灰度级为L级,那么图像的邻域平均灰度也分为L级。
    (2)设f(x,y)为像素点(x,y)的灰度值,g(x,y)为像素点(x,y)为中心的K*K的像
            素点集合的灰度平均值。令f(x,y)=i,g(x,y)=j,然后就形成了一个二元组(i,j)。
    (3)设二元组(i,j)出现的次数为fij,然后求出二元组对应的概率密度Pij,
            Pij=fij/N, i,j=1,2,…,L,其中N为图像像素点总数。
    (4)任意选取一个阈值向量(s,t)选取的阈值向量将图像的二维直方图划分成4个
            区域,B、C区域代表图像的前景和背景,A、D区域代表噪声点。

    (5)设C、B两个区域对应的概率分别为w1,w2,对应的均值矢量为u1,u2。整个图
            片所对应的均值矢量为uT。


    4.代码实现(opencv3):

        

    1. #include "opencv.hpp"
    2. #include "imgproc.hpp"
    3. #include "highgui.hpp"
    4. #include "iostream"
    5. #include "core.hpp"
    6. using namespace cv;
    7. using namespace std;
    8. int Otsu2D(Mat srcimage); //二维Otsu算法
    9. int main()
    10. {
    11. Mat srcimage, grayimage, dstimage;
    12. srcimage = imread("lena.jpg");
    13. namedWindow("原图", 0);
    14. imshow("原图", srcimage); //显示原图
    15. cvtColor(srcimage, grayimage, COLOR_RGB2GRAY); //得到灰度图
    16. double time0 = static_cast<double>(getTickCount()); //记录程序开始时间
    17. int thresholdValue = Otsu2D(grayimage); //调用二维Otsu函数
    18. time0 = ((double)getTickCount() - time0) / getTickFrequency();
    19. cout << "算法运行时间为:" << time0 << endl;
    20. cout << "Otsu阈值为:" << thresholdValue << endl;
    21. threshold(grayimage, dstimage, thresholdValue, 255, THRESH_BINARY); //将得到的阈值传入函数,得到分割效果图
    22. namedWindow("Otsu算法结果", 0);
    23. imshow("Otsu算法结果", dstimage);
    24. waitKey();
    25. return 0;
    26. }
    1. <code class="language-cpp">int Otsu2D(Mat srcimage)  
    2. {  
    3.     double Histogram[256][256];        //建立二维灰度直方图  
    4.     double TrMax = 0.0;                //用于存储矩阵的迹(矩阵对角线之和)  
    5.     int height = srcimage.rows;        //矩阵的行数  
    6.     int width = srcimage.cols;         //矩阵的列数  
    7.     int N = height*width;              //像素的总数  
    8.     int T;                             //最终阈值  
    9.     uchar *data = srcimage.data;  
    10.     for (int i = 0; i < 256; i++)  
    11.     {  
    12.         for (int j = 0; j < 256; j++)  
    13.         {  
    14.             Histogram[i][j] = 0;      //初始化变量  
    15.         }  
    16.     }  
    17.     for (int i = 0; i < height; i++)  
    18.     {  
    19.         for (int j = 0; j < width; j++)  
    20.         {  
    21.             int Data1 = data[i*srcimage.step + j];         //获取当前灰度值  
    22.             int Data2 = 0;                           //用于存放灰度的平均值  
    23.             for (int m = i - 1; m <= i + 1; m++)  
    24.             {  
    25.                 for (int n = j - 1; n <= j + 1; n++)  
    26.                 {  
    27.                     if ((m >= 0) && (m < height) && (n >= 0) && (n < width))  
    28.                         Data2 += data[m*srcimage.step + n];//邻域灰度值总和  
    29.                 }  
    30.             }  
    31.             Data2 = Data2 / 9;  
    32.             Histogram[Data1][Data2]++;                  //记录(i,j)的数量  
    33.         }  
    34.     }  
    35.     for (int i = 0; i < 256; i++)  
    36.         for (int j = 0; j < 256; j++)  
    37.             Histogram[i][j] /= N;     //归一化的每一个二元组的概率分布  
    38.   
    39.     double Fgi = 0.0;    //前景区域均值向量i分量  
    40.     double Fgj = 0.0;    //前景区域均值向量j分量  
    41.     double Bgi = 0.0;    //背景区域均值向量i分量  
    42.     double Bgj = 0.0;    //背景区域均值向量j分量  
    43.     double Pai = 0.0;    //全局均值向量i分量 panorama(全景)  
    44.     double Paj = 0.0;    //全局均值向量j分量  
    45.     double w0 = 0.0;     //前景区域联合概率密度  
    46.     double w1 = 0.0;     //背景区域联合概率密度  
    47.     double num1 = 0.0;   //遍历过程中前景区i分量的值  
    48.     double num2 = 0.0;   //遍历过程中前景区j分量的值  
    49.     double num3 = 0.0;   //遍历过程中背景区i分量的值  
    50.     double num4 = 0.0;   //遍历过程中背景区j分量的值  
    51.     int Threshold_s = 0; //阈值s  
    52.     int Threshold_t = 0; //阈值t  
    53.     double temp = 0.0;   //存储矩阵迹的最大值  
    54.     for(int i=0;i<256;i++)  
    55.     {  
    56.         for (int j = 0; j < 256; j++)  
    57.         {  
    58.             Pai += i*Histogram[i][j];   //全局均值向量i分量计算  
    59.             Paj += j*Histogram[i][j];   //全局均值向量j分量计算  
    60.         }  
    61.     }  
    62.     for (int i = 0; i < 256; i++)  
    63.     {  
    64.         for (int j = 0; j < 256; j++)  
    65.         {  
    66.             w0 += Histogram[i][j];        //前景的概率  
    67.             num1 += i*Histogram[i][j];    //遍历过程中前景区i分量的值  
    68.             num2 += j*Histogram[i][j];    //遍历过程中前景区j分量的值  
    69.   
    70.             w1 = 1 - w0;                  //背景的概率  
    71.             num3 = Pai - num1;            //遍历过程中背景区i分量的值  
    72.             num4 = Paj - num2;            //遍历过程中背景区j分量的值  
    73.   
    74.             Fgi = num1 / w0;                
    75.             Fgj = num2 / w1;  
    76.             Bgi = num3 / w0;  
    77.             Bgj = num4 / w1;  
    78.             TrMax = ((w0*Pai - num1)*(w0*Pai - num1) + (w0*Paj - num2)*(w0*Paj - num2)) / (w0*w1);  
    79.             if (TrMax > temp)  
    80.             {  
    81.                 temp = TrMax;  
    82.                 Threshold_s = i;  
    83.                 Threshold_t = j;  
    84.             }  
    85.         }  
    86.     }  
    87.     cout << Threshold_s << " " << Threshold_t << endl;  
    88.     T = Threshold_s;  
    89.     return T;  
    90. }  
    91. </code>  
    1. int Otsu2D(Mat srcimage)
    2. {
    3. double Histogram[256][256]; //建立二维灰度直方图
    4. double TrMax = 0.0; //用于存储矩阵的迹(矩阵对角线之和)
    5. int height = srcimage.rows; //矩阵的行数
    6. int width = srcimage.cols; //矩阵的列数
    7. int N = height*width; //像素的总数
    8. int T; //最终阈值
    9. uchar *data = srcimage.data;
    10. for (int i = 0; i < 256; i++)
    11. {
    12. for (int j = 0; j < 256; j++)
    13. {
    14. Histogram[i][j] = 0; //初始化变量
    15. }
    16. }
    17. for (int i = 0; i < height; i++)
    18. {
    19. for (int j = 0; j < width; j++)
    20. {
    21. int Data1 = data[i*srcimage.step + j]; //获取当前灰度值
    22. int Data2 = 0; //用于存放灰度的平均值
    23. for (int m = i - 1; m <= i + 1; m++)
    24. {
    25. for (int n = j - 1; n <= j + 1; n++)
    26. {
    27. if ((m >= 0) && (m < height) && (n >= 0) && (n < width))
    28. Data2 += data[m*srcimage.step + n];//邻域灰度值总和
    29. }
    30. }
    31. Data2 = Data2 / 9;
    32. Histogram[Data1][Data2]++; //记录(i,j)的数量
    33. }
    34. }
    35. for (int i = 0; i < 256; i++)
    36. for (int j = 0; j < 256; j++)
    37. Histogram[i][j] /= N; //归一化的每一个二元组的概率分布
    38. double Fgi = 0.0; //前景区域均值向量i分量
    39. double Fgj = 0.0; //前景区域均值向量j分量
    40. double Bgi = 0.0; //背景区域均值向量i分量
    41. double Bgj = 0.0; //背景区域均值向量j分量
    42. double Pai = 0.0; //全局均值向量i分量 panorama(全景)
    43. double Paj = 0.0; //全局均值向量j分量
    44. double w0 = 0.0; //前景区域联合概率密度
    45. double w1 = 0.0; //背景区域联合概率密度
    46. double num1 = 0.0; //遍历过程中前景区i分量的值
    47. double num2 = 0.0; //遍历过程中前景区j分量的值
    48. double num3 = 0.0; //遍历过程中背景区i分量的值
    49. double num4 = 0.0; //遍历过程中背景区j分量的值
    50. int Threshold_s = 0; //阈值s
    51. int Threshold_t = 0; //阈值t
    52. double temp = 0.0; //存储矩阵迹的最大值
    53. for(int i=0;i<256;i++)
    54. {
    55. for (int j = 0; j < 256; j++)
    56. {
    57. Pai += i*Histogram[i][j]; //全局均值向量i分量计算
    58. Paj += j*Histogram[i][j]; //全局均值向量j分量计算
    59. }
    60. }
    61. for (int i = 0; i < 256; i++)
    62. {
    63. for (int j = 0; j < 256; j++)
    64. {
    65. w0 += Histogram[i][j]; //前景的概率
    66. num1 += i*Histogram[i][j]; //遍历过程中前景区i分量的值
    67. num2 += j*Histogram[i][j]; //遍历过程中前景区j分量的值
    68. w1 = 1 - w0; //背景的概率
    69. num3 = Pai - num1; //遍历过程中背景区i分量的值
    70. num4 = Paj - num2; //遍历过程中背景区j分量的值
    71. Fgi = num1 / w0;
    72. Fgj = num2 / w1;
    73. Bgi = num3 / w0;
    74. Bgj = num4 / w1;
    75. TrMax = ((w0*Pai - num1)*(w0*Pai - num1) + (w0*Paj - num2)*(w0*Paj - num2)) / (w0*w1);
    76. if (TrMax > temp)
    77. {
    78. temp = TrMax;
    79. Threshold_s = i;
    80. Threshold_t = j;
    81. }
    82. }
    83. }
    84. cout << Threshold_s << " " << Threshold_t << endl;
    85. T = Threshold_s;
    86. return T;
    87. }

    5.运行结果(得到2个阈值,取其中一个即可):


    6.修改后的代码

    1. int Otsu2D(Mat srcimage)
    2. {
    3. double Histogram[256][256]; //建立二维灰度直方图
    4. double TrMax = 0.0; //用于存储矩阵的迹(矩阵对角线之和)
    5. int height = srcimage.rows; //矩阵的行数
    6. int width = srcimage.cols; //矩阵的列数
    7. int N = height*width; //像素的总数
    8. int T; //最终阈值
    9. uchar *data = srcimage.data;
    10. for (int i = 0; i < 256; i++)
    11. {
    12. for (int j = 0; j < 256; j++)
    13. {
    14. Histogram[i][j] = 0; //初始化变量
    15. }
    16. }
    17. for (int i = 0; i < height; i++)
    18. {
    19. for (int j = 0; j < width; j++)
    20. {
    21. int Data1 = data[i*srcimage.step + j]; //获取当前灰度值
    22. int Data2 = 0; //用于存放灰度的平均值
    23. for (int m = i - 1; m <= i + 1; m++)
    24. {
    25. for (int n = j - 1; n <= j + 1; n++)
    26. {
    27. if ((m >= 0) && (m < height) && (n >= 0) && (n < width))
    28. Data2 += data[m*srcimage.step + n];//邻域灰度值总和
    29. }
    30. }
    31. Data2 = Data2 / 9;
    32. Histogram[Data1][Data2]++; //记录(i,j)的数量
    33. }
    34. }
    35. for (int i = 0; i < 256; i++)
    36. {
    37. for (int j = 0; j < 256; j++)
    38. {
    39. Histogram[i][j] /= N; //归一化的每一个二元组的概率分布
    40. }
    41. }
    42. double S[256]; //统计前i行概率的数组
    43. double N1[256]; //统计遍历过程中前景区i分量的值
    44. double N2[256]; //统计遍历过程中前景区j分量的值
    45. S[0] = 0;
    46. N1[0] = 0;
    47. N2[0] = 0;
    48. for (int i = 1; i < 256; i++)
    49. {
    50. double x = 0,n1=0,n2=0;
    51. for (int j = 0; j < 256; j++)
    52. {
    53. x += Histogram[i - 1][j];
    54. n1 += ((i-1)*Histogram[i - 1][j]); //遍历过程中前景区i分量的值
    55. n2 += (j*Histogram[i - 1][j]); //遍历过程中前景区j分量的值
    56. }
    57. S[i] = x + S[i - 1];
    58. N1[i] = n1 + N1[i - 1];
    59. N2[i] = n2 + N2[i - 1];
    60. }
    61. double Pai = 0.0; //全局均值向量i分量 panorama(全景)
    62. double Paj = 0.0; //全局均值向量j分量
    63. int Threshold_s = 0; //阈值s
    64. int Threshold_t = 0; //阈值t
    65. int M = 0; //中间变量
    66. double temp = 0.0; //存储矩阵迹的最大值
    67. double num3 = 0.0; //遍历过程中背景区i分量的值
    68. double num4 = 0.0; //遍历过程中背景区j分量的值
    69. double Fgi = 0.0; //前景区域均值向量i分量
    70. double Fgj = 0.0; //前景区域均值向量j分量
    71. double Bgi = 0.0; //背景区域均值向量i分量
    72. double Bgj = 0.0; //背景区域均值向量j分量
    73. for(int i=0;i<256;i++)
    74. {
    75. for (int j = 0; j < 256; j++)
    76. {
    77. Pai += i*Histogram[i][j]; //全局均值向量i分量计算
    78. Paj += j*Histogram[i][j]; //全局均值向量j分量计算
    79. }
    80. }
    81. for (int i = 0; i < 256; i++)
    82. {
    83. double w0 = 0.0; //前景区域联合概率密度
    84. double w1 = 0.0; //背景区域联合概率密度
    85. double num1 = 0.0; //遍历过程中前景区i分量的值 与w0一样做相关处理
    86. double num2 = 0.0; //遍历过程中前景区j分量的值
    87. if (i >= 1)
    88. {
    89. w0 += S[i - 1];
    90. num1 += N1[i - 1];
    91. num2 += N2[i - 1];
    92. }
    93. for (int j = 0; j < 256; j++)
    94. {
    95. w0 += Histogram[i][j]; //前景的概率
    96. num1 += i*Histogram[i][j]; //遍历过程中前景区i分量的值
    97. num2 += j*Histogram[i][j]; //遍历过程中前景区j分量的值
    98. w1 = 1 - w0; //背景的概率
    99. num3 = Pai - num1; //遍历过程中背景区i分量的值
    100. num4 = Paj - num2; //遍历过程中背景区j分量的值
    101. }
    102. Fgi = num1 / w0;
    103. Fgj = num2 / w1;
    104. Bgi = num3 / w0;
    105. Bgj = num4 / w1;
    106. TrMax = ((w0*Pai - num1)*(w0*Pai - num1) + (w0*Paj - num2)*(w0*Paj - num2)) / (w0*w1);
    107. if (TrMax > temp)
    108. {
    109. temp = TrMax;
    110. Threshold_s = i;
    111. }
    112. }
    113. cout << Threshold_s <<endl;
    114. T = Threshold_s;
    115. return T;
    116. }

    展开全文
  • OTSU算法

    2015-11-09 22:35:58
    OTSU算法 OTSU算法以最佳门限将图像灰度直方图分割成两部分,使两部分类间方差取最大值,即分离性最大。 设图像灰度级 ,第 级象素 个,总象素,则第级灰度出现的概率为 。 设灰度门限值为 ,则图像像素按灰度级...

    OTSU算法

    OTSU算法以最佳门限将图像灰度直方图分割成两部分,使两部分类间方差取最大值,即分离性最大。

    设图像灰度级 ,第 级象素 个,总象素,则第级灰度出现的概率为

    设灰度门限值为 ,则图像像素按灰度级被分为两类:

    图像总平均灰度级:

    类的平均灰度级为: ,像素数为:

    类的平均灰度级为:   , 像素数为:

    两部分图像所占比例分别为:

          

    均值作处理:

         

    图像总均值可化为:  

    类间方差: 

    化为:     

    变化,使 最大的 即为所求之最佳门限。

    称为目标选择函数。

     

     

     

     

     

    My code:

           int i;          //图像灰度级 from 0--255

           double n[256];     //i级像素n[i]

           double Pro[256];   //i级灰度出现的概率

           double Tavegray = 0;   //图像总平均灰度级

           double avegray0[256];   //C0类的平均灰度级

           double avegray1[256];   //C1类的平均灰度级

           int N = 0;          //总像素

           int N0 = 0;         //C0类的像素数

           int N1 = 0;         //C1类的像素数

           double w0 = 0, w1 =0;     //分别为两部分图像所占比例

           double mean0 = 0, mean1 =0, mean = 0;  //C0C1 均值作处理,及图像总均值

           double var[256];        //类间方差

           int k = 0;      //灰度门限值

           double w[256];     //w[k]当门限值为k时,C0类图像所占的比值

           int Threshold;  //门限值

          

           //数组初始化

           for( i = 0; i <256; i++ )

           {

                  n[i] = 0;

                  Pro[i] = 0;

                  avegray0[i] = 0;

                  avegray1[i] = 0;

                  var[i] = 0;

                  w[i] = 0;

           }

           int I;

          

           LPSTRlpDibHdr = (LPSTR)::GlobalLock(hDib);

           LPBITMAPINFOHEADERlpbmi = (LPBITMAPINFOHEADER)lpDibHdr; //得到指向位图的信息头结构的指针

           LPSTRlpDibBits = lpDibHdr - lpbmi->biSize; //得到指向位图数据的指针

     

           for( int y = 0; y <lpbmi->biHeight; y++ )

           {

                  for( int x = 0; x < lpbmi->biWidth; x++ )

                  {

                         COLORREFc;

                         BYTE*pDate = (BYTE*)lpDibBits + (lpbmi->biHeight-y) * getpiont( hDib )+ ( x ) *3; //((lpbmi->biWidth*24+31)/32*4)

                         c= RGB( pDate[2], pDate[1], pDate[0]);

                  //     COLORREF cNew;

                  //     BYTE* pDateNew = (BYTE*)lpNewDibBits +(lpNewbmi->biHeight -y) * getpiont(hDib) + ( x ) * 3; //((lpbmi->biWidth*24+31)/32*4)

           //   cNew = RGB( pDateNew[2], pDateNew[1], pDateNew[0]);

                         I= ( pDate[2] + pDate[1] + pDate[0] ) / 3;

                         n[I]++;

                  }

           }

           //总象素

           for( i = 0; i <256; i++ )

           {

                  N += n[i];

           }

           //i级灰度出现的概率

           for( i = 0; i <256; i++ )

           {

                  Pro[i] = ( n[i] / N );

           }

           //图像总平均灰度级

           for( i = 1; i <=256; i++ )

           {

                  Tavegray += ( i * Pro[i - 1] );

           }

          

      for ( k = 0; k< 256; k++ )

      {

           for( i = 1; i < k+1; i++ )

           {

                  avegray0[k] += i*Pro[i];//C0类的平均灰度级

           //     N0 +=n[i];            //C0类像素数

                  w0 += Pro[i];

           }

           w[k]= w0;//C0部分图像所占比例

          

    //     avegray1[k]= Tavegray - avegray0[k];//C1类的平均灰度级和像素数

    //     N1= N - N0;                   //C1类像素数

           w1= 1 - w[k];    //C0部分图像所占比例

           if ( w[k] != 0 )  //C0均值作处理

           {

           mean0 = avegray0[k]/ w[k]; 

           }

           else mean0 = 0;

               

           mean1= ( Tavegray - avegray0[k] ) / ( 1 - w[k] );//C1均值作处理

           mean= ( w0 * mean0 ) + ( w1 * mean1 );

           var[k]= (( ( mean * w[k] ) - avegray0[k]) * ( ( mean * w[k] ) - avegray0[k])) / (w[k] * ( 1 - w[k]));

          

    }    

     

      int temp = 0;

      for( k = 0; k< 256; k++ )

      {

           if ( var[k] >temp)

           {

                  temp = var[k];

                  Threshold = k;

           }

    }

          

       

     

     

     

    OTSU 算法可以说是自适应计算单阈值(用来转换灰度图像为二值图像)的简单高效方法。下面的代码最早由Ryan Dibble提供,此后经过多人Joerg.Schulenburg,R.Z.Liu 等修改,补正。
     
    算法对输入的灰度图像的直方图进行分析,将直方图分成两个部分,使得两部分之间的距离最大。划分点就是求得的阈值。
      
     parameter:  *image             --- buffer for image
                         rows, cols        --- size of image
                         x0, y0, dx, dy   --- region of vector used for computing threshold
                         vvv                 --- debug option, is 0, no debug information outputed
     */
    /*======================================================================*/
    /*   OTSU global thresholdingroutine                                                */
    /*   takes a 2D unsigned char array pointer, number of rows,and        */
    /*   number of cols in the array. returns the value of thethreshold       */
    /*======================================================================*/
    int otsu (unsigned char *image, int rows, int cols, int x0, int y0, int dx, intdy, int vvv)
    {

      unsigned char*np;      // 图像指针
      int thresholdValue=1; // 阈值
      intihist[256];            // 图像直方图,256个点

      int i, j,k;          // various counters
      int n, n1, n2, gmin, gmax;
      double m1, m2, sum, csum, fmax, sb;

      // 对直方图置零...
      memset(ihist, 0, sizeof(ihist));

      gmin=255; gmax=0;
      // 生成直方图
      for (i = y0 + 1; i < y0 + dy - 1; i++) {
        np = &image[i*cols+x0+1];
        for (j = x0 + 1; j < x0 + dx - 1; j++) {
          ihist[*np]++;
          if(*np > gmax) gmax=*np;
          if(*np < gmin) gmin=*np;
          np++; /* next pixel */
        }
      }

      // set up everything
      sum = csum = 0.0;
      n = 0;

      for (k = 0; k <= 255; k++) {
        sum += (double) k * (double) ihist[k];    /*x*f(x) 质量矩*/
        n   +=ihist[k];                                        /*  f(x)    质量    */
      }

      if (!n) {
        // if n has no value, there is problems...
        fprintf (stderr, "NOT NORMAL thresholdValue =160\n");
        return (160);
      }

      // do the otsu global thresholdingmethod
      fmax = -1.0;
      n1 = 0;
      for (k = 0; k < 255; k++) {
        n1 += ihist[k];
        if (!n1) { continue; }
        n2 = n - n1;
        if (n2 == 0) { break; }
        csum += (double) k *ihist[k];
        m1 = csum / n1;
        m2 = (sum - csum) / n2;
        sb = (double) n1 *(double) n2 *(m1 - m2) * (m1 - m2);
        /* bbg: note: can be optimized. */
        if (sb > fmax) {
          fmax = sb;
          thresholdValue = k;
        }
      }

      // at this point we have ourthresholding value

      // debug code to displaythresholding values
      if ( vvv & 1 )
      fprintf(stderr,"# OTSU: thresholdValue = %d gmin=%dgmax=%d\n",
         thresholdValue, gmin, gmax);

      return(thresholdValue);
    }

    matlab的OTSU算法;

    function threshold=ostu(filename);

     

    x=imread('filename');

    %figure;

    %imshow(x);

    [m,n]=size(x);

    N=m*n;

    num=zeros(1,256);

    p=zeros(1,256);

     

    for i=1:m

    for j=1:n

    num(x(i,j)+1)=num(x(i,j)+1)+1;

    end

    end

     

    for i=0:255;

    p(i+1)=num(i+1)/N;

    end

     

    totalmean=0;

    for i=0:255;

    totalmean=totalmean+i*p(i+1);

    end

     

    maxvar=0;

     

    for k=0:255

    kk=k+1;

    zerosth=sum(p(1:kk));

     

    firsth=0;

    for h=0:k

    firsth=firsth+h*p(h+1);

    end

     

    var=totalmean*zerosth-firsth;

    var=var*var;

    var=var/(zerosth*(1-zerosth)+0.01);

    var=sqrt(var);

    if(var>maxvar)

    maxvar=var;

    point=k;

    end

     

    end

     

    threshold=point;

     


    网址:http://wenku.baidu.com/link?url=tOgi4aWIUlfxOSYWzXFRQnicPsTGDa048Yj9g_jySVmInlwnyhdCJp6qFVCQtYAqscDnnu7G4PCu-NDrw8uJ08Kt34gEptNX3GmcCWE-gE7

    展开全文
  • otsu算法

    万次阅读 2013-01-11 19:42:49
    otsu法(最大类间方差法,有时也... 所以 可以在二值化的时候 采用otsu算法来自动选取阈值进行二值化。otsu算法被认为是图像分割中阈值选取的最佳算法,计算简单,不受图像亮度和对比度的影响。因此,使类间方差最大的分

               otsu法(最大类间方差法,有时也称之为大津算法)使用的是聚类的思想,把图像的灰度数按灰度级分成2个部分,使得两个部分之间的灰度值差异最大,每个部分之间的灰度差异最小,通过方差的计算来寻找一个合适的灰度级别 来划分。  所以 可以在二值化的时候 采用otsu算法来自动选取阈值进行二值化。otsu算法被认为是图像分割中阈值选取的最佳算法,计算简单,不受图像亮度和对比度的影响。因此,使类间方差最大的分割意味着错分概率最小。


    设t为设定的阈值。

    wo: 分开后  前景像素点数占图像的比例

    uo:  分开后  前景像素点的平均灰度

    w1:分开后  被景像素点数占图像的比例

    u1:  分开后  被景像素点的平均灰度

    u=w0*u0 + w1*u1 :图像总平均灰度


    从L个灰度级遍历t,使得t为某个值的时候,前景和背景的方差最大, 则 这个 t  值便是我们要求得的阈值。

    其中,方差的计算公式如下:

    g=wo * (uo - u) * (uo - u) + w1 * (u1 - u) * (u1 - u)

    [             此公式计算量较大,可以采用:      g = wo * w1 * (uo - u1) * (uo - u1)                ]




    由于otsu算法是对图像的灰度级进行聚类,so  在执行otsu算法之前,需要计算该图像的灰度直方图。


    按照上面的解释参考代码如下:

    #include "stdafx.h"
    #include "stdio.h"
    #include "cv.h"
    #include "highgui.h"
    #include "Math.h"
    
    int Otsu(IplImage* src);
    
    int _tmain(int argc, _TCHAR* argv[])
    {
    	IplImage* img = cvLoadImage("c:\\aSa.jpg",0);
    	IplImage* dst = cvCreateImage(cvGetSize(img), 8, 1);
    	int threshold = Otsu(img);
    
    	cvThreshold(img, dst, threshold, 255, CV_THRESH_BINARY);
    
    
    	cvNamedWindow( "img", 1 );
    	cvShowImage("img", dst);
    
    
    	cvWaitKey(-1);
    
    	cvReleaseImage(&img);
    	cvReleaseImage(&dst);
    
    	cvDestroyWindow( "dst" );
    	return 0;
    }
    
    int Otsu(IplImage* src)  
    {  
    	int height=src->height;  
    	int width=src->width;      
    	long size = height * width; 
    
    	//histogram  
    	float histogram[256] = {0};  
    	for(int m=0; m < height; m++)
    	{  
    		unsigned char* p=(unsigned char*)src->imageData + src->widthStep * m;  
    		for(int n = 0; n < width; n++) 
    		{  
    			histogram[int(*p++)]++;  
    		}  
    	}  
    
    	int threshold;    
    	long sum0 = 0, sum1 = 0; //存储前景的灰度总和和背景灰度总和
    	long cnt0 = 0, cnt1 = 0; //前景的总个数和背景的总个数
    	double w0 = 0, w1 = 0; //前景和背景所占整幅图像的比例
    	double u0 = 0, u1 = 0;  //前景和背景的平均灰度
    	double variance = 0; //最大类间方差
    	int i, j;
    	double u = 0;
    	double maxVariance = 0;
    	for(i = 1; i < 256; i++) //一次遍历每个像素
    	{  
    		sum0 = 0;
    		sum1 = 0; 
    		cnt0 = 0;
    		cnt1 = 0;
    		w0 = 0;
    	    w1 = 0;
    		for(j = 0; j < i; j++)
    		{
    			cnt0 += histogram[j];
    			sum0 += j * histogram[j];
    		}
    
    		u0 = (double)sum0 /  cnt0; 
    		w0 = (double)cnt0 / size;
    
    		for(j = i ; j <= 255; j++)
    		{
    			cnt1 += histogram[j];
    			sum1 += j * histogram[j];
    		}
    
    		u1 = (double)sum1 / cnt1;
    		w1 = 1 - w0; // (double)cnt1 / size;
    
    		u = u0 * w0 + u1 * w1; //图像的平均灰度
    		printf("u = %f\n", u);
    		//variance =  w0 * pow((u0 - u), 2) + w1 * pow((u1 - u), 2);
    		variance =  w0 * w1 *  (u0 - u1) * (u0 - u1);
    		if(variance > maxVariance) 
    		{  
    			maxVariance = variance;  
    			threshold = i;  
    		} 
    	}  
    
    	printf("threshold = %d\n", threshold);
    	return threshold;  
    }  



    快哭了 把w1写成w0 ··害我debug 了好久~~总是不认真,脑袋浑浑噩噩的···这都看不出来。。。。委屈


    ==================

    对了,之前搜集的一个otsu的算法,代码如下(由于时间太久了,不知道出处了。。。膜拜大牛哈偷笑

    #include "stdafx.h"
    #include "stdio.h"
    #include "cv.h"
    #include "highgui.h"
    #include "Math.h"
    
    int Otsu(IplImage* src);
    
    int _tmain(int argc, _TCHAR* argv[])
    {
    	IplImage* img = cvLoadImage("c:\\aSa.jpg",0);
    	IplImage* dst = cvCreateImage(cvGetSize(img), 8, 1);
    	int threshold = Otsu(img);
    	printf("threshold = %d\n", threshold);
    	cvThreshold(img, dst, threshold, 255, CV_THRESH_BINARY);
    
    	cvNamedWindow( "img", 1 );
    	cvShowImage("img", dst);
    
    
    	cvWaitKey(-1);
    
    	cvReleaseImage(&img);
    	cvReleaseImage(&dst);
    	
    	cvDestroyWindow( "dst" );
    	return 0;
    }
    
    int Otsu(IplImage* src)  
    {  
    	int height=src->height;  
    	int width=src->width;      
    
    	//histogram  
    	float histogram[256] = {0};  
    	for(int i=0; i < height; i++)
    	{  
    		unsigned char* p=(unsigned char*)src->imageData + src->widthStep * i;  
    		for(int j = 0; j < width; j++) 
    		{  
    			histogram[*p++]++;  
    		}  
    	}  
    	//normalize histogram  
    	int size = height * width;  
    	for(int i = 0; i < 256; i++)
    	{  
    		histogram[i] = histogram[i] / size;  
    	}  
    
    	//average pixel value  
    	float avgValue=0;  
    	for(int i=0; i < 256; i++)
    	{  
    		avgValue += i * histogram[i];  //整幅图像的平均灰度
    	}   
    
    	int threshold;    
    	float maxVariance=0;  
    	float w = 0, u = 0;  
    	for(int i = 0; i < 256; i++) 
    	{  
    		w += histogram[i];  //假设当前灰度i为阈值, 0~i 灰度的像素(假设像素值在此范围的像素叫做前景像素) 所占整幅图像的比例
    		u += i * histogram[i];  // 灰度i 之前的像素(0~i)的平均灰度值: 前景像素的平均灰度值
    
    		float t = avgValue * w - u;  
    		float variance = t * t / (w * (1 - w) );  
    		if(variance > maxVariance) 
    		{  
    			maxVariance = variance;  
    			threshold = i;  
    		}  
    	}  
    
    	return threshold;  
    } 



    展开全文
  • 在MATLAB下实现OTSU算法,另外还有关于此算法的改进形式,对图像的阈值进行最佳的计算,提高二值化的效果
  • otsu算法仿真

    2016-05-22 14:10:38
    otsu算法仿真程序
  • otsu算法测试用图片,关于显微镜视图的二值划分,通过otsu算法可以将灰度图进行最佳二值分割。otsu算法测试用图片,关于显微镜视图的二值划分,通过otsu算法可以将灰度图进行最佳二值分割。
  • OTSU算法二维matlab实现

    2020-10-04 01:30:14
    OTSU算法 二维 matlab 代码脉络清晰 稍做修改可转化为C代码 OTSU算法 二维 matlab 代码脉络清晰 稍做修改可转化为C代码
  • Otsu算法原理及实现

    万次阅读 多人点赞 2017-10-21 13:04:29
    Otsu算法原理: Otsu算法实现: Otsu算法应用: 参考: https://en.wikipedia.org/wiki/Otsu%27s_method
  • OTSU算法matlab实现

    2017-08-11 17:26:04
    基于Matlab的OTSU算法实现。Matlab自带的全局阈值函数把阈值归一化了,这个代码没有归一化,直接输出全局阈值灰度值。
  • Otsu算法原理

    千次阅读 2019-08-26 23:16:26
    目录Otsu 算法Opencv API先给出原图直接二值化Ostu手动实现Ostu欢迎一起来参与leetcode刷题项目 Otsu 算法算法就是通过直方图的特性来自动选取阈值,尤其是有的图像的灰度直方图是双峰。 Opencv API 先给出原图 ...
  • C# OTSU算法

    2012-09-27 17:14:15
    这是一个用C#实现的OTSU算法,运行环境visual stdio 2008
  • Otsu算法实现

    2012-03-12 20:23:26
    图像分割中的OTSU算法,自己亲自试验,很有用的代码哦
  • OTSU算法是由日本学者OTSU于1979年提出的一种对图像进行二值化的高效算法。 利用阈值将原图像分成前景,背景两个图象。 前景:用n1,csum,m1来表示在当前阈值下的前景的点数,质量矩,平均灰度 背景:用n2, sum-...
  • 针对图像处理在追踪运动物体领域的应用,选择Otsu图像...最后在Opencv2.4.8和微软VS2010环境下对Otsu算法在该系统的应用进行测试,结果表明,Otsu算法能够快速高效的进行车辆检测及跟踪定位,具有较好的工程应用价值。
  • Otsu算法的研究及改进

    2010-07-16 17:59:50
    Otsu算法的研究及改进 Otsu算法Otsu算法
  • 二值化算法:Otsu算法、Bernsen算法、Niblack算法、循环阈值算法、迭代二值化算法等matlab代码,入门级
  • Otsu算法进行自动阈值分割
  • 基于移位协方差矩阵Otsu算法的信源数目估计
  • 设计并实现了一个基于OTSU算法的FPGA实时绕距测量系统。首先设计了视频图像灰度化的非浮点运算实现,然后详细讨论了OTSU算法的硬件实现方案,包括其原理、公式简化、流水线处理等。经过OTSU算法处理之后,接着通过...
  • Otsu算法理解

    2019-11-13 12:17:33
    Otsu算法理解 L:灰度级(0∼255为256个灰度级)L: 灰度级(0\sim255为256个灰度级)L:灰度级(0∼255为256个灰度级) ni:表示灰度级为i的像素个数n_i: 表示灰度级为i的像素个数ni​:表示灰度级为i的像素个数 MN:图像的像素...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 1,066
精华内容 426
关键字:

otsu算法