回転行列最適角度について

Question

目的:上画像の黄色部分の歯をAffine変換による回転により下画像のように縦向きに変換。
-90度から90度まで回転させ最適な角度を出力。

自分で考えたものとしては、2パターンあり、

1. 矩形で歯を囲みその中の歯以外のピクセルを合計し最小のものを最適とする。
2. 列ごとに歯部分のピクセルを数え最多のものを最適とする。

この2つだと実際に見てみると最適でないように感じました。
縦向きになったとどのように人が認知しているのか知りたいです。最適な角度を算出するアイデアをご教授願います。

添付フォルダ(tooth_array.zip内に3つnpyファイル、ipynbファイル)

#coding:utf-8
import numpy as np
import matplotlib.pyplot as plt

# アフィン変換で画像配列の回転
def rotate_affine(src, theta):
    # 元画像のサイズを取得
    h, w = src.shape[0], src.shape[1]

    # 出力画像用の配列生成（要素は全て0）
    dst = np.zeros((h,w))

    # degreeかｒradianに変換
    rd = np.radians(theta)

    # アフィン変換
    for y in range(0, h):
        for x in range(0, w):
            xi = int((x - int(w/2))*np.cos(rd) - (y - int(h/2))*np.sin(rd) + int(w/2))
            yi = int((x - int(w/2))*np.sin(rd) + (y - int(h/2))*np.cos(rd) + int(h/2))
            
            # 変換後の座標が範囲外でなければ出力画像配列に画素値を代入
            if yi < h - 1 and xi < w - 1 and yi > 0 and xi > 0:
                dst[y][x] = src[yi][xi]

    return dst

def pattern1(src):
    sum_cnt = 1000000
    best_theta = 0
    for theta in range(-90,90+1):
        dst = rotate_affine(src,theta)
        vertical_info = np.count_nonzero(dst ==255, axis=0)
        side_info = np.count_nonzero(dst ==255, axis=1)
        min_flag = False
        for i,num in enumerate(vertical_info):
            if num != 0 and min_flag == False:
                min_flag = True
                min_x = i
            if num == 0 and min_flag == True:
                max_x = i-1
                min_flag = False

        for i,num in enumerate(side_info):
            if num != 0 and min_flag == False:
                min_flag = True
                min_y = i
            if num == 0 and min_flag == True:
                max_y = i-1
                min_flag = False
        tooth_co_info = [min_y,min_x,max_y,max_x]
        ex_dst = dst[tooth_co_info[0]:tooth_co_info[2]+1,tooth_co_info[1]:tooth_co_info[3]+1]
        print("theta,sum:{},{}".format(theta,np.count_nonzero(ex_dst ==0, axis=0).sum()))
        if np.count_nonzero(ex_dst ==0, axis=0).sum() < sum_cnt:
            sum_cnt = np.count_nonzero(ex_dst ==0, axis=0).sum()
            best_theta = theta
    print("best")
    print(best_theta,sum_cnt)


def pattern2(src):
    max_cnt = 0
    best_theta = 0
    for theta in range(-90,90+1):
        dst = rotate_affine(src,theta)
        print("theta,count:{},{}".format(theta,np.count_nonzero(dst ==255, axis=0).max()))
        if np.count_nonzero(dst ==255, axis=0).max() > max_cnt:
            max_cnt = np.count_nonzero(dst ==255, axis=0).max()
            best_theta = theta
    print("best")
    print("theta,count:{},{}".format(theta,max_cnt))

if __name__ == '__main__':
    mask = np.load('mask.npy')
    bi_mask = np.where(mask==True,255,0)
    h, w = bi_mask.shape[0], bi_mask.shape[1]
    src = np.zeros((int((w**2 + h**2)**(1/2)), int((w**2 + h**2)**(1/2))))
    dx,dy = int((w**2 + h**2)**(1/2)/2) - int(w/2), int((w**2 + h**2)**(1/2)/2) - int(h/2)

    for y in range(0, h):#平行移動
        for x in range(0, w):
            src[y+dy][x+dx] = bi_mask[y][x]

    theta = 0
    dst = rotate_affine(src,theta)
    plt.imshow(dst)
    plt.show()
    
    pattern1(src) #切り替え必要
