2007. 6. 14.
石立 喬
Visual C++ 2005 Express Editionを用いた易しい画像処理（１５）
―――― ハフ（Hough）変換を用いて、直線と円弧を検出する ――――
画像処理の最終的な目的の一つに物体の認識がある。物体認識の前処理には、直線部分や円形部分を検出する方法が良く用いられる。直線や円弧を検出する方法の代表的なものとして、ハフ（Hough）変換が知られている。
概要
画像からエッジを検出し（必要があれば細線化し）、二値化された画像から直線や円弧を抽出するのにハフ変換を用いる。ハフ変換で円を検出する実例は、従来あまり紹介されていない。その理由の一つには、三次元の大量のメモリを必要とし、計算時間が長くかかることが考えられる。しかし、最近のコンピュータを用いると比較的簡単にできる。
一般に知られているハフ変換のアルゴリズムを改良した点は次の通りである。
１） 垂直に近い直線を逆変換すると無限大の係数が現れて逆変換が困難なので、ｘを変化させてｙを描くだけでなく、ｙを変化させてｘを描くことも加えた。
２） 検出すべき真の直線や円弧の近傍に発生する擬似的な直線や円弧が同時に検出されるのを避けるため、候補が一つ選ばれると、その近傍を将来検出される候補から除去した。
ハフ（Hough）変換とは
ハフ変換は、1962年にP.V.C. Houghによって米国特許が取られている古い技法であるが、最近のコンピュータの高速化とメモリ容量の増大により、見直されている。ハフ変換は、直線の検出や円弧の検出に用いられる。たとえば、直線の検出の場合、元になる直角座標上の点（x,y）を角度θと原点からの距離ρの二次元空間に変換し、角度θと距離ρ毎に、その個数（すなわち、ピクセルの数）をメモリ上に加算してゆく（Accumulator
Cellなどと呼ばれるものに…）。個数が最大になった組み合わせ（角度θと距離ρの）を元の直角座標に戻したのが、最も直線らしい点の集まり（最も多くの点が並んでいる）である。個数を下げてゆくと、次の候補が順次得られる。角度θと距離ρを細かく分けると、精度が上がるが、計算時間が長くなり、メモリ容量も増えるのが欠点である。円の検出の場合には、元になる直角座標上の点（i,j）を、円の中心点（x,y）と半径ｒの三次元空間に変換し、同様の処理を行う。
１）直線の検出
図１の左の図に示すように、点（x,y）を通る任意の直線は、その直線に交わる原点からの垂線の長さρと、その垂線の成す角度θで表すことができる。点（x,y）を通る直線は多数あり、点（x,y）が決まると、θとρの多数の組み合わせが得られる。もし、別の点が同じ直線状にあるとすると、二つの点は、θとρの同じ組み合わせを共有することになる。θとρの同じ組み合わせが多数あれば、元になる点は、同じ直線を共有すると言うよりは、「そこに直線がある」ことを示しているといえる。
図１の右の図は、θからρを求める式の原理を示す。
図１ 直線を検出するためのハフ変換式の原理
直角座標上で、原点からの距離がρでｘ座標軸と角度θをなす線分と垂直な直線上に点（x,y）が存在すると仮定すると、
の関係がある。上の式を用いて、直角座標上の点（x,y）を新しい二次元空間（θ,
ρ）上に変換すると、直角座標上の一点はθ-ρ空間上の一本の曲線に対応する。直角座標上の点が多数あると、θ-ρ空間上に多数の曲線が得られる。θ-ρ空間上でこれらの曲線が共有する点があれば、それは元のｘ-ｙ直角座標上では一本の線上に並ぶことになる。
逆変換は、
および
により、xとyを変化して直線を描く。y=f(x)の式は、sinθが0に近い場合にはyの間隔が広くなって適切でなく、x=f(y)の式は、cosθが0に近い付近ではxの間隔が広くなるので、両者を併用する。
２） 円弧の検出
半径がｒで、中心点が（x,y）の円弧の上に、座標上の点（i、j）が存在すると仮定すると、
の関係が成り立つ。上の式を用いて、直角座標上の点（i,j）を新しい三次元空間（x,y,r）上に変換すると、直角座標上の一点は三次元空間（x,y,r）の１枚の面に対応する。直角座標上の点が多数あると、x-y-r空間に多数の曲面が得られる。それらの曲面が共有する点があれば、それは元の直角座標上では、点（x,y）を中心にした半径rの一個の円弧上に並ぶことになる。
逆変換は、数式を用いるまでもなく、点（x,y）を中心に半径rの円を描けばよい。
直線を検出するハフ変換プログラムの構成
１）原画像を読み込んで、Bitmapクラスのbmap_srcとする。
２） 角度thetaは、半周のπラジアンを1024に分割する(THETA_RESOLUTION=1024)ので、から1023までを使用する。
３）距離rhoは、-500から499まで使えるようにする(RHO_RESOLUTION=1000)が、配列の添字にマイナスは使えないので、RHO_RESOLUTION/2(=500)を加算して、添字は0から999までとする。
４）角度thetaと距離rhoの組み合わせが発生する頻度を蓄積する直線検出用カウンタとして、counter[THETA_RESOLUTION,RHO_RESOLUTION]を用意する。
５）高速化のために、あらかじめ三角関数のテーブル（数表）を計算して用意する。
６）ハフ変換を実行する。具体的には、原画像bmap_srcをスキャンして、もし黒のピクセルがあれば、その点（i,j）を新しい二次元空間（θ,ρ）に変換しながら、counter[theta,rho]に頻度を累積する。
７) 以下のハフ逆変換は、LNUMBER_MAX(=15)で指定された回数（またはbreakで中断されるまで）だけ行う。
A） counter[theta,rho]の内容が最大になるtheta_maxとrho_maxと、その時のcounterの値counter_maxを求める。
B）theta_maxとrho_maxの近傍のcounter[theta,rho]をゼロにし、将来の候補から除外する。これは、疑似の直線が多数現れることを防ぐためである。除外の対象となる値を、theta_cut、rho_cutと呼ぶ。theta_cutがマイナスになる場合やTHETA_RESOLUTIONの範囲を超える場合には、補正を行い、rho_cutに対しては、符号の反転を行う。
C）theta_maxとrho_maxを、直角座標上の直線に変換して（xを変えながらyを描く、垂直および垂直に近い線は描けないので、theta_max!=0で事前に判定する）、Bitmap画像bmap_hに書き込む。
D）同様に、theta_maxとrho_maxを、直角座標上の直線に変換し（yを変えながらxを描く、水平および水平に近い線は描けないので、theta_max!=THETA_RESOLUTION/2で事前に判定する）、bmap_hに書き込む。
E）counter[theta,rho]のカウント数count_maxが一定値以下になれば、ノイズである恐れがあるので、breakで中止する。
８) Form1_Paint()で、ハフ逆変換結果が書き込まれたbmap_hを描画する。
円弧を検出するハフ変換プログラムの構成
１）原画像を読み込んで、Bitmapクラスのbmap_srcとする。
２） 円弧の中心と仮定する点(x,y)から、原画像の黒ピクセルが存在する位置までのX（横）方向の距離の最大をX_MAX(=80)、Y（縦）方向の距離の最大をY_MAX(=80)とする。このことは、これ以上大きい円弧は検出しないことを意味する。円弧の半径の最大値R_MAX(=115)は、X_MAXの二乗とY_MAXの二乗の平方根をとって設定しておく。
３）円弧の中心点(x,y)を表すxおよびy、半径radius、円弧を構成するピクセルの数countをまとめて、一つのユーザ定義クラスを作成し、その配列を用意する。
４）高速化のために、あらかじめ二乗和の平方根のテーブル（数表）を計算して用意する。
５）ハフ変換を実行する。原画像bmap_srcをスキャンして、もし黒のピクセルがあれば、その点（i,j）に対して、想定される円弧の中心点(x,y)をスキャンし、想定される半径radiusを求めて、空間（x,y,radius）に変換し、counter_r[x,y,radius]に頻度を累積する。
６) 以下のハフ逆変換は、RNUMBER_MAXで指定された回数（またはbreakで中断されるまで）だけ行う。
A） counter_r[x,y,radius]の内容が最大になるx_maxとy_maxとradius_max、および、その時のcounterの値counter_maxを求める。
B）ユーザ定義クラスCircleの配列circle1[n]に、上記データを設定する。
C）x_max、y_max、radius_maxの近傍のcounter_r[x,y,radius]をゼロにし、将来の候補から除外する。除外の対象となる値を、x_cut、y_cut、r_cutと呼ぶ。これは、疑似の直線が多数現れることを防ぐためである。
D）counter[theta,rho]のカウント数count_maxが一定値以下になれば、breakで中止し、その時までに検出された円弧の数nをnumber_rとして設定する。
７) Form1_Paint()で、number_rで与えられた数だけ、Circleクラスの配列circle1[n]にしたがって円を描く。
直線を検出するハフ変換プログラム
Bitmap^ bmap_src;
Bitmap^ bmap_h;
static int WIDTH=320,HEIGHT=240;
private: System::Void Form1_Load(System::Object^ sender, System::EventArgs^ e)
{
int i,j,k,n;
int x,y;
Color color1;
int counter_max;
//原画像を読み込む
bmap_src=gcnew Bitmap("A:/sample.gif");
int THETA_RESOLUTION=1024; //thetaの範囲は、0から1023まで
int RHO_RESOLUTION=1000; //rhoの範囲は、-500からまで499まで
int LNUMBER_MAX=15; //検出する直線の数の最大は15個
double pai=Math::PI/THETA_RESOLUTION;
int theta,rho;
int theta_max,rho_max,count;
int theta_cut,rho_cut;
//直線検出のためのカウンタを用意する
array<short,2>^ counter=gcnew array<short,2>(THETA_RESOLUTION,RHO_RESOLUTION);
//ハフ逆変換（直線）の結果を表示するbmap_hをインスタンス化する
bmap_h=gcnew Bitmap(WIDTH,HEIGHT);
//sinとcosのテーブルを用意する
array<double>^ sn=gcnew array<double>(THETA_RESOLUTION);
array<double>^ cs=gcnew array<double>(THETA_RESOLUTION);
for(i=0;i<THETA_RESOLUTION;i++){
sn[i]=Math::Sin(pai*i);
cs[i]=Math::Cos(pai*i);
}
//ハフ変換（直線）の実行
for(i=0;i<WIDTH;i++)
for(j=0;j<HEIGHT;j++){
color1=bmap_src->GetPixel(i,j);
if(color1.R<128){ //黒であった
for(theta=0;theta<THETA_RESOLUTION;theta++){
rho=(int)(i*cs[theta]+j*sn[theta]+0.5);
counter[theta,rho+500]++;
}
}
}
//ハフ逆変換（直線）の実行
for(n=0;n<LNUMBER_MAX;n++){
//counterが最大になるtheta_maxとrho_maxを求める
counter_max=0;
for(theta=0;theta<THETA_RESOLUTION;theta++)
for(rho=-RHO_RESOLUTION/2;rho<RHO_RESOLUTION/2;rho++)
if(counter[theta,rho+RHO_RESOLUTION/2]>counter_max){
counter_max=counter[theta,rho+RHO_RESOLUTION/2];
theta_max=theta;
rho_max=rho;
count=counter_max;
}
//counter[theta_max,rho_max]の近傍を0にする
for(i=-20;i<=20;i++)
for(j=-10;j<=10;j++){
theta_cut=theta_max+i;
rho_cut=rho_max+j;
if(theta_cut<0){
theta_cut+=THETA_RESOLUTION;
rho_cut=-rho_cut;
}
else if(theta_cut>THETA_RESOLUTION-1){
theta_cut-=THETA_RESOLUTION;
rho_cut=-rho_cut;
}
counter[theta_cut,rho_cut+RHO_RESOLUTION/2]=0;
//削除する
}
//逆ハフ変換（直線）した結果を表示する
if(theta_max!=0){ //垂直の線を除く
for(x=0;x<WIDTH;x++){
y=(int)((rho_max-x*cs[theta_max])/sn[theta_max]+0.5);
if(y>=HEIGHT || y<0) continue;
bmap_h->SetPixel(x,y,Color::Red);
}
}
if(theta_max!=THETA_RESOLUTION/2){ //水平の線を除く
for(y=0;y<HEIGHT;y++){
x=(int)((rho_max-y*sn[theta_max])/cs[theta_max]+0.5);
if(x>=WIDTH || x<0) continue;
bmap_h->SetPixel(x,y,Color::Red);
}
}
//直線を形成するピクセルが60個未満になったら表示しない
if(count<60) break;
}
}
private: System::Void Form1_Paint(System::Object^ sender, System::Windows::Forms::PaintEventArgs^
e) {
Graphics^ g=e->Graphics;
g->DrawImage(bmap_src,10,10,WIDTH,HEIGHT); //原画像を左に表示する
g->DrawImage(bmap_src,340,10,WIDTH,HEIGHT); //原画像を右にも表示する
g->DrawImage(bmap_h,340,10,WIDTH,HEIGHT); //逆ハフ変換の結果を原画像の上に表示する
}
円弧を検出するハフ変換プログラム
#pragma once
#include "Circle.h"
static int WIDTH=320,HEIGHT=240;
array<Circle^>^ circle1;
int number_r;
private: System::Void Form1_Load(System::Object^ sender, System::EventArgs^ e)
{
int X_MAX=80,Y_MAX=80,R_MAX=115; //検出する円弧の最大半径
int RNUMBER_MAX=10; //検出する円弧の数の最大は10個
int i,j,k,n;
int x,y,radius;
Color color1;
int counter_max;
int x_max,y_max,radius_max,count_max;
int x_cut,y_cut,r_cut;
//原画像を読み込む
bmap_src=gcnew Bitmap("A:/sample1.gif");
//円弧検出のためのカウンタを用意する
array<short,3>^ counter1=gcnew array<short,3>(WIDTH,HEIGHT,R_MAX);
//検出された円弧のデータを入れておくCircleクラスの配列を用意する
circle1=gcnew array<Circle^>(RNUMBER_MAX);
//斜線長（二乗和の平方根）のテーブルを用意する
array<int,2>^ diagonal=gcnew array<int,2>(X_MAX,Y_MAX);
for(i=0;i<X_MAX;i++)
for(j=0;j<Y_MAX;j++)
diagonal[i,j]=(int)(Math::Sqrt(i*i+j*j)+0.5);
//ハフ変換（円弧）の実行（かなり時間がかかる）
for(i=0;i<WIDTH;i++)
for(j=0;j<HEIGHT;j++){
color1=bmap_src->GetPixel(i,j);
if(color1.R<128){ //黒であった
for(x=0;x<WIDTH;x++)
for(y=0;y<HEIGHT;y++){
if(Math::Abs(x-i)>=X_MAX || Math::Abs(y-j)>=Y_MAX) continue;
radius=diagonal[Math::Abs(x-i),Math::Abs(y-j)];
counter_r[x,y,radius]++;
}
}
}
//ハフ逆変換（円弧）の実行
for(n=0;n<RNUMBER_MAX;n++){
//counterが最大になるx,y,radiusを求める
counter_max=0;
for(x=0;x<WIDTH;x++)
for(y=0;y<HEIGHT;y++)
for(radius=0;radius<R_MAX;radius++){
if(counter_r[x,y,radius]>counter_max){
counter_max=counter_r[x,y,radius];
x_max=x;
y_max=y;
radius_max=radius;
count_max=counter_max;
}
}
if(count_max<80){ //count_maxが80未満になったら、円弧として認識しない
number_r=n;
break;
}
circle1[n]=gcnew Circle(x_max,y_max,radius_max,count_max);
//counter_r[x_max,y_max,radius_max]の近傍を0にする
for(i=-5;i<=5;i++)
for(j=-5;j<=5;j++)
for(k=-5;k<=5;k++){
x_cut=x_max+i;
y_cut=y_max+j;
r_cut=radius_max+k;
if(x_cut<0 || x_cut>=WIDTH) continue;
if(y_cut<0 || y_cut>=HEIGHT) continue;
if(r_cut<0) continue;
counter_r[x_cut,y_cut,r_cut]=0; //削除する
}
}
}
private: System::Void Form1_Paint(System::Object^ sender, System::Windows::Forms::PaintEventArgs^
e) {
Graphics^ g=e->Graphics;
int x,y,r;
g->DrawImage(bmap_src,10,10,WIDTH,HEIGHT); //原画像を左に表示する
g->DrawImage(bmap_src,340,10,WIDTH,HEIGHT); //原画像を右にも表示する
for(int n=0;n<number_r;n++){
x=circle1[n]->X;
y=circle1[n]->Y;
r=circle1[n]->Radius;
g->DrawEllipse(Pens::Red,340+x-r,10+y-r,2*r,2*r);
}
}
Circleクラス
一応、クラスは作ってみたが、本格的なものではなく、手ごろに簡単に作ることを目的とした。たとえば、Circleクラスのメンバ変数は、本来public;
ではなく、private; にして、外部からの直接のアクセスを防ぐのが正攻法であるが、アクセス用のメンバ関数を記述するのが面倒なので、public;
のままにしてある。 X=x; の代わりに、this->X=x; と書けば、趣旨が明確になるが、これも省略した。クラスの作り方と、propertyキーワード（今回は使用しなかったが…）については、付録を参照のこと。
Circle.hの内容
#pragma once
ref class Circle
{
public:
int X;
int Y;
int Radius;
int Count;
public:
Circle(void);
Circle(int x,int y,int radius,int count);
};
Circle.cppの内容
#include "StdAfx.h"
#include "Circle.h"
Circle::Circle(void)
{
}
Circle::Circle(int x, int y, int radius, int count)
{
X=x;
Y=y;
Radius=radius;
Count=count;
}
直線を検出するハフ変換の結果
図２は、サンプル画像に対する直線の検出の結果を示す。完全に水平であるべき直線が、計算誤差により若干傾斜している。
図２ ハフ変換による直線の検出の結果
円弧を検出するハフ変換の結果
図３は、サンプル画像に対する円弧の検出の結果を示す。大きな円は半径がR_MAX以上なので、検出対象外である。小さい円も、構成するピクセルの数が規定値以下なので、対象外である。
図３ ハフ変換による円弧の検出の結果
付録1：クラスの作り方
折角のオブジェクト指向可能なプログラム開発環境があるのに、ユーザ定義クラスを使用しないのは、邪道と言うか、無謀というか、どちらにしても褒められた話ではない。しかし、実際には、ほとんど必要を感じさせられないし、かえって分かりにくくなる場合も多い。今回は、無理をして簡単なところからはじめてみたので、紹介する。
クラスの作り方を、クラス名Circleの例で説明すると、以下のようになる。
１） 統合開発環境のメニューから、「プロジェクト」→「クラスの追加」をクリックする。
２） 「クラスの追加」ダイアログボックスが開くので、すでに選択されている
カテゴリー欄------Visual C++ のC++
テンプレート欄----C++クラス
を確認して、「追加」をクリックする。
３） 「汎用C++クラスウイザード」ウインドウで、
クラス名------- Circle（入力する）
.hファイル----- (自動的に入る）
.cppファイル---（自動的に入る）
基本クラス----- 継承しなければ入力不要（空欄のまま）
アクセス------- publicになっている（一般的にはそのままで良い）
マネージ------- チェックを付ける（すでに付いているので、そのまま）
にして、「完了」をクリックする。
４） エディタウインドウが開いて、下記のように、Circle.h のスケルトン（コンストラクタのみ）が表示さ
れるので、適宜記述する。
#pragma once
ref class Circle
{
public:
Circle(void);
};
５） 「表示」→「ソリューションエクスプローラ」をクリックしてソリューションエクスプローラを表示さ
せ、Circle.cppをダブルクリックすると、下記のように、Circle.cppのスケルトン（コンストラクタの
み）が表示されるので、適宜記述する。
#include “StdAfx.h”
#include “Circle.h”
Circle::Circle(void)
{
}
６） Form1.hに、#include “Circle.h”を記述する。
付録2：propertyキーワード
下記のように、property int と宣言すると、get()やset()のメソッドが自動的に生成される。単にint
のみでは、Form1.hからこれらのメソッドを使用できない。
ref class Circle
{
public:
property int X;
property int Y;
property int Radius;
property int Count;
public:
Circle(void);
Circle(int x,int y,int radius,int count);
};
このように設定されたフィールド（メンバ変数）は、
circle1[n]->X::set(x_max);
circle1[n]->Y::set(y_max);
circle1[n]->Radius::set(radius_max);
circle1[n]->Count::set(count);
でForm1.hから設定できる。
また、下記により、Form1.hから受け取ることができる。
x=circle1[n]->X::get();
y=circle1[n]->Y::get();
r=circle1[n]->Radius::get();