¿Qué tiene de malo este código para la reconstrucción tomográfica por el método de Fourier?

He estado jugando con algoritmos de reconstrucción tomográfica recientemente. Ya tengo buenas implementaciones de trabajo de FBP, ART, un esquema iterativo tipo SIRT / SART e incluso uso de álgebra lineal recta (¡lento!). Esta pregunta no se trata de ninguna de esas técnicas ; las respuestas de la forma "por qué alguien lo haría de esa manera, aquí hay algo de código FBP" no son lo que estoy buscando.

Lo siguiente que quería hacer con este programa era " completar el conjunto " e implementar el llamado " método de reconstrucción de Fourier ". Comprendo que esto es básicamente que usted aplica una FFT 1D a las "exposiciones" del sinograma, las organiza como "radios radiales de una rueda" en el espacio de Fourier 2D (que esto es algo útil que se sigue directamente del teorema del corte central) , interpolar desde esos puntos a una cuadrícula regular en ese espacio 2D, y luego debería ser posible invertir la transformación de Fourier para recuperar el objetivo de escaneo original.

Suena simple, pero no he tenido mucha suerte de conseguir reconstrucciones que se parezcan al objetivo original.

El siguiente código de Python (numpy / SciPy / Matplotlib) es la expresión más concisa que pude encontrar de lo que estoy tratando de hacer. Cuando se ejecuta, muestra lo siguiente:

Figura 1: el objetivo

Figura 2: un sinograma del objetivo

Figura 3: las filas del sinograma FFT-ed

Figura 4: la fila superior es el espacio FFT 2D interpolado de las filas del sinograma del dominio de Fourier; la fila inferior es (para fines de comparación) la FFT 2D directa del objetivo. Este es el punto en el que empiezo a sospechar; las gráficas interpoladas a partir de las FFT de sinograma se parecen a las gráficas hechas directamente con 2D-FFT al objetivo ... y, sin embargo, diferentes.

Figura 5: la transformada inversa de Fourier de la Figura 4. Esperaba que esto fuera un poco más reconocible como objetivo de lo que realmente es.

¿Alguna idea de lo que estoy haciendo mal? No estoy seguro si mi comprensión de la reconstrucción del método de Fourier es fundamentalmente defectuosa, o si solo hay algún error en mi código.

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()

OK, finalmente lo resolví.

El truco básicamente se redujo a poner algunos fftshift/ ifftshifts en el lugar correcto para que la representación del espacio 2D de Fourier no fuera muy oscilatoria y condenada a ser imposible de interpolar con precisión. Al menos eso es lo que creo que lo arregló. La mayor parte de la comprensión limitada que tengo de la teoría de Fourier se basa en la formulación integral continua, y siempre encuentro el dominio discreto y las FFT un poco ... extravagantes.

Si bien encuentro que el código de Matlab es bastante críptico, debo dar crédito a esta implementación por al menos darme la confianza de que este algoritmo de reconstrucción puede expresarse de manera razonablemente compacta en este tipo de entorno.

Primero mostraré resultados, luego código:

Figura 1: un nuevo objetivo más complejo.

Figura 2: el sinograma (OK OK, es la transformación de radón) del objetivo.

Figura 3: las filas FFT-ed del sinograma (trazadas con DC en el centro).

Figura 4: el sinograma FFT-ed transformado en espacio FFT 2D (DC en el centro). El color es una función de valor absoluto.

Figura 4a: amplíe el centro del espacio FFT 2D solo para mostrar mejor la naturaleza radial de los datos del sinograma.

Figura 5: Fila superior: el espacio FFT 2D interpolado de las filas de sinograma FFT-ed dispuestas radialmente. Fila inferior: la apariencia esperada de simplemente 2D FFT-ing el objetivo.

Figura 5a: amplíe la región central de las subtramas en la Fig. 5 para mostrar que parecen estar bastante bien cualitativamente.

Figura 6: Prueba ácida: la FFT 2D inversa del espacio FFT interpolado recupera el objetivo. Lena todavía se ve bastante bien a pesar de todo lo que acabamos de hacerle pasar (probablemente porque hay suficientes "radios" de sinograma para cubrir el plano 2D FFT bastante denso; las cosas se ponen interesantes si reduce la cantidad de ángulos de exposición, por lo que esto ya no es cierto )

