Kittlerの方法でしきい値を求めるのを、エクセル2007とC#WPFで書いてみた

ヒストグラムを2分する時のしきい値はどこが最適なのかを求める
この場合なら左右対称なので真ん中

しきい値を4にして0～3と4～7に2分すれば丁度いい

それぞれの範囲の値から平均値と分散を計算して、それを使って正規分布をグラフにすると

こんな感じのができて、しきい値付近を拡大すると


青が0～3の範囲なので、3以上は誤差とみなして、その範囲は


この塗りつぶした範囲
4～7の赤の誤差に当たるのは4以下、この範囲も塗りつぶして表示してみると


この2つの誤差の合計が小さくなるところを探すのがKittlerの方法

式はこれで、Eが誤差でこの値が最小になったときのしきい値が求めるものになる
式をエクセルっぽく書くと

こうなって、記号の意味は

A範囲はしきい値の左側になるところ、B範囲は右側
logは

2を何乗すれば8になるのか、答えは3とか言うものらしくて、エクセルにもLOG関数が用意されていた

logの右側に小さく表記される数値は底と呼んで、これが10のときは特別に表記を省略することがあるらしい、これもエクセルには関数が用意されていた

LOG10関数、エクセルすごい
さっきのKittlerの式はどうやらこの底が省略表記されているものみたいなので、それで計算してみる
分散はVARP関数

さっきのヒストグラムのもとになる値の集団は

これで、これを並べて

ヒストグラムっぽくして、これで範囲の計算
範囲は、しきい値が2だった場合は

0～1がA範囲、2～7がB範囲

しきい値が3だった場合は

0～2がA範囲、3～7がB範囲ってことにしていく

結果

誤差が最小になったのはしきい値4のとき
左範囲がA範囲、右範囲がB範囲

Average関数で範囲の平均値

これは今回は使ってないわ、今気づいた

VARP関数で範囲の分散

Count関数で範囲の要素数、これを全要素数で割って範囲の割合

求めた分散と要素割合をLOG10関数で範囲の誤差量

IFで分岐しているのは、分散が0だったら誤差になりようがないのでエラー扱いにするために、誤差量を大きな値1にしている、これでこのしきい値は対象外になる

最後に左右の誤差量を合計して、その中から最小値になったところのしきい値が答えになる


別の値

いいと思う

値が偏っている場合

これは

一部でエラーになって、最小値もマイナスの値になっているけど、しきい値は6と判定、これは悪くないと思うけど、誤差がマイナスになるのはどうなんだろう？絶対値で考えればいいのかなあ？マイナスもエラー扱いにすると次の最小値になる閾値は5、これでもいいねえ

左右と中央にだけ値


左右対称だから期待するしきい値は4

いいねえ

ちょっと意地悪な0～5まで同じ個数のフラットな値

期待するしきい値は3


これもエラー出てるけど、最小値になる閾値は3、いいねえ

C#で書いてみる

//配列からヒストグラムの配列作成
private int[] MakeHistogram0to7(int[] vs)
{
    int[] histogram = new int[8];
    for (int i = 0; i < vs.Length; i++)
    {
        histogram[vs[i]]++;
    }
    return histogram;
}

テストなので配列の値は0～7のint型で決め打ちしているけど、実際の画像処理では256で決め打ちになる

配列を渡してヒストグラム作成したところ

/// <summary>
/// ヒストグラムの指定範囲のピクセルの個数
/// </summary>
/// <param name="histogram">ヒストグラム</param>
/// <param name="begin">範囲の始まり</param>
/// <param name="end">範囲の終わり(未満なので、100指定なら99まで計算する)</param>
/// <returns></returns>
private int HistogramCount(int[] histogram, int begin, int end)
{
    int count = 0;
    for (int i = begin; i < end; i++)
    {
        count += histogram[i];
    }
    return count;
}

ピクセルの個数ってあるけど要素数のこと
ヒストグラムを使ってカウント、比率の計算に使う

/// <summary>
/// ヒストグラムの指定範囲の分散を計算
/// </summary>
/// <param name="begin">範囲の始まり</param>
/// <param name="end">範囲の終わり(未満なので、100指定なら99まで計算する)</param>
/// <param name="count">範囲の画素数</param>
/// <param name="average">範囲の平均値</param>
/// <returns></returns>
private double HistogramVariance(int[] histogram, int begin, int end)
{
    long squareTotal = 0;
    long aveTotal = 0;
    long count = 0;//要素数
    for (long i = begin; i < end; i++)
    {
        squareTotal += i * i * histogram[i];//2乗の累計
        aveTotal += i * histogram[i];
        count += histogram[i];
    }
    //平均値
    double average = aveTotal / (double)count;
    //分散 = 2乗の平均 - 平均値の2乗
    return (squareTotal / (double)count) - (average * average);
}

