4dct

.pre-commit-config.yaml

@@ -15,7 +15,7 @@ repos:
   rev: v0.0.282
   hooks:
   - id: ruff
-    args: ["--fix", "--show-fixes"]
+    args: ["--fix", "--show-fixes","--extend-ignore=EXE001"]
 
 - repo: https://github.com/psf/black
   rev: 22.3.0

AutoscoperM/AutoscoperMLib/IO.py

@@ -74,6 +74,15 @@ def generateCameraCalibrationFile(camera: Camera, filename: str):
             + str(camera.vtkCamera.GetPosition()[2])
             + "]\n"
         )
+        f.write(
+            "view-up: ["
+            + str(camera.vtkCamera.GetViewUp()[0])
+            + ", "
+            + str(camera.vtkCamera.GetViewUp()[1])
+            + ", "
+            + str(camera.vtkCamera.GetViewUp()[2])
+            + "]\n"
+        )
         f.write("view-angle: " + str(camera.vtkCamera.GetViewAngle()) + "\n")
         f.write("image-width: " + str(camera.imageSize[0]) + "\n")
         f.write("image-height: " + str(camera.imageSize[1]) + "\n")
@@ -295,3 +304,35 @@ def writeTFMFile(filename: str, spacing: list[float], origin: list[float]):
     slicer.util.exportNode(transformNode, filename)
 
     slicer.mrmlScene.RemoveNode(transformNode)
+
+
+def writeTemporyFile(filename: str, data: vtk.vtkImageData):
+    """
+    Writes a temporary file to the slicer temp directory
+
+    :param filename: Output file name
+    :type filename: str
+
+    :param data: data
+    :type data: vtk.vtkImageData
+    """
+
+    slicerTempDirectory = slicer.app.temporaryPath
+
+    # write vtk image data as a vtk file
+    writer = vtk.vtkMetaImageWriter()
+    writer.SetFileName(os.path.join(slicerTempDirectory, filename))
+    writer.SetInputData(data)
+    writer.Write()
+
+
+def removeTemporyFile(filename: str):
+    """
+    Removes a temporary file from the slicer temp directory
+
+    :param filename: Output file name
+    :type filename: str
+    """
+
+    slicerTempDirectory = slicer.app.temporaryPath
+    os.remove(os.path.join(slicerTempDirectory, filename))

AutoscoperM/AutoscoperMLib/RadiographGeneration.py

@@ -11,9 +11,60 @@ def __init__(self) -> None:
         self.vtkCamera = vtk.vtkCamera()
         self.imageSize = [512, 512]
         self.id = -1
+        self.FrustumModel = None
 
+    def __str__(self) -> str:
+        dbgStr = f"Camera {self.id}\n"
+        dbgStr += f"Position: {self.vtkCamera.GetPosition()}\n"
+        dbgStr += f"Focal Point: {self.vtkCamera.GetFocalPoint()}\n"
+        dbgStr += f"View Angle: {self.vtkCamera.GetViewAngle()}\n"
+        dbgStr += f"Clipping Range: {self.vtkCamera.GetClippingRange()}\n"
+        dbgStr += f"View Up: {self.vtkCamera.GetViewUp()}\n"
+        dbgStr += f"Direction of Projection: {self.vtkCamera.GetDirectionOfProjection()}\n"
+        dbgStr += f"Distance: {self.vtkCamera.GetDistance()}\n"
+        dbgStr += f"Image Size: {self.imageSize}\n"
+        dbgStr += f"DID: {self.DID}\n"
+        dbgStr += "~" * 20 + "\n"
+        return dbgStr
 
-def generateNCameras(N: int, bounds: list[int], offset: int = 100, imageSize: tuple[int] = (512, 512)) -> list[Camera]:
+
+def _createFrustumModel(cam: Camera) -> None:
+    model = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLModelNode")
+    model.CreateDefaultDisplayNodes()
+    model.GetDisplayNode().SetColor(1, 0, 0)
+    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
+    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
 
