Rustで画像の二値化とフィルター処理を自作する

はじめに

前回はRust image crateを使用して、 Rustでの画像処理を試してみて自由度が高く結構楽しかったので、image crateには該当する機能はないけどもよく使用しそうな画像処理を作成しました。 今回作成した処理は以下になります。

作成したソースコードは https://gist.github.com/tk87s/ec42fd517697ed3507ac3ae7ec4a6715 にあります。

• thresholdを指定して二値化
• 大津の二値化法
• Gaussian(ガウシアン)フィルタ
• median(メディアン)フィルタ
• bilateral(バイラテラル)フィルタ

上2つは二値化に関するものになります。画像を白と黒の二種類で表すことを二値化といいます。 マスク処理等でよく使用し背景除去(URL)の際にも使用しました。 imageクレートのソースを探しましたが意外にも二値化に関する機能は存在しませんでした。

下3つはフィルタを使用した処理で、これらのフィルタは画像のノイズを除去するのに使用されます。 注目画素に対してその周囲の画素の情報をパラメータとして計算して、その値を注目画素の輝度に置き換える時に使用します。 フィルターの範囲は中心画素に対して対称であることが好ましいため、普通は一辺が奇数マス分になります。 辺の大きさによって取得したい特徴のスケールを選択することも可能です。 一般的にノイズは周囲の画素に対して外れ値であることが多いのでフィルタで修正しやすいことが予想できると思います。 Rust imageクレートにはGaussian blurのみ存在しますが、後の二つは存在しません。 またフィルタ処理にあたる image::imageops::filer3x3 が存在しますがこれは3x3カーネルのみになるので、 filer3x3 を任意のカーネルサイズで行えるように改造しました。

thresholdを指定して二値化

これは簡単で、まずグレースケール化を行い好きな閾値で輝度を0と255にすればよいです。

    
fn threshold(img: &GrayImage, thres: u8) -> GrayImage {
    let (w, h) = img.dimensions();
    let mut bin_img = image::ImageBuffer::new(w, h);
    for (x, y, pixel) in img.enumerate_pixels() {
        let image::Luma(v) = *pixel;
        bin_img.put_pixel(x, y, image::Luma(if v[0] > thres {[255]} else {[0]}));
    }
    bin_img
}
Copy

大津の二値化法

大津の二値化は輝度のヒストグラムを用いて、白と黒を適切に分離してくれる閾値を計算して二値化を行う手法です。 全体の輝度平均値と領域間の輝度平均値の差から計算される分散の比率が大きくなる(最大)になるように分離されます。 つまり以下の式が最大となる閾値を見つければよいです。

$w_0{w}_1{\displaystyle \frac{{({\mu }_0-{\mu }_1)}^2}{{(w_0+w_1)}^2}}$

$w_0$と$w_1$はピクセル数、${\mu }_0$と${\mu }_1$は輝度の平均

大津の二値化うまくいくケースも多いですが、ヒストグラムの形によってはうまくいかないことがあります。 例えばグラデーションが入った画像や、ヒストグラムの形が二峰性でないものは苦手な傾向があります。

    
fn otsu_threshold(img: &GrayImage) -> (u8, GrayImage) {
    let (w, h) = img.dimensions();
    let mut bin_img = image::ImageBuffer::new(w, h);
    let mut count = vec![0.0; 256];
    let mut sum = vec![0.0; 256];
    for &pixel in img.iter() {
        count[pixel as usize] += 1.0;
        sum[pixel as usize] += pixel as f64;
    }
    for i in 1..256 {
        count[i] += count[i-1];
        sum[i] += sum[i-1];
    }
 
    let mut thres = 0;
    let mut mx = 0.0;
    for i in 1..255 {
        let (mut ml, mut mr): (f64, f64) = (sum[i], sum[255] - sum[i]);
        let (cl, cr): (f64, f64) = (count[i], count[255] - count[i]);
        if cl == 0.0 || cr == 0.0 {continue;}
        ml /= cl;
        mr /= cr;
        let val = (cl * cr) / (cl + cr).powf(2.0) * (ml - mr).powf(2.0);
        if val > mx {
            thres = i as u8;
            mx = val;
        }
    }
    for (x, y, pixel) in img.enumerate_pixels() {
        let image::Luma(v) = *pixel;
        bin_img.put_pixel(x, y, image::Luma(if v[0] > thres {[255]} else {[0]}));
    }
    (thres, bin_img)
}
Copy

