画像処理とC++を用いたアナログメーターの自動読み取り　アルゴリズム実装

# はじめに
学生時代に研究室が企業からアナログメーターの遠隔読み取りができないかと依頼がありました。画像処理について研究していた私に白羽の矢が立ってしまいメインの研究と並行して簡易なアルゴリズムを考案し実装しました。
アルゴリズムの内容をレポート形式にまとめています。

# 開発背景
アナログメーターの自動読み取りは意外と人気のあるタスクです。
様々な企業もシステムを作成しています。
特に最近だとAIと画像処理の組み合わせが目立ちます。

[富士通九州システムズ アナログメーター認識ソフトウェア](https://www.fujitsu.com/jp/group/kyushu/solutions/business-technology/image-monitoring/analog-meter-recognizer/index.html)

[日立製作所 古いアナログメーターでも後付けで自動読み取り、点検効率化サービスを提供へ](https://monoist.atmarkit.co.jp/mn/articles/1801/29/news059.html)

[こちら](https://www.fujitsu.com/jp/group/kyushu/resources/news/press-releases/2011/0411.html)にシステム開発のモチベーションが記載されています。主にデータの収集や目視による作業コストの削減、さらには異常検知が狙いです。

「現状、データ送信機能のないアナログメーターが設置されている工場や施設では、定期的に作業員が巡回し、目視によるメーター数値の読み取りを行い、検針作業をおこなっています。異常値発生時には、現場で常時、メーター監視する必要があり作業員の負荷も大きくなります。 メーターの検針作業を自動化するために、データ送信機能付きのメーターに置き換える方法もありますが、装置を停止させる必要があるうえ、工事等のコストが大きいという問題があります。」

以下、上記のタスクについて考案した提案手法のロジックと実験結果です。文章は堅めに作成していたのでである口調になっていますがあまり気にしないでください。

# 提案システムの詳細
提案システムの詳細について図1に示す．また，各処理について図を用いて説明する．本システムでは圧力計の画像から針が示す数値を読み取ることを目的とする．

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/42f99b02-fd07-7da6-5d37-228785111d2d.png)
図1 フローチャート

## 詳細強化フィルタによる入力画像の強調
本システムはリアルタイムに変化する針が示すアナログ値を算出する必要がある．そのため，処理速度を考慮する必要がある．初めに図1(a)の入力画像に対してランツォス法を用いて0.5倍にダウンサンプリングする．想定する入力画像の解像度は1920x1080である．そして，この解像度に対してダウンサンプリングを行う．本システムで用いるランツォス法はバイキュービック法やバイリニア法よりも画素を決定するためのグリッドが大きく精度が高い特徴がある．次に，コントラストの強調とノイズの影響を軽減するために入力画像に対して詳細強化フィルタ(Detail Enhance Filter)を用いてエッジを強調する．詳細強化フィルタは鮮鋭化に用いられるフィルタである．入力画像に対して，詳細強化フィルタにより強調した結果を同図(b)に示す．同図(b)の結果より，メータ領域内部のコントラストが強調されていることが分かる．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/839b6958-3e10-548e-9ab0-aacae2acad05.png)
(a) 入力画像(1920x1080)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/9b7d103b-425b-a03b-aa31-0ee8000957d7.png)
(b) 詳細強化フィルタ(960x540)
図1 入力画像と詳細強化フィルタによる強調

## 確率的ハフ変換を用いた円領域検出
次に入力画像からメータのおおよその位置を特定するために確率的ハフ変換を用いて円領域の検出を行う．入力画像に対してグレイスケール変換とキャニー法による二値化を行う．確率的ハフ変換では複数の円が検出される．本システムではメータの内側の目盛りを通る円領域が必要なため，ハフ変換により最も円領域と判断された領域のみを検出する．確率的ハフ変換による円領域の検出結果を図2(a)に示す．そして，検出結果の円の中心と半径から同図(b)のマスク画像を生成する．同図(b)は背景とメータ領域を分離する際に用いるマスク画像である．なお，マスク画像の半径は検出された円の半径に半径の3分の1を加算している．円は目盛りの内側に検出されるため，マスクの半径を大きく取得することで外側の目盛りをマスク内に取得する．今後，このマスク画像を用いることでマスク領域内に限定して処理を行うことが可能となり，処理速度と精度を向上させることができる．また，確率的ハフ変換による円検出で検出される円はメータの目盛りの内側の必要がある．メータの設置の際は，図2のようにメータの上下がwebカメラに収まるように設置する．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/037f476a-9c09-24b8-b0b0-e0dad5d2833c.png)
(a) 確率的ハフ変換
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/bcbc6118-eb3f-3977-abc1-ff33fd1bf1f6.png)
(b) 背景分離マスク
図2 円検出と生成するマスク画像

