Diff of /cmac/aux_dicom.py [000000] .. [4cda31]

Switch to side-by-side view

--- a
+++ b/cmac/aux_dicom.py
@@ -0,0 +1,90 @@
+import os
+import glob
+import pydicom
+import numpy as np
+import nibabel as nib
+from dipy.align.reslice import reslice
+from dipy.align.imaffine import AffineMap
+
+def extract_cosines(ImageOrientationPatient):
+    row_cosine    = np.array(ImageOrientationPatient[:3])
+    column_cosine = np.array(ImageOrientationPatient[3:])
+    slice_cosine  = np.cross(row_cosine, column_cosine)
+    return row_cosine, column_cosine, slice_cosine
+
+def read_RefDs(dicom_folder_path):
+    cine_dicom_files = glob.glob(dicom_folder_path)
+    SliceLocations   = [pydicom.read_file(cine_dicom_file).SliceLocation 
+                        for cine_dicom_file in cine_dicom_files]
+    SliceOriginID    = np.array(SliceLocations, dtype=float).argmin()
+    
+    RefDs = pydicom.read_file(cine_dicom_files[SliceOriginID])
+    ImageOrientationPatient = RefDs.ImageOrientationPatient
+    ImagePositionPatient    = np.array(list(RefDs.ImagePositionPatient), dtype=float)
+    Zooms                   = np.array(list(RefDs.PixelSpacing) + [RefDs.SliceThickness], dtype=float)
+    
+    return ImageOrientationPatient, ImagePositionPatient, Zooms
+
+def read_affine_info(dicom_folder_path, viewer='slicer'):
+    
+    ImageOrientationPatient, ImagePositionPatient, Zooms = read_RefDs(dicom_folder_path)
+    
+    affine_axial        = np.diag(list(Zooms)+[1])
+    affine_axial[:3,3] += ImagePositionPatient
+    
+    row_cos, column_cos, slice_cos = extract_cosines(ImageOrientationPatient)
+    
+    ijk2ras = np.stack((row_cos, column_cos, slice_cos))
+    if viewer == "slicer":
+        ijk2ras = (ijk2ras*np.array([-1,-1,1])).T
+        ImagePositionPatient = ImagePositionPatient*np.array([-1,-1,1])
+
+    affine  = np.stack((ijk2ras[:,0]*Zooms[0],
+                        ijk2ras[:,1]*Zooms[1],
+                        ijk2ras[:,2]*Zooms[2],
+                        ImagePositionPatient), axis=1)
+    
+    return ijk2ras, np.vstack((affine,[[0,0,0,1]])), affine_axial   
+
+def reslice_3d(domain_nifti, codomain_nifti, mode='linear'):
+    
+    affine_map = AffineMap(affine=np.eye(4),
+                  domain_grid_shape=domain_nifti.shape, domain_grid2world=domain_nifti.affine, 
+                  codomain_grid_shape=codomain_nifti.shape, codomain_grid2world=codomain_nifti.affine)
+
+    return nib.Nifti1Image(affine_map.transform(codomain_nifti.get_fdata(),mode),domain_nifti.affine)
+
+def reslice_4d(domain_nifti, codomain_nifti, mode='linear'):
+    
+    reslice_arr = np.zeros(domain_nifti.shape[:3]+(codomain_nifti.shape[-1],))
+    affine_map  = AffineMap(affine=np.eye(4),
+                  domain_grid_shape=domain_nifti.shape[:3],     domain_grid2world=domain_nifti.affine, 
+                  codomain_grid_shape=codomain_nifti.shape[:3], codomain_grid2world=codomain_nifti.affine)
+    
+    for frame_id in range(reslice_arr.shape[-1]):
+        reslice_arr[:,:,:,frame_id] += affine_map.transform(codomain_nifti.get_fdata()[:,:,:,frame_id],mode)
+    
+    return nib.Nifti1Image(reslice_arr, domain_nifti.affine)
+
+
+def resample_nifti(subject_nifti, resolution=(1.5, 1.5), Nz=16, order=1, mode='nearest'):
+    """Resample a 3D or 4D (3D+time) cine-MRI nifti to a new in-plane `resolution` with `Nz` slices."""
+ 
+    data   = subject_nifti.get_fdata()
+    affine = subject_nifti.affine
+    zooms  = subject_nifti.header.get_zooms()[:3]
+    
+    new_zooms = (resolution[0], resolution[1], (zooms[2] * data.shape[2]) / Nz)
+     
+    data_resampled, affine_resampled = reslice(data, affine, zooms, new_zooms, order=order, mode=mode)
+    subject_nifti_resampled = nib.Nifti1Image(data_resampled, affine_resampled)
+     
+    x=subject_nifti_resampled.header.get_zooms()[:3]
+    y=new_zooms
+    if not np.allclose(x,y, rtol=1e-02):
+        print(subject_nifti.affine,affine)
+        print(zooms)
+        print(x,y)
+        warnings.warn('Output resolutions are different than expected!')
+         
+    return subject_nifti_resampled 
\ No newline at end of file