Merge branch 'master' of [email protected]:gschramm/pymirc.git

examples/fileio/view_nrrd_seg.py

@@ -0,0 +1,36 @@
+# small example on how to read nrrd files and segmentation as e.g.
+# produced by 3D slicer
+
+# to run it you need to install pynrrd
+
+import nrrd
+import numpy as np
+import pymirc.image_operations as pi
+import pymirc.viewer as pv
+
+def get_lps_affine_from_hdr(hdr):
+  aff = np.eye(4)
+  aff[:3,:3] = hdr['space directions']
+  aff[:3,-1] = hdr['space origin']
+
+  # convert the affine to LPS
+  if hdr['space'] == 'right-anterior-superior':
+    aff = np.diag([-1,-1,1,1]) @ aff
+  elif hdr['space'] == 'left-anterior-superior':
+    aff = np.diag([1,-1,1,1]) @ aff
+  elif hdr['space'] == 'right-posterior-superior':
+    aff = np.diag([-1,1,1,1]) @ aff
+
+  return aff
+
+ct, ct_hdr   = nrrd.read('CTChest.nrrd')
+seg, seg_hdr = nrrd.read('Segmentation preview.seg.nrrd')  
+
+ct_aff  = get_lps_affine_from_hdr(ct_hdr)
+seg_aff = get_lps_affine_from_hdr(seg_hdr) 
+
+ct_voxsize = np.sqrt((ct_aff**2).sum(0))[:-1]
+
+seg2 = pi.aff_transform(seg, np.linalg.inv(seg_aff) @ ct_aff, ct.shape, trilin = False)
+
+pv.ThreeAxisViewer([ct,seg2], voxsize=ct_voxsize)

examples/pet/nema_pet.py

@@ -1117,6 +1117,7 @@ def fit_WB_NEMA_sphere_profiles(vol,
                                 showprofiles = True,
                                 showrcs      = True,
                                 wm           = 'dist',
+                                nmax_spheres = 6,
                                 sameSignal   = False):
     """ Analyse the sphere profiles of a NEMA scan
 
@@ -1152,6 +1153,9 @@ def fit_WB_NEMA_sphere_profiles(vol,
     wm : str, optional
       the weighting method of the data (equal, dist, sqdist)
 
+    nmax_spheres:
+      maximum number of spheres to consider (default 6)
+
     Returns
     -------
     Dictionary
@@ -1169,7 +1173,9 @@ def fit_WB_NEMA_sphere_profiles(vol,
     slices  = NEMASubvols(vol, voxsizes, margin = margin)
 
     subvols = list()
-    for ss in slices: subvols.append(vol[ss])
+    for iss, ss in enumerate(slices): 
+      if iss < nmax_spheres:
+        subvols.append(vol[ss])
 
     if Rfix == None: Rfix = [None] * len(subvols)
     if len(Rfix) < len(subvols): Rfix = Rfix + [None] * (len(subvols) - len(Rfix))

pymirc/image_operations/__init__.py

@@ -8,3 +8,4 @@
 from .resample_img_cont           import resample_img_cont
 from .resample_cont_cont          import resample_cont_cont
 from .grad                        import grad, div, complex_grad, complex_div
+from .reorient                    import reorient_image_and_affine, flip_image_and_affine

pymirc/image_operations/reorient.py

@@ -0,0 +1,78 @@
+import nibabel as nib
+from nibabel.orientations import inv_ornt_aff
+
+def reorient_image_and_affine(img, aff):
+  """ Reorient an image and and affine such that the affine is approx. diagonal
+      and such that the elements on the main diagonal are positiv
+  
+  Parameters
+  ----------
+  img : 3D numpy array
+    containing the image
+
+  aff : 2D numpy array
+    (4,4) affine transformation matrix from image to anatomical coordinate system in
+    homogeneous coordinates
+
+  Returns
+  -------
+  a tuple with the reoriented image and the accordingly transformed affine
+
+  Note
+  ----
+  The reorientation uses nibabel's io_orientation() and apply_orientation()
+  """
+  ornt  = nib.io_orientation(aff)
+  img_t = nib.apply_orientation(img, ornt)
+  aff_t = aff.dot(inv_ornt_aff(ornt, img.shape))
+
+  return img_t, aff_t
+
+def flip_image_and_affine(img, aff, flip_dirs):
+  """ Flip an image along given axis and transform the affine accordingly
+
+  Parameters
+  ----------
+  img : numpy array (2D, 3D)
+    containing the image
+
+  aff : 2D numpy array
+    (4,4) affine transformation matrix from image to anatomical coordinate system in
+    homogeneous coordinates
+
+  flip_dirs : int or tuple of ints
+    containing the axis where the flip should  be applied
+
+  Returns
+  -------
+  a tuple with the flipped image and the accordingly transformed affine
+  """
+  if not isinstance(flip_dirs, tuple):
+    flip_dirs = (flip_dirs,)
+
+  # flip the image
+  img = np.flip(img, flip_dirs)
+
+  for flip_dir in flip_dirs:
+    # tmp is the diagonal of the flip affine (all ones, but -1 in the flip direction)
+    tmp            = np.ones(img.ndim + 1)
+    tmp[flip_dir]  = -1
+
+    # tmp2 is the (i,j,k,1) vector for the 0,0,0 voxel
+    tmp2           = np.zeros(img.ndim + 1)
+    tmp2[-1]       = 1
+
+    # tmp2 is the (i,j,k,1) vector for the last voxel along the flip direction
+    tmp3           = tmp2.copy()
+    tmp3[flip_dir] = dim[flip_dir] - 1
+
+    # this is the affine that does the flipping 
+    # the flipping is actual a shift to the center, followed by an inversion
+    # followed by in the inverse shift
+    flip_aff = np.diag(tmp)
+    flip_aff[flip_dir,-1] = ((aff @ tmp2) + (aff @ tmp3))[flip_dir]
+
+    # transform the affine
+    aff = flip_aff @ aff
+
+  return img, aff