Switch to Texas Tech University unrectified scans (#11)

Previously, the algorithm was designed to run on rectified scans from Texas Tech University. The rectification process consisted of a MATLAB script written by Jason Hill that would rotate the upper and lower volumes independently to better align. The script calculated the centroid and identified the spinal cord axis for better alignment.

This process proved to be extremely complex and unnecessary. In addition, significant portions of the upper and lower scans were cut off resulting in an entire vertebrae typically missing. This is non-ideal since Amy Givan's thesis is on quantifying the amount of VAT & SCAT on a vertebra-by-vertebra basis. Additionally, the typical verebrae cutoff was L3 or L2 and those are the locations, according to studies, where a significant measure of SCAT/VAT loss is recorded.

A script was written to move the unrectified TTU scans to where the rectified ones are. Thus, the actual loading of the scans is unchanged.

The config file was switched from XML to YAML for cleaner code. XML is quite confusing to write in Python. Updated configure window and segmentation algorithm to utilize the YAML config file instead.

Added automatic upper & lower scan alignment via the origin and spacing metadata in the NiFTi files. The images are aligned as close as possible by shifting the images around.

core/loadData.py

@@ -52,39 +52,69 @@ def _loadTexasTechDixonData(dataPath):
     niiFatLowerFilename = os.path.join(dataPath, 'fatLower.nii')
     niiWaterUpperFilename = os.path.join(dataPath, 'waterUpper.nii')
     niiWaterLowerFilename = os.path.join(dataPath, 'waterLower.nii')
-    configFilename = os.path.join(dataPath, 'config.xml')
+    configFilename = os.path.join(dataPath, 'config.yml')
 
-    if not (os.path.isfile(niiFatUpperFilename) and os.path.isfile(niiFatLowerFilename) and \
-            os.path.isfile(niiWaterUpperFilename) and os.path.isfile(niiWaterLowerFilename) and \
-            os.path.isfile(configFilename)):
+    if not (os.path.isfile(niiFatUpperFilename) and os.path.isfile(niiFatLowerFilename) and
+            os.path.isfile(niiWaterUpperFilename) and os.path.isfile(niiWaterLowerFilename)):
         raise Exception('Missing required files from source path folder.')
 
+    # Load the configuration file if it exists
+    if os.path.exists(configFilename):
+        with open(configFilename, 'r') as fh:
+            config = yaml.load(fh)
+    else:
+        # Otherwise create the config as an empty dictionary
+        config = {}
+
     # Load unrectified NIFTI files for the current dataPath
     niiFatUpper = nib.load(niiFatUpperFilename)
     niiFatLower = nib.load(niiFatLowerFilename)
     niiWaterUpper = nib.load(niiWaterUpperFilename)
     niiWaterLower = nib.load(niiWaterLowerFilename)
 
-    # Load config XML file
-    config = etree.parse(configFilename)
+    # Affine matrix, origin & spacing for upper image. Fat & water should have same info
+    upperAffineMatrixWT = niiFatUpper.header.get_best_affine()
+    upperAffineMatrix = np.array(upperAffineMatrixWT[:-1, :-1])
+    upperOrigin = np.array(upperAffineMatrixWT[:-1, -1])
+    spacing = np.linalg.norm(upperAffineMatrix, axis=1)
+
+    # -1 for flipped axes and 1 for non-flipped axes (RAS)
+    axesFlipped = np.where(np.diag(upperAffineMatrix) > 0, 1, -1)
 
-    # Get the root of the config XML file
-    configRoot = config.getroot()
+    # Odd hacks here to get the right amount of offset, this only applies to the listed subjects
+    # No idea why
+    # Toggle the axes flipped from 1 to -1 or -1 to 1. This way it will shift the opposite direction
+    if os.path.basename(dataPath) in ['Subject0003_Final', 'Subject0003_Initial']:
+        axesFlipped[0] *= -1
 
-    # Piece together upper and lower images for fat and water
-    # Retrieve the inferior and superior axial slice from config file for upper and lower images
-    imageUpperTag = configRoot.find('imageUpper')
-    imageLowerTag = configRoot.find('imageLower')
-    imageUpperInferiorSlice = int(imageUpperTag.attrib['inferiorSlice'])
-    imageUpperSuperiorSlice = int(imageUpperTag.attrib['superiorSlice'])
-    imageLowerInferiorSlice = int(imageLowerTag.attrib['inferiorSlice'])
-    imageLowerSuperiorSlice = int(imageLowerTag.attrib['superiorSlice'])
+    # Affine matrix, origin & spacing for lower image. Fat & water should have same info
+    lowerAffineMatrixWT = niiFatLower.header.get_best_affine()
+    lowerAffineMatrix = np.array(lowerAffineMatrixWT[:-1, :-1])
+    lowerOrigin = np.array(lowerAffineMatrixWT[:-1, -1])
 
     # Use inferior and superior axial slice to obtain the valid portion of the upper and lower fat and water images
-    fatUpperImage = niiFatUpper.get_data()[:, :, imageUpperInferiorSlice:imageUpperSuperiorSlice]
-    fatLowerImage = niiFatLower.get_data()[:, :, imageLowerInferiorSlice:imageLowerSuperiorSlice]
-    waterUpperImage = niiWaterUpper.get_data()[:, :, imageUpperInferiorSlice:imageUpperSuperiorSlice]
-    waterLowerImage = niiWaterLower.get_data()[:, :, imageLowerInferiorSlice:imageLowerSuperiorSlice]
+    fatUpperImage, fatLowerImage = niiFatUpper.get_data(), niiFatLower.get_data()
+    waterUpperImage, waterLowerImage = niiWaterUpper.get_data(), niiWaterLower.get_data()
+
+    # Take lower origin and subtract from the upper origin and divide by spacing
+    # For the Z dimension we want to add the fat lower shape to get the amount of slices that overlap
+    # Then convert to integer after rounding down to get the number of indices to shift to align
+    misalignedIndexAmount = np.floor((lowerOrigin - upperOrigin) / (spacing * axesFlipped) +
+                                     (0, 0, fatLowerImage.shape[2])).astype(int)
+
+    if misalignedIndexAmount[2] > 0:
+        fatLowerImage = fatLowerImage[:, :, :-misalignedIndexAmount[2]]
+        waterLowerImage = waterLowerImage[:, :, :-misalignedIndexAmount[2]]
+
+    if misalignedIndexAmount[1] != 0:
+        # Roll the fat upper image back the difference amount to align better to lower image
+        fatLowerImage = np.roll(fatLowerImage, misalignedIndexAmount[1], axis=1)
+        waterLowerImage = np.roll(waterLowerImage, misalignedIndexAmount[1], axis=1)
+
+    if misalignedIndexAmount[0] != 0:
+        # Roll the fat upper image back the difference amount to align better to lower image
+        fatLowerImage = np.roll(fatLowerImage, misalignedIndexAmount[0], axis=0)
+        waterLowerImage = np.roll(waterLowerImage, misalignedIndexAmount[0], axis=0)
 
     # Piece together the upper and lower parts of the fat and water images into one volume
     fatImage = np.concatenate((fatLowerImage.T, fatUpperImage.T), axis=0)
@@ -98,14 +128,11 @@ def _loadTexasTechDixonData(dataPath):
     constants.pathDir = dataPath
 
     # Create a NRRD header dictionary that will be used to save the intermediate debug NRRDs to view progress
-    constants.nrrdHeaderDict = {'space': 'right-anterior-superior'}
-    constants.nrrdHeaderDict['space directions'] = np.hstack((niiFatUpper.header['srow_x'][0:-1, None],
-                                                              niiFatUpper.header['srow_y'][0:-1, None],
-                                                              niiFatUpper.header['srow_z'][0:-1, None]))
-
-    constants.nrrdHeaderDict['space origin'] = (niiFatUpper.header['srow_x'][-1],
-                                                niiFatUpper.header['srow_y'][-1],
-                                                niiFatUpper.header['srow_z'][-1])
+    constants.nrrdHeaderDict = {
+        'space': 'right-anterior-superior',
+        'space directions': upperAffineMatrix,
+        'space origin': lowerOrigin
+    }
 
     return fatImage, waterImage, config
 

core/runSegmentation_TexasTechDixon.py

@@ -201,31 +201,19 @@ def runSegmentation(data):
     # The bias corrected fat and water images are going to be created in this directory regardless of debug constant
     os.makedirs(getDebugPath(''), exist_ok=True)
 
-    # Get the root of the config XML file
-    configRoot = config.getroot()
-
-    diaphragmSuperiorSlice = int(configRoot.find('diaphragm').attrib['superiorSlice'])
-    umbilicisTag = configRoot.find('umbilicis')
-    umbilicisInferiorSlice = int(umbilicisTag.attrib['inferiorSlice'])
-    umbilicisSuperiorSlice = int(umbilicisTag.attrib['superiorSlice'])
-    umbilicisLeft = int(umbilicisTag.attrib['left'])
-    umbilicisRight = int(umbilicisTag.attrib['right'])
-    umbilicisCoronal = int(umbilicisTag.attrib['coronal'])
-
-    # Load cardiac adipose tissue (CAT) tag and corresponding lines in axial plane
-    CATTag = configRoot.find('CAT')
-    CATAxial = []
-    CATPosterior = []
-    CATAnterior = []
-    # Foreach line in the CAT tag, append to the three arrays
-    for line in CATTag:
-        if line.tag != 'line':
-            print('Invalid tag for CAT, must be line')
-            continue
-
-        CATAxial.append(int(line.attrib['axial']))
-        CATPosterior.append(int(line.attrib['posterior']))
-        CATAnterior.append(int(line.attrib['anterior']))
+    # Load values from config dictionary
+    diaphragmAxial = config['diaphragmAxial']
+    umbilicis = config['umbilicis']
+    umbilicisInferior = umbilicis['inferior']
+    umbilicisSuperior = umbilicis['superior']
+    umbilicisLeft = umbilicis['left']
+    umbilicisRight = umbilicis['right']
+    umbilicisCoronal = umbilicis['coronal']
+
+    CATBounds = list([(x['axial'], x['posterior'], x['anterior']) for x in config['CATBounds']])
+    CATAxial = list([x[0] for x in CATBounds])
+    CATPosterior = list([x[1] for x in CATBounds])
+    CATAnterior = list([x[2] for x in CATBounds])
 
     # Convert three arrays to NumPy and get minimum/maximum axial slice
     # The min/max axial slice is used to determine the start and stopping point
@@ -245,9 +233,10 @@ def runSegmentation(data):
     # Perform bias correction on MRI images to remove inhomogeneity
     # If bias correction has been performed already, then load the saved data
     tic = time.perf_counter()
-    if os.path.exists(getDebugPath('fatImageBC.nrrd')) and os.path.exists(getDebugPath('waterImageBC.nrrd')):
-        fatImage, header = nrrd.read(getDebugPath('fatImageBC.nrrd'))
-        waterImage, header = nrrd.read(getDebugPath('waterImageBC.nrrd'))
+    if not constants.forceBiasCorrection and os.path.exists(getPath('fatImage.nrrd')) and os.path.exists(
+            getPath('waterImage.nrrd')):
+        fatImage, header = nrrd.read(getPath('fatImage.nrrd'))
+        waterImage, header = nrrd.read(getPath('waterImage.nrrd'))
 
         # Transpose image to get back into C-order indexing
         fatImage, waterImage = fatImage.T, waterImage.T
@@ -256,8 +245,8 @@ def runSegmentation(data):
         waterImage = correctBias(waterImage, shrinkFactor=constants.shrinkFactor, prefix='waterImageBiasCorrection')
 
         # If bias correction is performed, saved images to speed up algorithm in future runs
-        nrrd.write(getDebugPath('fatImageBC.nrrd'), fatImage.T, constants.nrrdHeaderDict)
-        nrrd.write(getDebugPath('waterImageBC.nrrd'), waterImage.T, constants.nrrdHeaderDict)
+        nrrd.write(getPath('fatImage.nrrd'), fatImage.T, constants.nrrdHeaderDict)
+        nrrd.write(getPath('waterImage.nrrd'), waterImage.T, constants.nrrdHeaderDict)
 
     toc = time.perf_counter()
     print('N4ITK bias field correction took %f seconds' % (toc - tic))
@@ -279,7 +268,7 @@ def runSegmentation(data):
     CAT = np.zeros(fatImage.shape, bool)
 
     for slice in range(0, fatImage.shape[0]):
-    # for slice in range(diaphragmSuperiorSlice, fatImage.shape[0]):
+        # for slice in range(diaphragmSuperiorSlice, fatImage.shape[0]):
         tic = time.perf_counter()
 
         fatImageSlice = fatImage[slice, :, :]
@@ -296,7 +285,7 @@ def runSegmentation(data):
         # Algorithm assumes that the skin is a closed contour and fully connects
         # This is a valid assumption but near the umbilicis, there is a discontinuity
         # so this draws a line near there to create a closed contour
-        if umbilicisInferiorSlice <= slice <= umbilicisSuperiorSlice:
+        if umbilicisInferior <= slice <= umbilicisSuperior:
             fatImageMask[umbilicisCoronal, umbilicisLeft:umbilicisRight] = True
 
         # Save fat and water masks for debugging
@@ -313,7 +302,7 @@ def runSegmentation(data):
         bodyMasks[slice, :, :] = bodyMask
 
         # Superior of diaphragm is divider between thoracic and abdominal region
-        if slice < diaphragmSuperiorSlice:
+        if slice < diaphragmAxial:
             fatVoidMask, abdominalMask, SCATSlice, VATSlice = \
                 segmentAbdomenSlice(slice, fatImageMask, waterImageMask, bodyMask)
 

core/runSegmentation_WashUDixon.py

@@ -115,7 +115,8 @@ def runSegmentation(data):
     # Perform bias correction on MRI images to remove inhomogeneity
     # If bias correction has been performed already, then load the saved data
     tic = time.perf_counter()
-    if os.path.exists(getDebugPath('fatImageBC.nrrd')) and os.path.exists(getDebugPath('waterImageBC.nrrd')):
+    if not constants.forceBiasCorrection and os.path.exists(getDebugPath('fatImageBC.nrrd')) and os.path.exists(
+            getDebugPath('waterImageBC.nrrd')):
         fatImage, header = nrrd.read(getDebugPath('fatImageBC.nrrd'))
         waterImage, header = nrrd.read(getDebugPath('waterImageBC.nrrd'))
 
@@ -255,7 +256,8 @@ def runSegmentation(data):
     if constants.debug:
         nrrd.write(getDebugPath('fatImageMask.nrrd'), skimage.img_as_ubyte(fatImageMasks).T, constants.nrrdHeaderDict,
                    compression_level=1)
-        nrrd.write(getDebugPath('waterImageMask.nrrd'), skimage.img_as_ubyte(waterImageMasks).T, constants.nrrdHeaderDict,
+        nrrd.write(getDebugPath('waterImageMask.nrrd'), skimage.img_as_ubyte(waterImageMasks).T,
+                   constants.nrrdHeaderDict,
                    compression_level=1)
         nrrd.write(getDebugPath('bodyMask.nrrd'), skimage.img_as_ubyte(bodyMasks).T, constants.nrrdHeaderDict,
                    compression_level=1)
@@ -266,7 +268,8 @@ def runSegmentation(data):
                    compression_level=1)
 
     # Save the results of adipose tissue segmentation and the original fat/water images
-    nrrd.write(getPath('fatImage.nrrd'), skimage.img_as_ubyte(fatImage).T, constants.nrrdHeaderDict, compression_level=1)
+    nrrd.write(getPath('fatImage.nrrd'), skimage.img_as_ubyte(fatImage).T, constants.nrrdHeaderDict,
+               compression_level=1)
     nrrd.write(getPath('waterImage.nrrd'), skimage.img_as_ubyte(waterImage).T, constants.nrrdHeaderDict,
                compression_level=1)
     nrrd.write(getPath('SCAT.nrrd'), skimage.img_as_ubyte(SCAT).T, constants.nrrdHeaderDict, compression_level=1)

core/runSegmentation_WashUUnknown.py

@@ -115,7 +115,7 @@ def runSegmentation(data):
     # Perform bias correction on MRI images to remove inhomogeneity
     # If bias correction has been performed already, then load the saved data
     tic = time.perf_counter()
-    if os.path.exists(getDebugPath('imageBC.nrrd')):
+    if not constants.forceBiasCorrection and os.path.exists(getDebugPath('imageBC.nrrd')):
         image, header = nrrd.read(getDebugPath('imageBC.nrrd'))
 
         # Transpose image to get back into C-order indexing