画像の拡大

例えば、下図のように2x2画素の画像を4x4の画像に拡大する場合、アフィン変換を使えばいいんでしょ!と、安易に考えていると、思わぬ落とし穴があったりもします。

大事なポイントとして、

●画像の座標の原点は左上の画素の中心が原点(0.0、0.0)となる。(例外もあります)

●アフィン変換の拡大縮小は原点を基準として拡大縮小される。

 

となります。

これを気にせず、ただ、アフィン変換で画像を2倍に拡大すると、左上の画素の中心を基点に画像が拡大されます。

これで、一見良さそうにも感じるのですが、拡大後の画像において、画素の中心が原点であることから、4x4画素の領域は下図の四角で示した領域であり、画像全体が左上に0.5画素ズレた状態になっていまいます。

 

 

アフィン変換で画像を拡大する時の変換前と変換後の状態は、以下のようになるのが正解です。

 

この変換をアフィン変換で実現するには以下のように行います。

 

変換前の状態

 

①画像全体を右下に(+0.5,+0.5)画素移動
$$\begin{pmatrix} 1 & 0 & 0.5 \\ 0 & 1 & 0.5 \\ 0 & 0 & 1 \end{pmatrix}$$
 

②画像を2倍に拡大

$$\begin{pmatrix} 2 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 1 \end{pmatrix}$$
 

③画像全体を左上に(-0.5,-0.5)画素移動

$$\begin{pmatrix} 1 & 0 & -0.5 \\ 0 & 1 & -0.5 \\ 0 & 0 & 1 \end{pmatrix}$$

 

となります。

この一連の変換をアフィン変換行列であらわすと

$$\begin{pmatrix} { x }^{ ‘ } \\ { y }^{ ‘ } \\ 1 \end{pmatrix}=\begin{pmatrix} 1 & 0 & -0.5 \\ 0 & 1 & -0.5 \\ 0 & 0 & 1 \end{pmatrix}\begin{pmatrix} 2 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 1 \end{pmatrix}\begin{pmatrix} 1 & 0 & 0.5 \\ 0 & 1 & 0.5 \\ 0 & 0 & 1 \end{pmatrix}\begin{pmatrix} x \\ y \\ 1 \end{pmatrix}\\ \\ $$

$$\begin{pmatrix} { x }^{ ‘ } \\ { y }^{ ‘ } \\ 1 \end{pmatrix}=\begin{pmatrix} 2 & 0 & 0.5 \\ 0 & 2 & 0.5 \\ 0 & 0 & 1 \end{pmatrix}\begin{pmatrix} x \\ y \\ 1 \end{pmatrix}\\ \\ $$

 

となり、単に拡大のアフィン変換行列だけを掛ければOKでは無いことが分かります。

 

ちなみに、C#で2x2画素の画像をPictureBoxのSizeModeプロパティをZoomにして、ImageプロパティにBitmapを設定すると、このようになります。

 

なんとなく、画像が左上にズレているようで、なんか怪しい!!

 

【関連記事】

【C#】画像の座標系

アフィン変換(平行移動、拡大縮小、回転、スキュー行列)

画素の補間(Nearest neighbor,Bilinear,Bicubic)の計算方法

画像の回転

 

画像処理のためのC#へ戻る

Canny edge detectionの処理アルゴリズム

Canny edge detectionは画像の輪郭部分を抽出するのに、よく用いられるのですが、詳細なアルゴリズムを理解しないまま使われている事も多いのではないでしょうか?(それ、私)

特に、ヒステリシスしきい値の部分などは、詳細な説明も少ないので、今回、まとめてみたいと思います。

 

まずは処理結果から。

 

処理前画像 Canny edge detection処理後

 

大まかな流れとしては

1.ノイズ除去

ガウシアンフィルタなどでノイズを除去します

 

2.輪郭抽出

ソーベルフィルタで輪郭を抽出します

 

3.非極大抑制

エッジの強さが極大となる部分以外を除外します。

詳細は後述

 

4.ヒステリシスしきい値処理

詳細は後述

 

完成!!

 

以上の4ステップで成り立っていますが、非極大抑制とヒステリシスしきい値処理がちょっとわかりづらいので、詳細は以下の通りです。

 

非極大抑制

ソーベルフィルタ後のエッジの強さを3D表示すると、こんな感じ↓になります。

 

この山の尾根のみを抽出するのが、非極大抑制となります。

 

具体的には、まず、ソーベルフィルタでエッジの強さを求めて、エッジの角度を求めます。

横方向

縦方向

エッジの強さ

エッジの角度(実際にはエッジの線に直行する向きが求まります)

 

エッジの向き(エッジの法線方向)が求まったら、

0°(-22.5°~22.5°)

45°(22.5°~47.5°)

90°(47.5°~112.5°)

135°(112.5°~157.5°)

の4つに分類します。

 

エッジの法線方向の3画素のエッジの強さを用いて、中央のエッジの強さが、残す2つの良さよりも大きければ、その部分が極大点ということになり、中央のエッジの強さが最大とならない部分を除外する処理が非極大抑制となります。

エッジの強さ 非極大抑制

 

この処理によりエッジの部分の細線化が行われます。

 

ヒステリシスしきい値処理

ヒステリシスしきい値の処理については、Webで探してもあまりわかりやすいのが無いような気がしますが、以下の通りです。

 

まず、非極大抑制処理で細線化された画像に対して、2つのしきい値を用いて二値化します。

2つのしきい値は、なんとなくエッジっぽい部分(しきい値「小」)と、確実にエッジな部分(しきい値「大」)となるしきい値を設定します。

 

非極大抑制処理画像

 

しきい値「小」 しきい値「大」

 

 

この2つの画像を用いて、しきい値「小」の中から、しきい値「大」につながっている部分のみを残します。

分かりやすいように、しきい値「小」の画像にしきい値「大」のエッジを赤くして重ね合わせて表示すると↓

 

 

この赤い線につながっていない部分を除外します。

 

この処理がまさにCanny edge detectionとなります。

 

(参考)

OpenCV 3.0.0-dev documentation Canny Edge Detection

 

関連記事

【OpenCV-Python】Canny(Canny edge detection)

 

画像処理アルゴリズムへ戻る

バイラテラルフィルタ