## 適応的二値化とラベリングによる目盛り領域検出
次に目盛りの領域を検出するためにラベリングによる領域分割を行う．そして，分割された複数の領域の横幅から目盛り領域を一意に検出する．詳細強化フィルタにより強調された画像に対して，ガウシアンフィルタを用いた適応的二値化を行う．適応的二値化の結果を図3(a)に示す．また，フィルタサイズは実験から51x51に設定している．次に，適応的二値化結果と背景を分離するマスク画像からマスク処理を行う．マスク処理後の結果を同図(b)に示す．これにより，背景のエッジを除去する．そして，同図(b)に対してラベリング処理を行う．最後にラベリングにより得られた複数の領域から最も領域の横幅が最大となる領域が目盛り領域となるため，最大横幅領域を検出する．ラベリングによる領域の検出結果を同図(c)に示す．

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/fefacf15-f610-b77d-5383-7e75b5261398.png)
(a) 適応的二値化
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/5203d222-32d5-8072-9fdd-9b7954595ef7.png)
(b) マスク処理
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/d44400ee-0e5e-3bea-c0d6-659063c7c25b.png)
(c) ラベリングによる最大幅領域
図3 適応的二値化とマスク処理とラベリングによる最大横幅領域の検出結果

本システムは針が目盛りと重なるメータを対象としている．そのため，対象の目盛り領域と針の領域は接続され，最大領域として検出される．今後，検出された領域から目盛りの始終点と針の芯線を求める．

## 細線化
次に，抽出された領域の結果から針の芯線を取得するために細線化を行う．図3(c)に対して細線化を行った結果を図4に示す．同図の結果から針の芯線と考えられる領域に画素が多く存在することが分かる．この特徴を用いることで針の芯線の角度を取得する．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/97de3984-9f94-13a5-46df-58b605fee199.png)
図4 細線化結果

同図の結果から針の芯線と考えられる領域に画素が多く存在することが分かる．この特徴を用いることで針の芯線の角度を取得する．

## 楕円近似
取得された目盛り領域はwebカメラから撮影されている．また，メータとwebカメラが必ずしも水平に設置できないため，領域が真円として取得できない場合がある．つまり，確率的ハフ変換によって検出された円の中心にずれが生じる場合がある．そのため，目盛り領域に対して楕円近似を行う．そして，メータの中心を楕円の中心に更新する．同図に対して楕円近似を行った結果を図5(a)に示す．同図の赤い点が得られた複数の楕円の中心を示している．また，円の中心を更新する際に用いる楕円は確率的ハフ変換で求められた円の中心とそれぞれの楕円の中心との距離を求め，最小距離である楕円の中心に変更する．同図では中央の楕円Aの中心座標がメータの中心座標に変更される．
![a.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/2b7ec96e-66fa-4059-3f0e-d7af0838a624.png)
(a)	楕円近似
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/59ff9991-1bb8-2d61-0522-af3c67193d41.png)
(b) ベクトルの例
図5 楕円近似結果とベクトルの例

## 針の芯線検出
次に，メータの針を検出する．図4の細線化結果から楕円の中心座標を持つ任意のベクトルを設定する．任意のベクトルを図5(b)のV1で示す．本システムでは座標値(0,0)と楕円の中心を任意のベクトルV1として設定している．そして，設定したベクトルと同図の最大領域に存在する各2値画素と楕円の中心から成るベクトルV2の角度を算出する．そして，算出された角度のヒストグラムを生成する．生成されるヒストグラムを図6に示す．同図のヒストグラムの横軸は0から360の角度を示しており，縦軸は頻度を示している． 細線化された結果は針の芯線の領域に画素を多く持つ．そのため，図6のヒストグラムで最大となる角度が針の芯線の角度となる．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/51a97592-a755-3d3d-ec0d-e4f265fd0f25.png)
図6 角度によるヒストグラム

