c#+ArcEngine实现OTSU方法计算栅格数据的二值化阈值并渲染

写作前面：关于本文方法的二值图片打开为全黑的问题：1、请在envi中打开即显示为黑白图；2、若在arcmap打开，请在符号系统中-设置根据值显示不同颜色(因为gis中默认最大value值与最小value值显示为相同颜色)。
OTSU方法是一种很好的阈值计算方法，在图像分割中经常用到，关于otsu原理网上有很多，这里直接上代码实现。
在计算完阈值后，添加栅格数据渲染方法显示分割效果。以此验证阈值计算的效果。

public static int Otsu_Binaryzation(IRasterLayer RasterLayer)
        {
            int OTSU = 0;      

        IRaster pRaster = RasterLayer.Raster;
        IRasterBandCollection pBandCol = pRaster as IRasterBandCollection;
        IRasterBand pBand = pBandCol.Item(0);//选择第一波段作为分级依据
        
        IRasterProps pRasterProps = pRaster as IRasterProps;            
        //获取图层的行列值   
        int Height = pRasterProps.Height;//当前栅格数据集的行数
        int Width = pRasterProps.Width;//当前栅格数据集的列数

        if (pBand.Histogram == null)
        {
            pBand.ComputeStatsAndHist();//创建直方图
        }

        IRasterStatistics pRasterStatistic = pBand.Statistics;
        double MaxValue = pRasterStatistic.Maximum;
        double MinValue = pRasterStatistic.Minimum;
        double MeanValue = pRasterStatistic.Mean;

        #region OTSU阈值计算
        IRasterHistogram pRasterHistogram = pBand.Histogram;
        double[] HistValues;
        HistValues = pRasterHistogram.Counts as double[];
        int HistValueCount = HistValues.GetUpperBound(0) + 1;
        double[] CountPixel = new double[HistValueCount];

        double w1 = 0, w2 = 0;
        double sum = HistValueCount;
        double mu1, mu2;
        double numerator = 0, denominator;
        double sigma;
        double tempmax = 0;

        int deta = (int)(0 - MinValue);
        //i为灰度级
        for (int i = (int)MinValue, j = 0; j < HistValueCount; i++,j++)
        //for (int i = 0; i <= HistValueCount; i++)//确保灰度值为0~~255
        {
            sum += i * HistValues[j];
        }
        for (int i = 0,j=(int)(MinValue+0.5); i < HistValueCount; i++,j++)
        {
            w1 += HistValues[i];
            numerator += j * HistValues[i];
            mu1 = numerator / w1;
            w2 = Width * Height - w1;
            mu2 = (sum - numerator) / w2;

            sigma = w1 * w2 * (mu1 - mu2) * (mu1 - mu2);

            //迭代
            if (sigma > tempmax)
            {
                tempmax = sigma;
                OTSU = j;
            }
        }
        //Class1.StrechDN(RasterLayer, OTSU,Convert .ToInt32(MinValue -1) , Convert.ToInt32(MaxValue+10));//赋值
        #endregion
        IRasterClassifyColorRampRenderer pRClassRend = new RasterClassifyColorRampRenderer() as IRasterClassifyColorRampRenderer;
        IRasterRenderer pRRend = pRClassRend as IRasterRenderer;

        pRRend.Raster = pRaster;
        pRRend.Update();
        pRClassRend.ClassCount = 2;//分为两类

        pRClassRend.set_Break(1, OTSU);
        pRClassRend.set_Break(2, MaxValue);

        IRgbColor pFromColor = new RgbColor() as IRgbColor;//设置颜色
        pFromColor.Red = 0; pFromColor.Green = 0; pFromColor.Blue = 0;//黑色

        IRgbColor pToColor = new RgbColor() as IRgbColor;//设置颜色
        pToColor.Red = 255; pToColor.Green = 255; pToColor.Blue = 255;//白色

        IPresetColorRamp colorRamp = new PresetColorRamp() as IPresetColorRamp;//明确定义每一种颜色

        colorRamp.Size = 2;//和前面的分级数对应
        colorRamp.set_PresetColor(1, pFromColor);//红色,可根据需要设定
        colorRamp.set_PresetColor(2, pToColor);//绿色

        IFillSymbol fillSymbol = new SimpleFillSymbol() as IFillSymbol;设置面填充符号    

        fillSymbol.Color = pFromColor;
        pRClassRend.set_Symbol(0, fillSymbol as ISymbol);

        fillSymbol.Color = pToColor;
        pRClassRend.set_Symbol(1, fillSymbol as ISymbol);

        RasterLayer.Renderer = pRRend;//渲染
        return OTSU;
    }

如果你觉得还可以点个赞吧，欢迎评论留言讨论。