From 026b2a46a91e8f3f7ed652a11eacbe5c1ab66262 Mon Sep 17 00:00:00 2001
From: steven <stevenxcao@gmail.com>
Date: Fri, 20 Jul 2018 16:11:09 -0700
Subject: [PATCH] Fixed bug in qsm recon

---
 functions.py |  4 ++++
 qsm_star.py  | 15 ++++++++++++---
 test.py      | 17 +++++++++++++++--
 3 files changed, 31 insertions(+), 5 deletions(-)

diff --git a/functions.py b/functions.py
index 1f38bbb..a8a836c 100644
--- a/functions.py
+++ b/functions.py
@@ -60,6 +60,10 @@ def left_pad(arr, axis):
     pad_size = tuple([(int(i == axis), 0) for i, a in enumerate(arr.shape)])
     return np.pad(arr, pad_size, 'constant')
 
+def right_pad(arr, axis):
+    pad_size = tuple([(0, int(i == axis)) for i, a in enumerate(arr.shape)])
+    return np.pad(arr, pad_size, 'constant')
+
 def unpad(arr, pad_size):
     pad_size = tuple(pad_size)
     while len(pad_size) < len(arr.shape):
diff --git a/qsm_star.py b/qsm_star.py
index 7db7bca..3c94478 100644
--- a/qsm_star.py
+++ b/qsm_star.py
@@ -3,10 +3,19 @@
 
 def qsm_star(phase, mask, voxel_size, TE, pad_size = (8,8,20), 
              B0 = 3, B0_dir = (0,0,1), tau = 1e-6, d2_thresh = .065):
-    d2 = calc_d2_matrix(phase.shape, voxel_size, B0_dir)
-    # sample = np.abs(d2) > d2_thresh
-    d2 = pad(d2, pad_size)
     phase, mask = pad(phase, pad_size), pad(mask, pad_size)
+    d2 = calc_d2_matrix(phase.shape, voxel_size, B0_dir)
+    
+#    sample = np.abs(d2) <= d2_thresh
+#    d2_thresh = np.zeros(d2.shape)
+#    d2_thresh[np.logical_and(sample, d2 >= 0)] = d2_thresh
+#    d2_thresh[np.logical_and(sample, d2 < 0)] = -d2_thresh
+#    phase_k = fftnc(phase)
+#    x_k = fftnc(ifftnc(phase_k / d2_thresh) * mask)
+#    x_k *= np.logical_not(sample).astype(int)
+#    x = ifftnc(x_k) * mask
+#    x = ifftnc(x * np.logical_not(sample)) * mask
+    
     A = lambda x: ifftnc(d2 * fftnc(x * mask))
     AT = A
     susc = sparsa(phase, np.zeros(phase.shape), A, AT, tau)
diff --git a/test.py b/test.py
index f77b917..57edb89 100644
--- a/test.py
+++ b/test.py
@@ -1,6 +1,6 @@
 import numpy as np
 import nibabel as nib
-# from matplotlib import pyplot as plt
+from matplotlib import pyplot as plt
 
 from read_dicom import read_dicom
 from create_mask import create_mask_qsm
@@ -17,4 +17,17 @@
 unwrapped = np.sum(unwrapped, axis = 3)
 tissue, mask2 = v_sharp(unwrapped, mask, voxel_size)
 susc = qsm_star(tissue, mask2, voxel_size, TE)
-# plt.gray()
\ No newline at end of file
+# plt.gray()
+
+def save(file, img):
+    img = np.rot90(img)
+    img = np.flip(img, 0)
+    image = nib.Nifti1Image(img, np.eye(4))
+    nib.save(image, file)
+    
+save('N041/tissue.nii', tissue)
+save('N041/susc.nii', susc)
+save('N041/mask.nii', mask.astype(int))
+save('N041/mask2.nii', mask2.astype(int))
+save('N041/phase.nii', phase)
+save('N041/unwrapped.nii', unwrapped)
\ No newline at end of file