円の中心座標を極座標変換を用いて求めたい

Question

円の中心座標を極座標変換を用いたちょっとめんどくさい方法で求めたいと思っています。
（わざわざ極座標変換なんかしなくても求まると思いますが、今回は極座標変換を使って求めたいです。）

まず上のような円が描かれた画像に対して、極座標変換を行います。極座標変換の基準となる座標を(X,Y)とすると、(X,Y)が円の中心からずれていれば変換後の画像は下のようなサインカーブっぽい線になると思います。

もし(X,Y)が円の中心と合っていれば、この下の画像のような縦に真っ直ぐな線が描かれると思います。

つまり変換後の線が真っ直ぐになるとき、(X,Y)は円の中心に近い値であると言えると思います。

【求めていること】
とある(X,Y)を基準に極座標変換したとき、変換後の線が最も真っ直ぐになるように(X,Y)を調整し、変換後の線が最も真っ直ぐになるような座標(Xc,Yc)を求めたい。

最終的に求めている出力は2つで、
・(Xc,Yc)の値
・(Xc,Yc)を基準に極座標変換した画像
です。

【皆様にお聞きしたいこと】
方法として(X,Y)の値を少しずつ変えていちいち代入して、最も最適な(X,Y)を探すのが良いのかなと思ったのですが、もし変換後の線の曲がり具合から直接(X,Y)を補正できるような方法をご存知の方がいればお聞きしたいです。
なければX,Yを代入する方法でやろうと思います。

Comments

[meg_]

質問者さんの言う「極座標変換」とはどのような計算のことでしょうか？

[mappys]

meg_様
ネットで調べて以下のコードで変換しています

img = cv2.imread('画像.png')

# キュービック補間 + 外れ値塗りつぶし + 極座標へリニアマッピング
flags = cv2.INTER_CUBIC + cv2.WARP_FILL_OUTLIERS + cv2.WARP_POLAR_LINEAR

# 引き数：画像, 変換後サイズ(幅W、高さH)、中心座標(X座標、Y座標)、半径R、変換フラグ
linear_polar_image = cv2.warpPolar(img, (W, H), (X,Y), R, flags)

[mappys]

とりあえずこれで極座標変換できてそうなのですが、正直変換フラグの中身の意味についてはイマイチ理解できてないです。

Answer 1

コーディングの話ではなく、極座標による円の中心の検出アルゴリズムについて回答します。

ある基準点を中心とする平面極座標系に円が描かれている場合、その基準点と円の位置関係は3通りとなります。
場合分けして考えてみましょう。

• 基準点が円の内部：

　基準点と円周上の点の距離の最大値MAXを与える偏角と最小値MINを与える偏角は180°ずれている
半径は(MAX+MIN)/2 であることと、ほぼ正しい中心がMAXを与える偏角方向にあることを使えば、ほぼ正しい中心の位置を推定することはたやすい。

• 基準点が円周上：

　基準点と円周上の点の距離の最小値は0である。
半径はMAX/2 であることと、ほぼ正しい中心がMAXを与える偏角方向にあることを使えば、ほぼ正しい中心の位置を推定することはたやすい。

• 基準点が円の外部：

　基準点と円周上の点の距離の最大値を与える偏角と最小値を与える偏角は等しい
半径は(MAX-MIN)/2 であることと、ほぼ正しい中心がMAXを与える偏角方向にあることを使えば、ほぼ正しい中心の位置を推定することはたやすい。

いずれにしてもほぼ正しい中心の位置は推定できるので、その周りを合わせた9点でやってみれば良いでしょう。
円周上にある場合は、内部あるいは外部とまとめて二通りとして処理しても問題ありません。

Comments

[mappys]

なるほど、すごくよくわかりました！ありがとうございます！

Answer 2

方法として(X,Y)の値を少しずつ変えていちいち代入して、最も最適な(X,Y)を探すのが良いのかなと思ったのですが

cv2.linearPolar() を使ってみましたが、円周の境界線が十分に太い事が前提条件になります(細いと変換後の線分に「隙間」が生じてしまうので適切な判定ができなくなるためです)。また、全面走査をしているので処理に時間がかかります(画像サイズ 478 x 464, Intel Core i5-8500T CPU@2.1GHz で 200 秒程度)。
なので、以下は参考程度に見ておいて下さい。

import cv2
import numpy as np

# read as grayscale
img = cv2.imread('circle.png', 0)
h, w = img.shape
r = h if w > h else w
center = []
# brute force ...
for x in range(w):
    for y in range(h):
        polar_image = cv2.linearPolar(img, (x, y), r, cv2.WARP_FILL_OUTLIERS)
        num = np.all(polar_image > 127, axis=0).sum()
        if num > 0:
            center.append((num, x, y))

# select line with maximum width
center = sorted(center, key=lambda x: x[0], reverse=True)
max_num = center[0][0]
center = [i for i in center if i[0] == max_num]

print(center)
polar_image = cv2.linearPolar(img, center[0][1:], r, cv2.WARP_FILL_OUTLIERS)
cv2.imshow("polor conversion", polar_image)
cv2.waitKey(0)
cv2.destroyAllWindows()

#
[(6, 241, 233)]

※ 左が元画像、右が極座標変換後。

Comments

[mappys]

melian様
回答ありがとうございます！すごく助かりました！
加えて質問してもよろしいでしょうか。書いていただいたコードで何をやっているか理解できていないところがございます。
#brute force内の

num = np.all(polar_image > 127, axis=0).sum()
if num > 0:

ここはどういった処理をしているのでしょうか。

[melian]

画像内に垂直線が何本あるのかを数えているのですが、127 というのは結構いい加減です(グレイスケール 0-255 の中央値)。
np.all(polar_image > 127, axis=0) では pixel を垂直方向で調べて、全ての pixel が127 以上の値を持っている場合に True が返ります(そうでなければ False)。sum() を実行することで True の個数(垂直線の本数)が判るので、 1 本以上あればリストに追加しています。

[mappys]

なるほど、理解しました！ありがとうございます！