Gaussian(ガウシアン)フィルタ

まずフィルタを適用するための実装を行います。 これはimage::imageops::filer3x3を任意の大きさのフィルターに適用できるように改造をしました。 注目画素に対して周囲の画素を表す相対位置の配列を作成し、それぞれの画素に対してカーネル内を走査して計算を行うだけです。

    
pub fn filter_nxn(image: &I, kernel: &[f32], ksize: usize) -> ImageBuffer>
    where
        I: GenericImageView,
        P: Pixel + 'static,
        S: Primitive + 'static,
    {
        assert!(kernel.len() == ksize*ksize);
    
        let r = (ksize / 2) as isize;
        let taps: Vec<(isize, isize)> = (-r..r+1).map(
            |x|
                (-r..r+1).map(|y| (x,y)).collect::>()
            )
            .flatten()
            .collect();
    
        let (width, height) = image.dimensions();
    
        let mut out = ImageBuffer::new(width, height);
    
        let max = S::DEFAULT_MAX_VALUE;
        let max: f32 = NumCast::from(max).unwrap();
    
        for y in 0..height {
            for x in 0..width {
                let mut t = vec![0.0; P::CHANNEL_COUNT as usize];
                let mut sum = 0.0;
                for (&k, &(a, b)) in kernel.iter().zip(taps.iter()) {
                    let x0 = x as isize + a;
                    let y0 = y as isize + b;
                    if x0 < 0 || y0 < 0 || x0 >= width as isize || y0>= height as isize {continue;}
    
                    let p = image.get_pixel(x0 as u32, y0 as u32);
    
                    let slice = p.channels();
    
                    let vec: Vec = slice.iter().map(|&x| NumCast::from(x).unwrap()).collect();
    
                    for (ti, &vi) in t.iter_mut().zip(vec.iter()) {
                        *ti += vi * k;
                    }
                    sum += k;
                }
                if sum == 0.0 {
                    sum = 1.0;
                }
                let t: P = *Pixel::from_slice(
                    &t.into_iter()
                    .map(|x| NumCast::from((x/sum).clamp(0.0, max)).unwrap())
                    .collect::>()
                );
                
                out.put_pixel(x, y, t);
            }
        }
        out
    }
Copy

ガウシアンフィルタの場合は以下のフィルタを適用すればよいです。 注目画素からの距離に対してガウス分布で重みづけを行います。 ガウス分布の分散が大きくなれば分布は広くなるので周囲の重みは大きくなりますのでより平滑化される傾向にあります。 上記コードはガウシアンフィルタにかかわらずSobelフィルターなど別のフィルタも適用可能です。

    
let kernel = [
    1.0/16.0, 1.0/8.0, 1.0/16.0, 
    1.0/8.0, 1.0/4.0, 1.0/8.0, 
    1.0/16.0, 1.0/8.0, 1.0/16.0, 
];
Copy

median(メディアン)フィルタ

メディアンフィルタはフィルタ内の全ての画素の中央値を適用します。 中央値は周囲の画素の輝度を配列に格納してソートし 列の長さ/2 のindexの値を取得すればよいです。 外れ値の影響を受けにくいので細かいノイズには強いですが、カーネルサイズが大きくなればぼやけやすいですのでエッジが少ない画像や部分に適用するのに向いていると思います。

    
for y in 0..height {
    for x in 0..width {
        let mut t = vec![vec![];  P::CHANNEL_COUNT as usize];
        for &(a, b) in taps.iter() {
            let x0 = x as isize + a;
            let y0 = y as isize + b;
            if x0 < 0 || y0 < 0 || x0 >= width as isize || y0>= height as isize {continue;}
 
            let p = image.get_pixel(x0 as u32, y0 as u32);
 
            let slice = p.channels();
 
            let vec: Vec = slice.iter().map(|&x| NumCast::from(x).unwrap()).collect();
 
            for (ti, &vi) in t.iter_mut().zip(vec.iter()) {
                ti.push(vi);
            }
        }
 
        for ti in t.iter_mut() {
            ti.sort_by(|a, b| a.partial_cmp(b).unwrap());
        }
        
        let t: P = *Pixel::from_slice(
            &t.into_iter()
            .map(|x| {
                    let n = x.len();
                    NumCast::from((x[n/2]).clamp(0.0, max)).unwrap()
                }
            )
            .collect::>()
        );
        
        out.put_pixel(x, y, t);
    }
}
Copy