forのカウンタでLong型を使っているのは、もしint型だと2乗の累計の掛け算のところで桁あふれするから
それにしてもヒストグラムを使うと計算が楽になるなあ

/// <summary>
/// ヒストグラムの指定範囲の誤差を返す、誤差 = 要素の比率 * Log10(分散 / 要素の比率)
/// </summary>
/// <param name="histogram"></param>
/// <param name="begin">範囲の開始点</param>
/// <param name="end">範囲の終わり(未満なので、100指定なら99まで計算する)</param>
/// <returns></returns>
private double KittlerSub(int[] histogram, int begin, int end)
{
    double varp = HistogramVariance(histogram, begin, end);//分散
    if (double.IsNaN(varp) || varp == 0)
        //分散が計算不能or0なら対象外になるように、大きな値(1.0)を返す
        return 1.0;
    else
    {
        int count = HistogramCount(histogram, 0, 7);//全要素数
        double ratio = HistogramCount(histogram, begin, end);
        ratio /= count;//画素数比率
        return ratio * Math.Log10(varp / ratio);
    }
}

今回の要
MathクラスにもLOGを計算するLog関数がある、しかもLog10関数もある！
範囲の要素数が0だった場合は分散が計算できなくてNaNっていう値で返ってくるのと、分散が0だった場合は対象外なので弾いている

//Kittlerの方法で閾値を求める、0～7までの値の配列専用
private int KittlerThreshold(int[] vs)
{
    int[] histogram = MakeHistogram0to7(vs);
    double aError, bError;
    double min = double.MaxValue;
    int threshold = 1;//閾値
    for (int i = 1; i < 8; i++)
    {
        aError = KittlerSub(histogram, 0, i);//A範囲(しきい値未満)の誤差
        bError = KittlerSub(histogram, i, 8);//B範囲(しきい値以上)の誤差
        var E = aError + bError;
        //誤差の総量が最小になるインデックスを閾値にする
        if (E < min)
        {
            min = E;
            threshold = i;
        }
    }
    return threshold;
}

これに配列を渡すとしきい値が返ってくる

今回のアプリ
20200413_Kittler2値化閾値.zip

github.com

ボタン押すだけ

一番上の数値列が、もとの配列の値
その値をヒストグラムっぽく表示して、その中に計算したしきい値を表示している
1～5番は決まった配列で

一番下のランダムは

毎回違う値の配列から計算





エクセルでテストした時と同じになったかな

---

コード全部
MainWindow.xaml

<Window x:Class="_20200413_Kittler2値化閾値.MainWindow"
        xmlns="http://schemas.microsoft.com/winfx/2006/xaml/presentation"
        xmlns:x="http://schemas.microsoft.com/winfx/2006/xaml"
        xmlns:d="http://schemas.microsoft.com/expression/blend/2008"
        xmlns:mc="http://schemas.openxmlformats.org/markup-compatibility/2006"
        xmlns:local="clr-namespace:_20200413_Kittler2値化閾値"
        mc:Ignorable="d"
        Title="MainWindow" Height="250" Width="500">
  <Grid>
    <StackPanel Orientation="Horizontal">
      <StackPanel Margin="10">
        <Button x:Name="Button1" Content="button1" Click="Button1_Click" Margin="0,4"/>
        <Button x:Name="Button2" Content="button2" Click="Button2_Click" Margin="0,4"/>
        <Button x:Name="Button3" Content="button3" Click="Button3_Click" Margin="0,4"/>
        <Button x:Name="Button4" Content="button4" Click="Button4_Click" Margin="0,4"/>
        <Button x:Name="Button5" Content="button5" Click="Button5_Click" Margin="0,4"/>
        <Button x:Name="ButtonRandom" Content="ランダム" Click="ButtonRandom_Click" Margin="0,4"/>
      </StackPanel>
      <StackPanel Margin="10">
        <TextBlock x:Name="TextBlock1" Text="text"/>
        <Separator Background="MediumAquamarine" Margin="0,10"/>
        <TextBlock x:Name="TextBlock2" Text="text"/>
        <StackPanel Orientation="Horizontal">
          <Separator Background="MediumAquamarine" Width="100"/>
          <TextBlock x:Name="TextBlockThreshold" Text="text" Margin="4,0"/>
          <Separator Background="MediumAquamarine" Width="100"/>
        </StackPanel>
        <TextBlock x:Name="TextBlock3" Text="text"/>
      </StackPanel>
    </StackPanel>
  </Grid>
