python: ラドン変換＋逆変換（ＣＴスキャン）

はじめに

今回はラドン変換とその逆変換を実装していく。ただし、数式や概念は理解したが少々説明しようと複雑になる。そこで、できるだけシンプルなpythonコードで実装することで、体感的に何をしているか理解できるような資料を目指す。

実装コード

今回用いるのは以前作成したアフィン変換と、今回新たに作成したラドン変換＋逆変換の処理である。

アフィン変換

まずアフィン変換は以下のように実装されている。ラドン変換のみ理解したい場合は飛ばしてもらってよい。アフィン変換についての詳細については以前投稿したPython: 自作アフィン変換（改）を見てください（宣伝）。

import cv2
import numpy as np
import matplotlib.pyplot as plt


def translation_matrix(tx, ty):
    return np.array([[1, 0, -tx],
                     [0, 1, -ty],
                     [0, 0, 1]])


def rotation_matrix(a):
    return np.array([[np.cos(a), -np.sin(a), 0],
                     [np.sin(a),  np.cos(a), 0],
                     [        0,          0, 1]])


def shear_matrix(mx, my):
    return np.array([[1, -mx, 0],
                     [-my, 1, 0],
                     [0,  0, 1]])    


def scaling_matrix(sx, sy):
    return np.array([[1/sx, 0, 0],
                     [0, 1/sy, 0],
                     [0,  0, 1]])


def affin(img, m):
    WID = np.max(img.shape)
    x = np.tile(np.linspace(-1, 1, WID).reshape(1, -1), (WID, 1))
    y = np.tile(np.linspace(-1, 1, WID).reshape(-1, 1), (1, WID))
    p = np.array([[x, y, np.ones(x.shape)]])
    dx, dy, _ = np.sum(p * m.reshape(*m.shape, 1, 1), axis=1)
    u = np.clip((dx + 1) * WID / 2, 0, WID-1).astype('i')
    v = np.clip((dy + 1) * WID / 2, 0, WID-1).astype('i')
    return img[v, u]

ラドン変換＋逆変換

次に以下が今回実装したラドン変換と逆変換の処理である。ごちゃごちゃと説明すると長くなるので、原理の説明についてはwikiに丸投げします。（汗）難しい数式が沢山出てきますが、アルゴリズムとしてシンプルに整理すると以下で実現できました。（自作なので本アルゴリズムについて参考になるものはないです）

def radon(img, n):
    M = [rotation_matrix(i * np.pi / (n - 1)) for i in range(n)]
    return np.array([np.sum(affin(img, m), axis=0) for m in M])

def dradon(radon_img):
    M = [rotation_matrix(i * np.pi / (len(radon_img) - 1)) for i in range(len(radon_img))]
    A = np.array([affin(np.tile(w, (len(w), 1)), m) for w, m in zip(radon_img, M)])
    return np.sum(A, axis=0)[::-1]

説明

今回はよくあるレナさんの画像にラドン変換をかけていく。ということでまずレナさんの画像を読み込みます。そして、一応確認しておく。実行するとこんな感じ。

img = cv2.resize(cv2.imread('lena.jpg', 0), (256, 256))
plt.imshow(img)
plt.show()

また、今回はアフィン変換による回転を用いるので、実際に回転させる処理をここで確認しておく。アフィン変換による回転は画像imgと回転行列mを渡すと実現できる。0~πまでをn分割してnパターンの回転行列Mを作る。あとは順番に回転行列を渡して変換して可視化していくと以下のようになる

n = 6
M = [rotation_matrix(i * np.pi / (n - 1)) for i in range(n)]
for m in M:
    plt.imshow(affin(img, m))
    plt.show()

ここからやっとラドン変換とかけていきます。ラドン変換では上のアフィン変換と同じように回転行列を作っていき、レナさんを回転させながら縦軸（axis=0）方向で総和を求めていきます。ラドン変換を下に表示しておきます（縦軸が角度です）。この処理はある方向からビームを照射しそれを観測することに対応しています。つまり、１８０度方向からぐるっとレナさんにビームを照射すると、各方向からの観測結果が以下の図のようになるということです。これのよいところは、中がどうなっているかわからなくても、ビームを外からぐるっと照射して観測すれば得られることです。

def radon(img, n):
    M = [rotation_matrix(i * np.pi / (n - 1)) for i in range(n)]
    return np.array([np.sum(affin(img, m), axis=0) for m in M])

radon_img = radon(img, 180)
plt.imshow(radon_img)

あとはこの結果から元のレナさんを復元できれば、１８０度方向からビームでぐるっと照射して得られた結果から、内部がどうなっているかわかるわけです。ということで、実際に復元してみましょう。やってることはシンプルですが、説明がややこしいので、アルゴリズムから理解してもらえるといいかと思います・・・(´・ω・｀)。（ざっというと、復元先の各ピクセルにはラドン変換の各角度に対して対応する値が必ず一つあり、それらを全て足した値がそのピクセルの値になるアルゴリズムになっています。）まぁ３行しかないので頑張って理解していただけると助かる・・・。

def dradon(radon_img):
    M = [rotation_matrix(i * np.pi / (len(radon_img) - 1)) for i in range(len(radon_img))]
    A = np.array([affin(np.tile(w, (len(w), 1)), m) for w, m in zip(radon_img, M)])
    return np.sum(A, axis=0)[::-1]

plt.imshow(dradon(radon_img))
plt.show()
plt.imshow(img)
plt.show()

復元画像は↓

元画像は↓

荒くラドン変換しているのでボヤっとしていますが復元できていることが確認できます。

まとめ

CTスキャンとかでよく使われるらしいのですが、ラドン変換はあまり実装出来ている例などが無く、またあっても複雑に実装しすぎていてわかりにくいのが多かった（特に逆変換）。私のは数行だし、理解するためには役立つのかなぁ～？(；・∀・)。まぁもしこれを見て分かりやすかったならgoodでも押しておいてください(´・ω・｀)