灰度化图像常见的阈值分割算法
1.OTSU算法。
关于最大类间方差法（otsu）的性能:
类间方差法对噪音和目标大小十分敏感，它仅对类间方差为单峰的图像产生较好的分割效果。
当目标与背景的大小比例悬殊时，类间方差准则函数可能呈现双峰或多峰，此时效果不好，但是类间方差法是用时最少的。
1 public int Otsu(Bitmap map) 2 { 3 int nThreshold = 0; 4 5 BitmapData mapData = map.LockBits(new Rectangle(0, 0, map.Width, map.Height), ImageLockMode.ReadWrite, map.PixelFormat); 6 double[] nHistogram = new double[256]; //灰度直方图 7 double[] dVariance = new double[256]; //类间方差 8 int N = map.Height * map.Width; //总像素数 9 for (int i = 0; i < 256; i++) 10 { 11 nHistogram[i] = 0.0; 12 dVariance[i] = 0.0; 13 } 14 unsafe 15 { 16 byte* ptr = (byte*)(mapData.Scan0); 17 for (int i = 0; i < map.Height; i++) 18 { 19 ptr = (byte*)mapData.Scan0 + i * mapData.Stride; 20 for (int j = 0; j < map.Width; j++) 21 { 22 nHistogram[(*ptr)]++; //建立直方图 23 ptr = ptr + 3; 24 } 25 } 26 } 27 map.UnlockBits(mapData); 28 double Pa = 0.0; //背景出现概率 29 double Pb = 0.0; //目标出现概率 30 double Wa = 0.0; //背景平均灰度值 31 double Wb = 0.0; //目标平均灰度值 32 double W0 = 0.0; //全局平均灰度值 33 double dData1 = 0.0, dData2 = 0.0; 34 for (int i = 0; i < 256; i++) //计算全局平均灰度 35 { 36 nHistogram[i] /= N; 37 W0 += i * nHistogram[i]; 38 } 39 for (int i = 0; i < 256; i++) //对每个灰度值计算类间方差 40 { 41 Pa += nHistogram[i]; 42 Pb = 1 - Pa; 43 dData1 += i * nHistogram[i]; 44 dData2 = W0 - dData1; 45 Wa = dData1 / Pa; 46 Wb = dData2 / Pb; 47 dVariance[i] = (Pa * Pb * Math.Pow((Wb - Wa), 2)); 48 } 49 //遍历每个方差，求取类间最大方差所对应的灰度值 50 double temp = 0.0; 51 for (int i = 0; i < 256; i++) 52 { 53 if (dVariance[i] > temp) 54 { 55 temp = dVariance[i]; 56 nThreshold = i; 57 } 58 } 59 return nThreshold; 60 }
2.一维交叉熵值法
1 public int LiLee(Bitmap map) 2 { 3 int nThreshold = 0; 4 5 BitmapData mapData = map.LockBits(new Rectangle(0, 0, map.Width, map.Height), ImageLockMode.ReadWrite, map.PixelFormat); 6 7 double[] dHistogram = new double[256]; //灰度直方图 8 double[] dEntropy = new double[256]; //每个像素的交叉熵 9 int N = map.Height * map.Width; //总像素数 10 for (int i = 0; i < 256; i++) 11 { 12 dHistogram[i] = 0.0; 13 dEntropy[i] = 0.0; 14 } 15 unsafe 16 { 17 byte* ptr = (byte*)(mapData.Scan0); 18 for (int i = 0; i < map.Height; i++) 19 { 20 ptr = (byte*)mapData.Scan0 + i * mapData.Stride; 21 for (int j = 0; j < map.Width; j++) 22 { 23 dHistogram[(*ptr)]++; //建立直方图 24 ptr = ptr + 3; 25 } 26 } 27 } 28 map.UnlockBits(mapData); 29 30 double Pa = 0.0; //区域1平均灰度值 31 double Pb = 0.0; //区域2平均灰度值 32 double P0 = 0.0; //全局平均灰度值 33 double Wa = 0.0; //第一部分熵 34 double Wb = 0.0; //第二部分的熵 35 double dData1 = 0.0, dData2 = 0.0; //中间值 36 double dData3 = 0.0, dData4 = 0.0; //中间值 37 38 for (int i = 0; i < 256; i++) //计算全局平均灰度 39 { 40 dHistogram[i] /= N; 41 P0 += i * dHistogram[i]; 42 } 43 44 for (int i = 0; i < 256; i++) 45 { 46 Wa = Wb = dData1 = dData2 = dData3 = dData4 = Pa = Pb = 0.0; 47 for (int j = 0; j < 256; j++) 48 { 49 if (j <= i) 50 { 51 dData1 += dHistogram[j]; 52 dData2 += j * dHistogram[j]; 53 } 54 else 55 { 56 dData3 += dHistogram[j]; 57 dData4 += j * dHistogram[j]; 58 } 59 } 60 Pa = dData2 / dData1; 61 Pb = dData4 / dData3; 62 for (int j = 0; j < 256; j++) 63 { 64 if (j <= i) 65 { 66 if ((Pa != 0) && (dHistogram[j] != 0)) 67 { 68 double d1 = Math.Log(dHistogram[j] / Pa); 69 Wa += j * dHistogram[j] * d1 / Math.Log(2); 70 } 71 } 72 else 73 { 74 if ((Pb != 0) && (dHistogram[j] != 0)) 75 { 76 double d2 = Math.Log(dHistogram[j] / Pb); 77 Wb += j * dHistogram[j] * d2 / Math.Log(2); 78 } 79 } 80 } 81 dEntropy[i] = Wa + Wb; 82 } 83 //遍历熵值，求取最小交叉熵所对应的灰度值 84 double temp = dEntropy[0]; 85 for (int i = 1; i < 256; i++) 86 { 87 if (dEntropy[i] < temp) 88 { 89 temp = dEntropy[i]; 90 nThreshold = i; 91 } 92 } 93 94 return nThreshold; 95 }
3.二维OTSU算法
1 public int TwoOtsu(Bitmap map) 2 { 3 int nThreshold = 0; 4 Bitmap mapClone = (Bitmap)map.Clone(); 5 BitmapData mapData = map.LockBits(new Rectangle(0, 0, map.Width, map.Height), ImageLockMode.ReadWrite, map.PixelFormat); 6 7 double[,] dHistogram = new double[256, 256]; //建立二维灰度直方图 8 double dTrMatrix = 0.0; //离散矩阵的迹 9 int N = map.Height * map.Width; //总像素数 10 for (int i = 0; i < 256; i++) 11 { 12 for (int j = 0; j < 256; j++) 13 dHistogram[i, j] = 0.0; //初始化变量 14 } 15 unsafe 16 { 17 byte* ptr = (byte*)mapData.Scan0; 18 for (int i = 0; i < map.Height; i++) 19 { 20 ptr = (byte*)mapData.Scan0 + i * mapData.Stride; 21 for (int j = 0; j < map.Width; j++) 22 { 23 int nData1 = *ptr; //当前的灰度值 24 int nData2 = 0; 25 int nData3 = 0; //注意9个值相加可能超过一个字节 26 for (int m = i - 1; m <= i + 1; m++) 27 { 28 for (int n = j - 1; n <= j + 1; n++) 29 { 30 if ((m >= 0) && (m < map.Height) && (n >= 0) && (n < map.Width)) 31 { 32 int a = mapClone.GetPixel(m, n).R; 33 nData3 += a;//当前的灰度值 34 }; 35 } 36 } 37 nData2 = nData3 / 9; //对于越界的索引值进行补零,邻域均值 38 dHistogram[nData1, nData2]++; 39 ptr = ptr + 3; 40 } 41 } 42 } 43 map.UnlockBits(mapData); 44 for (int i = 0; i < 256; i++) 45 for (int j = 0; j < 256; j++) 46 dHistogram[i, j] /= N; //得到归一化的概率分布 47 48 double Pai = 0.0; //目标区均值矢量i分量 49 double Paj = 0.0; //目标区均值矢量j分量 50 double Pbi = 0.0; //背景区均值矢量i分量 51 double Pbj = 0.0; //背景区均值矢量j分量 52 double Pti = 0.0; //全局均值矢量i分量 53 double Ptj = 0.0; //全局均值矢量j分量 54 double W0 = 0.0; //目标区的联合概率密度 55 double W1 = 0.0; //背景区的联合概率密度 56 double dData1 = 0.0; 57 double dData2 = 0.0; 58 double dData3 = 0.0; 59 double dData4 = 0.0; //中间变量 60 int nThreshold_s = 0; 61 int nThreshold_t = 0; 62 double temp = 0.0; //寻求最大值 63 for (int i = 0; i < 256; i++) 64 { 65 for (int j = 0; j < 256; j++) 66 { 67 Pti += i * dHistogram[i, j]; 68 Ptj += j * dHistogram[i, j]; 69 } 70 } 71 for (int i = 0; i < 256; i++) 72 { 73 for (int j = 0; j < 256; j++) 74 { 75 W0 += dHistogram[i, j]; 76 dData1 += i * dHistogram[i, j]; 77 dData2 += j * dHistogram[i, j]; 78 79 W1 = 1 - W0; 80 dData3 = Pti - dData1; 81 dData4 = Ptj - dData2; 82 /* W1=dData3=dData4=0.0; //对内循环的数据进行初始化 83 for(int s=i+1; s<256; s++) 84 { 85 for(int t=j+1; t<256; t++) 86 { 87 W1 += dHistogram[s][t]; 88 dData3 += s*dHistogram[s][t]; //方法2 89 dData4 += t*dHistogram[s][t]; //也可以增加循环进行计算 90 } 91 }*/ 92 93 Pai = dData1 / W0; 94 Paj = dData2 / W0; 95 Pbi = dData3 / W1; 96 Pbj = dData4 / W1; // 得到两个均值向量，用4个分量表示 97 dTrMatrix = ((W0 * Pti - dData1) * (W0 * Pti - dData1) + (W0 * Ptj - dData1) * (W0 * Ptj - dData2)) / (W0 * W1); 98 if (dTrMatrix > temp) 99 { 100 temp = dTrMatrix; 101 nThreshold_s = i; 102 nThreshold_t = j; 103 } 104 } 105 } 106 nThreshold = nThreshold_t; //返回结果中的灰度值 107 108 return nThreshold; 109 }
