Fixed bug in qsm recon

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):

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)

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()
+# 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)