ガウシアンフィルタなどのフィルタでは、ノイズをできるだけ除去しようとすると、輪郭もボケてしまうという欠点がありました。
この欠点を解決しようとした処理アルゴリズムがバイラテラルフィルタ(bilateral filter)です。

 

バイラテラルフィルタは処理前の画像データの配列をf(i, j)、処理後の画像データの配列をg(i, j)とすると

 

 

となります。

 

ただし、がカーネルのサイズ、σがガウシアンフィルタを制御、σが輝度差を制御しています。

と言われても、何だか式が難しくて良く分かりません。

 

でも、分母分子に出てくる最初のexpの部分はガウシアンフィルタで見たことがあるな~

という事に気が付けば、突破口が開けます。

 

2つ目のexpの部分が良く分からないので、とりあえず取っちゃってみて、

 

 

とすると、分母の部分がガウシアンフィルタと少し違うけど、Σの範囲が-W~Wなので、(2W+1)×(2W+1)のカーネルサイズのガウシアンフィルタになっています。

 

結局、分母はカーネルの値の合計なので、やっぱりガウシアンフィルタそのものだという事に気が付きます。

 

ここで、ガウシアンフィルタを使って、どうやれば輪郭をぼやかさず、ノイズだけを除去できるか?を考えると、カーネルの中心の輝度値と差の少ないところだけをガウシアンフィルタで平滑化すればいいのではないか?という発想が浮かびます。

 

例えば次のような画像において、

 

四角で囲まれた部分を拡大し注目すると、5x5のガウシアンフィルタの場合、中心の輝度値に近い部分を重み「1」、輝度差が大きい部分は重みが「0」になるようにして、

 

 

一般的な5x5のガウシアンフィルタの係数↓、

 

 

に重みをかけると、それぞれのカーネルの値は

 

カーネルの合計の169で割る カーネルの合計の209で割る

 

となります。

このようにして、場所、場所のカーネルの値を画像に合わせて変えて行くと、輪郭を残しつつノイズだけを除去できそうな感じがします。

 

しかし、輪郭付近の重みは「1」にするべきか?「0」にするべきか?少し悩みます。

そこで、カーネルの中心の輝度値との差に基づいて、重みを「0.3」や「0.8」のようなグレーゾーンを設けるために、輝度差を横軸にした正規分布を用いてみます。

 

正規分布のグラフはこんな感じ↓

 

になっていて、中心の0付近ほど値が大きく、外側(+方向、-方向)へ行くに従って値が小さくなります。

この性質を使ってカーネルの中心の輝度値との差「f(i, j)- f(i + m, j + n) 」を横軸にとった正規分布の式は

 

 

となり、これを重みに使うと、「1」「0」だった重みが、輝度差が小さいと重みが大きく、輝度差が大きいと重みが小さくなるように、なだらかに変化します。

 

この正規分布を重みとしたのが、まさにバイラテラルフィルタなのです。

つまりバイラテラルフィルタは

 

正規分布の重み付きガウシアンフィルタ

 

なのです。

 

ここで、最初に示したバイラテラルフィルタの式を見てみると、

分子は(2W+1)×(2W+1)サイズのカーネルの範囲内の輝度値に、ガウシアンフィルタの係数をかけ、さらにカーネルの中心との輝度差を用いた正規分布の値を掛け合わせた値、

分母はカーネルの合計値で、カーネルの合計値が1になるように調整しています。

 

Wはカーネルの大きさ、σは通常のガウシアンフィルタの係数と同じ。

σの値が、カーネルの中心の輝度値との差をどの程度許容するか?を制御している事が分かります。

そのため、σの値を大きくしていくと、ただのガウシアンフィルタの処理に近づき、輪郭もぼやけてしまいます。逆にσの値を小さくし過ぎると、ノイズ除去効果が弱くなります。

 

そこで、実際にはエッジを保持しつつノイズをできるだけ除去したい場合は、1回のバイラテラルフィルタでσ1、σの値を調整しようとするよりも、バイラテラルフィルタを何回か繰り返した方が効果的です。

 

原画 1回目 2回目
3回目 4回目 5回目

 

このように、バイラテラルフィルタはガウシアンフィルタのカーネルに輝度差に基づいて重みを付けている訳ですが、この、×××に基づいて重みを付ける処理という考え方は、処理を安定させるポイントになったりします。

 

例えばロバスト推定法なんかも近い感じ。

 

画像処理アルゴリズムへ戻る

疑似カラー(Pseudo-color)

三次元データやサーモグラフィのようにデータを画像にした場合には、モノクロで表示するよりも、色を付けて表示した方が見やすい場合があります。

 

モノクロ表示 疑似カラー表示

 

このグレーの色に疑似的に色を付ける方法を疑似カラー(Pseudo-color)と言います。この色の付け方は色相を使って青~緑~赤へと変化させてもいいのですが、もう少し簡単な方法を紹介します。
青~緑~赤~青へと色を変化させるには以下のようなパターンでR,G,Bの値を変化させます。

 

 

このパターンは青~緑~赤~青の色相で言うと一周分なので、一般的に用いられるのは青~緑~赤ぐらいまでなので、上図の0~240°部分を8Bitの輝度値に割り当て

 

 

のようにすると、このようなグレースケール↓に

 

疑似カラーを割りつけると、このように↓なります。

 

このようにクレースケールにカラーを割り当てることで、より画像を見やすくしています。モノクロ8Bitの画像データの場合は、画像データはそのままに、カラーテーブルを変更するだけで、疑似カラー表示する事が出来ます。

 

画像処理アルゴリズムへ戻る

 

色相、彩度、明度の計算方法

色相Hue:色合い)、彩度Saturation:鮮やかさ)、明度Brightness,Lightness,Intensity,Value:明るさ)については、以前、変換式には色相、彩度、明度ほかのページにまとめたのですが、実は訳も分からず公式だけをまとめていました。

 

で、なんだか気持ちが悪かったので色相、彩度、明度について、よ~く調べてみました。
私なりの理解ですが、以下にまとめました。

 

R、G、Bの色空間については、下図のようにRGBをXYZのように三次元座標で表すと、一辺の長さが255で表される立方体の範囲内で全ての色を表す事が出来ます。(R,G,B各8bitの場合)

 

この立方体を白(255、255、255)に位置から黒(0、0、0)の方向へ見て、R軸を右側に取ると、

 

 

