如何从4D图像中提取FFT的第一分量

2024-03-29 08:40:03 发布

您现在位置:Python中文网/ 问答频道 /正文

  • 我有4D(2D+沿z轴切片+时间帧)灰度图像,用于不同时刻的心脏跳动

  • 我确实喜欢沿时间轴进行傅里叶变换(分别针对每个切片),并分析基波(也称为H1分量,其中H代表希尔伯特空间),以便确定与ROI相对应的像素区域,该区域显示对心脏频率的最强响应

  • 我正使用python实现这一目的,我试图用下面的代码实现这一点,但我不确定这是否是正确的方法,因为我不知道如何确定截止频率以仅保持基波

This link to the image which I'm dealing with

import nibabel as nib
import numpy as np
import matplotlib.pyplot as plt

img = nib.load('patient057_4d.nii.gz')
f = np.fft.fft2(img)
#  Move the DC component of the FFT output to the center of the spectrum
fshift = np.fft.fftshift(f)
fshift_orig = fshift.copy()
# logarithmic transformation
 magnitude_spectrum = 20*np.log(np.abs(fshift))
# Create mask
rows, cols = img.shape
crow, ccol = int(rows/2), int(cols/2)
# Use mask to remove low frequency components
dist1 = 20
dist2 = 10
fshift[crow-dist1:crow+dist1, ccol-dist1:ccol+dist1] = 0
#fshift[crow-dist2:crow+dist2, ccol-dist2:ccol+dist2] = fshift_orig[crow-dist2:crow+dist2, ccol-dist2:ccol+dist2] 

# logarithmic transformation
magnitude_spectrum1 = 20*np.log(np.abs(fshift)) 
f_ishift = np.fft.ifftshift(fshift)
# inverse Fourier transform
img_back = np.fft.ifft2(f_ishift)
# get rid of imaginary part by abs
img_back = np.abs(img_back)
plt.figure(num = 'Im_Back')
plt.imshow(abs(fshift[:,:,2,2]).astype('uint8'),cmap='gray')
plt.show()


Tags: ofthetoimportfftimgasnp
1条回答
网友
1楼 · 发布于 2024-03-29 08:40:03
  • 解决方案是对每个切片分别进行三维傅里叶变换,然后只选择变换的第二个分量将其变换回空间,就是它
  • 这样做的好处是检测是否有东西沿着第三个轴移动(在我的例子中是时间)
for sl in range(img.shape[2]):
   #  -Fourier H1                    -
   # ff1[:, :, 1] H1 compnent 1, if 0 then DC
   ff1 = FFT.fftn(img[:,:,sl,:])
   fh = np.absolute(FFT.ifftn(ff1[:, :, 1])) 

   #  -Fourier H1                    -

相关问题 更多 >