#     pattern2(src)

Comments

[toast-uz]

興味深いけど、1と2をPythonコードで出して結果を表示しないと、プログラミングの話題にならないですね。PythonとかNumPyとかタグをつけるべきでもない。

[toast-uz]

コードの提示ありがとうございます。画像の中心でまわすことには、こだわらなくて、よいでしょうか？
（まわすだけでなく平行移動もして、より自然な見た目を追求してよいでしょうか？）
ちなみに、mainの方で平行移動されているようですが、これは何を基準に平行移動されていますか？

[amayan]

すべての歯に対して最適な角度を算出することができるなら大丈夫です。
元々歯の周りを囲む矩形のみしかなかったのですが、その場で回転させると歯が途切れて見えてしまいます。回転させるためndarrayのサイズを回転させる分大きくし、元々のndarrayの中心と拡張したndarrayの中心を一致させています。

Answer 1

scipyをフル活用することで全体的にすっきりと作ってみました。

• アイデアとしては、まずは歯の形を平面に置いた重心に中心を揃え、その後、歯の形を立ててバランスが取れる回転位置を求める、というものです。
• バランスが取れるかどうかは、歯全体を重さ密度一定の剛体としてみた場合に、垂直に立てた状態で、中心軸に対する慣性モーメントが最小になる、という定義をしました。
• 最初の重心、平行移動、回転は、全てscipy.ndimageを使い、モーメント演算だけは見当たらないので自作、そしてモーメント最小を求めるのにscipy.optimizeを使いました。
• なお、歯の画像は質問者様の最初の画像をキャプチャして作成しました（ですのでピクセルサイズは元とは異なります）。回転によって歯の画像が枠外にはみ出ない前提でいます。

以下がコードと実行結果の図です。-23.36度ほどの回転が最適である、と出ました。質問者様の結果と比較して、ほんの少しだけ右回転しており、とてもバランスが良いように見えます。

※ 当初は、左右で逆向きの回転モーメントを使ってましたが、慣性モーメントを使った方が安定した形になるとわかりました。回転モートメントだと横に伸ばした形でも左右対称であれば安定するのに対して、慣性モーメントだと横をなるべく小さくするような形が有利になります。フィギュアスケートの選手が腕をたたんで高速回転することと、原理は同じです。

import numpy as np
from scipy import ndimage, optimize
import matplotlib.pyplot as plt
import cv2

# 画像読み込みと2値化
img = cv2.imread('tooth.png', 0)
_, img = cv2.threshold(img, 100, 1, cv2.THRESH_BINARY)

# 平行移動して重心を中心にする
center = ndimage.center_of_mass(img)
img = ndimage.shift(img, np.array(img.shape)/2-np.array(center))

# 図形を立てた時の、中心軸に対する慣性モーメント
# mass:縦軸ピクセル数、radius:中心軸からの左右のピクセル距離
# 慣性モーメント = Σ mass * radius^2
def inertia(img):
    mass = img.sum(axis=0)
    radius = np.abs(np.arange(-len(mass)//2, len(mass)//2+1))
    if len(mass) != len(radius):
        radius = radius[radius != 0] - 0.5  # 中心位置の補正
    return (mass * radius * radius).sum()

# 回転させて慣性モーメントを測定する=最小化関数
def rotated_inertia(degree, img):
    return inertia(ndimage.rotate(img, degree, reshape=False))

# 最小化関数の最小値を求めることで、-90〜90度の回転時の最小慣性モーメントを得る
res = optimize.minimize_scalar(
    rotated_inertia, bounds=[-90,90], args=(img), method='Bounded')

# 最小値の時の回転角度を求める
print(res.x)

# 最小値の時の回転を実際に行う
img = ndimage.rotate(img, res.x, reshape=False)

plt.imshow(img)
plt.show()

結果

Comments

[amayan]

慣性モーメントかなりいいですね。
ありがとうございました。