mikgroup / sigpy

Python package for signal processing, with emphasis on iterative methods
BSD 3-Clause "New" or "Revised" License
304 stars 93 forks source link

Radial sample from cartesian K-Space data #107

Closed kingaza closed 2 years ago

kingaza commented 2 years ago

Describe the bug I try to interpolate the radial data from a cartesian k-space, and the coord data is generated by sigpy.mri.radial() But the reconstruction seems quite strange, and the sampled k-space value is not as expected

To Reproduce

import numpy as np
import matplotlib.pyplot as plt

import sigpy as sp
from sigpy.mri import samp
import sigpy.plot as pl

def fftc_2d(_data):
    data = _data.copy()
    data = np.fft.ifftshift(data)
    data = np.fft.fft2(data, norm="ortho")
    return np.fft.fftshift(data)    

def ifftc_2d(_data):
    data = _data.copy()
    data = np.fft.ifftshift(data)
    data = np.fft.ifft2(data, norm="ortho")
    return np.fft.fftshift(data)    

def show_image(image):
    plt.figure()
    plt.imshow(np.abs(image), cmap='gray')
    plt.show()

def show_kspace(kspace):
    plt.figure()
    plt.imshow(np.log(np.abs(kspace)), cmap='gray')
    plt.show()

# create shepp-logan image and k-space
image_shape = (160, 160)
shepp_image = sp.shepp_logan(image_shape)
show_image(shepp_image)
shepp_kspace = fftc_2d(shepp_image)
show_kspace(shepp_kspace)

# sample as radial style
spoke_num, sample_num = 60, 160
coord_shepp = samp.radial([spoke_num, sample_num, 2], image_shape)
pl.ScatterPlot(coord_shepp, title='Trajectory')

# if no offset to the coord, no image can be reconstructed
# BUT the k-space is very strange (the data in middle column shall have the largest magitude)
coord_offset = 80
kspace_nu = sp.interpolate(shepp_kspace, coord_shepp + coord_offset)
show_kspace(kspace_nu)

dcf_nu = (coord_shepp[..., 0]**2 + coord_shepp[..., 1]**2)**0.5
image_nu = sp.nufft_adjoint(kspace_nu * dcf_nu, coord_shepp)
show_image(image_nu)

Expected behavior A clear and concise description of what you expected to happen.

Screenshots If applicable, add screenshots to help explain your problem.

Desktop (please complete the following information):

kingaza commented 2 years ago

refer test_dcf