Cardiac Agatston scoring function for 3D Slicer

CardiacAgatstonScoring

# There are several variants of the metric. If you need to compute the metric slice by slice then
# you can use Mask volume effect to create a volume where all voxels are blanked out (set to -1000)
# except the calcifications in the selected vessel and compute the total score using this script.
# (based on Agatston (1990) - https://www.sciencedirect.com/science/article/pii/073510979090282T)
#
# Sample data set is available at:
# https://github.com/lassoan/PublicTestingData/releases/download/data/CardiacAgatstonScore.mrb

def computeAgatstonScore(volumeNode, minimumIntensityThreshold=130, minimumIslandSizeInMm2=1.0, verbose=False):
  import numpy as np
  import math
  import SimpleITK as sitk
  voxelArray = slicer.util.arrayFromVolume(volumeNode)
  areaOfPixelMm2 = volumeNode.GetSpacing()[0] * volumeNode.GetSpacing()[1]
  if verbose:
    print(f"Area of one pixel: {areaOfPixelMm2:.2f} mm2")
  sliceSpacing = volumeNode.GetSpacing()[2]
  minimumIslandSizeInPixels = int(round(minimumIslandSizeInMm2/areaOfPixelMm2))
  if verbose:
    print(f"Minimum island area: {minimumIslandSizeInPixels} voxels")
  numberOfSlices = voxelArray.shape[0]
  totalScore = 0
  for sliceIndex in range(numberOfSlices):
    voxelArraySlice = voxelArray[sliceIndex]
    maxIntensity = voxelArraySlice.max()
    if maxIntensity < minimumIntensityThreshold:
      continue
    # Get a thresholded image to analyse islands (connected components)
    # If island size less than minimum size then it will be discarded.
    thresholdedVoxelArraySlice = voxelArraySlice>minimumIntensityThreshold
    sliceImage = sitk.GetImageFromArray(voxelArraySlice)
    thresholdedSliceImage = sitk.GetImageFromArray(thresholdedVoxelArraySlice.astype(int))
    labelImage = sitk.ConnectedComponentImageFilter().Execute(thresholdedSliceImage)
    stats = sitk.LabelStatisticsImageFilter()
    stats.Execute(sliceImage, labelImage)
    sliceScore = 0
    if verbose:
      numberOfNonZeroVoxels = 0
      numberOfIslands = 0
      minWeightFactor = 1e6
      maxWeightFactor = -1
      minIslandSize = 1e6
      maxIslandSize = -1
    for labelIndex in stats.GetLabels():
      if labelIndex == 0:
        continue
      islandSizeVoxels = stats.GetCount(labelIndex)
      if islandSizeVoxels < minimumIslandSizeInPixels:
        continue
      maxIntensity = stats.GetMaximum(labelIndex)
      weightFactor = math.floor(maxIntensity/100)
      if weightFactor > 4:
        weightFactor = 4.0
      numberOfNonZeroVoxelsInIsland = islandSizeVoxels
      sliceScore += numberOfNonZeroVoxelsInIsland * areaOfPixelMm2 * weightFactor
      if verbose:
        if minWeightFactor > weightFactor:
            minWeightFactor = weightFactor
        if maxWeightFactor < weightFactor:
            maxWeightFactor = weightFactor
        if minIslandSize > islandSizeVoxels:
            minIslandSize = islandSizeVoxels
        if maxIslandSize < islandSizeVoxels:
            maxIslandSize = islandSizeVoxels
        numberOfNonZeroVoxels += numberOfNonZeroVoxelsInIsland
        numberOfIslands += 1
    totalScore += sliceScore
    if (sliceScore > 0) and verbose:
      print(f"Slice {sliceIndex} score: {sliceScore:.1f} (from {numberOfIslands} islands of size {minIslandSize} to {maxIslandSize},"
        f" in total {numberOfNonZeroVoxels} voxels, weight factor: {minWeightFactor} to {maxWeightFactor})")
      slicer.app.processEvents()
  normalizedTotalScore = totalScore * sliceSpacing / 3.0
  if verbose:
    print(f"Total Score without slice space scaling: {totalScore}. After normalizing from {sliceSpacing} to 3.0mm spacing: {normalizedTotalScore}.")
  return normalizedTotalScore

print("Total Agatston score: {0}".format(computeAgatstonScore(getNode('Panoramix-cropped masked'), verbose=True)))