From 3ca27534e867f7f5d649793d6791840732b4aaea Mon Sep 17 00:00:00 2001
From: Tim <jbt5jf@virginia.edu>
Date: Mon, 22 Apr 2024 00:51:54 -0400
Subject: [PATCH] minor edits and add recon test file

---
 brkraw/lib/recon.py     | 35 ++++++++++++++++---------------
 tests/recon_api_test.py | 46 +++++++++++++++++++++++++++++++++++++++++
 2 files changed, 64 insertions(+), 17 deletions(-)
 create mode 100644 tests/recon_api_test.py

diff --git a/brkraw/lib/recon.py b/brkraw/lib/recon.py
index 144f72c..9efec11 100644
--- a/brkraw/lib/recon.py
+++ b/brkraw/lib/recon.py
@@ -24,6 +24,21 @@
 SUPPORTED_PROTOCOLS = ['rare','localizer' ,'gre', 'msme',      
                        'mge','dess', 'fisp', 'flash']
 
+def reconstruction(scanobj,process='image', **kwargs):
+    # Ensure Scans are Image Based
+    ACQ_dim_desc = [scanobj.acqp.get('ACQ_dim_desc')] if isinstance(scanobj.acqp.get('ACQ_dim_desc'), str) else scanobj.acqp.get('ACQ_dim_desc')
+    if 'Spectroscopic' in ACQ_dim_desc:
+        warnings.warn('Scan is spectroscopic')
+        process = 'readout' 
+    
+    # Reconstruction Processing
+    recoObj = Reconstruction(scanobj)
+    if process == 'readout':
+        return recoObj.sort_fid()
+    elif process == 'kspace':
+        return recoObj.process_kspace()
+    return recoObj.reconstruct(rms=kwargs['rms'] if 'rms' in kwargs.keys() else True) 
+
 class Reconstruction:
     def __init__(self, scanobj:'ScanObj', reco_id:'int'=1) -> None: 
         self.acqp       = scanobj.acqp
@@ -177,9 +192,9 @@ def process_kspace(self):
                                                               map_index[(NI+1)*(NR+1)-1])[:,:,:,np.newaxis],
                                                  [1,1,1,self.NRecs])
         
-        # Zeropad KSPACE
-        newdata_dims=[1, 1, 1]
+        # Zeropad KSPACE 
         RECO_ft_size = self.reco.get('RECO_ft_size')
+        newdata_dims=[1, 1, 1]
         newdata_dims[0:len(RECO_ft_size)] = RECO_ft_size
         newdata = np.zeros(shape=newdata_dims+[self.NRecs, self.NI, self.NR], dtype=complex)
         for NR in range(self.NR):
@@ -205,18 +220,4 @@ def reconstruct(self, kspace=None, rms=True):
         if rms:
             image = np.sqrt(np.mean(np.square(np.abs(image)), axis=3))
         return image
-    
-def reconstruction(scanobj,process='image', **kwargs):
-    # Ensure Scans are Image Based
-    ACQ_dim_desc = [scanobj.acqp.get('ACQ_dim_desc')] if isinstance(scanobj.acqp.get('ACQ_dim_desc'), str) else scanobj.acqp.get('ACQ_dim_desc')
-    if 'Spectroscopic' in ACQ_dim_desc:
-        warnings.warn('Scan is spectroscopic')
-        process = 'readout' 
-    
-    # Reconstruction Processing
-    recoObj = Reconstruction(scanobj)
-    if process == 'readout':
-        return recoObj.sort_fid()
-    elif process == 'kspace':
-        return recoObj.process_kspace()
-    return recoObj.reconstruct(rms=kwargs['rms'] if 'rms' in kwargs.keys() else True) 
\ No newline at end of file
+    
\ No newline at end of file
diff --git a/tests/recon_api_test.py b/tests/recon_api_test.py
new file mode 100644
index 0000000..169cde4
--- /dev/null
+++ b/tests/recon_api_test.py
@@ -0,0 +1,46 @@
+from brkraw.app import tonifti as tonii
+from brkraw.lib.recon import reconstruction
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+import nibabel as nib
+
+fiddata_path = 'YOUR.PvDatasets'
+studyobj = tonii.brkraw.BrkrawToNifti(fiddata_path)
+scan_id  = 10 # REPLACE WITH YOUR SCAN NUMBER
+reco_id  = 1
+scanobj  = studyobj.get_scan(scan_id)
+
+reconobj = np.abs(np.squeeze(reconstruction(scanobj, rms=True)))
+affine   = studyobj.get_affine(scan_id, reco_id)
+dataobj  = studyobj.get_dataobj(scan_id, reco_id)
+
+print(reconobj.shape)
+plt.subplot(1,2,1)
+if len(reconobj.shape) == 2:
+    plt.imshow(reconobj[:,:].T)
+    plt.title(f"{scan_id},{scanobj.acqp['ACQ_scan_name']}")
+elif len(reconobj.shape) == 3:
+    plt.imshow(reconobj[:,:,reconobj.shape[2]//2].T)
+    plt.title(f"{scan_id},{scanobj.acqp['ACQ_scan_name']}")
+elif len(reconobj.shape) == 4:
+    plt.imshow(reconobj[:,:,reconobj.shape[2]//2,0].T)
+    plt.title(f"{scan_id},{scanobj.acqp['ACQ_scan_name']}")  
+plt.subplot(1,2,2)
+if len(reconobj.shape) == 2:
+    plt.imshow(dataobj[:,:].T)
+    plt.title(f"{scan_id},dataobj")
+elif len(reconobj.shape) == 3:
+    plt.imshow(dataobj[:,:,dataobj.shape[2]//2].T)
+    plt.title(f"{scan_id},dataobj")
+elif len(reconobj.shape) == 4:
+    plt.imshow(dataobj[:,:,dataobj.shape[2]//2,0].T)
+    plt.title(f"{scan_id},dataobj")  
+plt.show()
+
+print(reconobj.shape, dataobj.shape) 
+assert np.prod(dataobj.shape) == np.prod(reconobj.shape), "Shape mismatched"
+
+niiobj = nib.Nifti1Image(reconobj/np.max(np.abs(reconobj)), affine)
+niiobj.to_filename('reconfile.nii.gz')