I have 4D( 2D + slices along z axis + time frames) gray-scale image for the heart beating on different moments.
I do like to take Fourier Transform along time axis(for each slice separately), and analyze the fundamental Harmonic (also called H1 component, where H stands for Hilbert Space) so I can determine pixel regions corresponding to ROI which show strongest response to cardiac frequency.
I'm using python for this purpose, and I tried to do that with the following code, but I'm not sure that this is the correct way to do it, because I don't know how to determine the cut-frequency to keep only the fundamental Harmonic.
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()
```