のように、正六角形となります。
この時、の方向を0°として、反時計回りにの位置が120°240°色相(0~360°(2π))を定めます。
彩度は一番外側の六角形に対して、どの割合の位置に配されているかを0~1.0で表したものが彩度となります。

 

詳細は後述しますが、色相彩度はカメラやパソコンなどの性能評価(使いやすさ、価格、処理速度など)を表す時に用いるレーダーチャート(クモの巣グラフ)もどきみたいな物?!と思うと、自分の中で少し整理ができました。

 

さらに、この六角形の高さ方法に明度を割り振ると、HSV(六角錐モデル)やHLS(双六角錐モデル)となります。
それぞれの違いは明度の定義が異なり、R、G、Bの最大輝度値をImax、最小輝度値をIminとしたときに

明度V = Imax

としたものがHSV

明度L = ( Imax + Imin ) / 2

としたものがHLSとなり、明度の値は0~1.0で表されます。

これを立体で表すと

 

HSV(六角錐モデル)

HLS(双六角錐モデル)

となります。
このHSV、HLSともに、六角錐の斜面の部分が彩度が1.0となります。

 

以下、色相、彩度、明度の詳細な計算方法です。

 

HSVの計算方法

はじめにR、G、Bの輝度値の範囲を0~255から0~1.0となるように変換します。
(R、G、Bのそれぞれの値を255で割ります。)

 

【色相Hの求め方】
下図のように、0°方向にR、120°方向にG、240°方向にBだけ進み、最後の点の位置のR軸に対する角度が色相となります。

 

 

この最後の点の座標は中心を(x、y) = (0、0) とすると、R、G、Bの方向のなす角度から

 

 

となり、xとyより色相Hが求まります。

 

 

ただし、アークタンジェントの計算が出来ない場合など、この方法とは別に、近似的に求める方法もあります。(こちらの方が一般的)

 

下図を見ても分かる?ようにR、G、Bの成分の比を比べ、

Rが最大の場合、色相は-60°(300°)~60° (R方向の0°±60°)
Gが最大の場合、色相は 60°~180° (G方向の120°±60°)
Bが最大の場合、色相は 180°~300° (B方向の240°±60°)

の範囲内に色相は収まります。

 

 

以下、Rの値が最大の場合を例に取って説明したいと思います。

 

下図のように、2つの矢印の長さが分かれば、その矢印の比で角度60°を分割することで、角度(色相)を近似することが出来ます。

 

 

R、G、Bの大きさがR≧G≧Bの場合

 

色相H = 60° × (G – B) / (R – B)

R、G、Bの大きさがR≧B≧Gの場合

 

色相H = 60° × (G – B) / (R – G)

 

となります。
ただし、この場合、色相の値が負となるので、

 

色相H = 60° × (G – B) / (R – G) + 360°

 

とします。

 

と、なる理屈を理解するのに苦労しました...
図中に書いてある黄色い正三角形がポイント!
正三角形なので、三辺の長さが等しい分けで。
R以外のGやBが最大となる場合も理屈は同じです。
120°づつ回転させて考えてみると分かります。

 

この式を一般的に書くと、R、G、Bの成分のうち、最大の成分をImax、最小の成分をIminとすると

ImaxがRのとき

ImaxがGのとき

ImaxがBのとき

となります。

 

【明度Vの求め方】

明度は、もともとHSVの定義よりR、G、Bの成分のうち、最大の成分をImaxとすると

明度V = Imax

 

とします。
明度Vの範囲は0~1.0となります。

 

【彩度Sの求め方】
R、G、Bの成分のうち、最大の成分をImax、最小の成分をIminとすると

 

彩度S = (Imax – Imin) / Imax

 

となります。
彩度Sの範囲は0~1.0となります。

 

 

HLSの計算方法

【色相Hの求め方】
色相HはHSVの色相Hの求め方と同じです。

 

【明度Lの求め方】
明度Lは、もともとHLSの定義よりR、G、Bの成分のうち、最大の成分をImaxとすると

 

明度L = ( Imax + Imin ) / 2

 

とします。
明度Lの範囲は0~1.0となります。

 

【彩度Sの求め方】
R、G、Bの成分のうち、最大の成分をImax、最小の成分をIminとすると

 

L≦0.5のとき

彩度S = (Imax – Imin) / (Imax + Imin)

L>0.5のとき

彩度S = (Imax – Imin) / (2 – Imax – Imin)

 

となります。
彩度Sの範囲は0~1.0となります。

 

以下、補足説明です。
彩度Sは下図の外側の六角形に対して、内側の六角形の大きさの割合で求められます。

 

 

この六角形の大きさはR,G,Bの輝度値が最大となる軸上で考えると比較的分かりやすいと思います。
今回はRの値が最大となる場合とします。

 

HLSは双六角錐モデルであるため、明度Lが0.5以下の場合、外側の六角形の大きさ(上図のE’の位置)は

 

E’ = Imax + Imin

 

で求まります。

 

 

明度Lが0.5より大きい場合、外側の六角形の大きさ(上図のE’の位置)は

 

E’ = 2 – Imax – Imin

 

となります。

 

 

以上のことから、最初の彩度Sの式が求まります。

 

色相、彩度、明度を使った色判別時の注意点

色相および彩度を用いると、画像の明るさ(明度)が変動しても似た色の領域を抽出する事が可能となりますが、彩度の値が小さい場合、つまりR、G、Bの値がそれぞれ近い場合は色相の値が不安定になります。

 

例えば、
(R、G、B) = (121、120、120)の場合、 色相H = 0°
(R、G、B) = (120、121、120)の場合、 色相H = 120°
(R、G、B) = (120、120、121)の場合、 色相H = 240°

 

と、ほんの少しのR、G、Bの値の違いでも色相の値は大きく異なります。

 

また、色相Hは角度で表されるので、例えば1°も359°も値こそ離れていますが、どちらも0°±1°の範囲内で角度的には近いので、色相Hの値で単純に二値化処理することで色の領域を抽出する場合は注意して下さい。

 

画像処理アルゴリズムへ戻る

 

色相、彩度、明度の公式

カラーの画像処理をする時には、これら色相などの知識は必須となります。
Windows標準で付いてくるペイントで、色の作成の表示をすると、雰囲気が分かると思います。

 

 