Aquí está el código; muestra las tramas en menos de 15 segundos en SciPy de 64 bits de Debian / Wheezy en un 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()

Actualización 17-02-2013: si te ha interesado lo suficiente como para atravesar ese lote, puedes encontrar más resultados del programa de autoaprendizaje del que forma parte en este póster . El cuerpo del código en este repositorio también puede ser de interés (aunque tenga en cuenta que el código no es tan simple como el anterior). Puedo intentar volver a empaquetarlo como un "cuaderno" de IPython en algún momento.

No sé exactamente dónde está el problema, pero el teorema de corte significa que estos dos casos especiales deberían ser ciertos:

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

Por lo tanto, siga su código e intente encontrar el punto donde estos dejan de ser equivalentes, trabajando hacia adelante desde el sinograma y hacia atrás desde el FFT 2D generado.

Esto no se ve bien:

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]])

¿La FFT 2D que está generando gira 90 grados de lo que debería ser?

Sugeriría trabajar con magnitud y fase en lugar de real e imaginario, para que pueda ver más fácilmente lo que está sucediendo:

(Las esquinas blancas no se hacen log(abs(0)), no son un problema)

Creo que la razón teórica real de por qué la primera solución no hizo el trabajo viene del hecho de que las rotaciones se hacen con respecto a los centros de las imágenes, lo que induce un desplazamiento de [S/2, S/2], lo que significa que cada una de las filas de su sinogramno es de 0a S, sino más bien de -S/2a S/2. En su ejemplo, el desplazamiento es en realidad offset = np.floor(S/2.). Tenga en cuenta que esto funciona para Spares o impares, y es equivalente a lo que hizo en su código S/2(aunque ser más explícito evita problemas, cuando Ses un float, por ejemplo).

Supongo que los retrasos de fase que este cambio introduce en la transformación de Fourier (FT) están en el origen de lo que hablas en tu segundo mensaje: las fases están en mal estado, y uno necesita compensar ese cambio para poder aplicar la inversión de la transformada de radón. Uno necesita profundizar más en esa teoría para estar seguro de lo que se necesita exactamente para que el inverso funcione como se espera.

Para compensar ese desplazamiento, puede usar fftshift como lo hizo (lo que coloca el centro de cada fila al principio, y dado que usar DFT en realidad corresponde al cálculo de la Transformada de Fourier de una señal S-periódica, termina con el material correcto ), o compensar explícitamente este efecto en la transformada compleja de Fourier, al calcular el sinogramFT. En la práctica, en lugar de:

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

puede eliminar ifftshifty multiplicar cada fila por un vector correctivo:

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)

Esto proviene de las propiedades de transformación de Fourier, cuando se considera un cambio de tiempo (consulte la página de Wikipedia de FT para el "teorema de cambio", y solicite el cambio igual a - offset- porque colocamos la imagen de nuevo alrededor del centro).

Del mismo modo, puede aplicar la misma estrategia a la reconstrucción y reemplazarla fftshiftpor la corrección de las fases, en ambas dimensiones, pero en la otra dirección (compensación):

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))
        )
    )

Bueno, esto no mejora su solución, sino que arroja otra luz sobre los aspectos teóricos de su pregunta. ¡Espero que ayude!

Además, no me gusta mucho usarlo fftshiftporque tiende a confundirse con la forma en que fftse calcula. En este caso, sin embargo, debe colocar el centro del FT en el centro de la imagen antes de la interpolación para obtener fft2(o al menos tener cuidado al configurarlo r, para que pueda hacerlo completamente fftshiftgratis), y fftshiftrealmente es útil. allí. Sin embargo, prefiero mantener el uso de esa función para fines de visualización, y no dentro del "núcleo" de cómputo. :-)

Atentamente,

Jean-Louis

PD: ¿has intentado reconstruir la imagen sin recortar el círculo? eso da un efecto de desenfoque bastante bueno en las esquinas, sería bueno tener esa función en programas como Instagram, ¿no?