Что не так с этим кодом для томографической реконструкции методом Фурье?

19

Недавно я играл с алгоритмами томографической реконструкции. У меня уже есть хорошие рабочие реализации FBP, ART, SIRT / SART-подобная итерационная схема и даже использование прямой линейной алгебры (медленно!). Этот вопрос не о какой-либо из этих техник ; ответы на вопрос «почему кто-то так поступил, вместо этого вот код FBP» - это не то, что я ищу.

Следующее, что я хотел сделать с этой программой, это « завершить набор » и реализовать так называемый « метод восстановления Фурье ». Мое понимание этого в основном заключается в том, что вы применяете 1D БПФ к синограммным «экспозициям», размещаете их как радиальные «спицы колеса» в 2D-пространстве Фурье (что это полезно сделать, следует непосредственно из теоремы о центральном срезе) интерполировать из этих точек в регулярную сетку в этом двумерном пространстве, и тогда должна быть возможность обратного преобразования Фурье для восстановления исходной цели сканирования.

Звучит просто, но мне не повезло, что я получил какие-либо реконструкции, которые выглядят как оригинальная цель.

Приведенный ниже код на языке Python (numpy / SciPy / Matplotlib) является наиболее кратким выражением из того, что я пытаюсь сделать. При запуске отображается следующее:

Рисунок 1: цель рисунок 1

Рисунок 2: синограмма цели fig2

Рисунок 3: FFT-ed строки синограммы fig3

Фиг.4: верхний ряд - это двумерное пространство БПФ, интерполированное из строк синограммы области Фурье; нижний ряд (для сравнения) представляет собой прямое 2D БПФ цели. В этот момент я начинаю подозревать; графики, интерполированные из БПФ с синограммой, выглядят аналогично графикам, полученным путем непосредственного 2D-БПФ цели ... и все же отличаются. fig4

Рисунок 5: обратное преобразование Фурье, показанное на рисунке 4. Я надеялся, что это будет более узнаваемой целью, чем на самом деле. fig5

Есть идеи, что я делаю не так? Не уверен, что мое понимание реконструкции метода Фурье в корне неверно, или в моем коде есть какая-то ошибка.

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

import scipy.interpolate
import scipy.fftpack
import scipy.ndimage.interpolation

S=256  # Size of target, and resolution of Fourier space
A=359  # Number of sinogram exposures

# Construct a simple test target
target=np.zeros((S,S))
target[S/3:2*S/3,S/3:2*S/3]=0.5
target[120:136,100:116]=1.0

plt.figure()
plt.title("Target")
plt.imshow(target)

# Project the sinogram
sinogram=np.array([
        np.sum(
            scipy.ndimage.interpolation.rotate(
                target,a,order=1,reshape=False,mode='constant',cval=0.0
                )
            ,axis=1
            ) for a in xrange(A)
        ])

plt.figure()
plt.title("Sinogram")
plt.imshow(sinogram)

# Fourier transform the rows of the sinogram
sinogram_fft_rows=scipy.fftpack.fftshift(
    scipy.fftpack.fft(sinogram),
    axes=1
    )

plt.figure()
plt.subplot(121)
plt.title("Sinogram rows FFT (real)")
plt.imshow(np.real(np.real(sinogram_fft_rows)),vmin=-50,vmax=50)
plt.subplot(122)
plt.title("Sinogram rows FFT (imag)")
plt.imshow(np.real(np.imag(sinogram_fft_rows)),vmin=-50,vmax=50)

# Coordinates of sinogram FFT-ed rows' samples in 2D FFT space
a=(2.0*math.pi/A)*np.arange(A)
r=np.arange(S)-S/2
r,a=np.meshgrid(r,a)
r=r.flatten()
a=a.flatten()
srcx=(S/2)+r*np.cos(a)
srcy=(S/2)+r*np.sin(a)

# Coordinates of regular grid in 2D FFT space
dstx,dsty=np.meshgrid(np.arange(S),np.arange(S))
dstx=dstx.flatten()
dsty=dsty.flatten()

# Let the central slice theorem work its magic!
# Interpolate the 2D Fourier space grid from the transformed sinogram rows
fft2_real=scipy.interpolate.griddata(
    (srcy,srcx),
    np.real(sinogram_fft_rows).flatten(),
    (dsty,dstx),
    method='cubic',
    fill_value=0.0
    ).reshape((S,S))
fft2_imag=scipy.interpolate.griddata(
    (srcy,srcx),
    np.imag(sinogram_fft_rows).flatten(),
    (dsty,dstx),
    method='cubic',
    fill_value=0.0
    ).reshape((S,S))