基本的に以下の色相、彩度、明度を用いて色を表すのですが、変換式にいくつかの種類があります。

 

色相(Hue)

色合いを表します。
赤や緑、青などに色を0~360°(0~2π)の角度を用いて表します。

彩度(Saturation)

鮮やかさを表します。
と一般的に言われるのですが、鮮やかさ?と言われても、いまいちピンと来ませんが、下記に示した式から見ても分かるように、R,G,Bの値にどれだけ開きがあるか?を示しています。
このことは逆にいうと、R,G,Bの値に開きが無い場合は、グレーに近い事から、彩度は如何にグレーっぽく無いか?、という事から、どれだけ純色(赤、緑、青、黄、シアン、紫など)に近いか?を表しています。

明度(Brightness,Lightness,Intensity,Value)

色の明るさを表します。

 

HSV変換

6角錐モデルとも言います。

 

【RGB⇒HSV変換】

Imax = Max(R,G,B)
Imin = Min(R,G,B)
とすると

 

R = Imaxのとき

H = 60×(G – B) / (Imax – Imin)

G = Imaxのとき

H = 60×(B – R) / (Imax – Imin) + 120

B = Imaxのとき

H = 60×(R – G) / (Imax – Imin) + 240


S = (Imax – Imin) / Imax


V = Imax


【HSV⇒RGB変換】

h = floor(H / 60)     floor()は切り捨て処理
P = V × (1 – S)
Q = V × (1 – S × (H / 60 – h))
T = V × (1 – S × (1 – H / 60 + h))
とすると

 

h= 0のとき

R = V, G = T, B = P

h= 1のとき

R = Q, G = V, B = P

h= 2のとき

R = P, G = V, B = T

h= 3のとき

R = P, G = Q, B = V

h= 4のとき

R = T, G = P, B = V

h= 5のとき

R = V, G = P, B = Q


HLS変換

双6角錐モデルとも言います。

 

【RGB⇒HLS変換】

Imax = Max(R,G,B)
Imin = Min(R,G,B)
とすると

 

R = Imaxのとき

H = 60×(G – B) / (Imax – Imin)

G = Imaxのとき

H = 60×(B – R) / (Imax – Imin) + 120

R = Imaxのとき

H = 60×(R – G) / (Imax – Imin) + 240


L = (Imax + Imin) / 2

L ≦ 0.5のとき

S = (Imax – Imin) / (Imax + Imin)

L > 0.5のとき

S = (Imax – Imin) / (2 – Imax – Imin)


【HLS⇒RGB変換】

h < 0のとき

h’ = h  + 360

h ≧ 360のとき

h’ = h  – 360

その他

h’ = h

L ≦ 0.5のとき

M2 = L × (1 + S)

L > 0.5のとき

M2 = L + S – L × S


M1 = 2 × L – M2


h’ < 60のとき

X = M1 + (M2 – M1) × h’ / 60

60 ≦ h’ < 180のとき

X = M2

180 ≦ h’ < 240のとき

X = M1 + (M2 – M1) × (240 – h’ ) / 60

240 ≦ h’ ≦ 360のとき

X = M1

 

とすると

 

R = X ただし、h = H + 120とする

G = X ただし、h = Hとする

B = X ただし、h = H – 120とする

 

カラー変換用関数