@@ -29,6 +80,9 @@ def generateNCameras(N: int, bounds: list[int], offset: int = 100, imageSize: tu
     :param imageSize: Image size. Defaults to [512,512].
     :type imageSize: list[int]
 
+    :param camDebugMode: Whether or not to show the cameras in the scene. Defaults to False.
+    :type camDebugMode: bool
+
     :return: List of cameras
     :rtype: list[Camera]
     """
@@ -67,18 +121,83 @@ def generateNCameras(N: int, bounds: list[int], offset: int = 100, imageSize: tu
         camera = Camera()
         camera.vtkCamera.SetPosition(points.GetPoint(i))
         camera.vtkCamera.SetFocalPoint(center)
-        camera.vtkCamera.SetViewAngle(25)
-        camera.vtkCamera.SetClippingRange(0.1, 1000)
+        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
+    :type fiduaicalNode: slicer.vtkMRMLMarkupsFiducialNode
+    :param volumeBounds: Bounds of the volume
+    :type volumeBounds: list[int]
+    :param clippingRange: Clipping range
+    :type clippingRange: tuple[int]
+    :param viewAngle: View angle
+    :type viewAngle: int
+    :param imageSize: Image size. Defaults to [512,512].
+    :type imageSize: list[int]
+    :param cameraDebug: Whether or not to show the cameras in the scene. Defaults to False.
+    :type cameraDebug: bool
+    :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 generateVRG(
     camera: Camera,
-    volumeNode: slicer.vtkMRMLVolumeNode,
+    volumeImageData: vtk.vtkImageData,
     outputFileName: str,
     width: int,
     height: int,
@@ -89,8 +208,8 @@ def generateVRG(
     :param camera: Camera
     :type camera: Camera
 
-    :param volumeNode: Volume node
-    :type volumeNode: slicer.vtkMRMLVolumeNode
+    :param volumeImageData: Volume image data
+    :type volumeImageData: vtk.vtkImageData
 
     :param outputFileName: Output file name
     :type outputFileName: str
@@ -105,19 +224,19 @@ def generateVRG(
     # create the renderer
     renderer = vtk.vtkRenderer()
     renderer.SetBackground(1, 1, 1)
-    renderer.SetUseDepthPeeling(1)
-    renderer.SetMaximumNumberOfPeels(100)
-    renderer.SetOcclusionRatio(0.1)
 
     # create the render window
     renderWindow = vtk.vtkRenderWindow()
-    renderWindow.SetOffScreenRendering(1)
+    # renderWindow.SetOffScreenRendering(1)
     renderWindow.SetSize(width, height)
     renderWindow.AddRenderer(renderer)
 
+    renderWindowInteractor = vtk.vtkRenderWindowInteractor()
+    renderWindowInteractor.SetRenderWindow(renderWindow)
+
     # create the volume mapper
     volumeMapper = vtk.vtkGPUVolumeRayCastMapper()
-    volumeMapper.SetInputData(volumeNode.GetImageData())
+    volumeMapper.SetInputData(volumeImageData)
     volumeMapper.SetBlendModeToComposite()
 
     # create the volume property
@@ -129,8 +248,6 @@ def generateVRG(
     volumeProperty.SetSpecular(0.2)
 
     # create a piecewise function for scalar opacity
-    # first point is X: 300, O: 0.00
-    # second point is X: 950, O: 0.20
     opacityTransferFunction = vtk.vtkPiecewiseFunction()
     opacityTransferFunction.AddPoint(-10000, 0.05)
     opacityTransferFunction.AddPoint(0, 0.00)
@@ -145,11 +262,12 @@ def generateVRG(
     # add the volume to the renderer
     renderer.AddVolume(volume)
     renderer.SetActiveCamera(camera.vtkCamera)
-    renderer.ResetCamera()
 
     # render the image
     renderWindow.Render()
 
+    renderWindowInteractor.Start()
+
     # save the image
     writer = vtk.vtkTIFFWriter()
     writer.SetFileName(outputFileName)
@@ -178,6 +296,7 @@ def _calculateDataIntensityDensity(camera: Camera, whiteRadiographFName: str) ->
     :param whiteRadiographFName: White radiograph file name
     :type whiteRadiographFName: str
     """
+
     import numpy as np
     import SimpleITK as sitk
 
@@ -255,13 +374,34 @@ def optimizeCameras(
         def progressCallback(_x):
             return None
 
+    # Parallel calls to cliModule
+    cliModule = slicer.modules.calculatedataintensitydensity
+    cliNodes = []
     for i in range(len(cameras)):
         camera = cameras[i]
         vrgFName = glob.glob(os.path.join(cameraDir, f"cam{camera.id}", "*.tif"))[0]
-        _calculateDataIntensityDensity(camera, vrgFName)
-        progress = ((i + 1) / len(cameras)) * 29 + 42
+        # _calculateDataIntensityDensity(camera, vrgFName)
+        cliNode = slicer.cli.run(
+            cliModule,
+            None,
+            {"whiteRadiographFName": vrgFName},
+            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)) * 50 + 40
         progressCallback(progress)
+        slicer.mrmlScene.RemoveNode(cliNodes[i])
 
-    cameras.sort(key=lambda x: x.DID, reverse=True)
+    cameras.sort(key=lambda x: x.DID)
 
     return cameras[:nOptimizedCameras]