bilateral(バイラテラル)フィルタ

上記のメディアンフィルタはエッジには弱いですが、bilateralフィルタはエッジを保持することができます。 大まかに式で表すと以下のようになります。

• $w_1=gaussian((\mathrm{カ}\mathrm{ー}\mathrm{ネ}\mathrm{ル}\mathrm{内}\mathrm{画}\mathrm{素}\mathrm{と}\mathrm{注}\mathrm{目}\mathrm{画}\mathrm{素}\mathrm{か}\mathrm{ら}\mathrm{の}\mathrm{距}\mathrm{離}{)}^2,sigma)$
• $w_2=gaussian((\mathrm{注}\mathrm{目}\mathrm{画}\mathrm{素}\mathrm{と}\mathrm{カ}\mathrm{ー}\mathrm{ネ}\mathrm{ル}\mathrm{内}\mathrm{画}\mathrm{素}\mathrm{と}\mathrm{の}\mathrm{輝}\mathrm{度}\mathrm{の}\mathrm{差}{)}^2,sigma)$
• $I={\displaystyle \frac{\sum (\mathrm{カ}\mathrm{ー}\mathrm{ネ}\mathrm{ル}\mathrm{内}\mathrm{画}\mathrm{素}\mathrm{の}\mathrm{輝}\mathrm{度})\cdot w_1{w}_2}{\sum w_1{w}_2}}$

式を見ると注目画素との距離と同時に注目画素との輝度の差を重みとして考慮しています。 つまり注目画素と似た画素を重視することで周辺の特徴が損なわれるような可能性が減るため、 ぼやけを抑止しエッジを保持することにつながります。

    
for y in 0..height {
    for x in 0..width {
        let mut t = vec![0.0;  P::CHANNEL_COUNT as usize];
        let mut sum = vec![0.0;  P::CHANNEL_COUNT as usize];
        let cp = image.get_pixel(x as u32, y as u32);
        let slice = cp.channels();
        let cv: Vec = slice.iter().map(|&x| NumCast::from(x).unwrap()).collect();
 
        for &(a, b) in taps.iter() {
            let x0 = x as isize + a;
            let y0 = y as isize + b;
            if x0 < 0 || y0 < 0 || x0 >= width as isize || y0>= height as isize {continue;}
 
            let gau = gaussian((a*a + b*b) as f32, 1.0);
            let p = image.get_pixel(x0 as u32, y0 as u32);
 
            let slice = p.channels();
 
            let offset: Vec = slice.iter().map(|&x| NumCast::from(x).unwrap()).collect();
 
            for i in 0..P::CHANNEL_COUNT as usize {
                let rw = gau * gaussian((offset[i] - cv[i]).powf(2.0), sigma);
                t[i] += offset[i] * rw;
                sum[i] += rw;
            }
        }
        
        let t: P = *Pixel::from_slice(
            &t.into_iter().zip(sum.into_iter())
            .map(|(x, s)| {
                    NumCast::from((x/s).clamp(0.0, max)).unwrap()
                }
            )
            .collect::>()
        );
        
        out.put_pixel(x, y, t);
    }
Copy

画像の比較

以下の画像はノイズを与えた元の画像に対して、それぞれフィルタ処理を施したものになります。 左から元画像、ガウシアンフィルタ、メディアンフィルタ、バイラテラルフィルタの結果になります。

上記のようなゴマ塩ノイズとよばれる細かいノイズはメディアンフィルタとの相性が良いです。 ただぼやけやすいのでエッジ部分を避けノイズの領域だけを絞って処理をするなどの工夫をするといいと思います。 その他の二つはどちらかといえば平滑化、均すような処理に近いですので、 ガウシアンノイズのような輝度が上下に少しづつブレるようなノイズには強い反面、欠陥のような輝度差の大きいノイズは消しきれないことがあります。 カーネルサイズやガウス分布のシグマ(分散)によって、ぼやけ方やノイズの除去のされ方は変わりますのでいろいろ試してみるといいかと思います。 普通はカーネルサイズや分散が大きくなればぼやけやすいですが、ノイズも隠しやすくなる傾向にあります。

まとめ

Rustを使って比較的よく使用する処理を実装しました。 imageクレートのおかげで躓いたところがなく快適に作成することができました。

Reference

1. https://crates.io/crates/image
2. https://en.wikipedia.org/wiki/Otsu%27s_method