【Win32APIの場合】 VBの表記例
‘HLS変換(Windows 2000以降、またはInternet Explorer 5.0がインストールされてある環境。(SHLWAPI.DLL Version 5.00以上)
‘h (色相)
‘赤(0)、黄(40)、緑(80)、シアン(120)、青(160)、マゼンダ(200)の順に定義0~239まで設定可
‘L (明度)
‘色の明るさをあらわす。0~240まで設定可。0が黒、240が白になる。
‘s (彩度)
‘0~240まで設定可。240が純色になる。
Public Declare Sub ColorRGBToHLS Lib “SHLWAPI.DLL” _

(ByVal clrRGB As Long, _
pwHue As Integer, _
pwLuminance As Integer, _
pwSaturation As Integer)

Public Declare Function ColorHLSToRGB Lib “SHLWAPI.DLL” _

(ByVal wHue As Integer, _
ByVal wLuminance As Integer, _
ByVal wSaturation As Integer) As Long

 

最初に紹介したペイントの色の作成では、この関数と同様の変換(設定値を含めて同じ)をしています。

 

【.NET Frameworkの場合】
System.Drawing.Color構造体にて

GetHueメソッド      HSBのH(色相)を取得
GetSaturationメソッド   HSBのS(彩度)を取得
GetBrightnessメソッド   HSBのB(明度)を取得

RGB⇒HSB変換はなし?

 

【OpenCVの場合】
cvCvtColor関数にて、以下の変換に対応
XYZ, YCrCb(YCC), HSV, HLS, L*a*b, L*u*v
その他 同じcvCvtColor関数で Bayer変換に対応

 

カラー画像処理例

そのカラー画像を色相、彩度、明度に分解し、それぞれの値(主に色相)でフィルタリング処理(バンドパスフィルタ)を行い、カラー画像に逆変換する事により、特定の色だけを抽出する事が可能になります。これにより、色の位置や個数などの検査をする事が可能となります。

処理前 処理後

ただし、彩度の値の小さな色(白やグレーに近い色)は彩度の値が不安的になりがちなので、カラー画像処理には不向きです。

 

画像処理アルゴリズムへ戻る

 

外周画像の処理

平滑化フィルタやメディアンフィルタなどの注目画素の周辺画素を用いた画像フィルタ処理では、画像の外周部分が下図のように画像の外側を参照してしまうため、処理ができなくなります。

 

 

この外周部分を処理する方法はいくつかあるのですが、代表的な方法を紹介します。
以下、5×5サイズのカーネルを用いた場合の処理を例にとって紹介します。

 

■外周部分の輝度値を画像の外側にコピーして補間する方法

5×5サイズのカーネルの場合、画像の外側に2画素分、画像の輝度値を参照
してしまうので、この2画素分の輝度値を画像の外周部分の輝度値をコピーして
輝度値を参照します。

 

おそらく?この手法が一般的だと思います。

 

■外周部分を中心にして対称の位置にある輝度値を外側にコピーして補間する方法

例えば、カーネルが座標(-1、-2)の画素を参照する場合は座標(1、2)の画素の輝度値を参照するようにします。

 

 

他にも、単純に外周部分は処理をしないで黒(輝度値=0)でマスクしたり、カーネルが画像の外側を参照する場合にカーネルを形を変えて画像の外側を参照しないようにしたりする方法などもあります。

 

画像処理アルゴリズムへ戻る

 

画像の回転

画像を回転する場合、任意点周りの回転移動でも紹介したように回転行列を使って、例えば、画像の中心周りに画像を回転させると、下図のように回転後の画像が虫食い状態になってしまいます。

 

回転前の画像 回転後の画像

 

こうならないようにするためには、回転前の画像の座標を回転行列を使って回転後の座標を計算するのではなく、回転後の座標が回転前の画像のどの座標を参照しているのかを計算し、画像を変換します。

つまり、画素を置きに行くのではなく、拾いに行くようにします。

 

【回転前の画像を回転行列を使って変換】

 

【回転後の座標が回転前のどの座標を参照しているかを計算して変換】

 

実際の変換処理は以下のように行います。

 

回転前の画像の座標を(x,y)、回転後の画像の座標を(X,Y)、画像の中心座標を(Cx,Cy)とすると単純に画像の中心周りに座標を回転すると以下のような行列で表されます。

 

 

この行列を回転前の画像の座標(x、y)に関して解けばいいので、行列の式の両辺にそれぞれ逆行列をかければいいので、以下のような手順で行列を解いていきます。

 

これで、回転前の座標(x、y)に関して解くことができます。
でも、逆行列を解くのは面倒と思う事なかれ。
(+Cx,+Cy)の平行移動の逆行列は(-Cx、-Cy)方向への移動と同じ、+θ方向への回転の逆行列は-θ方向への回転と同じなので、

 

 

とすれば、逆行列を解くことなく回転後の座標(x、y)に関して座標を解くことができます。
このようにして画像の回転処理を行うと、このようになります。

 

回転前の画像 回転後の画像

 

でも、回転後の画像はなんかギザギザしてしまっていますが、これは回転前の画像の座標を計算するとほとんどの場合、画素の画素の間の座標となってしまいますが、上図の例ではこの座標を四捨五入して輝度値を参照しているためで、bilinearやbicubicなどの補間を使うと少しは滑らかな画像となります。

今回は画像の回転について紹介していますが、画像の拡大縮小についても同様の考え方で処理することができます。

 

画像処理アルゴリズムへ戻る

 

関連記事

アフィン変換(平行移動、拡大縮小、回転、スキュー行列)

任意点周りの回転移動(アフィン変換)

画素の補間(Nearest neighbor,Bilinear,Bicubic)の計算方法

アフィン変換(平行移動、拡大縮小、回転、スキュー行列)

画像の拡大縮小、回転、平行移動などを行列を使って座標を変換する事をアフィン変換と呼びます。

X,Y座標の二次元データをアフィン変換するには、変換前の座標を(x, y)、変換後の座標を(x’,y’)とすると回転や拡大縮小用の2行2列の行列と、平行移動用に2行1列の行列を使って

$$\left(\begin{array}{c}x^{‘}\\ y^{‘}\end{array}\right)=
\left(\begin{array}{c}a & b\\ c & d\end{array}\right)
\left(\begin{array}{c}x\\ y\end{array}\right)
+ \left(\begin{array}{c}T_{x}\\ T_{y}\end{array}\right)$$

のように表現される場合もありますが、回転、拡大縮小、平行移動を1つの3x3の行列にまとめて

$$\left(\begin{array}{c}x^{‘}\\ y^{‘} \\ 1\end{array}\right)=
\left(\begin{array}{c}a & b & c\\ d & e & f\\ 0 & 0 & 1\end{array}\right)
\left(\begin{array}{c}x\\ y\\ 1\end{array}\right)$$

と3x3の行列で表現する場合もあります。

この表現を同次座標系と呼びます。

同次座標系では一見、3行目は無駄なようにも見えるのですが、この意味の無いような1行を追加する事で、平行移動も同じ行列の積で表現でき、逆行列を使う事で、アフィン変換後の座標からアフィン変換前の座標も簡単に求める事ができるようになります。

私は3行3列の行列を用いた同次座標系のアフィン変換しかしていない、というか断然おススメなので、同次座標系で説明したいと思います。

後半でアフィン変換の実用例を示しているので、その部分で、同次座標系の恩恵を感じてもらえると嬉しいです。

 

変換前の画像を以下のようにすると、

各種変換は以下の通りとなります。

 

拡大縮小

X軸方向の拡大率をSx、Y軸方向の拡大率をSyとすると拡大縮小のアフィン変換は

 

 

と表されます。

 

例)X軸方向に2倍

 

例)Y軸方向に2倍

 

例)X軸、Y軸方向に2倍

 

例)Y軸方向に-1倍

このように、ある軸(上記の例ではX軸)に対して反転する処理の事を鏡映と呼びます。

 

平行移動

X軸方向にTx、Y軸方向にTyだけ移動するアフィン変換は

 

 

のように表されます。

 

 

回転

原点を中心に反時計回りにθ°回転する時のアフィン変換は

 

 

のように表されます。

 

 

スキュー(せん断)

四角形の画像を平行四辺形に変形する処理をスキューまたはせん断といいます。
このアフィン変換は

 

 

アフィン変換の実用方法

画像処理で使われるアフィン変換は、下図のようにピクチャボックスなどの左上が原点[座標が(0, 0)]で右方向が+X、下方向が+Y、時計周りが+θ方向となる場合が多いかと思います。

また、平行移動、拡大縮小、回転などの行列を紹介しましたが、これらの行列を1回だけで処理する事はまれで、それぞれの行列を組み合わせてアフィン変換を行います。
さらに、拡大縮小、回転、スキューのアフィン変換行列は、あくまでも原点を基点として、変換される事に注意が必要です。

例えば、画像が原点の位置に無い場合、画像をX,Y方向に2倍の大きさにしようとしたとき、単に拡大縮小のアフィン変換行列で座標を計算すると、表示位置も2倍、原点の位置から離れた位置に移動します。

これらの事を踏まえ、画像を幅方向に2倍、高さ方向に3倍し、画像を反時計方向に90°回転、さらに、画像の左上の座標が(50, 50)から(30,180)へ移動するアフィン変換を例題にとって考えたいと思います。

