Eu tenho brincado com algoritmos de reconstrução tomográfica recentemente. Eu já tenho boas implementações de trabalho de FBP, ART, um esquema iterativo semelhante ao SIRT / SART e até mesmo usando álgebra linear direta (lenta!). Esta pergunta não é sobre nenhuma dessas técnicas ; respostas do formulário "por que alguém faria dessa maneira, aqui está um código FBP?" não é o que estou procurando.
A próxima coisa que eu queria fazer com este programa era " completar o conjunto " e implementar o chamado " método de reconstrução de Fourier ". Meu entendimento disso é basicamente que você aplica uma FFT 1D às "exposições" do sinograma, organiza-as como "raios de uma roda" radiais no espaço 2D de Fourier (que é algo útil a seguir segue diretamente do teorema da fatia central) , interpole a partir desses pontos para uma grade regular nesse espaço 2D e, em seguida, será possível inverter a transformação de Fourier para recuperar o destino da verificação original.
Parece simples, mas não tive muita sorte em fazer reconstruções parecidas com o alvo original.
O código Python (numpy / SciPy / Matplotlib) abaixo é sobre a expressão mais concisa que eu poderia ter do que estou tentando fazer. Quando executado, ele exibe o seguinte:
Figura 1: o alvo
Figura 2: um sinograma do alvo
Figura 3: as linhas de sinograma da FFT-ed
Figura 4: a linha superior é o espaço FFT 2D interpolado das linhas de sinograma do domínio Fourier; a linha inferior é (para fins de comparação) a FFT 2D direta do alvo. Este é o ponto em que estou começando a suspeitar; as plotagens interpoladas das FFTs do sinograma são semelhantes às plotagens feitas diretamente por 2D-FFTs no alvo ... e ainda assim diferentes.
Figura 5: a transformação inversa de Fourier da Figura 4. Eu esperava que isso fosse um pouco mais reconhecível como o alvo do que realmente é.
Alguma idéia do que estou fazendo de errado? Não tenho certeza se meu entendimento da reconstrução do método de Fourier é fundamentalmente falho ou se há algum erro no meu 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()
Respostas:
OK, eu finalmente resolvi isso.
O truque basicamente foi colocar alguns
fftshift
/ifftshift
s no lugar certo, de modo que a representação do espaço 2D Fourier não fosse extremamente oscilatória e condenada a ser impossível interpolar com precisão. Pelo menos é o que acho que consertou. A maior parte do entendimento limitado que tenho da teoria de Fourier se baseia na formulação integral contínua, e sempre acho o domínio discreto e as FFTs um pouco ... peculiares.Embora eu ache o código do matlab bastante enigmático, tenho que creditar essa implementação por pelo menos me dar a confiança de que esse algoritmo de reconstrução pode ser expresso de maneira razoavelmente compacta nesse tipo de ambiente.
Primeiro, mostrarei os resultados e depois codificarei:
Figura 1: um novo alvo, mais complexo.
Figura 2: o sinograma (OK OK, é a transformação Radon) do alvo.
Figura 3: as linhas do sinograma com FFT-ed (plotadas com CD no centro).
Figura 4: o sinograma FFT-ed transformado em espaço 2D FFT (CD no centro). A cor é uma função do valor absoluto.
Figura 4a: Aumente o zoom no centro do espaço FFT 2D apenas para mostrar melhor a natureza radial dos dados do sinograma.
Figura 5: Linha superior: o espaço 2D da FFT interpolado das linhas de sinograma dispostas radialmente na FFT-ed. Linha inferior: a aparência esperada da simples FFT 2D no alvo.
Figura 5a: Aumente o zoom na região central das subparcelas na Fig5 para mostrar que estas parecem estar em concordância qualitativa.
Figura 6: Teste de ácido: FFT 2D inverso do espaço interpolado da FFT recupera o alvo. Lena ainda parece muito boa, apesar de tudo que a colocamos (provavelmente porque existem "raios" de sinograma suficientes para cobrir o plano FFT 2D de maneira bastante densa; as coisas ficam interessantes se você reduzir o número de ângulos de exposição, para que isso não seja mais verdade) )
Aqui está o código; traz os gráficos em menos de 15s no SciPy de 64 bits do Debian / Wheezy em um i7.
Atualização 17-02-2013: se você estiver interessado o suficiente para percorrer esse lote, poderá encontrar mais resultados do programa de auto-estudo do qual fazia parte, na forma deste pôster . O corpo do código neste repositório também pode ser interessante (embora observe que o código não é tão otimizado quanto o descrito acima). Posso tentar reembalá-lo como um "notebook" do IPython em algum momento.
fonte
Não sei exatamente onde está o problema, mas o teorema da fatia significa que esses dois casos especiais devem ser verdadeiros:
Portanto, siga seu código e tente encontrar o ponto em que eles param de ser equivalentes, trabalhando para frente a partir do sinograma e para trás a partir da FFT 2D gerada.
Isso não parece certo:
O FFT 2D que você está gerando é girado 90 graus do que deveria ser?
Sugiro trabalhar com magnitude e fase, em vez de real e imaginário, para que você possa ver mais facilmente o que está acontecendo:
(Os cantos brancos estão desinformados
log(abs(0))
, não são um problema)fonte
Acredito que a razão teórica real porque a primeira solução não funcionou vem do fato de que as rotações são feitas em relação aos centros das imagens, induzindo um deslocamento de
[S/2, S/2]
, o que significa que cada uma das linhas da suasinogram
não é a partir0
deS
, mas sim de-S/2
paraS/2
. No seu exemplo, o deslocamento é realmenteoffset = np.floor(S/2.)
. Observe que isso funciona paraS
pares ou ímpares e é equivalente ao que você fez no seu códigoS/2
(embora seja mais explícito evite problemas, quandoS
é umfloat
, por exemplo).Meu palpite é que os atrasos de fase que essa mudança introduz na transformada de Fourier (FT) estão na origem do que você fala em sua segunda mensagem: as fases estão confusas e é necessário compensar essa mudança para poder aplique a inversão da transformação Radon. É preciso cavar mais nessa teoria para ter certeza do que é exatamente necessário para que o inverso funcione conforme o esperado.
Para compensar esse deslocamento, você pode usar o fftshift como fez (o que coloca o centro de cada linha no início e, como o uso da DFT realmente corresponde ao cálculo da Transformada de Fourier de um sinal periódico S, você acaba com as coisas certas ) ou compense explicitamente esse efeito na complexa transformação de Fourier, ao calcular o
sinogram
FT. Na prática, em vez de:você pode remover
ifftshift
e multiplicar cada linha por um vetor corretivo:Isso vem das propriedades da transformação de Fourier, ao considerar um deslocamento no tempo (verifique a página da Wikipédia do FT para o "teorema do deslocamento" e aplique o deslocamento igual a
- offset
- porque colocamos a imagem de volta no centro).Da mesma forma, você pode aplicar a mesma estratégia à reconstrução e substituir
fftshift
pela correção das fases, em ambas as dimensões, mas na outra direção (compensação de retorno):Bem, isso não melhora sua solução, mas lança outra luz sobre os aspectos teóricos da sua pergunta. Espero que ajude!
Além disso, não gosto muito de usá-
fftshift
lo porque ele tende a mexer com a maneira comofft
é calculado. Nesse caso, no entanto, você precisa colocar o centro do FT no centro da imagem antes da interpolaçãofft2
(ou pelo menos tenha cuidado ao definirr
- para poder torná-lo completamentefftshift
livre!), E éfftshift
realmente útil há. No entanto, prefiro manter o uso dessa função para fins de visualização, e não dentro do "núcleo" da computação. :-)Cumprimentos,
Jean-Louis
PS: você tentou reconstruir a imagem sem cortar o círculo? que dá um efeito de desfoque bem legal nos cantos, seria bom ter esse recurso em programas como o Instagram, não é?
fonte