plt.figure()
plt.suptitle("FFT2 space")
plt.subplot(221)
plt.title("Recon (real)")
plt.imshow(fft2_real,vmin=-10,vmax=10)
plt.subplot(222)
plt.title("Recon (imag)")
plt.imshow(fft2_imag,vmin=-10,vmax=10)

# Show 2D FFT of target, just for comparison
expected_fft2=scipy.fftpack.fftshift(scipy.fftpack.fft2(target))

plt.subplot(223)
plt.title("Expected (real)")
plt.imshow(np.real(expected_fft2),vmin=-10,vmax=10)
plt.subplot(224)
plt.title("Expected (imag)")
plt.imshow(np.imag(expected_fft2),vmin=-10,vmax=10)

# Transform from 2D Fourier space back to a reconstruction of the target
fft2=scipy.fftpack.ifftshift(fft2_real+1.0j*fft2_imag)
recon=np.real(scipy.fftpack.ifft2(fft2))

plt.figure()
plt.title("Reconstruction")
plt.imshow(recon,vmin=0.0,vmax=1.0)

plt.show()
timday
источник
... потому что здесь есть код для этого Вещи, которые должны быть в центре, находятся по краям, а вещи, которые должны быть по краям, находятся в центре, как будто сдвиг фазы на 90 градусов где-то не должен быть?
эндолит
1
Код, который вы связали, предназначен для метода отфильтрованной обратной проекции (FBP). Который основан на той же математике центрального среза, но никогда явно не пытается построить двухмерное изображение области Фурье. Вы можете посмотреть на подавление низких частот фильтром FBP как на компенсацию более высокой плотности "спиц" центрального среза в середине. В методе реконструкции Фурье, который я пытаюсь реализовать, это просто проявляется как более высокая плотность точек для интерполяции. Я свободно признаю, что пытаюсь внедрить немного использованную технику, и в литературе ее освещение ограничено,
Timday
Ой, да, ты прав. Вот версия в C . Я немного просмотрел и опубликовал некоторые вещи. Я посмотрю больше позже.
Эндолит

Ответы:

15

ОК, я наконец-то взломал это.

Трюк в основном сводился к тому, чтобы поместить some fftshift/ ifftshifts в правильное место, чтобы представление 2D в пространстве Фурье не было дико колебательным и было обречено на невозможность точной интерполяции. По крайней мере, это то, что я думаю исправил. Большая часть моего ограниченного понимания теории Фурье основана на непрерывной интегральной формулировке, и я всегда нахожу дискретную область и БПФ немного ... странными.

Хотя я нахожу код Matlab довольно загадочным, я должен отдать должное этой реализации, по крайней мере, чтобы дать мне уверенность, что этот алгоритм реконструкции может быть достаточно компактно выражен в такой среде.

Сначала я покажу результаты, затем код:

Рисунок 1: новая, более сложная цель. Рисунок 1

Рисунок 2: синограмма (ОК, ОК, преобразование Радона) цели. fig2

Рисунок 3: FFT-ряды синограммы (нанесены с DC в центре). Рис3

Рисунок 4: синограмма с БПФ, преобразованная в пространство 2D БПФ (постоянный ток в центре). Цвет является функцией абсолютного значения. fig4

Рисунок 4a: Увеличьте центр 2D-пространства БПФ, чтобы лучше показать радиальную природу данных синограммы. Fig4a

Рисунок 5: Верхний ряд: пространство 2D БПФ, интерполированное из радиально расположенных строк синограммы с БПФ. Нижний ряд: ожидаемое появление от просто 2D БПФ-цели.
Fig5

Рис. 5а. Увеличьте центральную область вспомогательных участков на рис. 5, чтобы показать, что они выглядят довольно качественно в хорошем соответствии. Fig5a

Рисунок 6: Кислотный тест: обратное 2D БПФ интерполированного пространства БПФ восстанавливает цель. Лена все еще выглядит неплохо, несмотря на все, что мы только что с ней справились (возможно, потому, что «спиц» синограмм достаточно, чтобы достаточно плотно покрыть плоскость 2D БПФ; все становится интересным, если вы уменьшите количество углов экспонирования, так что это уже не так). ). введите описание изображения здесь

Вот код; выводит графики менее чем за 15 секунд на 64-битной SciPy Debian / Wheezy на i7.

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

import scipy.interpolate
import scipy.fftpack
import scipy.misc
import scipy.ndimage.interpolation

S=256 # Size of target, and resolution of Fourier space
N=259 # Number of sinogram exposures (odd number avoids redundant direct opposites)

V=100 # Range on fft plots

# Convenience function
def sqr(x): return x*x

# Return the angle of the i-th (of 0-to-N-1) sinogram exposure in radians.
def angle(i): return (math.pi*i)/N