この変換は、一発で変換するアフィン変換行列を考えるのではなく、平行移動拡大縮小回転に分けて、アフィン変換の順番を考えます。


拡大縮小と回転が原点を基点とするため、画像を原点の位置へ移動するため、x方向に-50、y方向に-50の平行移動します。

平行移動のアフィン変換行列は

$$\begin{pmatrix} 1 & 0 & -50 \\ 0 & 1 & -50 \\ 0 & 0 & 1 \end{pmatrix}$$


拡大縮小か回転のどちらからでも構いませんが、x方向に2倍、y方向に3倍の拡大縮小を行います。

拡大縮小のアフィン変換行列は

$$\begin{pmatrix} 2 & 0 & 0 \\ 0 & 3 & 0 \\ 0 & 0 & 1 \end{pmatrix}$$


反時計周りに90°回転(ー90°回転)を行います。

回転のアフィン変換行列は

$$\begin{pmatrix} cos(-90) & -sin(-90) & 0 \\ sin(-90) & cos(-90) & 0 \\ 0 & 0 & 1 \end{pmatrix}$$


最後に目的の位置へ移動するのに、x方向に+30、y方向に+180の平行移動します。

平行移動のアフィン変換行列は

$$\begin{pmatrix} 1 & 0 & 30 \\ 0 & 1 & 180 \\ 0 & 0 & 1 \end{pmatrix}$$


このようにアフィン変換を平行移動、拡大縮小、回転に分解して、変換の手順を考える事が大事です。
今回の場合は

変換前 → 平行移動 → 拡大縮小 → 回転 → 平行移動 → 変換後

としています。

画像に対してアフィン変換を行う場合、考え方としては、平行移動や拡大縮小などに分解して考えますが、画像データを毎回変換するのではなく、アフィン変換行列をまとめて計算し、一発でアフィン変換処理を行います。

具体的には、今回のアフィン変換処理を行列で表すと

$$\begin{pmatrix} { x }^{ ‘ } \\ { y }^{ ‘ } \\ 1 \end{pmatrix}=\begin{pmatrix} 1 & 0 & 30 \\ 0 & 1 & 180 \\ 0 & 0 & 1 \end{pmatrix}\begin{pmatrix} cos(-90) & -sin(-90) & 0 \\ sin(-90) & cos(-90) & 0 \\ 0 & 0 & 1 \end{pmatrix}\begin{pmatrix} 2 & 0 & 0 \\ 0 & 3 & 0 \\ 0 & 0 & 1 \end{pmatrix}\begin{pmatrix} 1 & 0 & -50 \\ 0 & 1 & -50 \\ 0 & 0 & 1 \end{pmatrix}\begin{pmatrix} x \\ y \\ 1 \end{pmatrix}$$

行列部分をまとめて計算し

$$\begin{pmatrix} { x }^{ ‘ } \\ { y }^{ ‘ } \\ 1 \end{pmatrix}=\begin{pmatrix} 0 & 3 & -120 \\ -2 & 0 & 280 \\ 0 & 0 & 1 \end{pmatrix}\begin{pmatrix} x \\ y \\ 1 \end{pmatrix}$$

となり、アフィン変換行列を一発で処理する事ができます。

さらに、変換後の座標(x’, y’)から、変換前の座標(x, y)を求める場合があるのですが、その時は、行列の左側からアフィン変換行列の逆行列を掛けて、

$$\begin{pmatrix} 0 & 3 & -120 \\ -2 & 0 & 280 \\ 0 & 0 & 1 \end{pmatrix}^{ -1 }\begin{pmatrix} { x }^{ ‘ } \\ { y }^{ ‘ } \\ 1 \end{pmatrix}=\begin{pmatrix} 0 & 3 & -120 \\ -2 & 0 & 280 \\ 0 & 0 & 1 \end{pmatrix}^{ -1 }\begin{pmatrix} 0 & 3 & -120 \\ -2 & 0 & 280 \\ 0 & 0 & 1 \end{pmatrix}\begin{pmatrix} x \\ y \\ 1 \end{pmatrix}$$

$$\begin{pmatrix} x \\ y \\ 1 \end{pmatrix}=\begin{pmatrix} 0 & 3 & -120 \\ -2 & 0 & 280 \\ 0 & 0 & 1 \end{pmatrix}^{ -1 }\begin{pmatrix} { x }^{ ‘ } \\ { y }^{ ‘ } \\ 1 \end{pmatrix}$$

$$\begin{pmatrix} x \\ y \\ 1 \end{pmatrix}=\begin{pmatrix} 0 & -0.5 & 140 \\ 0.333 & 0 & 40 \\ 0 & 0 & 1 \end{pmatrix}\begin{pmatrix} { x }^{ ‘ } \\ { y }^{ ‘ } \\ 1 \end{pmatrix}$$

として、変換前の座標を求める事ができます。

まとめ

  • アフィン変換は平行移動、拡大縮小、回転、スキューに分けて、変換の順番を考える
  • 実際のアフィン変換は、アフィン変換行列をまとめて計算し、一発で処理を行う
  • 変換前の座標は逆行列を使うと、求めることができる

この順番の考え方は、ほとんどの場合、画像を原点へ平行移動し、拡大縮小、回転を行ってから、目的の位置へ移動すると、アフィン変換行列が求まりますが、他にも、特定の点を基準に拡大縮小や回転を行う場合は、その点を原点へ移動すればアフィン変換行列が求まります。

(参考)

 

注意点

アフィン変換では任意の3×3(2×3)の行列で表す事ができるので、任意形状に変換できそうにも思えるのですが、四角形が平行四辺形にまでは変形できるものの、台形には変形できないのでご注意願います。
この台形に変形できる処理は射影変換(ホモグラフィ)と呼びます。

アフィン変換は今回の説明のように、画像を移動、変形させるための手法として説明されますが、もう少し汎用的に座標変換として捉えると応用範囲が広がります。

例えば、データのグラフを表示する時に、横方向、縦方向に拡大/縮小した時に、表示するデータの範囲を求める場合などに応用すると、少し便利です。

 

画像処理アルゴリズムへ戻る

 

関連記事

実際にOpenCVを使って行うアフィン変換については、こちらのページ↓にまとめました。

【OpenCV-Python】アフィン変換(同次座標系を使用)

 