</Window>

MainWindow.xaml.cs

using System;
using System.Windows;


//0～7の値の配列を使ってKittlerの方法を使って閾値を求める
namespace _20200413_Kittler2値化閾値
{

    public partial class MainWindow : Window
    {
        public MainWindow()
        {
            InitializeComponent();
            this.Title = this.ToString();
        }

        //得られた閾値とヒストグラムをアプリに表示
        private void Hyouzi(int[] vs, int threshold)
        {
            TextBlockThreshold.Text = $"閾値＝{threshold}";
            string str;
            str = string.Join(", ", vs);
            TextBlock1.Text = str;
            int[] histogram = MakeHistogram0to7(vs);

            TextBlock2.Text = GetStr(histogram, 0, threshold);
            TextBlock3.Text = GetStr(histogram, threshold, 8);
        }
      //表示する文字列作成
        private string GetStr(int[] histogram, int begin, int end)
        {
            string str = "";
            for (int i = begin; i < end; i++)
            {
                for (int j = 0; j < histogram[i]; j++)
                {
                    str += i + ", ";
                }
                str += Environment.NewLine;
            }
            str = str.Remove(str.Length - 2);
            return str;
        }

        //Kittlerの方法で閾値を求める、0～7までの値の配列専用
        private int KittlerThreshold(int[] vs)
        {
            int[] histogram = MakeHistogram0to7(vs);
            double aError, bError;
            double min = double.MaxValue;
            int threshold = 1;//閾値
            for (int i = 1; i < 8; i++)
            {
                aError = KittlerSub(histogram, 0, i);//A範囲(しきい値未満)の誤差
                bError = KittlerSub(histogram, i, 8);//B範囲(しきい値以上)の誤差
                var E = aError + bError;
                //誤差の総量が最小になるインデックスを閾値にする
                if (E < min)
                {
                    min = E;
                    threshold = i;
                }
            }
            return threshold;
        }

        /// <summary>
        /// ヒストグラムの指定範囲の誤差を返す、誤差 = 要素の比率 * Log10(分散 / 要素の比率)
        /// </summary>
        /// <param name="histogram"></param>
        /// <param name="begin">範囲の開始点</param>
        /// <param name="end">範囲の終わり(未満なので、100指定なら99まで計算する)</param>
        /// <returns></returns>
        private double KittlerSub(int[] histogram, int begin, int end)
        {
            double varp = HistogramVariance(histogram, begin, end);//分散
            if (double.IsNaN(varp) || varp == 0)
                //分散が計算不能or0なら対象外になるように、大きな値(1.0)を返す
                return 1.0;
            else
            {
                int count = HistogramCount(histogram, 0, 7);//全要素数
                double ratio = HistogramCount(histogram, begin, end);
                ratio /= count;//画素数比率
                return ratio * Math.Log10(varp / ratio);
            }
        }

        /// <summary>
        /// ヒストグラムの指定範囲の分散を計算
        /// </summary>
        /// <param name="begin">範囲の始まり</param>
        /// <param name="end">範囲の終わり(未満なので、100指定なら99まで計算する)</param>
        /// <param name="count">範囲の画素数</param>
        /// <param name="average">範囲の平均値</param>
        /// <returns></returns>
        private double HistogramVariance(int[] histogram, int begin, int end)
        {
            long squareTotal = 0;
            long aveTotal = 0;
            long count = 0;//要素数
            for (long i = begin; i < end; i++)
            {
                squareTotal += i * i * histogram[i];//2乗の累計
                aveTotal += i * histogram[i];
                count += histogram[i];
            }
            //平均値
            double average = aveTotal / (double)count;
            //分散 = 2乗の平均 - 平均値の2乗
            return (squareTotal / (double)count) - (average * average);
        }

        /// <summary>
        /// ヒストグラムの指定範囲のピクセルの個数
        /// </summary>
        /// <param name="histogram">ヒストグラム</param>
        /// <param name="begin">範囲の始まり</param>
        /// <param name="end">範囲の終わり(未満なので、100指定なら99まで計算する)</param>
        /// <returns></returns>
        private int HistogramCount(int[] histogram, int begin, int end)
        {
            int count = 0;
            for (int i = begin; i < end; i++)
            {
                count += histogram[i];
            }
            return count;
        }