# Prepare a target image
x,y=np.meshgrid(np.arange(S)-S/2,np.arange(S)-S/2)
mask=(sqr(x)+sqr(y)<=sqr(S/2-10))
target=np.where(
    mask,
    scipy.misc.imresize(
        scipy.misc.lena(),
        (S,S),
        interp='cubic'
        ),
    np.zeros((S,S))
    )/255.0

plt.figure()
plt.title("Target")
plt.imshow(target)
plt.gray()

# Project the sinogram (ie calculate Radon transform)
sinogram=np.array([
        np.sum(
            scipy.ndimage.interpolation.rotate(
                target,
                np.rad2deg(angle(i)), # NB rotate takes degrees argument
                order=3,
                reshape=False,
                mode='constant',
                cval=0.0
                )
            ,axis=0
            ) for i in xrange(N)
        ])

plt.figure()
plt.title("Sinogram")
plt.imshow(sinogram)
plt.jet()

# Fourier transform the rows of the sinogram, move the DC component to the row's centre
sinogram_fft_rows=scipy.fftpack.fftshift(
    scipy.fftpack.fft(
        scipy.fftpack.ifftshift(
            sinogram,
            axes=1
            )
        ),
    axes=1
    )

plt.figure()
plt.subplot(121)
plt.title("Sinogram rows FFT (real)")
plt.imshow(np.real(sinogram_fft_rows),vmin=-V,vmax=V)
plt.subplot(122)
plt.title("Sinogram rows FFT (imag)")
plt.imshow(np.imag(sinogram_fft_rows),vmin=-V,vmax=V)

# Coordinates of sinogram FFT-ed rows' samples in 2D FFT space
a=np.array([angle(i) for i in xrange(N)])
r=np.arange(S)-S/2
r,a=np.meshgrid(r,a)
r=r.flatten()
a=a.flatten()
srcx=(S/2)+r*np.cos(a)
srcy=(S/2)+r*np.sin(a)

# Coordinates of regular grid in 2D FFT space
dstx,dsty=np.meshgrid(np.arange(S),np.arange(S))
dstx=dstx.flatten()
dsty=dsty.flatten()

plt.figure()
plt.title("Sinogram samples in 2D FFT (abs)")
plt.scatter(
    srcx,
    srcy,
    c=np.absolute(sinogram_fft_rows.flatten()),
    marker='.',
    edgecolor='none',
    vmin=-V,
    vmax=V
    )

# Let the central slice theorem work its magic!
# Interpolate the 2D Fourier space grid from the transformed sinogram rows
fft2=scipy.interpolate.griddata(
    (srcy,srcx),
    sinogram_fft_rows.flatten(),
    (dsty,dstx),
    method='cubic',
    fill_value=0.0
    ).reshape((S,S))

plt.figure()
plt.suptitle("FFT2 space")
plt.subplot(221)
plt.title("Recon (real)")
plt.imshow(np.real(fft2),vmin=-V,vmax=V)
plt.subplot(222)
plt.title("Recon (imag)")
plt.imshow(np.imag(fft2),vmin=-V,vmax=V)

# Show 2D FFT of target, just for comparison
expected_fft2=scipy.fftpack.fftshift(
    scipy.fftpack.fft2(
        scipy.fftpack.ifftshift(
            target
            )
        )
    )

plt.subplot(223)
plt.title("Expected (real)")
plt.imshow(np.real(expected_fft2),vmin=-V,vmax=V)
plt.subplot(224)
plt.title("Expected (imag)")
plt.imshow(np.imag(expected_fft2),vmin=-V,vmax=V)

# Transform from 2D Fourier space back to a reconstruction of the target
recon=np.real(
    scipy.fftpack.fftshift(
        scipy.fftpack.ifft2(
            scipy.fftpack.ifftshift(fft2)
            )
        )
    )

plt.figure()
plt.title("Reconstruction")
plt.imshow(recon,vmin=0.0,vmax=1.0)
plt.gray()

plt.show()

Обновление 2013-02-17: Если вы были достаточно заинтересованы, чтобы пройтись по этому участку, вы можете найти дополнительную информацию о программе самообучения, частью которой она была, в виде этого плаката . Тело кода в этом репозитории также может представлять интерес (хотя обратите внимание, что код не так упрощен, как приведенный выше). Я могу попробовать и упаковать его как «ноутбук» IPython в какой-то момент.

timday
источник
3

Я не знаю точно, где проблема, но теорема среза означает, что эти два частных случая должны быть верными:

fft2(target)[0] = fft(sinogram[270])
fft2(target)[:,0] = fft(sinogram[0])

