Quel est le problème avec ce code pour la reconstruction tomographique par la méthode de Fourier?

J'ai récemment joué avec des algorithmes de reconstruction tomographique. J'ai déjà de belles implémentations de travail de FBP, ART, un schéma itératif de type SIRT / SART et même en utilisant l'algèbre linéaire droite (lente!). Cette question ne concerne aucune de ces techniques ; les réponses du formulaire "pourquoi quelqu'un le ferait-il de cette façon, voici du code FBP à la place" ne sont pas ce que je recherche.

La prochaine chose que je voulais faire avec ce programme était de « compléter l'ensemble » et de mettre en œuvre la soi-disant « méthode de reconstruction de Fourier ». Ma compréhension de cela est essentiellement que vous appliquez une FFT 1D aux "expositions" sinogrammes, que vous les organisez comme des "rayons d'une roue" radiaux dans l'espace de Fourier 2D (que c'est une chose utile à faire découle directement du théorème de la tranche centrale) , interpoler à partir de ces points vers une grille régulière dans cet espace 2D, puis il devrait être possible d'inverser la transformée de Fourier pour récupérer la cible de balayage d'origine.

Cela semble simple, mais je n'ai pas eu beaucoup de chance d'obtenir des reconstructions qui ressemblent à la cible d'origine.

Le code Python (numpy / SciPy / Matplotlib) ci-dessous est à propos de l'expression la plus concise que j'ai pu trouver de ce que j'essaie de faire. Lorsqu'il est exécuté, il affiche les éléments suivants:

Figure 1: la cible

Figure 2: un sinogramme de la cible

Figure 3: les rangées de sinogrammes FFT-ed

Figure 4: la rangée du haut est l'espace FFT 2D interpolé à partir des rangées de sinogrammes du domaine de Fourier; la ligne du bas est (à des fins de comparaison) la FFT 2D directe de la cible. C'est le moment où je commence à devenir suspect; les tracés interpolés à partir des FFT sinogrammes ressemblent aux tracés réalisés en FFT 2D directement sur la cible ... et pourtant différents.

Figure 5: la transformée de Fourier inverse de la figure 4. J'aurais espéré que ce serait un peu plus reconnaissable comme cible qu'elle ne l'est réellement.

Des idées sur ce que je fais mal? Je ne sais pas si ma compréhension de la reconstruction de la méthode de Fourier est fondamentalement erronée, ou il y a juste un bug dans mon code.

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, je l'ai finalement craqué.

L'astuce consistait essentiellement à mettre certains fftshift/ ifftshifts au bon endroit afin que la représentation de l'espace de Fourier 2D ne soit pas extrêmement oscillatoire et vouée à être impossible à interpoler avec précision. C'est du moins ce que je pense avoir corrigé. La majeure partie de ma compréhension limitée de la théorie de Fourier est basée sur la formulation intégrale continue, et je trouve toujours le domaine discret et les FFT un peu ... excentriques.

Bien que je trouve le code matlab plutôt cryptique, je dois créditer cette implémentation pour au moins me donner la confiance que cet algorithme de reconstruction peut être exprimé de manière assez compacte dans ce type d'environnement.

Je vais d'abord montrer les résultats, puis coder:

Figure 1: une nouvelle cible plus complexe.