        //未使用
        /// <summary>
        /// ヒストグラムの指定範囲の平均値
        /// </summary>
        /// <param name="histogram">ヒストグラム</param>
        /// <param name="begin">範囲開始点</param>
        /// <param name="end">範囲終了点(未満なので、100指定なら99まで計算する)</param>
        /// <returns></returns>
        private double HistogramAverage(int[] histogram, int begin, int end)
        {
            long total = 0;
            long count = 0;
            for (int i = begin; i < end; i++)
            {
                total += i * histogram[i];
                count += histogram[i];
            }
            return total / (double)count;
        }


        //配列からヒストグラムの配列作成
        private int[] MakeHistogram0to7(int[] vs)
        {
            int[] histogram = new int[8];
            for (int i = 0; i < vs.Length; i++)
            {
                histogram[vs[i]]++;
            }
            return histogram;
        }


        private void Button1_Click(object sender, RoutedEventArgs e)
        {
            //0～7までの値の配列
            var vs = new int[] {
                3,  3,  3,  3,  3,
                3,  3,  3,  3,  3,
                5,  5,  5,  5,  5,
                4,  4,  4,  4,  4,
                2,  2,  2,  6,  1,
                5,  7,  0,  6,  1,};
            Hyouzi(vs, KittlerThreshold(vs));
        }

        private void Button2_Click(object sender, RoutedEventArgs e)
        {
            int[] vs = new int[] {
                0,0,0,0,0,
                0,0,0,0,0,
                4,4,4,4,4,
                3,3,3,3,3,
                7,7,7,7,7,
                7,7,7,7,7 };
            Hyouzi(vs, KittlerThreshold(vs));
        }

        private void Button3_Click(object sender, RoutedEventArgs e)
        {
            var vs = new int[] {
            7,  7,  7,  7,  7,
            7,  7,  7,  7,  7,
            7,  7,  7,  7,  7,
            6,  1,  7,  6,  4,
            5,  5,  4,  6,  5,
            3,  6,  2,  4,  3,
            };
            Hyouzi(vs, KittlerThreshold(vs));
        }

        private void Button4_Click(object sender, RoutedEventArgs e)
        {
            var vs = new int[] {
            2,  2,  2,  2,  0,
            2,  2,  2,  2,  7,
            5,  5,  5,  5,  1,
            5,  5,  5,  5,  6,
            3,  3,  3,  1,  1,
            4,  4,  4,  6,  6,
            };
            Hyouzi(vs, KittlerThreshold(vs));
        }

        private void Button5_Click(object sender, RoutedEventArgs e)
        {
            var vs = new int[] {
            0,  0,  0,  0,  0,
            1,  1,  1,  1,  1,
            2,  2,  2,  2,  2,
            3,  3,  3,  3,  3,
            4,  4,  4,  4,  4,
            5,  5,  5,  5,  5,
            };
            Hyouzi(vs, KittlerThreshold(vs));
        }

        private void ButtonRandom_Click(object sender, RoutedEventArgs e)
        {
            var r = new Random();
            var vs = new int[30];
            for (int i = 0; i < vs.Length; i++)
            {
                vs[i] = r.Next(8);
            }
            Hyouzi(vs, KittlerThreshold(vs));
        }
    }
}

---

参照したところ
www.thothchildren.com
図を見てわかりやすいのはここ

広範囲文字認識のための動画像を用いた2値画像生成手法
https://www.researchgate.net/publication/49940273_guangfantongwenzirenshinotamenodonghuaxiangwoyongita2zhihuaxiangshengchengshoufa
ここのリンクにあるPDFファイル

www.jstage.jst.go.jp ここのリンクにあるPDFファイル
読んでもわかんないけどグラフだけ見て、大津の2値化とは違うんだなあってのがわかる

kenyu.red

Kittlerの方法はググっても難しいのしか出てこない、解る人には式があるんだからそれでいいんだろうけどねえ、僕クラスになると見かた読み方がわからん
LOGは高校で習ったって言うけど、全く憶えがないんだよねえ「俺のログには何も無いな」
ラジアン、分散、行列を教えてもらってないのは文系高校＋年代で確定していて、もしかしたらLOGもそうなのかも(ログ疑惑)
あと、分散を標準偏差と勘違いしていたせいで期待するしきい値が出てこなくて3日くらい詰まっていた
今回のもこれであっているのわからん、いつもどおりだな👈😺ヨシ

---

関連記事
次の記事は2日後
gogowaten.hatenablog.com 2値化画像アプリにKittlerの方法を追加した

前回のWPF記事は4日前
gogowaten.hatenablog.com

大津の2値化を試したのは1年4ヶ月前
gogowaten.hatenablog.com

1年前
gogowaten.hatenablog.com 今回の記事書いてる途中で前にもLOG関数使ったことある？なあと思って検索したら1年前にあった、すっかり忘れてたわ