関数radonとiradonの中身が知りたいです

Question

実現したいこと

「該当のソースコード」はこちらにあります。

11行目の関数radon、40行目の関数iradonのソースがどうなっているかを知ることができるかどうか、という質問です。radonやiradonは等間隔のデータに対して適用できるものです。
不等間隔のデータに対して、ラドン変換やラドン逆変換をやりたいときは自分でコードを書かなければいけないと思うのですが、その参考に等間隔のときのソースコードを読みたいと思っています。

該当のソースコード

import numpy as np
import matplotlib.pyplot as plt
from skimage.data import shepp_logan_phantom
from skimage.transform import radon, rescale

image = shepp_logan_phantom()

fig, ax2 = plt.subplots(1, 1, figsize=(4, 4))

theta = np.linspace(0., 180., max(image.shape), endpoint=False)
sinogram = radon(image, theta=theta)
dx, dy = 0.5 * 180.0 / max(image.shape), 0.5 / sinogram.shape[0]
ax2.imshow(sinogram, cmap=plt.cm.Greys_r,
           extent=(-dx, 180.0 + dx, -dy, sinogram.shape[0] + dy),
           aspect='auto')
fig.tight_layout()
plt.show()

#######################################
# Reconstruction with the Filtered Back Projection (FBP)

from skimage.transform.radon_transform import _get_fourier_filter

filters = ['ramp', 'shepp-logan', 'cosine', 'hamming', 'hann']

for ix, f in enumerate(filters):
    response = _get_fourier_filter(2000, f)
    plt.plot(response, label=f)

plt.xlim([0, 1000])
plt.xlabel('frequency')
plt.legend()
plt.show()

#####################################
# Applying the inverse radon transformation with the 'ramp' filter, we get:

from skimage.transform import iradon

reconstruction_fbp = iradon(sinogram, theta=theta, filter_name='ramp')
error = reconstruction_fbp - image

imkwargs = dict(vmin=-0.2, vmax=0.2)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5),
                               sharex=True, sharey=True)
ax1.set_title("Reconstruction\nFiltered back projection")
ax1.imshow(reconstruction_fbp, cmap=plt.cm.Greys_r)
ax2.set_title("Reconstruction error\nFiltered back projection")
ax2.imshow(reconstruction_fbp - image, cmap=plt.cm.Greys_r, **imkwargs)
plt.show()

試したこと

ネットでいろいろ検索しましたが、ソースは見つけられませんでした。

Comments

[melian]

scikit-image の API リファレンスにソースコードへのリンクがあります。

skimage.transform — skimage 0.23.2 documentation
https://scikit-image.org/docs/stable/api/skimage.transform.html

radon
https://scikit-image.org/docs/stable/api/skimage.transform.html#skimage.transform.radon
https://github.com/scikit-image/scikit-image/blob/v0.23.2/skimage/transform/radon_transform.py#L16-L114

iradon
https://scikit-image.org/docs/stable/api/skimage.transform.html#skimage.transform.iradon
https://github.com/scikit-image/scikit-image/blob/v0.23.2/skimage/transform/radon_transform.py#L187-L323

[yyicp]

>melian様
ありがとうございました。
githubの方は見つけられていなかったので助かりました。

Answer 1

不等間隔のデータに対して、ラドン変換やラドン逆変換をやりたいときは自分でコードを書かなければいけないと思う

skimage.transformのradon, iradon関数についての質問です。
この関数にはthetaという引数があります。
この引数を指定することで、角度が不等間隔のデータに対して、ラドン変換やラドン逆変換を行うことができると解釈してよいでしょうか？

Geminiの回答

skimage.transformのradon、iradon関数におけるtheta引数について

はい、その解釈は正しいです。skimage.transformモジュールのradon関数とiradon関数におけるtheta引数は、角度が不等間隔のデータに対してラドン変換とラドン逆変換を行うために使用できます。

不等間隔角度におけるラドン変換

radon関数は、画像配列とオプションのtheta配列を引数として受け取ります。theta配列は、ラドン変換を計算する角度を指定します。thetaが提供されない場合、変換は0度と180度の間で等間隔の角度で計算されます。しかし、theta配列を提供すると、変換は指定された角度で計算されます。

以下のコードは、0度、45度、90度の角度で画像のラドン変換を計算する例です。

import numpy as np
from skimage.transform import radon

img = np.zeros((10, 10))
img[2:4, 2:4] = 1
img[7:9, 2:4] = 1

theta = [0, np.pi/4, np.pi/2]
sinogram = radon(img, theta=theta)

不等間隔角度におけるラドン逆変換

iradon関数は、シノグラム配列とオプションのtheta配列を引数として受け取ります。theta配列は、シノグラムが計算された角度を指定します。thetaが提供されない場合、シノグラムは0度と180度の間で等間隔の角度で計算されたと仮定されます。しかし、theta配列を提供すると、逆変換は指定された角度を使用して実行されます。

以下のコードは、0度、45度、90度の角度で計算されたシノグラムから画像を再構築する例です。

import numpy as np
from skimage.transform import iradon

sinogram = ...  # 前の例からシノグラムを読み込む

theta = [0, np.pi/4, np.pi/2]
img_rec = iradon(sinogram, theta=theta)

不等間隔角度を使用する利点

ラドン変換とラドン逆変換に不等間隔角度を使用する利点はいくつかあります。

• 精度向上: 場合によっては、不等間隔角度を使用することで、再構築の精度を向上させることができます。これは、角度の選択を再構築する特定のオブジェクトに合わせることができるためです。
• 計算コストの削減: 特定の角度でのみオブジェクトの再構築に関心がある場合は、thetaを使用してその角度のみを指定できます。これにより、再構築の計算コストを削減できます。

まとめ

skimage.transformモジュールのradon関数とiradon関数におけるtheta引数は、角度が不等間隔のデータに対してラドン変換とラドン逆変換を実行するために使用できます。これは、再構築の精度を向上させたり、計算コストを削減したりするために役立ちます。

Comments

[yyicp]

どうもありがとうございました。

例えば、
https://github.com/scikit-image/scikit-image/blob/main/skimage/transform/radon_transform.py
にも"evenly spaced"と書いてあったのでradonは等間隔のときにしか使えないと思っていました。
いろいろ調べてみても、不等間隔での計算例を見つけられなかったので、
theta = [0, np.pi/4, np.pi/2]（これは等間隔ですが）
のように与えるのが正しいのか要確認ですね。
（正しければこの後の作業がとても楽になるので正しいことを祈っています）

[yyicp]

よく読んでいなくて、勘違いしていました。

Radon変換について、「該当のソースコード」で10行目でthetaの配列を定義して、11行目でradonの引数でtheta=thetaとしているので、これで正しそうです。（等間隔でなくても良さそうです）

Radon逆変換について、上記URLでiradonの記述を見ると

theta : array, optional
Reconstruction angles (in degrees). Default: m angles evenly spaced
between 0 and 180 (if the shape of `radon_image` is (N, M)).

再構成角度 (度単位)。デフォルト: m 個の等間隔の角度
0 から 180 の間 (「radon_image」の形状が (N, M) の場合)。

となっていたので、デフォルトが等間隔なだけで、自分で指定すれば等間隔でなくても良さそうです。

Answer 2

from skimage.transform import radon, rescale
とあるのでskimageで検索するとscikit-imageが正式名称っぽいので今度はそれで検索するとGitHubのリポジトリが出てくるのでリポジトリ内でradonを検索するとこちらが出てきました

https://github.com/scikit-image/scikit-image/blob/main/skimage/transform/radon_transform.py

Comments

[yyicp]

どうもありがとうございました。