Поворот одной линии в двухмерном дискретном изображении действительно сложен. Вы всегда получаете грубое приближение, интерполяция не очень помогает.
Процесс, который вы намереваетесь следовать (я добавил фильтрацию):
- Для каждой проекции в преобразовании Радона:
- примените БПФ
- примените фильтр клина
- запишите это как строку, хотя источник двумерного комплексного изображения
- повернуть это изображение для соответствия ориентации проекции
- накапливает результат в выходном изображении в частотной области путем суммирования
- . Примените 2D IFFT к изображению в частотной области, чтобы получить восстановленное изображение
Поскольку мы знаем, что операция IFFT коммутирует с суммированием, можно переместить операцию IFFT в l oop:
- Для каждая проекция в преобразовании Радона:
- применить БПФ
- применить фильтр клина
- записать это в виде строки, хотя источник двумерного комплексного изображения
- поверните это изображение, чтобы соответствовать t Ориентация проекции
- Применить 2D IFFT
- , чтобы накапливать результат в выходном изображении в пространственной области путем суммирования
Также вращения и операции IFFT коммутируют, поэтому вышеприведенное идентично:
- Для каждого проецирования в преобразовании Радона:
- применить FFT
- применить фильтр клина
- запишите это как строку, хотя начало 2D сложного изображения
- Примените 2D IFFT
- поверните это изображение, чтобы соответствовать ориентации проекции
- накапливать результат в выходном изображении в пространственной области суммированием
В этом последнем случае мы вращаем изображение в пространственной области, которое является плавным; это не одиночная линия, нарисованная в противном случае пустым изображением, это полностью ограниченная полосой функция, которая может быть правильно интерполирована. В этом случае результат поворота намного лучше.
Этот последний процесс почти такой же, как и алгоритм обратной проекции. Кроме того, мы можем понять, что двумерное IFFT изображения с одной строкой данных через начало координат (все остальное изображение равно нулю) - это то же самое, что взять 1-мерное IFFT и воспроизвести его во всех строках изображения. Это экономит немного вычислений.
Вот код. Первый метод будет (несколько исправлений из кода OP, но выходные данные по-прежнему не распознаются!):
import numpy as np
import matplotlib.pyplot as plt
from skimage.io import imread
from skimage.data import shepp_logan_phantom
from skimage.transform import radon, rescale
from skimage.transform import iradon
from skimage.transform import rotate
import cv2
x = shepp_logan_phantom()
x = cv2.resize(x, (128,128), interpolation = cv2.INTER_AREA)
theta = np.linspace(0, 180, len(x), endpoint=False)
R = radon(x, theta)
filter = np.abs(np.fft.fftfreq(128))
fullFft2D = np.zeros((128,128)).astype(np.complex128)
for i in range(len(theta)):
temp_ = np.zeros((128,128)).astype(np.complex128)
temp_[64,:] = np.fft.fftshift(filter * np.fft.fft(R[:,i]))
fft_real = rotate(np.real(temp_), theta[i], order=0, center=(64,64))
fft_imag = rotate(np.imag(temp_), theta[i], order=0, center=(64,64))
fullFft2D += fft_real + 1j*fft_imag
y = np.fft.ifft2(np.fft.ifftshift(fullFft2D))
plt.imshow(np.fft.fftshift(y.real)); plt.show()
Исправления включают в себя: (1) Источник в частотной области fftshift ed размера 128 в 64, а не 63. (2) Вращение явно выполняется вокруг начала координат. (3) У ОП была опечатка: fft_real + 1j* fft_real. (4) Добавлена фильтрация клина. (5) Не включая 180 градусов в преобразовании Радона (так как он идентичен 0 градусам). (6) Использование действительной части IFFT, а не абсолютного значения.
При вычислении через частотную область, если вы ожидаете реальный результат, но получаете нетривиальный (тривиальный == почти нулевой) мнимый компонент, что-то не так. В приведенном выше коде мнимая составляющая нетривиальна. Это результат ротации данных, которые не могут быть правильно интерполированы. Вращение просто разрушает изменения успеха.
Последний метод будет:
y = np.zeros((128,128))
for i in range(len(theta)):
tmp_ = np.zeros((128,128)).astype(np.complex128)
tmp_[0,:] = filter * np.fft.fft(R[:,i])
y += rotate(np.fft.ifft2(tmp_).real, -theta[i], center=(64,64))
plt.imshow(y); plt.show()
Этот код несколько упрощен, потому что нам не нужно использовать fftshift, мы можем написать прямая в начале координат, как ожидается от БПФ (строка 0). Полученный результат правильно воспроизводит фантом.