python - 如何从 4D 图像中提取 FFT 的第一个分量
问题描述
我有4D(2D +沿z轴的切片+时间帧)灰度图像,用于不同时刻的心脏跳动。
我确实喜欢沿时间轴进行傅里叶变换(分别针对每个切片),并分析基本谐波(也称为 H1 分量,其中 H 代表希尔伯特空间),因此我可以确定对应于 ROI 的像素区域,这些区域对心脏的反应最强频率。
我为此目的使用python,并尝试使用以下代码执行此操作,但我不确定这是正确的方法,因为我不知道如何确定截止频率只保留基本谐波。
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()
解决方案
- 解决方案是对每个切片分别进行傅里叶变换 3D,然后仅选择变换的第二个分量将其变换回空间空间,仅此而已。
- 这样做的好处是检测是否有东西沿着第三轴移动(在我的例子中是时间)。
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-----------------------------------------
推荐阅读
- python - 字典中存在键时出现KeyError
- vba - 交叉表查询访问
- javascript - 如何在 Javascript 中通过 datepicker 过滤 JSON 数据?
- opencl - 英特尔 OpenCL 是否有任何 C++ 支持?
- python - PyQt5 中的多线程应用程序中的分段错误(错误)
- python - 将python中的类对象与指定条件进行比较
- python - 翻译 django 模型选择
- android - RenderFlex 溢出 Flutter
- javascript - Promise 在每个对象数组中处于待处理状态
- c - 尝试通过 C 中的结构实现堆栈,但得到以下代码的运行时错误。谁能解释并指出出了什么问题?