## 目盛りの始終点検出
次に目盛りの始終点間の角度，開始点と検出した針の角度を求めるため，メータの始終点を算出する．検出された楕円の中心を用いて領域を分割する．前提条件としてメータは始点が左下，終点が右下であり，ほぼ水平に設置されているとする．楕円の中心からx方向，y方向に各2分割する．すなわち4分割することになる．そして，下側の左右それぞれの領域から始終点を求める．図7に 4分割した際の，下側の左右それぞれの結果を示す．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/b86976ca-c483-23ac-62ef-0dde0c3c0690.png)
(a) 始点側の領域分割結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/dd5448e9-2eef-5bef-f1f2-dec6a2d494ab.png)
(b) 終点側の領域分割結果
図7 目盛り領域の始終点分割結果

同図の(a)，(b)それぞれの結果に対して，y座標の値が最大となる画素を探索する．そして，その画素をメータの始終点とする．これにより，始終点と楕円の中心を用いて始点と終点の角度から目盛り間隔ごとの数値を求め，始点と針の角度から最終的なアナログ値を算出する．算出されたアナログ値の結果とプロット結果を図8に示す．現在，メータの目盛り間の総数と終点の最大値は手動で設定している．なお，今回の実験では一目盛りの間隔は5.3度(始点から終点までの角度/目盛り間隔の総数)となる．

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/987845eb-f144-dddb-7e7b-62b544f249fc.png)
図8 アナログ値の検出結果(0.64)
図8の結果はアナログ値0.638782である．また，図9に示すように圧力計は圧力の関係から目盛りの始点から一目盛り目までの目盛り間隔が一定ではない．そのため，本システムではユーザが手動で設定する．本実験で用いたメータでは，0から0.02の目盛りの間隔を通常の目盛りの間隔の0.3倍と設定している．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/1b888808-693c-3940-e5af-87c50da1ea62.png)
図9 目盛り間隔の割合

## アナログ値の算出結果
　条件が同様で針が変化するフレームに対して実験を行った．実験結果を図10に示す．同図(a),(c),(e),(g)はそれぞれの入力フレームであり，同図(b),(d),(f),(h)は同図(a),(c),(e),(g)と対になるそれぞれの算出したアナログ値と楕円の中心，始終点，針の芯線のプロット結果である．

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/31af40e1-1cdc-dcb8-15c1-a42a5185cc94.png)
(a) 入力フレーム1
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/0ba470b4-6dd6-7557-aafc-70d287899524.png)
(b) アナログ値(0.79)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/e65cf371-92fe-1596-01bf-d16c13621a23.png)
(c) 入力フレーム2
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/56fbb4e1-845a-da0c-c643-d018b98a2dfd.png)
(d) アナログ値(0.74)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/c5888d0e-879d-1aba-1802-9b93c965c5a4.png)
(e) 入力フレーム3
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/06a8727d-b422-a89f-0991-2b6642da3f2a.png)
(f) アナログ値(0.66)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/c87f5f03-f812-f252-b5d7-d41fcfa15434.png)
(g) 入力フレーム4
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/1da9b452-0bfe-7915-b1c9-0ac122a129e1.png)
(h) アナログ値(0.64)
図10 アナログ値算出結果

## 処理速度について
本手法による処理速度について示す．実験環境を表1に示す．また，処理速度を表2に示す．

表1 実験環境
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/523662f1-96d9-6ff0-5745-f6ec8b52ad98.png)

表2 処理速度の測定結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135147/23ca6e42-0b0d-d191-5d19-d1921bd59f6f.png)
同表より3FPS程度で処理が可能であることが分かる．

# 実装
大雑把ではありますが興味があればご確認ください。
https://github.com/SE4MLESS/pressure_gaugesA2D

# 終わりに
ご一読ありがとうございました。時間が限られていた際の実装とアルゴリズムの検討だったので詰めが甘い部分が多かったです。
工場では意外と需要のあるアナログメーターの自動読み取りタスクのアルゴリズム実装でした。

LGTMやアルゴリズムの改善・疑問点などありましたら、是非コメントお願い致します。