灰度共生矩阵GLCM的计算

首先贴一个链接：http://www.fxyqpx.org/ZGGX/html/gx20150311.htm#zz
各种概念性的东西论文里都有了，但是我在网上查找C#代码的时候，老是找不全，要么就是有些博客写的有问题。

这是计算图像各个特征参数的函数：

        /// <summary>
        /// 
        /// </summary>
        /// <param name="bmp">传入一张图片</param>
        /// <param name="radius">需要计算的区域的半径</param>
        /// <param name="cenX">圆心横坐标</param>
        /// <param name="cenY">圆心纵坐标</param>
        /// <param name="JJ">能量</param>
        /// <param name="HH">熵</param>
        /// <param name="GG">对比度</param>
        /// <param name="QQ">同质性</param>
        /// <param name="CC">相关性</param>
         public void getEntropy(Bitmap bmp, int radius, int cenX, int cenY, out double[] JJ, out double[] HH, out double[] GG, out double[] QQ, out double[] CC
         {
            double[, ,] p = getP(bmp, radius, cenX, cenY);//得到共生矩阵

            double J = 0;//能量，J值越大，纹理越粗，反之越细
            double H = 0;//熵，H值越大，纹理较多
            double G = 0;//对比度，G值越大，纹理越清晰
            double Q = 0;//同质性，Q值越大，图像纹理不同区域间变化小，局部非常均匀
            JJ = new double[3];
            HH = new double[3];
            GG = new double[3];
            QQ = new double[3];
            CC = new double[3];//相关性
            double[] ui = new double[3];
            double[] uj = new double[3];
            double[] si = new double[3];
            double[] sj = new double[3];
            for (int i = 0; i < 3; i++)
            {
                ui[i] = 0;
                uj[i] = 0;
                si[i] = 0;
                sj[i] = 0;
                CC[i] = 0;
            }
            for (int n = 0; n < 3; n++)
            {
                for (int ii = 0; ii < 16; ii++)
                    for (int jj = 0; jj < 16; jj++)
                    {
                        if (p[ii, jj, n] != 0)
                        {
                            J += p[ii, jj, n] * p[ii, jj, n];
                            H -= p[ii, jj, n] * Math.Log(p[ii, jj, n]);
                            G += (ii - jj) * (ii - jj) * p[ii, jj, n];
                            Q += p[ii, jj, n] / (1 + (ii - jj) * (ii - jj));
                        }
                        ui[n] = ii * p[ii, jj, n] + ui[n];
                        uj[n] = jj * p[ii, jj, n] + uj[n];
                    }
                JJ[n] = J;
                HH[n] = H;
                GG[n] = G;
                QQ[n] = Q;
            }
            for (int n = 0; n < 3; n++)
            {
                for (int ii = 0; ii < 16; ii++)
                    for (int jj = 0; jj < 16; jj++)
                    {
                        si[n] = (ii - ui[n]) * (ii - ui[n]) * p[ii, jj, n] + si[n];
                        sj[n] = (jj - uj[n]) * (ii - uj[n]) * p[ii, jj, n] + sj[n];
                        CC[n] = ii * jj * p[ii, jj, n] + CC[n];
                    }
                CC[n] = (CC[n] - ui[n] * uj[n]) / (si[n] * sj[n]);
            }
        }

下面是计算灰度共生矩阵的代码

        public double[,,] getP(Bitmap bmp, int radius, int cenX, int cenY)//计算共生矩阵，灰度级弄成16
        {
            int[,] g = new int[2*radius,2*radius];
            //转换灰度级
            for(int w=cenX-radius;w<cenX+radius;w++)
                for (int h = cenY - radius; h < cenY + radius; h++)
                {
                    int Gray = (30 * bmp.GetPixel(w, h).R + 59 * bmp.GetPixel(w, h).G + 11 * bmp.GetPixel(w, h).B)/100;//计算灰度值
                    g[w - cenX + radius, h - cenY + radius] = Gray / 16;
                }
            double[, ,] p = new double[16, 16, 3];//三个方向的，向右一步，向下一步，向右下一步
            for (int n = 0; n < 3; n++)//都转换成概率
                for (int i = 0; i < 16; i++)
                    for (int j = 0; j < 16; j++)
                    {
                        p[i, j, n] = 0;
                    }
            for (int m = 0; m < 16; m++)
                for (int n = 0; n < 16;n++ )
                    for (int w = 0; w < 2 * radius; w++)
                    {
                        for (int h = 0; h < 2 * radius; h++)
                        {
                            if (g[w, h] == m)
                            {
                                if (h < 2 * radius - 1)
                                    if (g[w, h + 1] == n)
                                        p[m, n, 0] = p[m, n, 0] + 1;
                                if (h < 2 * radius - 1 && w < 2 * radius - 1)
                                    if (g[w + 1, h + 1] == n)
                                        p[m, n, 1] = p[m, n, 1] + 1;
                                if (w < 2 * radius - 1)
                                    if (g[w + 1, h] == n)
                                        p[m, n, 2] = p[m, n, 2] + 1;

                            }
                        }

                    }

            for (int n = 0; n < 3; n++)//都转换成概率
            {
                double sum = 0;
                for (int i = 0; i < 16; i++)
                    for (int j = 0; j < 16; j++)
                    {
                        if (p[i, j, n] != 0)
                            sum += p[i, j, n];
                    }
                for (int i = 0; i < 16; i++)
                    for (int j = 0; j < 16; j++)
                    {
                        if (p[i, j, n] != 0)
                            p[i, j, n]=p[i,j,n]/sum;
                    }
            }
            return p;
        }