CTスキャンの原理でString Artを描く

はじめに

突然ですが、皆さんは以下の画像のような作品をYouTubeなどで見かけたことはありますでしょうか？

有名な舌を出したAlbert Einsteinですね。
しかし、よく見ると円周上の点同士を結んだ線のみで描かれていることが分かります。
YouTubeで「String Art」や「Thread Art」などと調べると、実際に糸を板に打ち付けた釘へとくくりつけて、このような作品を作成している様子を見ることができます。

初めて見たときにはとても人間業とは思えず、世界にはとんだ天才も居るものだと思いましたが、その実、計算により糸を張る手順を得ていたのでした。

今回は、そんなString ArtをCTの原理を使って描いていきたいと思います。

追記

GitHubに簡単なGUI付きのプログラムを上げました。

CTの原理

CTスキャンやMRIでは体の断面画像を得ることができますが、どのような計算をしているのでしょうか？
その答えはRadon変換（ラドン変換）と呼ばれる計算です。

ラドン変換

ラドン変換とは、CTで説明すると、物質をある角度から透過光（放射線）で見たときの明るさの分布を180°分（透過のため残り180°は同じ結果となる）集める処理に相当します。

上のGIFの右側のようなある角度から見た明るさ分布を縦軸、その角度を横軸に取ったグラフをサイノグラムあるいはシノグラムと呼びます。

逆ラドン変換

シノグラムから画像を復元する操作をそのまま逆ラドン変換といいますが、どのような複雑な計算をするのでしょうか…？
気になって調べると直交座標と極座標、実空間と周波数空間を行ったり来たりするいかにも難しそうな計算式が出てきますが、実際の復元操作自体は極めて単純です。

ある角度から見た分布をその角度方向に引き伸ばした画像を用意します。

左0°、右30°

これを全ての角度について足し合わせます。

ね、簡単でしょう？

実はこれだけでは、上のようにやや不鮮明な復元像しか得られないので、実際にはここにちょっとした工夫を施します。
先ほど座標空間を行ったり来たり云々と書いていましたが、ここの変数変換の関係で周波数空間にて$|\omega|$をヤコビアンとしてかけ合わせるのが適切（だそう）です。

かなり鮮明な復元像が得られましたね。

実用上は$|\omega|$以外にも様々なフィルターを周波数空間でかけ合わせることがあるみたいです。

String Art

String Artについて日本語の情報はほとんど見当たりませんが、英語では幾らか情報が公開されており、それらのほぼ全ては貪欲法によるアプローチです。
実のところ、今回紹介する方法よりも貪欲法のほうがString Artには適しているのですが（両方試した肌感覚上）、他の人がやっていない＋CTと芸術が結びつく面白さからラドン変換を応用した方法を紹介できればと思います。
その他のアプローチとしては、線形代数的に最適化を試みている以下の動画なども面白いです（結局貪欲法に落ち着いていましたが笑）。

ラドン変換をStringArtに応用する

逆ラドン変換ではある角度方向の濃淡を、その方向に引き伸ばした画像を合成していくのでした。
この操作、0,1で2値表現したら糸の有無を表現できるのではないでしょうか？
…というのがこの手法のモチベーションです。

例えば、次のように円周上に等間隔に60個のノードを置いたとします。これらノードから垂直方向に対応するノード同士を繋ぎます。

このノード同士のつながりをエッジと呼ぶとします。
全てのエッジはこの図形を360/60度ずつ30回回転させた、いずれかの場合に含まれます。

この30個のパターンのそれぞれが、先程の画像

に対応するということです。

非等間隔なラドン変換

このアイデアを実装するに当たり、まず考えなければならないのがラドン変換です。
上の画像から分かるようにエッジの間隔は等間隔ではないので、等間隔に並ぶ画像のピクセルからなんとかして対応する位置の値を得なければいけません。
そこでnp.interpを使用して内挿することとしました。

画像の準備

img = plt.imread(r"画像のパス（正方形）")
w = 600
cripped_img = cv2.resize(img, (w, w))
mask = ((np.mgrid[-1:1:w*1j, -1:1:w*1j]**2).sum(0)**.5) > 1
cripped_img[mask] = 0
gray_img = cv2.cvtColor(cripped_img, cv2.COLOR_BGR2GRAY).astype(float)
plt.imshow(gray_img, cmap="gray")