マイクロソフトの.NETやDirectXで扱うアフィン変換行列は、行と列が逆になり、行列を掛ける順番も逆(右側から掛ける)になります。
この仕様については、下記ページ↓にまとめました。

【C#.NET】マイクロソフト仕様のアフィン変換

円形度

輪郭追跡を行うと周囲長を求めることができますが、この周囲長を用いた代表的な特徴量の円形度を紹介します。

 

円形度とは円らしさを表す値で値が1となる時、もっとも円に近くなります。
定義は面積(画素数)をS周囲長Lとすると、

 

円形度 = 4πS/L2

 

となります。

 

【円形度の例】

面積S = 96
周囲長L = 34.627円形度 = 1.006
面積S = 87
周囲長L = 33.456円形度 = 0.977

なぜそうなるのか?理屈は簡単で、
円の半径とすると、円の面積Sπr周囲長L2πrより、

 

円形度 = 4πS / L2 = 4π(πr2) / (2πr)2 = 1

 

となり、円形度4πS/L2に近いときもっとも円らしいということになります。

 

画像処理アルゴリズムへ戻る

関連記事

【OpenCV-Python】円形度

輪郭追跡処理アルゴリズム

二値化された画像において、各連結部分の境界部分を求める事を輪郭追跡といいます。
この輪郭追跡を行う事で周囲長や連結領域の高さなどを求めることが出来ます。

 

以下、この輪郭追跡の処理アルゴリズムを紹介します。

 

まず、左上からラスタスキャンを行い、白の部分を見つけます。

 

 

最初の白の部分(上図の赤丸の部分)をスタート地点として、反時計回りに輪郭の部分を検出します。
最初の白の部分は左上からラスタスキャンしたため、スタート地点の右上、上、左上、左には白の部分が無い事が確定しているので、最初に左下の部分から反時計回りに白の部分を検索し、白の部分が見つかったら次の輪郭へ検索を続けます。
最後にスタート地点に戻ったら処理は終了です。

 

この時、輪郭の向きを下図のように定義しておくと、

 

最初の図の輪郭は

 

0→0→1→0→2→3→4→3→2→5→5→6→4→0→6→5

 

の順で変化しています。(この番号はチェインコードといいます。)
ここにちょっとした規則があります。
例えば、の次に5、1の次に6、・・・などは絶対に来ません。(なぜかは、ちょっと考えると分かります。)
一般的に表現すると、
一つ前の輪郭の向きから時計回りに3つ分の向きには輪郭が存在しない。
という事になります。
この性質を利用して、輪郭を検索するときは、
一つ前の輪郭の向きの値をVold
輪郭の検索を開始する向きの値を Vnew
とすると、
Vnew = (Vold + 6) % 8; (%は8で割った時の余りの計算)
となるVnew の向きから輪郭を検索すれば効率的に輪郭を検索することが可能となります。

 

また、輪郭の周囲長を求める時、縦、横方向の輪郭は長さ、斜め方向の輪郭は√2 として
長さを足していくと、周囲長が求まります。

 

そして、ここからがあまり書かれていない注意点!
スタート地点の白の画素のパターンが以下のパターンとなる場合、スタート地点を2回通る事になるので、ご注意下さい。

※グレーの部分は白でも黒でもよい部分

 

画像処理アルゴリズムへ戻る

 

ラベリング

二値化画像処理された画像において、白の部分(または黒の部分)が連続した画素に同じ番号を割り振る処理を
ラベリングと言います。

 

通常、同じ番号ごとの面積(画素数)や幅、高さなどの特徴量を求めて欠陥検査や分類処理などに用いられます。

 

ラベリングには、二値化された画像の縦、横方向に連続している部分を同じラベルにする4連結と、
縦、横、斜め方向に連続している部分を同じラベルにする8連結との2種類の処理があります。
(4近傍、8近傍という場合もあります。)

【二値化画像】
【4連結】 【8連結】

以下、8連結の場合において、ラベリング処理アルゴリズムを紹介します。

 

まず、画像全ての画素のラベル番号を0(ゼロ)で初期化しておき、ラベリングで番号を割り付けるための
ラベリング番号のルックアップテーブルを用意しておきます。(テーブルの使い方の詳細は後ほど)

 

そして、画像の左上からラスタスキャンを行い、画素の色がの位置を検索します。

 

白の画素の左上、上、右上、左の画素のラベル番号を参照し、全て0(ゼロ)の場合は、最後に割り振った番号+1のラベル番号を割り振ります。
もし、参照した画素のラベル番号が複数存在した場合は、最小の番号を割り振ります。

 

最初に番号を割り振った様子↓

 

この処理を左上からラスタスキャンして続けていきます。

 

新しい番号が割り振られるとき↓

 

新しい番号が割り振られるとき↓ その2

 

参照した画素のラベル番号が複数存在した場合、最小の番号を割り振ります。
このとき、用いなかったラベル番号(下図の例では3)のルックアップテーブルの番号を最小の番号に書き換えます。

 

同様に、処理を続けていきます↓

 

全ての白の画素に番号が割り振れれた様子↓

 

ただし、よ~く見ると、ラベル番号の2や3など、連続していない部分が存在している事に気がつきます。
そこで、ルックアップテーブルの出番です。

 

処理の途中でルックアップテーブルの番号を変更した番号を
2→1、3→1、6→5
に変更すると、連続した領域は同じ番号になります。

 

このまま番号を変更しても良いのですが、ラベル番号が虫食いになるのも気持ちが悪いので、以下のようにルックアップテーブルの番号を詰め直してから、ルックアップテーブルを参照し、全画素のラベル番号を修正すると、以下のようになります。

 

基本的なアルゴリズムは上記、紹介したアルゴリズムで良いのですが、ラベリング処理を行う画像のパターンによってはうまくラベル番号が割り振られない場合があります。
そこは、何とかして(ちょっと説明が難しい部分)うまく割り振られるようにチャレンジしてみて下さい。

 

参考までに、ラベリング処理をすると失敗しやすい特徴的な画像を添付しておきます。↓

◆Level.1

 

◆Level.2

 

◆Level.Max!?

 

ラベリング処理のアルゴリズムには、この手法の他にも輪郭線を追跡し、閉じた領域を1つのラベルとする手法や、ランレングス的に求める方法もあります。

 