Figure 2: le sinogramme (OK OK, c'est la transformée de Radon) de la cible.

Figure 3: les lignes FFT-ed du sinogramme (tracé avec DC au centre).

Figure 4: le sinogramme FFT-ed transformé en espace FFT 2D (DC au centre). La couleur est fonction de la valeur absolue.

Figure 4a: Zoomez au centre de l'espace FFT 2D juste pour mieux montrer la nature radiale des données du sinogramme.

Figure 5: rangée du haut: l'espace FFT 2D interpolé à partir des rangées de sinogrammes FFT éditées radialement. Rangée du bas: l'apparence attendue d'une FFT simplement 2D sur la cible.

Figure 5a: Zoomez sur la région centrale des sous-parcelles de la Fig5 pour montrer qu'elles semblent être en assez bon accord qualitativement.

Figure 6: Test d'acide: FFT 2D inverse de l'espace FFT interpolé récupère la cible. Lena a toujours l'air bien malgré tout ce que nous venons de lui faire comprendre (probablement parce qu'il y a suffisamment de "rayons" sinogrammes pour couvrir le plan FFT 2D de manière assez dense; les choses deviennent intéressantes si vous réduisez le nombre d'angles d'exposition, ce n'est donc plus vrai ).

Voici le code; fait apparaître les tracés en moins de 15 secondes sur Debian / Wheezy 64 bits SciPy sur 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()

Mise à jour du 17/02/2013: Si vous avez été suffisamment intéressé pour parcourir ce lot, vous trouverez plus de résultats du programme d'autoformation dont il faisait partie sous la forme de cette affiche . Le corps de code dans ce référentiel peut également être intéressant (bien que le code ne soit pas aussi rationalisé que ci-dessus). Je peux essayer de le reconditionner en tant que "bloc-notes" IPython à un moment donné.

Je ne sais pas exactement où est le problème, mais le théorème de la tranche signifie que ces deux cas spéciaux devraient être vrais:

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

Suivez donc votre code et essayez de trouver le point où ceux-ci cessent d'être équivalents, en travaillant à partir du sinogramme et en arrière à partir de la FFT 2D générée.

Cela ne semble pas correct:

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 vous générez est tournée de 90 degrés par rapport à ce qu'elle devrait être?

Je suggère de travailler avec l'amplitude et la phase plutôt que réel et imaginaire, afin que vous puissiez voir plus facilement ce qui se passe:

(Les coins blancs sont -inf de faire log(abs(0)), ils ne sont pas un problème)

Je crois que la véritable raison théorique pour laquelle la première solution n'a pas fonctionné vient du fait que les rotations sont effectuées par rapport aux centres des images, induisant un décalage de [S/2, S/2], ce qui signifie que chacune des lignes de votre sinogramn'est pas de 0à S, mais plutôt de -S/2à S/2. Dans votre exemple, le décalage est en fait offset = np.floor(S/2.). Notez que cela fonctionne pour Spair ou impair, et est équivalent à ce que vous avez fait dans votre code S/2(bien qu'être plus explicite évite les problèmes, quand Sest un float, par exemple).

Je suppose que les retards de phase que ce décalage introduit dans la transformée de Fourier (FT) sont à l'origine de ce dont vous parlez dans votre deuxième message: les phases sont gâchées, et il faut compenser ce décalage afin de pouvoir appliquer l'inversion de la transformée de Radon. Il faut approfondir cette théorie pour être sûr de ce qui est exactement nécessaire pour que l'inverse fonctionne comme prévu.

Pour compenser ce décalage, vous pouvez soit utiliser fftshift comme vous l'avez fait (ce qui place le centre de chaque ligne au début, et puisque l'utilisation de la DFT correspond en fait au calcul de la transformée de Fourier d'un signal S-périodique, vous vous retrouvez avec le bon truc ), ou compenser explicitement cet effet dans la transformée de Fourier complexe, lors du calcul du sinogramFT. En pratique, au lieu de:

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

vous pouvez supprimer le ifftshiftet multiplier chaque ligne par un vecteur correctif:

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)

Cela vient des propriétés de la transformée de Fourier, quand on considère un décalage temporel (consultez la page wikipedia FT pour le "théorème de décalage", et appliquez le décalage égal à - offset- parce que nous remettons l'image autour du centre).

De même, vous pouvez appliquer la même stratégie à la reconstruction, et la remplacer fftshiftpar la correction des phases, dans les deux dimensions, mais dans l'autre sens (compensation arrière):

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

Eh bien, cela n'améliore pas votre solution, mais éclaire plutôt les aspects théoriques de votre question. J'espère que cela pourra aider!

De plus, je n'aime pas tellement utiliser fftshiftcar il a tendance à jouer avec la façon dont le fftest calculé. Dans ce cas, cependant, vous devez placer le centre du FT au centre de l'image avant l'interpolation pour obtenir fft2(ou du moins soyez prudent lors du réglage r- afin que vous puissiez le rendre complètement- fftshiftgratuit!), Et fftshiftest vraiment pratique Là. Je préfère cependant conserver l'utilisation de cette fonction à des fins de visualisation, et non au sein du "noyau" de calcul. :-)

Meilleures salutations,

Jean-Louis

PS: avez-vous essayé de reconstruire l'image sans recadrer le cercle? cela donne un effet de flou assez cool dans les coins, ce serait bien d'avoir une telle fonctionnalité dans des programmes comme Instagram, non?