n = 360 # 円周の分割数（偶数）
w = gray_img.shape[0]
cos = np.cos(np.r_[:np.pi:(n//2+1)*1j])[1:-1]
sino = np.zeros((n//2, n//2-1))

for i in range(n//2):
    trans = cv2.getRotationMatrix2D((w//2, w//2), i*360/n, 1.)
    s = cv2.warpAffine(gray_img, trans, (w, w))
    sino[i] = np.interp(cos, np.r_[-1:1:w*1j], s.mean(0))

黒線が画像を回転して縦方向に平均を計算した値で、赤点線がエッジの位置で内挿した値を示します。

こうして得られた非等間隔シノグラムに先述の$|\omega|$フィルターを周波数空間でかけて得られたものが以下となります。

fsino = np.fft.fft(sino, axis=1)
fsino *= 2 * abs(np.fft.fftfreq(n//2-1))
sino = np.fft.ifft(fsino, axis=1).real
plt.pcolormesh(sino.T)
plt.colorbar()

シノグラムの2値化

得られたシノグラムは実数となっています。
しかし、表現したいのはその角度に線を引くか否かの2択ですので、シノグラムを上手いこと2値化してあげる必要があります。

ここでは2種類の2値化手段を提案します。
1つ目はシンプルに閾値を設けて2値化する方法、もう一つは誤差拡散法を用いる方法です。

誤差拡散処理をするプログラム

J. F. Jarvis, C. N. Judice, and W. H. Ninke of Bell Labsらのカーネルで拡散率を半分にしています。深い理由はありません。

norm_sino = (sino - sino.min()) / (sino.max() - sino.min())
h, w = norm_sino.shape
threshold = .85
karnel = np.array([[0, 0, 0, 7, 5], [3, 5, 7, 5, 3], [1, 3, 5, 3, 1]]) / 48 / 2
bi = (norm_sino > threshold) * 1.
error = norm_sino - bi
ed_sino = np.zeros(norm_sino.shape)
for i in range(h):
    for j in range(w):
        ed_sino[i, j] = (bi[i, j]+error[i, j]) > threshold
        e = (norm_sino[i, j] - (bi[i, j]+error[i, j] > threshold)) * karnel
        for di, dj in [[0, 1], [0, 2], [1, -2], [1, -1], [1, 0], [1, 1], [1, 2], [2, -2], [2, -1], [2, 0], [2, 1], [2, 2]]:
            if 0<=i+di<h and 0<=j+dj<w:
                error[i+di, j+dj] += e[di, dj+2]

以下はそれぞれ閾値と誤差拡散法でだいたい同じ場所が残るように2値化したシノグラムです。

誤差拡散法を使用すると高密度になりすぎないように調節できます。

2値化したシノグラムからString Artを描く

シノグラムの各列はある角度に伸びるエッジを描くか否かを端から順に示しているので、それに従って線を引いていきます。

以下は実装例です。
matplotlibで数千本の線を描くと重くなりがちなのでpyqtgraphを活用しています。
pyqtgraphはGraphItemというノードとエッジから成る"グラフ"を描く機能が提供されているため、それを使用します。
また、Jupytetでも実行するたびにシャットダウンしないように既存のインスタンスを確認するようにしてあります。

import pyqtgraph as pg
from pyqtgraph.Qt import QtCore, QtGui
from PyQt5 import QtWidgets
import numpy as np

app = QtWidgets.QApplication.instance() or QtWidgets.QApplication([])

c = 0 # 0: 黒地に白い糸, 1: 白地に黒い糸
win = pg.GraphicsLayoutWidget(show=True, title="String Art")
win.resize(800, 800)
win.setBackground(['k', 'w'][c])
plot = win.addPlot()

item = pg.GraphItem()
plot.addItem(item)
pos = np.c_[np.cos(np.r_[:2*np.pi:2*np.pi/n]), np.sin(np.r_[:2*np.pi:2*np.pi/n])]
e = ed_sino < .5 if c else ed_sino > .5 # ブーリアン型に変更

adj = []
for i in range(n//2):
    idx = np.roll(np.r_[:n], i)
    adj += np.c_[idx[1:n//2], idx[n//2+1:][::-1]][e[i]].tolist()
adj = np.array(adj)
item.setData(pos=pos, adj=adj, pen=pg.mkPen(['w', 'k'][c], width=.2), symbolBrush=['w', 'k'][c])

plot.setAspectLocked(True)

app.exec()

閾値による2値化

誤差拡散法による2値化

無事、CTの原理を応用してString Artを描くことができました。

おまけ -貪欲法によるString Art-

貪欲法のアルゴリズムはとても単純で、

• 初期位置のノードをランダムに選択する
• 現在地から全てのノードに向かう線を考えて、線が通過するピクセルの値の平均値が最も小さい（黒線の場合）ノードに移動する
• 通過したビクセルを少し明るくする（値を足す）
• 同様に次の移動先を選定する
• 以上の繰り返し

これだけです。

貪欲法の利点として、

• 1本ずつ足していくので、任意の本数、任意のタイミングで更新を打ち切ることが可能
• 順にノードを移動していくので、実際に糸を使用して作品として制作することを考えた場合に、糸が連続した1連の手順が得られる
• 移動先の選定処理を工夫することで、「絶対に真白いピクセルは通らない」などの挙動のアレンジが可能
• 処理の過程が可視化できるので見ていて楽しい

などがあります。

実装上は、「通過するピクセル」のインデックスをあらかじめ計算して保持しておくと計算を効率化・高速化できます。

「通過するピクセル」インデックス計算

n = 360
w = gray_img.shape[0]
cos, sin = np.cos(np.r_[:2*np.pi:2*np.pi/n]), np.sin(np.r_[:2*np.pi:2*np.pi/n])
lin = np.c_[:1:2j*w]
con = []
for i in tqdm(range(n)):
    a = (((cos - cos[i])*lin + cos[i] + 1) / 2 * w).astype(int).clip(0, w-1).tolist()
    b = (((sin - sin[i])*lin + sin[i] + 1) / 2 * w).astype(int).clip(0, w-1).tolist()
    c = np.dstack([a, b])
    con.append({j: np.unique(k, axis=1) for j, k in enumerate(c.transpose(1, 2, 0)) if i!=j})

貪欲法版String Art

import numpy as np
import pyqtgraph as pg
from pyqtgraph.Qt import QtCore, QtGui
from PyQt5 import QtWidgets

c = 0
temperature = 30 # 増やすと全体に満遍なく線を張ってくれるが、増やしすぎるとコントラストが下がる

nail = np.random.randint(0, n)
im_ = gray_img.astype(float)

app = QtWidgets.QApplication.instance()
if app is None:
    app = QtWidgets.QApplication([])
win = pg.GraphicsLayoutWidget(show=True)
win.setBackground(['k', 'w'][c])
view = win.addViewBox()
view.setAspectLocked(True)
pos = np.c_[np.sin(np.r_[:2*np.pi:2*np.pi/n]), -np.cos(np.r_[:2*np.pi:2*np.pi/n])]
adj = np.array([], int).reshape(0, 2)
item = pg.GraphItem()
view.addItem(item)
item.setData(pos=pos, adj=adj, pen=pg.mkPen(['w', 'k'][c], width=.2), symbolBrush=['w', 'k'][c])

def update():
    global nail, adj
    _nail = (max, min)[c](con[nail], key=lambda x: im_[con[nail][x][0], con[nail][x][1]].mean())
    a, b = con[nail][_nail]
    im_[a, b] -= temperature * (-1)**c
    adj = np.vstack([adj, [nail, _nail]])
    item.setData(pos=pos, adj=adj, pen=pg.mkPen(['w', 'k'][c], width=.2), symbolBrush=['w', 'k'][c])
    nail = _nail

timer = QtCore.QTimer()
timer.timeout.connect(update)
timer.start(1)
is_paused = False

def toggle_pause():
    global is_paused
    if is_paused:
        timer.start(1)
    else:
        timer.stop()
    is_paused = not is_paused

def keyPressEvent(event):
    if event.key() == QtCore.Qt.Key_Space:
        toggle_pause()

win.keyPressEvent = keyPressEvent

app.exec_()

…ｯ楽しい...！

おわりに

String Artの計算について解説してみました。
そろそろ同級生の第一次結婚ブームが来る時期との噂ですので、相手次第では作って贈るかもしれません笑

ご興味がありましたら、是非皆さんも作ってみてください。

800本台という少ない本数で綺麗に出来たお気に入りです。