Merge pull request #64 from jonbmartin/tseg-bug-fix

SENSE time-segmented off-resonance correction bug fix
mikgroup · Oct 16, 2020 · 6e15434 · 6e15434
2 parents a125b77 + aa117d8
commit 6e15434
Show file tree

Hide file tree

Showing 2 changed files with 29 additions and 7 deletions.
diff --git a/sigpy/mri/util.py b/sigpy/mri/util.py
@@ -66,7 +66,7 @@ def tseg_off_res_b_ct(b0, bins, lseg, dt, T):
     """
 
     # create time vector
-    t = np.linspace(0, T, T/dt)
+    t = np.linspace(0, T, np.int(T/dt))
     hist_wt, bin_edges = np.histogram(np.imag(2j * np.pi * np.concatenate(b0)),
                                       bins)
 

diff --git a/tests/mri/test_linop.py b/tests/mri/test_linop.py
@@ -32,9 +32,6 @@ def test_sense_model(self):
         img = sp.randn(img_shape, dtype=np.complex)
         mps = sp.randn(mps_shape, dtype=np.complex)
 
-        mask = np.zeros(img_shape)
-        mask[::2, ::2] = 1.0
-
         A = linop.Sense(mps)
 
         check_linop_adjoint(A, dtype=np.complex)
@@ -49,9 +46,6 @@ def test_sense_model_batch(self):
         img = sp.randn(img_shape, dtype=np.complex)
         mps = sp.randn(mps_shape, dtype=np.complex)
 
-        mask = np.zeros(img_shape, dtype=np.complex)
-        mask[::2, ::2] = 1.0
-
         A = linop.Sense(mps, coil_batch_size=1)
         check_linop_adjoint(A, dtype=np.complex)
         npt.assert_allclose(sp.fft(img * mps, axes=[-1, -2]),
@@ -73,6 +67,34 @@ def test_noncart_sense_model(self):
         npt.assert_allclose(sp.fft(img * mps, axes=[-1, -2]).ravel(),
                             (A * img).ravel(), atol=0.1, rtol=0.1)
 
+    def test_sense_tseg_off_res_model(self):
+        img_shape = [16, 16]
+        mps_shape = [8, 16, 16]
+
+        img = sp.randn(img_shape, dtype=np.complex)
+        mps = sp.randn(mps_shape, dtype=np.complex)
+
+        y, x = np.mgrid[:16, :16]
+        coord = np.stack([np.ravel(y - 8), np.ravel(x - 8)], axis=1)
+        coord = coord.astype(np.float)
+
+        d = np.sqrt(x * x + y * y)
+        sigma, mu, a = 2, 0.25, 400
+        b0 = a * np.exp(-((d - mu) ** 2 / (2.0 * sigma ** 2)))
+        tseg = {"b0": b0, "dt": 4e-6, "lseg": 1, "n_bins": 10}
+
+        F = sp.linop.NUFFT(mps_shape, coord)
+        b, ct = sp.mri.util.tseg_off_res_b_ct(b0=b0, bins=10, lseg=1, dt=4e-6,
+                                              T=coord.shape[0] * 4e-6)
+        B1 = sp.linop.Multiply(F.oshape, b.T)
+        Ct1 = sp.linop.Multiply(img_shape, ct.reshape(img_shape))
+        S = sp.linop.Multiply(img_shape, mps)
+
+        A = linop.Sense(mps, coord=coord, tseg=tseg)
+
+        check_linop_adjoint(A, dtype=np.complex)
+        npt.assert_allclose(B1 * F * S * Ct1 * img, A * img)
+
     def test_noncart_sense_model_batch(self):
         img_shape = [16, 16]
         mps_shape = [8, 16, 16]