Hough変換との接点
38年前のはなし。
学生から社会人になって職場に配属されたわけですが、１年先輩は画像処理チップ設計のため長期出張中だという。その画像処理というのが今回テーマのHough変換です。
当時はリアルタイムにHough変換できるシステムが世の中にない時代。
クルマを制御するには、ビデオレート程度のリアルタイム性が必要ってことで、専用チップを開発していたわけです。１年目の新人にそんなことを任せる、とんでもない所に配属されたと感じました。
Hough変換については、直線は極座標系で、原点からの垂線の距離と角度で一意に特定されるので、この方法を使って画像中の直線成分を抽出する仕組みと記憶していました。しかし、これまで具体的に計算したり試してみたりしたことはありません。
今回は、制御工学ではありませんが、「Hough変換の基本的な仕組み」についてExcelを使って確認してみるという試みです。
Hough変換の基本的な仕組み
Hough変換は初期は画像の中の直線成分の抽出に使われていましたが、現在では一般化され、円や図形の抽出にも広く使われています。
このブログでは、ある平面上にある複数の点が、同一直線上にあることをいかにして検出するかについて、その基本的な仕組みにを説明します。
極座標への変換
直線は極座標系で原点からの垂線の距離と角度で一意に特定されます。直線は原点からの垂線の距離 R と垂線の角度 θ を使うと、次のように表せます。
\[R=x\cos\theta+y\sin\theta \]
\(R\)：原点から直線までの垂線距離
\(\theta\) ：その垂線の角度
(R,θ)パラメータ空間（Hough空間）の生成
極座標では直線はRとθで一意に特定できるのですが、それを複数の点が直線状に並んでいるかの判断にどうやって使うのでしょうか。
複数の点が直線状に並んでいれば、それらの点は共通の（ R, θ ）を持ちますが、その直線の（ R, θ ）をいかにして知ることができるでしょうか。
具体的には、各点について、θを0～180degスキャンし、Rを計算してみるという方法になります。ある平面上にある3点について、θを0～180degスキャンしてRを計算した結果が以下のグラフです。
150degを過ぎた辺りに交点があります。
この交点になった（ R, θ ）で特定される直線上に3点が並んでいることを表しています。
投票による直線の抽出
平面上に点が多くなるとその組み合わせの数だけ交点ができてしまいます。多くの点が並んでいる直線は、その交点で交わる数も多くなるので、その交点で交わる数をカウントします。
先程の平面上にある3点について、θと求めたRの組み合わせ（ R, θ ）に投票した結果が下記のグラフです。
このようにすると3点のR値が交差している（ R, θ ）を特定でき、点列の直線成分を抽出することができます。
EXCELによるHough変換の確認
ダウンロード
こちら→リンク からExcelファイルをダウンロードできます。
ファイルの説明
シート【HoughTransform】に40×40セルのHough変換対象領域と、（ R, θ ）パラメータ空間グラフ、投票結果のグラフがあります。
40×40セルのHough変換対象領域には、1～40の値を格納したセルを並べています。
上部のスイッチ【Hough変換実行】を押すと、40×40セルの対象領域のHough変換が実行され、投票結果の上位10番までの（ R, θ ）で特定された直線が描画されます。
スイッチ【ラインクリア】を押すと計算結果と描画されたラインを消すことができます。
1～40の値を格納したセルの配置を変更することで、セルの配置に応じた直線成分の抽出結果を確認することができます。
VBAの説明
スイッチ【Hough変換実行】を押すと以下を実行します。
対象領域40×40セルの左下(1,1)から全セルについて以下を実行
- 1～40の値を格納したセルがあった場合には、θを0～180deg間2deg毎にRを計算
\(R=x\cos\theta+y\sin\theta \) \(x,y\)はセルのNo. (1～40)です。 - 結果を（R,θ）パラメータ空間の表に格納
- パラメータ空間の（R,θ）の該当投票値に1を加算
- 投票値が上位10番の範囲であれば上位10番を更新する
対象領域40×40セルの処理が完了した後、投票結果の上位10番までの（ R, θ ）で
特定された直線を描画する。
この処理はVBAで実装しています。VBAのコードはこんな感じです。
'////////////////////////////////////////////////// ' ' Hough変換処理 ' ' シートHoughTransformの40×40の対象セル領域について ' 1～40の値が格納されたセルの並びの直線成分を抽出する ' '////////////////////////////////////////////////// Const sizeofX As Integer = 40 Const sizeofY As Integer = 40 Const offsetX As Integer = 1 Const offsetY As Integer = 9 Const Pi As Double = 3.14159265358979 Dim angleRadian As Double Dim Radius As Double Dim Vote As Integer Dim AxisDeg As Integer Dim AxisElement As Integer Dim AxisRadius As Integer Dim VoteMinimum As Integer Dim NoMinimum As Integer Type voteTyp Vote As Integer AxisDeg As Integer AxisRadius As Integer Radius As Double angleRadian As Double End Type Dim VoteTop10(9) As voteTyp Type indexCrosspointTyp x As Integer y As Integer End Type Dim indexCrosspoint(3) As indexCrosspointTyp Dim ws As Worksheet Dim shp As Shape Dim rng As Range '////////////////////////////////////////////////// ' Hough変換処理 '////////////////////////////////////////////////// Sub HoughTransform() Dim indexX, indexY As Integer Dim elementNo As Integer Dim angleDegree As Integer Dim intRadius As Integer 'MsgBox "HoughTransform" '////////////////////////////////////////////////// ' 描画ラインをクリア '////////////////////////////////////////////////// Call Module1.ClearLines VoteMinimum = 0 NoMinimum = 0 For indexY = 1 To sizeofY For indexX = 1 To sizeofX elementNo = Cells(sizeofY + offsetY - indexY + 1, indexX + offsetX) ' MsgBox elementNo If elementNo <> 0 Then AxisElement = elementNo + offsetY '////////////////////////////////////////////////// ' 原点からの垂線の角度θを2degごとに更新 '////////////////////////////////////////////////// For angleDegree = 0 To 178 Step 2 AxisDeg = (angleDegree / 2) + 55 ' BC列 -> 55 angleRadian = angleDegree * Pi / 180 ' degをradに変換 '////////////////////////////////////////////////// ' 垂線距離 Radius 算出 ' R = x cosθ + y sinθ '////////////////////////////////////////////////// Radius = (indexX * Cos(angleRadian)) + (indexY * Sin(angleRadian)) Cells(AxisElement, AxisDeg) = Radius '////////////////////////////////////////////////// ' 投票 ～交点数の抽出 '////////////////////////////////////////////////// intRadius = CInt(Radius) AxisRadius = 114 - intRadius ' 114行 Radius = 0 Vote = Cells(AxisRadius, AxisDeg) Vote = Vote + 1 Cells(AxisRadius, AxisDeg) = Vote '////////////////////////////////////////////////// ' 投票Top10を更新 '////////////////////////////////////////////////// Call Module1.recordTop10 Next angleDegree End If Next indexX Next indexY '////////////////////////////////////////////////// ' 検出したラインを描画 '////////////////////////////////////////////////// Call Module1.calCrosspoint ' DEBUG ' For i = 0 To 9 ' Debug.Print VoteTop10(i).Vote ' Next i End Sub '////////////////////////////////////////////////// ' 描画ラインをクリアする '////////////////////////////////////////////////// Sub ClearLines() '////////////////////////////////////////////////// '対象のシートを設定 '////////////////////////////////////////////////// Set ws = ThisWorkbook.Sheets("HoughTransform") '////////////////////////////////////////////////// ' 計算結果をクリア '////////////////////////////////////////////////// ws.Range("BC10:EN49").Clear ws.Range("BC54:EN174").Clear '////////////////////////////////////////////////// ' 投票TOP10をクリア '////////////////////////////////////////////////// For i = 1 To 9 VoteTop10(i).Vote = 0 Next i '////////////////////////////////////////////////// '対象の範囲を設定 '////////////////////////////////////////////////// Set rng = ws.Range("B10:AO49") '////////////////////////////////////////////////// 'シート内の全ての図形 '////////////////////////////////////////////////// For Each shp In ws.Shapes ' 図形の位置が範囲内にあるか確認 If Not Intersect(shp.TopLeftCell, rng) Is Nothing Then shp.Delete ' 図形を削除 End If Next shp End Sub '////////////////////////////////////////////////// ' 投票値Top10を記録・更新する '////////////////////////////////////////////////// Sub recordTop10() ' 初期化 flgComplete = 0 If (Vote > VoteMinimum) Then '////////////////////////////////////////////////// ' 同一要素の更新確認 '////////////////////////////////////////////////// For i = 0 To 9 If (VoteTop10(i).AxisDeg = AxisDeg And VoteTop10(i).AxisRadius = AxisRadius) Then ' 同一要素 VoteTop10(i).Vote = Vote VoteTop10(i).AxisDeg = AxisDeg VoteTop10(i).AxisRadius = AxisRadius VoteTop10(i).Radius = Radius VoteTop10(i).angleRadian = angleRadian flgComplete = 1 Exit For End If Next i '////////////////////////////////////////////////// ' 未更新の場合、最小値のエリアに記録する '////////////////////////////////////////////////// If (flgComplete = 0) Then VoteTop10(NoMinimum).Vote = Vote VoteTop10(NoMinimum).AxisDeg = AxisDeg VoteTop10(NoMinimum).AxisRadius = AxisRadius VoteTop10(NoMinimum).Radius = Radius VoteTop10(NoMinimum).angleRadian = angleRadian End If '////////////////////////////////////////////////// ' 最小値の更新 '////////////////////////////////////////////////// VoteMinimum = VoteTop10(0).Vote NoMinimum = 0 For i = 1 To 9 If (VoteTop10(i).Vote < VoteMinimum) Then VoteMinimum = VoteTop10(i).Vote NoMinimum = i End If Next i End If End Sub '////////////////////////////////////////////////// ' 検出したラインを描画するため、x=1, x=40, y=1, y=40 ' の交差ポイント(セル)を求める '////////////////////////////////////////////////// Sub calCrosspoint() Dim r, theta As Double Dim x, y As Integer Dim nCrosspoint As Integer For i = 0 To 9 r = VoteTop10(i).Radius theta = VoteTop10(i).angleRadian nCrosspoint = 0 If (Cos(theta) <> 0) Then '////////////////////////////////////////////////// ' 2degごとの計算のため三角関数の最小値は 0以外は 0.0348 ' xは領域サイズが40までであればIntegerでもオーバーフローしない '////////////////////////////////////////////////// '////////////////////////////////////////////////// ' Crosspoint y=1 '////////////////////////////////////////////////// x = CInt((r - Sin(theta)) / Cos(theta)) ' Debug.Print "Crosspoint y=1 x=" & x If (x >= 1 And x <= sizeofX) Then indexCrosspoint(nCrosspoint).x = x indexCrosspoint(nCrosspoint).y = 1 ' Crosspoint y=1 nCrosspoint = nCrosspoint + 1 End If '////////////////////////////////////////////////// ' Crosspoint y=40 '////////////////////////////////////////////////// x = CInt((r - sizeofY * Sin(theta)) / Cos(theta)) ' Debug.Print "Crosspoint y=40 x=" & x If (x >= 1 And x <= sizeofX) Then indexCrosspoint(nCrosspoint).x = x indexCrosspoint(nCrosspoint).y = sizeofY ' Crosspoint y=40 nCrosspoint = nCrosspoint + 1 End If End If If (Sin(theta) <> 0) Then '////////////////////////////////////////////////// ' Crosspoint x=1 '////////////////////////////////////////////////// y = CInt((r - Cos(theta)) / Sin(theta)) ' Debug.Print "Crosspoint x=1 y=" & y If (y >= 1 And y <= sizeofY) Then indexCrosspoint(nCrosspoint).x = 1 ' Crosspoint x=1 indexCrosspoint(nCrosspoint).y = y nCrosspoint = nCrosspoint + 1 End If '////////////////////////////////////////////////// ' Crosspoint x=40 '////////////////////////////////////////////////// y = CInt((r - sizeofX * Cos(theta)) / Sin(theta)) ' Debug.Print "Crosspoint x=40 y=" & y If (y >= 1 And y <= sizeofY) Then indexCrosspoint(nCrosspoint).x = sizeofX ' Crosspoint x=40 indexCrosspoint(nCrosspoint).y = y nCrosspoint = nCrosspoint + 1 End If End If '////////////////////////////////////////////////// ' 検出したラインを描画 ' 同一の交差セルを2回選択する場合があるためチェック '////////////////////////////////////////////////// If (nCrosspoint >= 2) Then 'Call Module1.DisplayLine(1, 1, 40, 40) 'Call Module1.DisplayLine(indexCrosspoint(0).x, indexCrosspoint(0).y, indexCrosspoint(1).x, indexCrosspoint(1).y) If (indexCrosspoint(0).x <> indexCrosspoint(1).x Or indexCrosspoint(0).y <> indexCrosspoint(1).y) Then Call Module1.DisplayLine(indexCrosspoint(0).x, indexCrosspoint(0).y, indexCrosspoint(1).x, indexCrosspoint(1).y) ElseIf (nCrosspoint > 2) Then Call Module1.DisplayLine(indexCrosspoint(0).x, indexCrosspoint(0).y, indexCrosspoint(2).x, indexCrosspoint(2).y) End If End If Next i End Sub '////////////////////////////////////////////////// ' 検出したラインを描画 '////////////////////////////////////////////////// Sub DisplayLine(ByVal x1, y1, x2, y2 As Integer) Dim px1, px2, py1, py2 As Integer Dim indexX1, indexY1, indexX2, indexY2 As Integer indexX1 = x1 + offsetX '始点の横セルNo indexY1 = sizeofY + offsetY + 1 - y1 '始点の縦セルNo indexX2 = x2 + offsetX '終点の横セルNo indexY2 = sizeofY + offsetY + 1 - y2 '終点の縦セルNo px1 = Cells(indexY1, indexX1).Left + Cells(indexY1, indexX1).Width / 2 '始点の横ポイント py1 = Cells(indexY1, indexX1).Top + Cells(indexY1, indexX1).Height / 2 '始点の縦ポイント px2 = Cells(indexY2, indexX2).Left + Cells(indexY2, indexX2).Width / 2 '終点の横ポイント py2 = Cells(indexY2, indexX2).Top + Cells(indexY2, indexX2).Height / 2 '終点の縦ポイント ' px1 = Cells(10, 2).Left + Cells(10, 2).Width / 2 '始点の横ポイント ' py1 = Cells(10, 2).Top + Cells(10, 2).Height / 2 '始点の縦ポイント ' ' px2 = Cells(49, 41).Left + Cells(49, 41).Width / 2 '終点の横ポイント ' py2 = Cells(49, 41).Top + Cells(49, 41).Height / 2 '終点の縦ポイント ' ' ws.Shapes.AddLine px1, py1, px2, py2 With ws.Shapes.AddLine(px1, py1, px2, py2).Line 'ラインを描画 .ForeColor.RGB = vbBlue End With End Sub
おわりに
「Hough変換の基本的な仕組み」についてExcelで確認してみました。
Hough変換の基本的なアイデアは、1962年に Paul Hough によって特許取得されたものですが、現在では多くのクルマに搭載され運転支援に使われている技術になっています。
最後までお読みいただき、ありがとうございました。
0 件のコメント:
コメントを投稿