OtsuThreshold.cs

using IdApp.Cv.Statistics;
 
namespace IdApp.Cv.Transformations.Thresholds
{
    public static partial class TransformationOperations
    {
        public static float OtsuThreshold(this Matrix<float> M)
        {
            return M.OtsuThreshold(256);
        }
 
        public static float OtsuThreshold(this Matrix<float> M, int Buckets)
        {
            return M.OtsuThreshold(Buckets, 0, 1f);
        }
 
        public static float OtsuThreshold(this Matrix<float> M, int Buckets, float Min, float Max)
        {
            int[] Histogram = M.Histogram(Buckets, Min, Max);
            int NrPixels = M.Width * M.Height;
            int i, j, c = Buckets;
 
            long SumB = 0;
            long WSumHistogram = 0;
            float t, mF, wF, wB = 0;
            float BestT = 0;
            float MaxT = 0;
 
            for (i = 1; i < c; i++)
                WSumHistogram += i * Histogram[i];
 
            for (i = 0; i < c; i++)
            {
                wF = NrPixels - wB;
                if (wB > 0 && wF > 0)
                {
                    mF = (WSumHistogram - SumB) / wF;
                    t = (SumB / wB) - mF;
                    t = wB * wF * t * t;
                    if (t >= MaxT)
                    {
                        BestT = i;
                        MaxT = t;
                    }
                }
 
                j = Histogram[i];
                wB += j;
                SumB += i * j;
            }
 
            return BestT * (Max - Min) / Buckets + Min;
        }
 
        public static int OtsuThreshold(this Matrix<int> M)
        {
            return M.OtsuThreshold(256);
        }
 
        public static int OtsuThreshold(this Matrix<int> M, int Buckets)
        {
            return M.OtsuThreshold(Buckets, 0, 0x01000000);
        }
 
        public static int OtsuThreshold(this Matrix<int> M, int Buckets, int Min, int Max)
        {
            int[] Histogram = M.Histogram(Buckets, Min, Max);
            int NrPixels = M.Width * M.Height;
            int i, j, c = Buckets;
 
            long SumB = 0;
            long WSumHistogram = 0;
            float t, mF, wF, wB = 0;
            float BestT = 0;
            float MaxT = 0;
 
            for (i = 1; i < c; i++)
                WSumHistogram += i * Histogram[i];
 
            for (i = 0; i < c; i++)
            {
                wF = NrPixels - wB;
                if (wB > 0 && wF > 0)
                {
                    mF = (WSumHistogram - SumB) / wF;
                    t = (SumB / wB) - mF;
                    t = wB * wF * t * t;
                    if (t >= MaxT)
                    {
                        BestT = i;
                        MaxT = t;
                    }
                }
 
                j = Histogram[i];
                wB += j;
                SumB += i * j;
            }
 
            return (int)(BestT * (Max - Min) / Buckets + Min + 0.5);
        }
    }
}