diff --git a/AutoscoperM/AutoscoperM.py b/AutoscoperM/AutoscoperM.py index 71e63ff..4a82b84 100644 --- a/AutoscoperM/AutoscoperM.py +++ b/AutoscoperM/AutoscoperM.py @@ -18,7 +18,7 @@ ) from slicer.util import VTKObservationMixin -from AutoscoperMLib import IO, RadiographGeneration, SubVolumeExtraction +from AutoscoperMLib import IO, SubVolumeExtraction # # AutoscoperM @@ -36,10 +36,6 @@ def __init__(self, parent): self.parent.categories = [ "Tracking", ] - self.parent.dependencies = [ - "CalculateDataIntensityDensity", - "VirtualRadiographGeneration", - ] self.parent.contributors = [ "Anthony Lombardi (Kitware)", "Amy M Morton (Brown University)", @@ -216,8 +212,6 @@ def setup(self): # Pre-processing Library Buttons self.ui.tiffGenButton.connect("clicked(bool)", self.onGeneratePartialVolumes) - self.ui.vrgGenButton.connect("clicked(bool)", self.onGenerateVRG) - self.ui.manualVRGGenButton.connect("clicked(bool)", self.onManualVRGGen) self.ui.configGenButton.connect("clicked(bool)", self.onGenerateConfig) self.ui.segmentationButton.connect("clicked(bool)", self.onSegmentation) @@ -228,11 +222,6 @@ def setup(self): os.path.join(slicer.mrmlScene.GetCacheManager().GetRemoteCacheDirectory(), "AutoscoperM-Pre-Processing") ) - # Dynamic camera frustum functions - self.ui.mVRG_markupSelector.connect("currentNodeChanged(vtkMRMLNode*)", self.onMarkupNodeChanged) - self.ui.mVRG_ClippingRangeSlider.connect("valuesChanged(double,double)", self.updateClippingRange) - self.ui.mVRG_viewAngleSpin.connect("valueChanged(int)", self.updateViewAngle) - # Make sure parameter node is initialized (needed for module reload) self.initializeParameterNode() @@ -279,8 +268,6 @@ def initializeParameterNode(self): # so that when the scene is saved and reloaded, these settings are restored. self.setParameterNode(self.logic.getParameterNode()) - 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 @@ -470,151 +457,6 @@ def onGeneratePartialVolumes(self): self.onLoadPV() # onLoadPV has a call to the "success" display, remove the one here so the user doesn't get two. - def onGenerateVRG(self): - """ - This function optimizes the camera positions for a given volume and then - generates a VRG file for each optimized camera. - """ - - with slicer.util.tryWithErrorDisplay("Failed to compute results", waitCursor=True): - self.updateProgressBar(0) - - # Set up and validate inputs - volumeNode = self.ui.volumeSelector.currentNode() - mainOutputDir = self.ui.mainOutputSelector.currentPath - width = self.ui.vrgRes_width.value - height = self.ui.vrgRes_height.value - 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( - volumeNode=volumeNode, - mainOutputDir=mainOutputDir, - width=width, - height=height, - nPossibleCameras=nPossibleCameras, - nOptimizedCameras=nOptimizedCameras, - tmpDir=tmpDir, - tfmPath=tfmPath, - cameraSubDir=cameraSubDir, - vrgSubDir=vrgSubDir, - ): - raise ValueError("Invalid inputs") - return - if not self.logic.validatePaths(mainOutputDir=mainOutputDir): - raise ValueError("Invalid paths") - return - if nPossibleCameras < nOptimizedCameras: - raise Exception("Failed to generate VRG: more optimized cameras than possible cameras") - return - - # center the volume - self.logic.createPathsIfNotExists(os.path.join(mainOutputDir, tfmPath)) - if not self.logic.IsSequenceVolume(volumeNode): - volumeNode.AddCenteringTransform() - tfmNode = slicer.util.getNode(f"{volumeNode.GetName()} centering transform") - volumeNode.HardenTransform() - volumeNode.SetAndObserveTransformNodeID(None) - tfmPath = os.path.join(mainOutputDir, tfmPath, "Origin2Dicom.tfm") - tfmNode.Inverse() - slicer.util.saveNode(tfmNode, tfmPath) - slicer.mrmlScene.RemoveNode(tfmNode) - else: - # Just export the first frame - currentNode, _ = self.logic.getItemInSequence(volumeNode, 0) - currentNode.AddCenteringTransform() - tfmNode = currentNode.GetParentTransformNode() - currentNode.HardenTransform() - currentNode.SetAndObserveTransformNodeID(None) - tfmPath = os.path.join(mainOutputDir, tfmPath, "Origin2Dicom.tfm") - tfmNode.Inverse() - slicer.util.saveNode(tfmNode, tfmPath) - slicer.mrmlScene.RemoveNode(tfmNode) - - # Harden and remove the transform from the sequence - for i in range(1, volumeNode.GetNumberOfDataNodes()): - currentNode, _ = self.logic.getItemInSequence(volumeNode, i) - currentNode.AddCenteringTransform() - tfmNode = currentNode.GetParentTransformNode() - currentNode.HardenTransform() - currentNode.SetAndObserveTransformNodeID(None) - slicer.mrmlScene.RemoveNode(tfmNode) - - numFrames = 1 - currentNode = volumeNode - curName = volumeNode.GetName() - if self.logic.IsSequenceVolume(currentNode): - numFrames = volumeNode.GetNumberOfDataNodes() - currentNode, curName = self.logic.getItemInSequence(volumeNode, 0) - bounds = [0] * 6 - currentNode.GetBounds(bounds) - - # Generate all possible camera positions - camOffset = self.ui.camOffSetSpin.value - cameras = RadiographGeneration.generateNCameras( - nPossibleCameras, bounds, camOffset, [width, height], self.ui.camDebugCheckbox.isChecked() - ) - - # TODO: Validate that both the tmp directory and the camera calibration directory are empty before starting - - cameraDir = os.path.join(mainOutputDir, cameraSubDir) - if not os.path.exists(cameraDir): - os.mkdir(cameraDir) - for cam in cameras: # Generate Camera Calib files - camFName = os.path.join(cameraDir, f"cam{cam.id}.json") - IO.generateCameraCalibrationFile(cam, camFName) - - self.updateProgressBar(10) - - # Generate initial VRG for each camera - for i in range(numFrames): - filename = self.logic.cleanFilename(curName, i) if not self.ui.idxOnly.isChecked() else i - self.logic.generateVRGForCameras( - os.path.join(mainOutputDir, cameraSubDir), - currentNode, - os.path.join(mainOutputDir, tmpDir), - [width, height], - filename=filename, - ) - progress = (i + 1) / numFrames * 40 + 10 - self.updateProgressBar(progress) - - if self.logic.IsSequenceVolume(volumeNode): - currentNode, curName = self.logic.getNextItemInSequence(volumeNode) - - # Optimize the camera positions - bestCameras = RadiographGeneration.optimizeCameras( - cameras, os.path.join(mainOutputDir, tmpDir), nOptimizedCameras, progressCallback=self.updateProgressBar - ) - - shutil.rmtree(os.path.join(mainOutputDir, cameraSubDir)) - - # Move the optimized VRGs to the final directory and generate the camera calibration files - self.logic.generateCameraCalibrationFiles( - bestCameras, - os.path.join(mainOutputDir, tmpDir), - os.path.join(mainOutputDir, vrgSubDir), - os.path.join(mainOutputDir, cameraSubDir), - progressCallback=self.updateProgressBar, - ) - - # Clean Up - if self.ui.removeVrgTmp.isChecked(): - shutil.rmtree(os.path.join(mainOutputDir, tmpDir)) - - slicer.util.messageBox("Success!") - - if not self.logic.IsSequenceVolume(volumeNode): - firstNode = volumeNode - else: - firstNode, _ = self.logic.getItemInSequence(volumeNode, 0) - firstNode.SetAndObserveTransformNodeID(tfmNode.GetID()) - firstNode.HardenTransform() - slicer.mrmlScene.RemoveNode(tfmNode) - def onGenerateConfig(self): """ Generates a complete config file (including all partial volumes, radiographs, @@ -792,113 +634,6 @@ def onLoadPV(self): slicer.util.messageBox("Success!") - def onManualVRGGen(self): - with slicer.util.tryWithErrorDisplay("Failed to compute results.", waitCursor=True): - markupsNode = self.ui.mVRG_markupSelector.currentNode() - volumeNode = self.ui.volumeSelector.currentNode() - mainOutputDir = self.ui.mainOutputSelector.currentPath - viewAngle = self.ui.mVRG_viewAngleSpin.value - clippingRange = ( - self.ui.mVRG_ClippingRangeSlider.minimumValue, - self.ui.mVRG_ClippingRangeSlider.maximumValue, - ) - width = self.ui.vrgRes_width.value - height = self.ui.vrgRes_height.value - vrgDir = self.ui.vrgSubDir.text - cameraDir = self.ui.cameraSubDir.text - if not self.logic.validateInputs( - markupsNode=markupsNode, - volumeNode=volumeNode, - mainOutputDir=mainOutputDir, - viewAngle=viewAngle, - clippingRange=clippingRange, - width=width, - height=height, - vrgDir=vrgDir, - cameraDir=cameraDir, - ): - raise Exception("Failed to generate VRG: invalid inputs") - return - if not self.logic.validatePaths(mainOutputDir=mainOutputDir): - raise Exception("Failed to generate VRG: invalid output directory") - return - self.logic.createPathsIfNotExists( - os.path.join(mainOutputDir, vrgDir), os.path.join(mainOutputDir, cameraDir) - ) - - if self.logic.vrgManualCameras is None: - self.onMarkupNodeChanged(markupsNode) # create the cameras - - if not self.logic.IsVolumeCentered(volumeNode): - logging.warning("Volume is not centered at the origin. This may cause issues with Autoscoper.") - - for cam in self.logic.vrgManualCameras: - IO.generateCameraCalibrationFile(cam, os.path.join(mainOutputDir, cameraDir, f"cam{cam.id}.json")) - - self.updateProgressBar(0) - - numFrames = 1 - currentNode = volumeNode - curName = currentNode.GetName() - if self.logic.IsSequenceVolume(currentNode): - numFrames = volumeNode.GetNumberOfDataNodes() - currentNode, curName = self.logic.getItemInSequence(volumeNode, 0) - - for i in range(numFrames): - filename = self.logic.cleanFilename(curName, i) - self.logic.generateVRGForCameras( - os.path.join(mainOutputDir, cameraDir), - currentNode, - os.path.join(mainOutputDir, vrgDir), - [width, height], - filename=filename, - ) - - progress = ((i + 1) / numFrames) * 100 - self.updateProgressBar(progress) - - if self.logic.IsSequenceVolume(volumeNode): - currentNode, curName = self.logic.getNextItemInSequence(volumeNode) - - self.updateProgressBar(100) - slicer.util.messageBox("Success!") - - def onMarkupNodeChanged(self, node): - if node is None: - if self.logic.vrgManualCameras is not None: - # clean up - for cam in self.logic.vrgManualCameras: - slicer.mrmlScene.RemoveNode(cam.FrustumModel) - self.logic.vrgManualCameras = None - return - if self.logic.vrgManualCameras is not None: - # clean up - for cam in self.logic.vrgManualCameras: - slicer.mrmlScene.RemoveNode(cam.FrustumModel) - self.logic.vrgManualCameras = None - # get the volume nodes - volumeNode = self.ui.volumeSelector.currentNode() - self.logic.validateInputs(volumeNode=volumeNode) - bounds = self.logic.GetRASBounds(volumeNode) - self.logic.vrgManualCameras = RadiographGeneration.generateCamerasFromMarkups( - node, - bounds, - (self.ui.mVRG_ClippingRangeSlider.minimumValue, self.ui.mVRG_ClippingRangeSlider.maximumValue), - self.ui.mVRG_viewAngleSpin.value, - [self.ui.vrgRes_width.value, self.ui.vrgRes_height.value], - True, - ) - - def updateClippingRange(self, min, max): - for cam in self.logic.vrgManualCameras: - cam.vtkCamera.SetClippingRange(min, max) - RadiographGeneration._updateFrustumModel(cam) - - def updateViewAngle(self, value): - for cam in self.logic.vrgManualCameras: - cam.vtkCamera.SetViewAngle(value) - RadiographGeneration._updateFrustumModel(cam) - # # AutoscoperMLogic @@ -924,7 +659,6 @@ def __init__(self): self.AutoscoperProcess = qt.QProcess() self.AutoscoperProcess.setProcessChannelMode(qt.QProcess.ForwardedChannels) self.AutoscoperSocket = None - self.vrgManualCameras = None @staticmethod def IsSequenceVolume(node: Union[slicer.vtkMRMLNode, None]) -> bool: @@ -1300,136 +1034,6 @@ def extractSubVolumeForVRG( return newVolumeImageData, bounds - def generateVRGForCameras( - self, - cameraDir: str, - volumeNode: slicer.vtkMRMLVolumeNode, - outputDir: str, - size: list[int], - filename: str, - ) -> None: - """ - Generates VRG files for each camera in the cameras list - - :param cameraDir: Directory containing the camera JSON files - :param volumeNode: volume node - :param outputDir: output directory - :param size: size of the VRG - :param filename: filename of the VRG - """ - self.createPathsIfNotExists(outputDir) - # Apply a thresh of 0 to the volume to remove air from the volume - thresholdScalarVolume = slicer.modules.thresholdscalarvolume - parameters = { - "InputVolume": volumeNode.GetID(), - "OutputVolume": volumeNode.GetID(), - "ThresholdValue": 0, - "ThresholdType": "Below", - "Lower": 0, - } - slicer.cli.runSync(thresholdScalarVolume, None, parameters) - - # write a temporary volume to disk - volumeFileName = "AutoscoperM_VRG_GEN_TEMP.mhd" - IO.writeTemporyFile(volumeFileName, self.convertNodeToData(volumeNode)) - - # Execute CLI for each camera - cliModule = slicer.modules.virtualradiographgeneration - parameters = { - "inputVolumeFileName": os.path.join(slicer.app.temporaryPath, volumeFileName), - "cameraDir": cameraDir, - "radiographMainOutDir": outputDir, - "outputFileName": f"{filename}.tif", - "outputWidth": size[0], - "outputHeight": size[1], - } - cliNode = slicer.cli.run(cliModule, None, parameters) # run asynchronously - - # Note: CLI nodes are currently not executed in parallel. See https://github.com/Slicer/Slicer/pull/6723 - # This just allows the UI to remain responsive while the CLI nodes are running for now. - - # Wait for all the CLI nodes to finish - while cliNode.GetStatusString() != "Completed": - slicer.app.processEvents() - if cliNode.GetStatus() & cliNode.ErrorsMask: - # error - errorText = cliNode.GetErrorText() - slicer.mrmlScene.RemoveNode(cliNode) - raise ValueError("CLI execution failed: " + errorText) - slicer.mrmlScene.RemoveNode(cliNode) - - def generateCameraCalibrationFiles( - self, - bestCameras: list[RadiographGeneration.Camera], - tmpDir: str, - finalDir: str, - calibDir: str, - progressCallback: Optional[callable] = None, - ) -> None: - """ - Copies the optimized VRGs from the temporary directory to the final directory - and generates the camera calibration files - - :param bestCameras: list of optimized cameras - :param tmpDir: temporary directory - :param finalDir: final directory - :param calibDir: calibration directory - :param progressCallback: progress callback, defaults to None - """ - self.validatePaths(tmpDir=tmpDir) - self.createPathsIfNotExists(finalDir, calibDir) - if not progressCallback: - logging.warning( - "[AutoscoperM.logic.generateCameraCalibrationFiles] " - "No progress callback provided, progress bar will not be updated" - ) - - def progressCallback(x): - return x - - for idx, cam in enumerate(bestCameras): - IO.generateCameraCalibrationFile(cam, os.path.join(calibDir, f"cam{cam.id}.json")) - cameraDir = os.path.join(finalDir, f"cam{cam.id}") - self.createPathsIfNotExists(cameraDir) - # Copy all tif files from the tmp to the final directory - for file in glob.glob(os.path.join(tmpDir, f"cam{cam.id}", "*.tif")): - shutil.copy(file, cameraDir) - - progress = ((idx + 1) / len(bestCameras)) * 10 + 90 - progressCallback(progress) - - @staticmethod - def convertNodeToData(volumeNode: slicer.vtkMRMLVolumeNode) -> vtk.vtkImageData: - """ - Converts a volume node to a vtkImageData object - """ - imageData = vtk.vtkImageData() - imageData.DeepCopy(volumeNode.GetImageData()) - imageData.SetSpacing(volumeNode.GetSpacing()) - origin = list(volumeNode.GetOrigin()) - imageData.SetOrigin(origin) - - mat = vtk.vtkMatrix4x4() - volumeNode.GetIJKToRASMatrix(mat) - if mat.GetElement(0, 0) < 0 and mat.GetElement(1, 1) < 0: - origin[0:2] = [x * -1 for x in origin[0:2]] - imageData.SetOrigin(origin) - - # Ensure we are in the correct orientation (RAS vs LPS) - imageReslice = vtk.vtkImageReslice() - imageReslice.SetInputData(imageData) - - axes = vtk.vtkMatrix4x4() - axes.Identity() - axes.SetElement(0, 0, -1) - axes.SetElement(1, 1, -1) - - imageReslice.SetResliceAxes(axes) - imageReslice.Update() - imageData = imageReslice.GetOutput() - - return imageData - @staticmethod def getItemInSequence(sequenceNode: slicer.vtkMRMLSequenceNode, idx: int) -> slicer.vtkMRMLNode: """ @@ -1492,21 +1096,6 @@ def GetVolumeSpacing(node: Union[slicer.vtkMRMLVolumeNode, slicer.vtkMRMLSequenc return AutoscoperMLogic.getItemInSequence(node, 0)[0].GetSpacing() return node.GetSpacing() - @staticmethod - def GetRASBounds(node: Union[slicer.vtkMRMLVolumeNode, slicer.vtkMRMLSequenceNode]) -> list[float]: - bounds = [0] * 6 - if AutoscoperMLogic.IsSequenceVolume(node): - AutoscoperMLogic.getItemInSequence(node, 0)[0].GetRASBounds(bounds) - else: - node.GetRASBounds(bounds) - return bounds - - @staticmethod - def IsVolumeCentered(node: Union[slicer.vtkMRMLVolumeNode, slicer.vtkMRMLSequenceNode]) -> bool: - if AutoscoperMLogic.IsSequenceVolume(node): - return AutoscoperMLogic.getItemInSequence(node, 0)[0].IsCentered() - return node.IsCentered() - @staticmethod def loadTransformFromFile(transformFileName: str) -> slicer.vtkMRMLLinearTransformNode: return slicer.util.loadNodeFromFile(transformFileName) diff --git a/AutoscoperM/AutoscoperMLib/IO.py b/AutoscoperM/AutoscoperMLib/IO.py index b8ccd52..43c561e 100644 --- a/AutoscoperM/AutoscoperMLib/IO.py +++ b/AutoscoperM/AutoscoperMLib/IO.py @@ -7,8 +7,6 @@ import slicer import vtk -from .RadiographGeneration import Camera - def loadSegmentation(segmentationNode: slicer.vtkMRMLSegmentationNode, filename: str): """ @@ -35,34 +33,6 @@ def loadSegmentation(segmentationNode: slicer.vtkMRMLSegmentationNode, filename: return None -def generateCameraCalibrationFile(camera: Camera, filename: str): - """ - Generates a VTK camera calibration json file from the given camera. - - :param camera: Camera - :param filename: Output file name - """ - import json - - contents = {} - contents["@schema"] = "https://autoscoperm.slicer.org/vtkCamera-schema-1.0.json" - contents["version"] = 1.0 - contents["focal-point"] = camera.vtkCamera.GetFocalPoint() - contents["camera-position"] = camera.vtkCamera.GetPosition() - contents["view-up"] = camera.vtkCamera.GetViewUp() - contents["view-angle"] = camera.vtkCamera.GetViewAngle() - contents["image-width"] = camera.imageSize[0] - contents["image-height"] = camera.imageSize[1] - # The clipping-range field is not used by Autoscoper, it is used to communicate - # information between AutoscoperM and VirtualRadiographGeneration modules within Slicer. - contents["clipping-range"] = camera.vtkCamera.GetClippingRange() - - contents_json = json.dumps(contents, indent=4) - - with open(filename, "w+") as f: - f.write(contents_json) - - def generateConfigFile( mainDirectory: str, subDirectories: list[str], @@ -265,7 +235,7 @@ def writeTRA(fileName: str, transforms: list[vtk.vtkMatrix4x4]) -> None: traFile.write(",".join(row) + "\n") -def writeTemporyFile(filename: str, data: vtk.vtkImageData) -> str: +def writeTemporaryFile(filename: str, data: vtk.vtkImageData) -> str: """ Writes a temporary file to the slicer temp directory @@ -285,7 +255,7 @@ def writeTemporyFile(filename: str, data: vtk.vtkImageData) -> str: return writer.GetFileName() -def removeTemporyFile(filename: str): +def removeTemporaryFile(filename: str): """ Removes a temporary file from the slicer temp directory diff --git a/AutoscoperM/AutoscoperMLib/RadiographGeneration.py b/AutoscoperM/AutoscoperMLib/RadiographGeneration.py deleted file mode 100644 index 86566a6..0000000 --- a/AutoscoperM/AutoscoperMLib/RadiographGeneration.py +++ /dev/null @@ -1,239 +0,0 @@ -import math -from typing import Optional - -import slicer -import vtk - - -class Camera: - def __init__(self) -> None: - self.DID = 0 - self.vtkCamera = vtk.vtkCamera() - self.imageSize = [512, 512] - self.id = -1 - self.FrustumModel = None - - def __str__(self) -> str: - return "\n".join( - [ - f"Camera {self.id}", - f"Position: {self.vtkCamera.GetPosition()}", - f"Focal Point: {self.vtkCamera.GetFocalPoint()}", - f"View Angle: {self.vtkCamera.GetViewAngle()}", - f"Clipping Range: {self.vtkCamera.GetClippingRange()}", - f"View Up: {self.vtkCamera.GetViewUp()}", - f"Direction of Projection: {self.vtkCamera.GetDirectionOfProjection()}", - f"Distance: {self.vtkCamera.GetDistance()}", - f"Image Size: {self.imageSize}", - f"DID: {self.DID}", - "~" * 20, - ] - ) - - -def _createFrustumModel(cam: Camera) -> None: - model = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLModelNode") - model.CreateDefaultDisplayNodes() - model.GetDisplayNode().SetColor(1, 0, 0) # Red - model.GetDisplayNode().SetOpacity(0.3) - model.GetDisplayNode().SetVisibility(True) - - model.SetName(f"cam{cam.id}-frustum") - - cam.FrustumModel = model - - _updateFrustumModel(cam) - - -def _updateFrustumModel(cam: Camera) -> None: - if cam.FrustumModel is None: - _createFrustumModel(cam) - return - # The equations of the six planes of the frustum in the order: left, right, bottom, top, far, near - # Given as A, B, C, D where Ax + By + Cz + D = 0 for each plane - planesArray = [0] * 24 - aspectRatio = cam.vtkCamera.GetExplicitAspectRatio() - - cam.vtkCamera.GetFrustumPlanes(aspectRatio, planesArray) - - planes = vtk.vtkPlanes() - planes.SetFrustumPlanes(planesArray) - - hull = vtk.vtkHull() - hull.SetPlanes(planes) - pd = vtk.vtkPolyData() - hull.GenerateHull(pd, [-1000, 1000, -1000, 1000, -1000, 1000]) - - cam.FrustumModel.SetAndObservePolyData(pd) - - -def generateNCameras( - N: int, bounds: list[int], offset: int = 100, imageSize: tuple[int] = (512, 512), camDebugMode: bool = False -) -> list[Camera]: - """ - Generate N cameras - - :param N: Number of cameras to generate - :param bounds: Bounds of the volume - :param offset: Offset from the volume. Defaults to 100. - :param imageSize: Image size. Defaults to [512,512]. - :param camDebugMode: Whether or not to show the cameras in the scene. Defaults to False. - - :return: List of cameras - """ - # find the center of the bounds - center = [(bounds[0] + bounds[1]) / 2, (bounds[2] + bounds[3]) / 2, (bounds[4] + bounds[5]) / 2] - - # find the largest dimension of the bounds - largestDimension = max([bounds[1] - bounds[0], bounds[3] - bounds[2], bounds[5] - bounds[4]]) - - # find the distance from the center to the bounds - r = largestDimension / 2 + offset - - points = vtk.vtkPoints() - points.SetNumberOfPoints(N) - points.SetDataTypeToDouble() - points.Allocate(N) - - # use the spherical fibonacci algorithm to generate the points - goldenRatio = (1 + math.sqrt(5)) / 2 - i = range(0, N) - theta = [2 * math.pi * i / goldenRatio for i in i] - phi = [math.acos(1 - 2 * (i + 0.5) / N) for i in i] - x = [math.cos(theta[i]) * math.sin(phi[i]) for i in i] - y = [math.sin(theta[i]) * math.sin(phi[i]) for i in i] - z = [math.cos(phi[i]) for i in i] - # scale the points to the radius - x = [r * x[i] for i in i] - y = [r * y[i] for i in i] - z = [r * z[i] for i in i] - for px, py, pz in zip(x, y, z): - points.InsertNextPoint(px + center[0], py + center[1], pz + center[2]) - - # create the cameras - cameras = [] - for i in range(N): - camera = Camera() - camera.vtkCamera.SetPosition(points.GetPoint(i)) - camera.vtkCamera.SetFocalPoint(center) - camera.vtkCamera.SetViewAngle(30) - camera.vtkCamera.SetClippingRange(0.1, r + largestDimension) - camera.vtkCamera.SetViewUp(0, 1, 0) - camera.id = i - camera.imageSize = imageSize - cameras.append(camera) - - if camDebugMode: - # Visualize the cameras - markups = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsFiducialNode") - for cam in cameras: - print(cam) - # add a point to the markups node - markups.AddControlPoint( - cam.vtkCamera.GetPosition()[0], - cam.vtkCamera.GetPosition()[1], - cam.vtkCamera.GetPosition()[2], - f"cam{cam.id}", - ) - # lock the point - markups.SetNthControlPointLocked(markups.GetNumberOfControlPoints() - 1, True) - - _createFrustumModel(cam) - - return cameras - - -def generateCamerasFromMarkups( - fiduaicalNode: slicer.vtkMRMLMarkupsFiducialNode, - volumeBounds: list[int], - clippingRange: tuple[int], - viewAngle: int, - imageSize: tuple[int] = (512, 512), - cameraDebug: bool = False, -) -> list[Camera]: - """ - Generate cameras from a markups fiducial node - - :param fiduaicalNode: Markups fiducial node - :param volumeBounds: Bounds of the volume - :param clippingRange: Clipping range - :param viewAngle: View angle - :param imageSize: Image size. Defaults to [512,512]. - :param cameraDebug: Whether or not to show the cameras in the scene. Defaults to False. - - :return: List of cameras - """ - center = [ - (volumeBounds[0] + volumeBounds[1]) / 2, - (volumeBounds[2] + volumeBounds[3]) / 2, - (volumeBounds[4] + volumeBounds[5]) / 2, - ] - n = fiduaicalNode.GetNumberOfControlPoints() - cameras = [] - for i in range(n): - camera = Camera() - camera.vtkCamera.SetPosition(fiduaicalNode.GetNthControlPointPosition(i)) - camera.vtkCamera.SetFocalPoint(center) - camera.vtkCamera.SetViewAngle(viewAngle) - camera.vtkCamera.SetClippingRange(clippingRange[0], clippingRange[1]) - camera.vtkCamera.SetViewUp(0, 1, 0) - camera.id = fiduaicalNode.GetNthControlPointLabel(i) - camera.imageSize = imageSize - if cameraDebug: - _createFrustumModel(camera) - cameras.append(camera) - return cameras - - -def optimizeCameras( - cameras: list[Camera], - cameraDir: str, - nOptimizedCameras: int, - progressCallback: Optional[callable] = None, -) -> list[Camera]: - """ - Optimize the cameras by finding the N cameras with the best data intensity density. - - :param cameras: Cameras - :param cameraDir: Camera directory - :param nOptimizedCameras: Number of optimized cameras to find - - :return: Optimized cameras - """ - import os - - if not progressCallback: - - def progressCallback(_x): - return None - - # Parallel calls to cliModule - cliModule = slicer.modules.calculatedataintensitydensity - cliNodes = [] - for i in range(len(cameras)): - camera = cameras[i] - vrgDirName = os.path.join(cameraDir, f"cam{camera.id}") - cliNode = slicer.cli.run( - cliModule, - None, - {"whiteRadiographDirName": vrgDirName}, - wait_for_completion=False, - ) - cliNodes.append(cliNode) - - for i in range(len(cameras)): - while cliNodes[i].GetStatusString() != "Completed": - slicer.app.processEvents() - if cliNode.GetStatus() & cliNode.ErrorsMask: - # error - errorText = cliNode.GetErrorText() - slicer.mrmlScene.RemoveNode(cliNode) - raise ValueError("CLI execution failed: " + errorText) - cameras[i].DID = float(cliNodes[i].GetOutputText()) # cliNodes[i].GetParameterAsString("dataIntensityDensity") - progress = ((i + 1) / len(cameras)) * 40 + 50 - progressCallback(progress) - slicer.mrmlScene.RemoveNode(cliNodes[i]) - - cameras.sort(key=lambda x: x.DID) - - return cameras[:nOptimizedCameras] diff --git a/AutoscoperM/CMakeLists.txt b/AutoscoperM/CMakeLists.txt index 9af705f..f752650 100644 --- a/AutoscoperM/CMakeLists.txt +++ b/AutoscoperM/CMakeLists.txt @@ -6,7 +6,6 @@ set(MODULE_PYTHON_SCRIPTS ${MODULE_NAME}.py ${MODULE_NAME}Lib/__init__.py ${MODULE_NAME}Lib/IO.py - ${MODULE_NAME}Lib/RadiographGeneration.py ${MODULE_NAME}Lib/SubVolumeExtraction.py ) diff --git a/AutoscoperM/Resources/UI/AutoscoperM.ui b/AutoscoperM/Resources/UI/AutoscoperM.ui index 19e9e66..761cc9d 100644 --- a/AutoscoperM/Resources/UI/AutoscoperM.ui +++ b/AutoscoperM/Resources/UI/AutoscoperM.ui @@ -218,36 +218,6 @@ - - - - true - - - Delete Temporary VRG Files - - - true - - - true - - - - - - - Only use indices for radiograph filename - - - - - - - VRG Temp Subdirectory: - - - @@ -286,16 +256,6 @@ - - - - Camera Debug Mode - - - false - - - @@ -337,13 +297,6 @@ - - - - VRG-Temp - - - @@ -523,217 +476,6 @@ - - - - VRG Generation - Manual Camera Placement - - - false - - - - - - Camera Positions - - - - - - - 0.100000000000000 - - - 2000.000000000000000 - - - - - - - Clipping Range - - - - - - - View Angle - - - - - - - true - - - - vtkMRMLMarkupsFiducialNode - - - - - - - - - - - - - - Generate VRGs from Markups - - - - - - - 0.100000000000000 - - - 2000.000000000000000 - - - 0.100000000000000 - - - 0.100000000000000 - - - 300.000000000000000 - - - Qt::Horizontal - - - - - - - 360 - - - 30 - - - - - - - 0.100000000000000 - - - 2000.000000000000000 - - - 300.000000000000000 - - - - - - - - - - true - - - VRG Generation - Automatic Camera Placement - - - false - - - true - - - - - - 0 - - - 1000 - - - 400 - - - - - - - Camera Offset: - - - - - - - # of Optimized Cameras: - - - - - - - 10 - - - 50 - - - - - - - # of Possible Cameras: - - - - - - - false - - - 0 - - - 1000 - - - 400 - - - Qt::Horizontal - - - - - - - Generate VRGs - - - - - - - 2 - - - 2 - - - - - - @@ -996,85 +738,5 @@ - - camOffSetSpin - valueChanged(int) - camOffSetSlider - setValue(int) - - - 953 - 667 - - - 553 - 667 - - - - - camOffSetSlider - valueChanged(int) - camOffSetSpin - setValue(int) - - - 553 - 667 - - - 953 - 667 - - - - - AutoscoperM - mrmlSceneChanged(vtkMRMLScene*) - mVRG_markupSelector - setMRMLScene(vtkMRMLScene*) - - - 508 - 429 - - - 507 - 409 - - - - - mVRG_ClippingRangeSlider - minimumValueChanged(double) - mVRG_clippingRangeMinBox - setValue(double) - - - 567 - 499 - - - 189 - 499 - - - - - mVRG_ClippingRangeSlider - maximumValueChanged(double) - mVRG_clippingRangeMaxBox - setValue(double) - - - 567 - 499 - - - 944 - 499 - - - diff --git a/CMakeLists.txt b/CMakeLists.txt index 22c84bd..d5efaa5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -42,8 +42,6 @@ endif() # Extension modules add_subdirectory(AutoscoperM) add_subdirectory(TrackingEvaluation) -add_subdirectory(CalculateDataIntensityDensity) -add_subdirectory(VirtualRadiographGeneration) add_subdirectory(Hierarchical3DRegistration) ## NEXT_MODULE diff --git a/CalculateDataIntensityDensity/CMakeLists.txt b/CalculateDataIntensityDensity/CMakeLists.txt deleted file mode 100644 index 1588932..0000000 --- a/CalculateDataIntensityDensity/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -#----------------------------------------------------------------------------- -set(MODULE_NAME CalculateDataIntensityDensity) - -SlicerMacroBuildScriptedCLI( - NAME ${MODULE_NAME} - ) diff --git a/CalculateDataIntensityDensity/CalculateDataIntensityDensity.py b/CalculateDataIntensityDensity/CalculateDataIntensityDensity.py deleted file mode 100644 index ca123ee..0000000 --- a/CalculateDataIntensityDensity/CalculateDataIntensityDensity.py +++ /dev/null @@ -1,95 +0,0 @@ -#!/usr/bin/env python-real - -import concurrent.futures as cf -import glob -import os -import sys - -import numpy as np -import SimpleITK as sitk - - -def calcDID(whiteRadiographFName): - MEAN_COMPARISON = 185 - # Read in the white radiograph - whiteRadiograph = sitk.ReadImage(whiteRadiographFName) - - # Superpixel Segmentation - slicImageFilter = sitk.SLICImageFilter() - slicImageFilter.SetSuperGridSize([85, 85, 85]) - 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)) - - -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`. - - :param whiteRadiographFName: White radiograph file name - - return Data intensity density - """ - 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}") - - if len(whiteRadiographFiles) > 1: # Avoid overhead in the single 3DCT case - dids = [] - with cf.ThreadPoolExecutor() as executor: - futures = [executor.submit(calcDID, wrFName) for wrFName in whiteRadiographFiles] - for future in cf.as_completed(futures): - dids.append(future.result()) - return np.mean(dids) - return calcDID(whiteRadiographFiles[0]) - - -if __name__ == "__main__": - expected_args = [ - "whiteRadiographFileName", - # Value reported on standard output - "DID", - ] - expected_args = [f"<{arg}>" for arg in expected_args] - if len(sys.argv) < len(expected_args): - print(f"Usage: {sys.argv[0]} {' '.join(expected_args)}") - sys.exit(1) - print(main(sys.argv[1])) diff --git a/CalculateDataIntensityDensity/CalculateDataIntensityDensity.xml b/CalculateDataIntensityDensity/CalculateDataIntensityDensity.xml deleted file mode 100644 index eaaac5e..0000000 --- a/CalculateDataIntensityDensity/CalculateDataIntensityDensity.xml +++ /dev/null @@ -1,41 +0,0 @@ - - - Tracking.Advanced - 0 - CalculateDataIntensityDensity - - 0.1.0. - https://autoscoperm.slicer.org/ - - Anthony Lombardi (Kitware) - Amy M Morton (Brown University) - Bardiya Akhbari (Brown University) - Beatriz Paniagua (Kitware) - Jean-Christophe Fillion-Robin (Kitware) - - - - whiteRadiographDirName - - 1 - - input - - - DID - - 2 - - output - - - outliers - - 3 - - output - - - - - diff --git a/Hierarchical3DRegistration/Hierarchical3DRegistration.py b/Hierarchical3DRegistration/Hierarchical3DRegistration.py index 9e03a04..5317e2f 100644 --- a/Hierarchical3DRegistration/Hierarchical3DRegistration.py +++ b/Hierarchical3DRegistration/Hierarchical3DRegistration.py @@ -30,10 +30,6 @@ def __init__(self, parent): self.parent.categories = [ "Tracking", ] - self.parent.dependencies = [ - "CalculateDataIntensityDensity", - "VirtualRadiographGeneration", - ] self.parent.contributors = [ "Anthony Lombardi (Kitware)", "Amy M Morton (Brown University)", diff --git a/VirtualRadiographGeneration/CMakeLists.txt b/VirtualRadiographGeneration/CMakeLists.txt deleted file mode 100644 index 02872d6..0000000 --- a/VirtualRadiographGeneration/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -#----------------------------------------------------------------------------- -set(MODULE_NAME VirtualRadiographGeneration) - -SlicerMacroBuildScriptedCLI( - NAME ${MODULE_NAME} - ) diff --git a/VirtualRadiographGeneration/VirtualRadiographGeneration.py b/VirtualRadiographGeneration/VirtualRadiographGeneration.py deleted file mode 100644 index 721155c..0000000 --- a/VirtualRadiographGeneration/VirtualRadiographGeneration.py +++ /dev/null @@ -1,168 +0,0 @@ -#!/usr/bin/env python-real - -import concurrent.futures as cf -import glob -import json -import os -import sys - -import vtk - - -def generateVRG( - camera: vtk.vtkCamera, - volumeImageData: vtk.vtkImageData, - outputFileName: str, - width: int, - height: int, -) -> None: - """ - Generate a virtual radiograph from the given camera and volume node - - :param camera: Camera - :param volumeImageData: Volume image data - :param outputFileName: Output file name - :param width: Width of the output image - :param height: Height of the output image - """ - - # find the min and max scalar values - hist = vtk.vtkImageHistogramStatistics() - hist.SetInputData(volumeImageData) - hist.Update() - minVal = hist.GetMinimum() - maxVal = hist.GetMaximum() - - # create the renderer - renderer = vtk.vtkRenderer() - renderer.SetBackground(1, 1, 1) # Set background to white - renderer.SetUseDepthPeeling(True) - - # create the render window - renderWindow = vtk.vtkRenderWindow() - renderWindow.SetOffScreenRendering(1) - renderWindow.SetSize(width, height) - renderWindow.AddRenderer(renderer) - - # create the volume mapper - volumeMapper = vtk.vtkGPUVolumeRayCastMapper() - volumeMapper.SetInputData(volumeImageData) - volumeMapper.SetBlendModeToComposite() - volumeMapper.SetBlendModeToComposite() - volumeMapper.SetUseJittering(False) - - # Set the transfer functions for opacity, gradient and color - opacityTransferFunction = vtk.vtkPiecewiseFunction() # From the Slicer CT XRay preset - opacityTransferFunction.AddPoint(minVal, 0.0) - opacityTransferFunction.AddPoint(1500, 0.05) - opacityTransferFunction.AddPoint(maxVal, 0.05) - - gradTransferFunction = vtk.vtkPiecewiseFunction() # From the Slicer CT XRay preset - gradTransferFunction.AddPoint(0, 1) - gradTransferFunction.AddPoint(255, 1) - - colorTransferFunction = vtk.vtkColorTransferFunction() - colorTransferFunction.AddRGBPoint(maxVal, 1, 1, 1) - colorTransferFunction.AddRGBPoint(minVal, 0, 0, 0) - - volumeProperty = vtk.vtkVolumeProperty() - volumeProperty.SetInterpolationTypeToLinear() - volumeProperty.ShadeOff() - volumeProperty.SetScalarOpacity(opacityTransferFunction) - volumeProperty.SetGradientOpacity(gradTransferFunction) - volumeProperty.SetColor(colorTransferFunction) - - # create the volume - volume = vtk.vtkVolume() - volume.SetMapper(volumeMapper) - volume.SetProperty(volumeProperty) - - # add the volume to the renderer - renderer.AddVolume(volume) - renderer.SetActiveCamera(camera) - - # render the image - renderWindow.Render() - - # save the image - writer = vtk.vtkTIFFWriter() - writer.SetFileName(outputFileName) - - windowToImageFilter = vtk.vtkWindowToImageFilter() - windowToImageFilter.SetInput(renderWindow) - windowToImageFilter.SetScale(1) - windowToImageFilter.SetInputBufferTypeToRGB() - - # convert the image to grayscale - luminance = vtk.vtkImageLuminance() - luminance.SetInputConnection(windowToImageFilter.GetOutputPort()) - - writer.SetInputConnection(luminance.GetOutputPort()) - writer.Write() - - -def _createVTKCameras(cameraDir: str, radiographMainDir: str): - """ - Generates a vtkCamera object from the given parameters - """ - cameras = {} - for camFName in glob.glob(os.path.join(cameraDir, "*.json")): - with open(camFName) as f: - camJSON = json.load(f) - - cam = vtk.vtkCamera() - cam.SetPosition(camJSON["camera-position"]) - cam.SetFocalPoint(camJSON["focal-point"]) - cam.SetViewUp(camJSON["view-up"]) - cam.SetViewAngle(camJSON["view-angle"]) - cam.SetClippingRange(camJSON["clipping-range"]) - - cameraSubDirName = os.path.basename(camFName).split(".")[0] - cameraDirName = os.path.join(radiographMainDir, cameraSubDirName) - if not os.path.exists(cameraDirName): - os.mkdir(cameraDirName) - - cameras[cam] = cameraDirName - - return cameras - - -if __name__ == "__main__": - expected_args = [ - "inputVolumeFileName", - "cameraDir", - "radiographMainOutDir", - "outputFileName", - "outputWidth", - "outputHeight", - ] - expected_args = [f"<{arg}>" for arg in expected_args] - if len(sys.argv[1:]) != len(expected_args): - print(f"Usage: {sys.argv[0]} {' '.join(expected_args)}") - sys.exit(1) - - volumeData = sys.argv[1] - cameraDir = sys.argv[2] - radiographMainOutDir = sys.argv[3] - outputFileName = sys.argv[4] - outputWidth = int(sys.argv[5]) - outputHeight = int(sys.argv[6]) - - # create the camera - cameras = _createVTKCameras(cameraDir, radiographMainOutDir) - - # Read the mhd file - reader = vtk.vtkMetaImageReader() - reader.SetFileName(volumeData) - reader.Update() - - # generate the virtual radiograph - with cf.ThreadPoolExecutor() as executor: - futures = [ - executor.submit( - generateVRG, cam, reader.GetOutput(), os.path.join(camDir, outputFileName), outputWidth, outputHeight - ) - for cam, camDir in cameras.items() - ] - for future in cf.as_completed(futures): - future.result() diff --git a/VirtualRadiographGeneration/VirtualRadiographGeneration.xml b/VirtualRadiographGeneration/VirtualRadiographGeneration.xml deleted file mode 100644 index 79a2209..0000000 --- a/VirtualRadiographGeneration/VirtualRadiographGeneration.xml +++ /dev/null @@ -1,62 +0,0 @@ - - - Tracking.Advanced - 0 - VirtualRadiographGeneration - - 0.1.0. - https://autoscoperm.slicer.org/ - - Anthony Lombardi (Kitware) - Amy M Morton (Brown University) - Bardiya Akhbari (Brown University) - Beatriz Paniagua (Kitware) - Jean-Christophe Fillion-Robin (Kitware) - - - - inputVolumeFileName - - 0 - - input - - - cameraDir - - 2 - - input - - - radiographMainOutDir - - 3 - - input - - - outputFileName - - 4 - - input - - - outputWidth - - 5 - - input - - - outputHeight - - 6 - - input - - - - -