diff --git a/AutoscoperM/AutoscoperM.py b/AutoscoperM/AutoscoperM.py index 8807eb1..2bfdcb0 100644 --- a/AutoscoperM/AutoscoperM.py +++ b/AutoscoperM/AutoscoperM.py @@ -2,6 +2,7 @@ import glob import logging import os +import re import shutil import time import zipfile @@ -221,6 +222,7 @@ def setup(self): self.ui.segmentationButton.connect("clicked(bool)", self.onSegmentation) self.ui.loadPVButton.connect("clicked(bool)", self.onLoadPV) + self.ui.volumeSelector.connect("currentNodeChanged(vtkMRMLNode*)", self.logic.check4D) # Default output directory self.ui.mainOutputSelector.setCurrentPath( @@ -278,6 +280,10 @@ def initializeParameterNode(self): # so that when the scene is saved and reloaded, these settings are restored. self.setParameterNode(self.logic.getParameterNode()) + if self.ui.volumeSelector.currentNode() is not None: + self.logic.check4D(self.ui.volumeSelector.currentNode()) + if self.ui.mVRG_markupSelector.currentNode() is not None: + self.onMarkupNodeChanged(self.ui.mVRG_markupSelector.currentNode()) # Select default input nodes if nothing is selected yet to save a few clicks for the user # NA @@ -469,6 +475,7 @@ def onGenerateVRG(self): nPossibleCameras = self.ui.posCamSpin.value nOptimizedCameras = self.ui.optCamSpin.value tmpDir = self.ui.vrgTempDir.text + tfmPath = self.ui.tfmSubDir.text cameraSubDir = self.ui.cameraSubDir.text vrgSubDir = self.ui.vrgSubDir.text if not self.logic.validateInputs( @@ -479,6 +486,7 @@ def onGenerateVRG(self): nPossibleCameras=nPossibleCameras, nOptimizedCameras=nOptimizedCameras, tmpDir=tmpDir, + tfmPath=tfmPath, cameraSubDir=cameraSubDir, vrgSubDir=vrgSubDir, ): @@ -491,8 +499,17 @@ def onGenerateVRG(self): logging.error("Failed to generate VRG: more optimized cameras than possible cameras") return + # center the volume + self.logic.createPathsIfNotExists(os.path.join(mainOutputDir, tfmPath)) + self.logic.centerVolume(volumeNode, os.path.join(mainOutputDir, tfmPath)) + + numFrames = 1 + currentNode = volumeNode + if self.logic.is_4d: + numFrames = volumeNode.GetNumberOfDataNodes() + currentNode = self.logic.getItemInSequence(volumeNode, 0) bounds = [0] * 6 - volumeNode.GetBounds(bounds) + currentNode.GetBounds(bounds) # Generate all possible camera positions camOffset = self.ui.camOffSetSpin.value @@ -503,14 +520,18 @@ def onGenerateVRG(self): self.updateProgressBar(10) # Generate initial VRG for each camera - self.logic.generateVRGForCameras( - cameras, - volumeNode, - os.path.join(mainOutputDir, tmpDir), - width, - height, - progressCallback=self.updateProgressBar, - ) + for i in range(numFrames): + filename = self.logic.cleanFilename(currentNode.GetName(), i) + self.logic.generateVRGForCameras( + cameras, + currentNode, + os.path.join(mainOutputDir, tmpDir), + [width, height], + filename=filename, + progressCallback=self.updateProgressBar, + ) + + currentNode = self.logic.getNextItemInSequence(volumeNode) if self.logic.is_4d else currentNode # Optimize the camera positions bestCameras = RadiographGeneration.optimizeCameras( @@ -696,14 +717,36 @@ def onManualVRGGen(self): if self.logic.vrgManualCameras is None: self.onMarkupNodeChanged(markupsNode) # create the cameras - self.logic.generateVRGForCameras( - self.logic.vrgManualCameras, - volumeNode, - os.path.join(mainOutputDir, vrgDir), - width, - height, - progressCallback=self.updateProgressBar, - ) + # Check if the volume is centered at the origin + bounds = [0] * 6 + if self.logic.is_4d: + # get the bounds of the first frame + volumeNode.GetNthDataNode(0).GetRASBounds(bounds) + else: + volumeNode.GetRASBounds(bounds) + + center = [(bounds[0] + bounds[1]) / 2, (bounds[2] + bounds[3]) / 2, (bounds[4] + bounds[5]) / 2] + center = [round(x) for x in center] + if center != [0, 0, 0]: + logging.warning("Volume is not centered at the origin. This may cause issues with Autoscoper.") + + numFrames = 1 + currentNode = volumeNode + if self.logic.is_4d: + numFrames = volumeNode.GetNumberOfDataNodes() + currentNode = self.logic.getItemInSequence(volumeNode, 0) + + for i in range(numFrames): + filename = self.logic.cleanFilename(currentNode.GetName(), i) + self.logic.generateVRGForCameras( + self.logic.vrgManualCameras, + currentNode, + os.path.join(mainOutputDir, vrgDir), + [width, height], + filename=filename, + progressCallback=self.updateProgressBar, + ) + currentNode = self.logic.getNextItemInSequence(volumeNode) if self.logic.is_4d else currentNode self.updateProgressBar(100) @@ -726,6 +769,7 @@ def onMarkupNodeChanged(self, node): # get the volume nodes volumeNode = self.ui.volumeSelector.currentNode() self.logic.validateInputs(volumeNode=volumeNode) + volumeNode = self.logic.getItemInSequence(volumeNode, 0) if self.logic.is_4d else volumeNode bounds = [0] * 6 volumeNode.GetBounds(bounds) self.logic.vrgManualCameras = RadiographGeneration.generateCamerasFromMarkups( @@ -773,6 +817,7 @@ def __init__(self): self.AutoscoperProcess.setProcessChannelMode(qt.QProcess.ForwardedChannels) self.AutoscoperSocket = None self.vrgManualCameras = None + self.is_4d = False def setDefaultParameters(self, parameterNode): """ @@ -1047,9 +1092,9 @@ def generateVRGForCameras( cameras: list[RadiographGeneration.Camera], volumeNode: slicer.vtkMRMLVolumeNode, outputDir: str, - width: int, - height: int, - progressCallback: Optional[callable] = None, + size: list[int], + filename: str, + progressCallback=None, ) -> None: """ Generates VRG files for each camera in the cameras list @@ -1057,8 +1102,11 @@ def generateVRGForCameras( :param cameras: list of cameras :param volumeNode: volume node :param outputDir: output directory - :param width: width of the radiographs - :param height: height of the radiographs + :type outputDir: str + :param size: size of the VRG + :type size: list[int] + :param filename: filename of the VRG + :type filename: str :param progressCallback: progress callback, defaults to None """ self.createPathsIfNotExists(outputDir) @@ -1091,6 +1139,7 @@ def progressCallback(x): cliModule = slicer.modules.virtualradiographgeneration cliNodes = [] for cam in cameras: + logging.info(f"Generating {filename}.tif for {cam.id}") cameraDir = os.path.join(outputDir, f"cam{cam.id}") self.createPathsIfNotExists(cameraDir) camera = cam.vtkCamera @@ -1101,9 +1150,9 @@ def progressCallback(x): "cameraViewUp": [camera.GetViewUp()[0], camera.GetViewUp()[1], camera.GetViewUp()[2]], "cameraViewAngle": camera.GetViewAngle(), "clippingRange": [camera.GetClippingRange()[0], camera.GetClippingRange()[1]], - "outputWidth": width, - "outputHeight": height, - "outputFileName": os.path.join(cameraDir, "1.tif"), + "width": size[0], + "height": size[1], + "outputFileName": os.path.join(cameraDir, f"{filename}.tif"), } cliNode = slicer.cli.run(cliModule, None, parameters) # run asynchronously cliNodes.append(cliNode) @@ -1194,3 +1243,118 @@ def convertNodeToData(self, volumeNode: slicer.vtkMRMLVolumeNode) -> vtk.vtkImag imageData = imageReslice.GetOutput() return imageData + + def check4D(self, node: slicer.vtkMRMLNode) -> bool: + """ + Checks if the volume is 4D + """ + self.is_4d = type(node) == slicer.vtkMRMLSequenceNode + + def centerVolume(self, volumeNode: slicer.vtkMRMLVolumeNode, transformPath: str) -> None: + """ + A requirement for Autoscoper is that the center of the volume is at the origin. + This method will center the volume and save the transform to the transformPath + + :param volumeNode: volume node + :type volumeNode: slicer.vtkMRMLVolumeNode + :param transformPath: path to save the transform to + :type transformPath: str + :return: None + """ + + # Get the bounds of the volume + bounds = [0] * 6 + if self.is_4d: + volumeNode.GetNthDataNode(0).GetRASBounds(bounds) + else: + volumeNode.GetRASBounds(bounds) + + # Get the center of the volume + center = [0] * 3 + for i in range(3): + center[i] = (bounds[i * 2] + bounds[i * 2 + 1]) / 2 + + center_rounded = [round(x) for x in center] # don't want to move the volume if its off by a small amount + if center_rounded == [0, 0, 0]: + return # Already centered + # Create a transform node + transformNode = slicer.vtkMRMLTransformNode() + transformNode.SetName("CenteringTransform") + slicer.mrmlScene.AddNode(transformNode) + + # Get the transform matrix + matrix = vtk.vtkMatrix4x4() + + # Move the center of the volume to the origin + matrix.SetElement(0, 3, -center[0]) + matrix.SetElement(1, 3, -center[1]) + matrix.SetElement(2, 3, -center[2]) + + # Set the transform matrix + transformNode.SetMatrixTransformToParent(matrix) + + # Apply the transform to the volume + num_frames = 1 + curVol = volumeNode + if self.is_4d: + num_frames = volumeNode.GetNumberOfDataNodes() + for i in range(num_frames): + if self.is_4d: + curVol = self.getItemInSequence(volumeNode, i) + curVol.SetAndObserveTransformNodeID(transformNode.GetID()) + + # Harden the transform + slicer.modules.transforms.logic().hardenTransform(curVol) + curVol.SetAndObserveTransformNodeID(None) + + # # Invert and save the transform + matrix.Invert() + transformNode.SetMatrixTransformToParent(matrix) + slicer.util.exportNode(transformNode, os.path.join(transformPath, "Origin2DICOMCenter.tfm")) + + def getItemInSequence(self, sequenceNode: slicer.vtkMRMLSequenceNode, idx: int) -> slicer.vtkMRMLNode: + """ + Returns the item at the specified index in the sequence node + + :param sequenceNode: sequence node + :type sequenceNode: slicer.vtkMRMLSequenceNode + :param idx: index + :type idx: int + :return: item at the specified index + :rtype: slicer.vtkMRMLNode + """ + if type(sequenceNode) != slicer.vtkMRMLSequenceNode: + logging.error("[AutoscoperM.logic.getItemInSequence] sequenceNode must be a sequence node") + return None + + if idx >= sequenceNode.GetNumberOfDataNodes(): + logging.error(f"[AutoscoperM.logic.getItemInSequence] index {idx} is out of range") + return None + + browserNode = slicer.modules.sequences.logic().GetFirstBrowserNodeForSequenceNode(sequenceNode) + browserNode.SetSelectedItemNumber(idx) + return browserNode.GetProxyNode(sequenceNode) + + def getNextItemInSequence(self, sequenceNode: slicer.vtkMRMLSequenceNode) -> slicer.vtkMRMLNode: + """ + Returns the next item in the sequence + + :param sequenceNode: sequence node + :type sequenceNode: slicer.vtkMRMLSequenceNode + :return: next item in the sequence + :rtype: slicer.vtkMRMLNode + """ + if type(sequenceNode) != slicer.vtkMRMLSequenceNode: + logging.error("[AutoscoperM.logic.getNextItemInSequence] sequenceNode must be a sequence node") + return None + + browserNode = slicer.modules.sequences.logic().GetFirstBrowserNodeForSequenceNode(sequenceNode) + browserNode.SelectNextItem() + return browserNode.GetProxyNode(sequenceNode) + + def cleanFilename(self, volumeName: str, index: Optional[int] = None) -> str: + filename = ( + re.sub(r"\s+", "_", f"{volumeName}_{index}") if index is not None else re.sub(r"\s+", "_", f"{volumeName}") + ) # Remove spaces + filename = re.sub(r"[^\w]", "", filename) # Remove non alphanumeric characters + return re.sub(r"__+", "_", filename) # Remove double or more underscores diff --git a/AutoscoperM/AutoscoperMLib/RadiographGeneration.py b/AutoscoperM/AutoscoperMLib/RadiographGeneration.py index 6404110..a05ef09 100644 --- a/AutoscoperM/AutoscoperMLib/RadiographGeneration.py +++ b/AutoscoperM/AutoscoperMLib/RadiographGeneration.py @@ -221,7 +221,6 @@ def optimizeCameras( :return: Optimized cameras :rtype: list[Camera] """ - import glob import os if not progressCallback: @@ -234,11 +233,11 @@ def progressCallback(_x): cliNodes = [] for i in range(len(cameras)): camera = cameras[i] - vrgFName = glob.glob(os.path.join(cameraDir, f"cam{camera.id}", "*.tif"))[0] + vrgDirName = os.path.join(cameraDir, f"cam{camera.id}") cliNode = slicer.cli.run( cliModule, None, - {"whiteRadiographFileName": vrgFName}, + {"whiteRadiographDirName": vrgDirName}, wait_for_completion=False, ) cliNodes.append(cliNode) diff --git a/AutoscoperM/Resources/UI/AutoscoperM.ui b/AutoscoperM/Resources/UI/AutoscoperM.ui index 9cbb0ff..59c91ba 100644 --- a/AutoscoperM/Resources/UI/AutoscoperM.ui +++ b/AutoscoperM/Resources/UI/AutoscoperM.ui @@ -17,7 +17,7 @@ - 0 + 1 @@ -151,6 +151,7 @@ vtkMRMLScalarVolumeNode + vtkMRMLSequenceNode diff --git a/CalculateDataIntensityDensity/CalculateDataIntensityDensity.py b/CalculateDataIntensityDensity/CalculateDataIntensityDensity.py index 3733533..3d9668b 100644 --- a/CalculateDataIntensityDensity/CalculateDataIntensityDensity.py +++ b/CalculateDataIntensityDensity/CalculateDataIntensityDensity.py @@ -1,10 +1,11 @@ #!/usr/bin/env python-real +import glob import os import sys -def main(whiteRadiographFName: str) -> float: +def main(whiteRadiographDirName: str) -> float: """ Calculates the data intensity density of the given camera on its corresponding white radiograph. Internal function used by :func:`optimizeCameras`. @@ -17,53 +18,71 @@ def main(whiteRadiographFName: str) -> float: import numpy as np import SimpleITK as sitk - MEAN_COMPARISON = 170 # 255 / 3 * 2 - - # Read in the white radiograph - if not isinstance(whiteRadiographFName, str): - raise TypeError(f"whiteRadiographFName must be a string, not {type(whiteRadiographFName)}") - if not os.path.isfile(whiteRadiographFName): - raise FileNotFoundError(f"File {whiteRadiographFName} not found.") - whiteRadiograph = sitk.ReadImage(whiteRadiographFName) - - # Superpixel Segmentation - slicImageFilter = sitk.SLICImageFilter() - slicImageFilter.SetSuperGridSize([15, 15, 15]) # smaller grid size = finer grid overall default is [50,50,50] - labelImage = slicImageFilter.Execute(whiteRadiograph) - - # Get the mean pixel value for each label - labelStatsFilter = sitk.LabelStatisticsImageFilter() - labelStatsFilter.Execute(whiteRadiograph, labelImage) - N = labelStatsFilter.GetNumberOfLabels() - labelMeanColors = np.zeros((N, 1)) - labelWidth, labelHeight = labelImage.GetSize() - labels = list(labelStatsFilter.GetLabels()) - labels.sort() - for labelIdx, labelValue in enumerate(labels): - labelMeanColors[labelIdx, 0] = labelStatsFilter.GetMean(labelValue) - - # Create a binary label from the labelImage where all '1' are labels whose meanColor are < MEAN_COMPARISON - labelShapeFilter = sitk.LabelShapeStatisticsImageFilter() - labelShapeFilter.Execute(labelImage) - binaryLabels = np.zeros((labelWidth, labelHeight)) - for labelIdx, labelValue in enumerate(labels): - if labelValue == 0: - continue - if labelMeanColors[labelIdx, 0] < MEAN_COMPARISON: - pixels = list(labelShapeFilter.GetIndexes(labelValue)) - for j in range(0, len(pixels), 2): - y = pixels[j] - x = pixels[j + 1] - binaryLabels[x, y] = 1 - - # Calculate the Data Intensity Density - # Largest Region based off of https://discourse.itk.org/t/simpleitk-extract-largest-connected-component-from-binary-image/4958/2 - binaryImage = sitk.Cast(sitk.GetImageFromArray(binaryLabels), sitk.sitkUInt8) - componentImage = sitk.ConnectedComponent(binaryImage) - sortedComponentImage = sitk.RelabelComponent(componentImage, sortByObjectSize=True) - largest = sortedComponentImage == 1 - - return np.sum(sitk.GetArrayFromImage(largest)) + MEAN_COMPARISON = 185 + + whiteRadiographFiles = glob.glob(os.path.join(whiteRadiographDirName, "*.tif")) + + if not isinstance(whiteRadiographDirName, str): + raise TypeError(f"whiteRadiographDirName must be a string, not {type(whiteRadiographDirName)}") + if not os.path.isdir(whiteRadiographDirName): + raise FileNotFoundError(f"Directory {whiteRadiographDirName} not found.") + if len(whiteRadiographFiles) == 0: + raise FileNotFoundError(f"No white radiographs found in {whiteRadiographDirName}") + + dids = [] + for whiteRadiographFName in whiteRadiographFiles: + # Read in the white radiograph + whiteRadiograph = sitk.ReadImage(whiteRadiographFName) + + # Superpixel Segmentation + slicImageFilter = sitk.SLICImageFilter() + slicImageFilter.SetSuperGridSize([15, 15, 15]) # smaller grid size = finer grid overall default is [50,50,50] + labelImage = slicImageFilter.Execute(whiteRadiograph) + + # Get the mean pixel value for each label + labelStatsFilter = sitk.LabelStatisticsImageFilter() + labelStatsFilter.Execute(whiteRadiograph, labelImage) + N = labelStatsFilter.GetNumberOfLabels() + labelMeanColors = np.zeros((N, 1)) + labelWidth, labelHeight = labelImage.GetSize() + labels = list(labelStatsFilter.GetLabels()) + labels.sort() + for labelIdx, labelValue in enumerate(labels): + labelMeanColors[labelIdx, 0] = labelStatsFilter.GetMean(labelValue) + + # Create a binary label from the labelImage where all '1' are labels whose meanColor are < MEAN_COMPARISON + labelShapeFilter = sitk.LabelShapeStatisticsImageFilter() + labelShapeFilter.Execute(labelImage) + binaryLabels = np.zeros((labelWidth, labelHeight)) + for labelIdx, labelValue in enumerate(labels): + if labelValue == 0: + continue + if labelMeanColors[labelIdx, 0] < MEAN_COMPARISON: + pixels = list(labelShapeFilter.GetIndexes(labelValue)) + for j in range(0, len(pixels), 2): + y = pixels[j] + x = pixels[j + 1] + binaryLabels[x, y] = 1 + + # Calculate the Data Intensity Density + # Largest Region based off of https://discourse.itk.org/t/simpleitk-extract-largest-connected-component-from-binary-image/4958/2 + binaryImage = sitk.Cast(sitk.GetImageFromArray(binaryLabels), sitk.sitkUInt8) + componentImage = sitk.ConnectedComponent(binaryImage) + sortedComponentImage = sitk.RelabelComponent(componentImage, sortByObjectSize=True) + largest = sortedComponentImage == 1 + + dids.append(np.sum(sitk.GetArrayFromImage(largest))) + + # Calculate the average DID and check for any statistical outliers using interquartile range + dids = np.array(dids) + dids.sort() + q1, q3 = np.percentile(dids, [25, 75]) + iqr = q3 - q1 + lower_bound = q1 - (1.5 * iqr) + upper_bound = q3 + (1.5 * iqr) + # flag any outliers by saving the indices of the outliers + outliers = np.where((dids < lower_bound) | (dids > upper_bound)) + return np.mean(dids), outliers if __name__ == "__main__": @@ -76,4 +95,4 @@ def main(whiteRadiographFName: str) -> float: if len(sys.argv) < len(expected_args): print(f"Usage: {sys.argv[0]} {' '.join(expected_args)}") sys.exit(1) - print(main(sys.argv[1])) + print(main(sys.argv[1])[0]) diff --git a/CalculateDataIntensityDensity/CalculateDataIntensityDensity.xml b/CalculateDataIntensityDensity/CalculateDataIntensityDensity.xml index ff53988..eaaac5e 100644 --- a/CalculateDataIntensityDensity/CalculateDataIntensityDensity.xml +++ b/CalculateDataIntensityDensity/CalculateDataIntensityDensity.xml @@ -15,10 +15,10 @@ - whiteRadiographFileName - + whiteRadiographDirName + 1 - + input @@ -28,6 +28,13 @@ output + + outliers + + 3 + + output +