Comments (4)
Thank you very much for your reply. I replicated your results with PULSAR_recon.py
using my own phantom data and realized that PULSAR for matlab is producing the same results. I got confused about the results using basic_grappa.py
because I forgot to add Gaussian noise to the k-space like I was doing in my matlab script. And thank you for the tip how to improve the image quality.
from pygrappa.
Hi @estab1 , thanks for reporting! I'll have more time to look at this a little later today, but I think I've run into relative scaling issues in the past, so normalizing both the reference and reconstruction before comparing may help. I've not used PULSAR before, but I can take a look and see if there are any obvious differences in reconstruction technique -- is the demo recon.m
script a good place to start?
from pygrappa.
Thanks for your quick reply. recon.m
contains all the parameters needed for the reconstruction. And grappa.m
contains the algorithm itself. So you probably need both to understand the reconstruction technique.
from pygrappa.
So I took a look at that recon.m
and grappa.m
and I think these are similar parameters to what's happening there:
- R=2
- 32 ACS lines
- 90% of FE lines used
The following modified basic_grappa.py
script attempts to replicate these:
PULSAR_recon.py
'''Replicate PULSAR recon.'''
import numpy as np
import matplotlib.pyplot as plt
from phantominator import shepp_logan
from pygrappa import grappa
if __name__ == '__main__':
# Generate fake sensitivity maps: mps
N = 128
ncoils = 4
xx = np.linspace(0, 1, N)
x, y = np.meshgrid(xx, xx)
mps = np.zeros((N, N, ncoils))
mps[..., 0] = x**2
mps[..., 1] = 1 - x**2
mps[..., 2] = y**2
mps[..., 3] = 1 - y**2
# generate 4 coil phantom
ph = shepp_logan(N)
imspace = ph[..., None]*mps
imspace = imspace.astype('complex')
ax = (0, 1)
kspace = np.fft.fftshift(np.fft.fft2(
np.fft.ifftshift(imspace, axes=ax), axes=ax), axes=ax)
# crop 32x90% window from the center of k-space for calibration
pd = 16
ctr = int(N/2)
calib = kspace[ctr-pd:ctr+pd, int(.05*N):-int(.05*N), :].copy()
# calibrate a kernel
kernel_size = (7, 7)
# undersample by a factor of 2 in ky
kspace[:, 1::2, ...] = 0
# reconstruct:
res = grappa(
kspace, calib, kernel_size, coil_axis=-1, lamda=0.01)
# replace calib
#res[ctr-pd:ctr+pd, int(.05*N):-int(.05*N), :] = calib
# SOS recon
res = np.sqrt(np.sum(np.abs(np.fft.fftshift(np.fft.ifft2(
np.fft.ifftshift(res, axes=ax), axes=ax), axes=ax))**2, axis=-1))
# SOS truth
truth = np.sqrt(np.sum(np.abs(imspace)**2, axis=-1))
residual = np.abs(res - truth)
plt.title('Residual | recon - truth |')
plt.imshow(residual, vmin=0, vmax=residual.flatten().max(), cmap='gray')
plt.colorbar()
plt.show()
The results I'm getting seem alright -- I was not able to run PULSAR because I don't have easy access to MATLAB. Do you know if it runs in Octave? A couple notes:
- The unmodified
basic_grappa.py
is reconstructing a different scenario: R=4 (Rx=Ry=2) with a smaller calibration region (20x20) - you are comparing the coil reconstructions to truth -- this in general will not look good because the SNR of any single coil is worse than the combined reconstruction. I recommend doing SOS on the reconstructed coil images for a good comparison
- try tuning the kernel size: larger kernels run slower but produce less residual, I found
(7, 7)
to be a good compromise - replacing the calibration region into your reconstructed k-space will help (if the calibration region came from your collected data, that is)
- The quality of the coil sensitivities might also be in question -- I generated 4 very simple non-complex coil sensitivities. You might try experimenting with more/better coils
EDIT: screen shots for posterity
Cranking kernel size and replacing calibration region:
from pygrappa.
Related Issues (20)
- pip installation fails in Mac HOT 8
- bad result in vcgrappa when R=3 HOT 6
- Benchmarking/testing with real data
- API references missing cgrappa
- cgrappa breaks on 7T real data
- CI error for python 3.5
- ENH: multidimensional data support for cgsense
- ValueError when using complex64 with cgrappa HOT 10
- ENH: revamp cgrappa HOT 1
- ERROR when installing pygrappa in a singularity image HOT 7
- Precision on the kernel size HOT 13
- igrappa doesn't return coil dimension
- Regularly undersampled GRAPPA HOT 1
- BUG: radialgrappaop should use lamda0, not lamda
- Slow performance of mdgrappa HOT 7
- skimage.util.pad not available HOT 2
- Keeping calibration weights for reuse? HOT 5
- Defination of kernel sizes not agree with most papers HOT 2
- Local utils package cannot be imported properly HOT 3
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from pygrappa.