Так что следуйте вашему коду и попытайтесь найти точку, где они перестают быть эквивалентными, работая вперед от синограммы и обратно от сгенерированного 2D БПФ.

Это не выглядит правильно:

In [47]: angle(expected_fft2[127:130,127:130])
Out[47]: 
array([[-0.07101021,  3.11754929,  0.02299738],
       [ 3.09818784,  0.        , -3.09818784],
       [-0.02299738, -3.11754929,  0.07101021]])

In [48]: fft2_ = fft2_real+1.0j*fft2_imag

In [49]: angle(fft2_[127:130,127:130])
Out[49]: 
array([[ 3.13164353, -3.11056554,  3.11906449],
       [ 3.11754929,  0.        , -3.11754929],
       [ 3.11519503,  3.11056604, -2.61816765]])

2D БПФ, который вы генерируете, поворачивается на 90 градусов от того, что должно быть?

Я бы посоветовал работать с амплитудой и фазой, а не с реальными и воображаемыми, чтобы вам было легче видеть, что происходит:

введите описание изображения здесь

(Белые углы -инф от выполнения log(abs(0)), они не проблема)

эндолиты
источник
2

Я считаю , что фактическая теоретической причина , почему первым решение не работа исходит из того , что севооборот сделаны в отношении центров изображений, вызывая смещение [S/2, S/2], что означает , что каждый из строк из вашего sinogramне от 0до S, а точнее из -S/2к S/2. В вашем примере смещение на самом деле offset = np.floor(S/2.). Обратите внимание, что это работает для Sчетного или нечетного и эквивалентно тому, что вы сделали в своем коде S/2(хотя, будучи более явным, избегает проблем, например, когда Sесть float).

Я предполагаю, что фазовые задержки, которые этот сдвиг вводит в преобразовании Фурье (FT), являются источником того, о чем вы говорите во втором сообщении: фазы перепутаны, и нужно скомпенсировать этот сдвиг, чтобы иметь возможность применить инверсию преобразования Радона. Нужно углубиться в эту теорию, чтобы быть уверенным в том, что именно нужно для того, чтобы обратное сработало так, как ожидалось.

Чтобы компенсировать это смещение, вы можете либо использовать fftshift, как вы это делали (который помещает центр каждой строки в начало, и поскольку использование DFT фактически соответствует вычислению преобразования Фурье для S-периодического сигнала, вы получаете правильные данные ) или явно компенсировать этот эффект в комплексном преобразовании Фурье при вычислении sinogramFT. На практике вместо:

sinogram_fft_rows=scipy.fftpack.fftshift(
    scipy.fftpack.fft(
        scipy.fftpack.ifftshift(
            sinogram,
            axes=1
            )
        ),
    axes=1
    )

Вы можете удалить ifftshiftи умножить каждую строку на корректирующий вектор:

offset = np.floor(S/2.)
sinogram_fft_rows = scipy.fftpack.fftshift(
    scipy.fftpack.fft(sinogram, axis=1)
    * (np.exp(1j * 2.* np.pi * np.arange(S) * offset / S)),
    axes=1)

Это происходит из свойств преобразования Фурье при рассмотрении сдвига во времени (проверьте «теорему о сдвиге» на странице Википедии FT и примените сдвиг, равный - offset- потому что мы поместили изображение обратно вокруг центра).

Аналогично, вы можете применить ту же стратегию к реконструкции и заменить ее fftshiftкоррекцией фаз в обоих измерениях, но в другом направлении (компенсация назад):

recon=np.real(
    scipy.fftpack.ifft2(
        scipy.fftpack.ifftshift(fft2)
        *  np.outer(np.exp(- 1j * 2.* np.pi * np.arange(S) * offset / S),
                    np.exp(- 1j * 2.* np.pi * np.arange(S) * offset / S))
        )
    )

Что ж, это не улучшает ваше решение, а скорее проливает новый свет на теоретические аспекты вашего вопроса. Надеюсь, это поможет!

Кроме того, я не очень люблю использовать, fftshiftпотому что это имеет тенденцию возиться с тем, как fftвычисляется. В этом случае, однако, вам нужно поместить центр FT в центр изображения до интерполяции, чтобы получить fft2(или, по крайней мере, быть осторожным при настройке r- чтобы вы могли сделать это полностью- fftshiftбесплатно!), И это fftshiftдействительно удобно. там. Однако я предпочитаю использовать эту функцию в целях визуализации, а не в «ядре» вычислений. :-)

С наилучшими пожеланиями,

Жан-Луи

PS: вы пытались восстановить изображение, не обрезая круг? это дает довольно крутой эффект размытия углов, было бы неплохо иметь такую ​​функцию в таких программах, как Instagram, не так ли?

Жан-Луи Дюррио
источник