画像処理アルゴリズムへ戻る

 

細線化

二値化された画像において、線の中心1画素分だけを残すように線を細くする処理を細線化といいます。細線化のアルゴリズムにはHilditch、田村の方法、Zhang Suenなど、いろいろあるのですが、ここでは田村の方法について説明します。
細線化は黒の線を細くする場合と白の線を細くする場合が考えられますが、以下では白の線を細くするのを前提として説明しますので、黒の線を細くする場合は白と黒を読み替えて下さい。

 

細線化処理前 細線化処理後

 

ざっくり言うと、収縮処理をすると線幅が細くなりますが、線の長さが短くなったり、細い線は消えてしまったりするので、線の端点や交点、線幅が1となった画素などを除去しないように条件を付けながら収縮処理を行います。

 

処理をフローチャートを以下に示します。

 

 

画像全体にラスタスキャンを行い、注目画素の3×3画素の並びがパターン1に該当する場合は中心の画素を除去(白から黒に変える)します。ただし、除去しないパターンに該当する場合は除去を行いません。

パターン1に該当する画素が1画素も無ければその時点で終了です。

次にパターン1の時と同様に3×3画素の並びがパターン2に該当する場合は中心の画素を除去(白から黒に変える)します。ただし、除去しないパターンに該当する場合は除去を行いません。

パターン2に該当する画素が1画素も無ければその時点で終了です。そうでない場合はパターン1の処理に戻り、除去する画素がなくなるまでパターン1、パターン2の除去処理を繰り返します。

 

 

 

画像処理アルゴリズムへ戻る

 

膨張・収縮・オープニング・クロージング

膨張・収縮処理では一般的に二値化された白黒の画像に対して処理が行われ、

注目画素の周辺に1画素でも白い画素があれば白に置き換える処理を膨張(Dilation)
逆に周辺に1画素でも黒い画素があれば黒に置き換える処理を収縮(Erosion)といいます。

 

【元画像】

 

【膨張処理】

 

【収縮処理】

二値化された画像ではなく、グレー画像に対して処理を行う場合は、膨張の場合、注目画素の近傍の最大輝度値を注目画素の輝度値に置き換えます。
収縮の場合は最小輝度値に置き換えることでグレー画像に対して処理を行います。
この処理を二値化画像と区別するために、それぞれ最大値フィルタ最小値フィルタと呼ぶ場合もあります。

 

【処理例】

元画像 膨張 膨張
膨張 収縮 収縮
収縮 収縮 収縮
収縮 膨張 膨張
膨張 連続処理

膨張・収縮処理は単独で処理を行う事はまれで、膨張・収縮を繰り返し処理を行う場合が多くあります。
とくに、同じ回数分だけ膨張して収縮する処理をクロージング(Closing)

同じ回数分だけ収縮して膨張する処理をオープニング(Opening)とよびます。

 

元画像 クロージング
元画像 オープニング

 

この処理例を見てもわかるように、オープニング、クロージング処理では小さいパターン細いパターンの除去を行います。

 

さらに、クロージングした画像から元画像を差し引いた処理をブラックハット(Black-Hat)
元画像からオープニングした画像を差し引いた処理をトップハット(Top-Hat)といいます。

 

元画像 ブラックハット
元画像 トップハット

 

これらの処理を見ても分かるように、膨張・収縮処理を用いることで、欠け、断線、ホコリやゴミといった欠陥の検査や、細い線状のパターンの文字や線などの抽出への応用が考えられます。

 

画像処理アルゴリズムへ戻る

 

アンシャープマスキング(鮮鋭化フィルタ)

アンシャープマスキング【unsharp masking】とは、画像のボヤけた輪郭を強調するフィルタ処理です。

アンシャープマスキング処理前 アンシャープマスキング処理後

 

【処理アルゴリズム】

元画像データ
元画像データの平滑化処理を行い、元画像から平滑化データを差し引きます。
差し引いた分を元画像に上乗せします

 

この処理により、平滑化処理でボヤけてしまう分だけ、逆に輪郭をくっきりさせる事ができます。

 

この処理をカーネルで表現すると

となります。
kは任意レートで、値が大きいほど、輪郭が鮮明になる効果が高くなります。
(後半の元画像から平滑化画像を引いている部分の処理はラプラシアンフィルタと呼ばれます。)

 

今回は3×3のマスクサイズを例にとって紹介しましたが、マスクサイズを大きくするとマスクサイズよりも小さなパターンを強調する効果を得る事もできます。

 

 

画像処理アルゴリズムへ戻る

 

エッジ抽出(Sobel,Prewittなど)

下図のように画像の輪郭(エッジ)を抽出する方法を紹介します。

エッジ抽出前 エッジ抽出後

上図はソーベルフィルタの例です。

 

微分フィルタ

画像の輝度値に対して、隣り合う画素の輝度差が大きいほど、画像のエッジだと考えることができ、
画像の輝度値データ 画像輝度値データに対して

 

横方向の差分(偏微分)以下のカーネルで求めることができる。

0 0 0
0 -1 1
0 0 0

 

縦方向の差分(偏微分)下のは以下のカーネルで求めることができる。

0 1 0
0 -1 0
0 0 0

または

0 0 0
0 -1 0
0 1 0

※画像データの並びがボトムアップかボトムダウンかにより使い分けます。

 

このとき、エッジの強さは

 

エッジの傾きは

として求めることができる。

 

プリューウィットフィルタ(Prewitt filter)

微分フィルタではノイズに弱いため、横方向(または縦方向)のエッジを計算してから
縦方向(または横方向)に平均化処理を行う手法です。カーネルには以下とおりです。

 

横方向の差分

-1 0 1
-1 0 1
-1 0 1

 

縦方向の差分

1 1 1
0 0 0
-1 -1 -1

または

-1 -1 -1
0 0 0
1 1 1

 

ソーベルフィルタ(Sobel filter)

プリューウィットフィルタではノイズを除去するのに平滑化処理を行っていましたが、ガウシアン平滑化処理を行ったのがソーベルフィルタとなります。カーネルには以下とおりです。

 

横方向の差分

-1 0 1
-2 0 2
-1 0 1

 

縦方向の差分

1 2 1
0 0 0
-1 -2 -1

または

-1 -2 -1
0 0 0
1 2 1

 

画像処理アルゴリズムへ戻る