Merge pull request #62 from arcadelab/dev

Dev

README.md

@@ -82,7 +82,8 @@ in the base directory. Then do `cd docs` and `make html` to build the static sit
 The following minimal example loads a CT volume from a NifTi `.nii.gz` file and simulates an X-ray projection:
 
 ```python
-from deepdrr import geo, Volume, MobileCArm, Projector
+from deepdrr import geo, Volume, MobileCArm
+from deepdrr.projector import Projector # separate import for CUDA init
 import matplotlib.pyplot as plt
 
 volume = Volume.from_nifti('/path/to/ct_image.nii.gz')
@@ -147,6 +148,41 @@ This capability has not been tested in version 1.0. For tool insertion, we recom
 2. The density of the tool needs to be provided via hard coding in the file 'load_dicom_tool.py' (line 127). The pose of the tool/implant with respect to the CT volume requires manual setup. We provide one example origin setting at line 23-24.
 3. The tool/implant will supersede the anatomy defined by the CT volume intensities. To this end, we sample the CT materials and densities at the location of the tool in the tool volume, and subtract them from the anatomy forward projections in detector domain (to enable different resolutions of CT and tool volume). Further information can be found in the IJCARS article.
 
+### Using DeepDRR Simultaneously with PyTorch
+
+Some issues may arise when using DeepDRR at the same time as PyTorch due to conflicts between pycuda's CUDA initialization and PyTorch CUDA initialization. The best workaround we know of is to first initialize the PyCUDA context (by importing `deepdrr.projector`) and then run your model on a dummy batch before creating a `Projector` object. For mysterious reasons (likely involving overlapping GPU resources and the retrograde of Mercury), this seems to work.
+
+```Python
+import torch
+from torch import nn
+from torchvision import models
+
+import deepdrr
+from deepdrr.projector import Projector # initializes PyCUDA
+
+# Before creating a Projector, run backprop to initialize PyTorch
+criterion = nn.CrossEntropyLoss()
+model = models.resnet50() # Your model here
+model.cuda()
+optimizer = optim.SGD(model.parameters(), lr=0.001, momentum=0.9)
+optimizer.zero_grad()
+x = torch.ones((32, 3, 224, 224), dtype=torch.float32).cuda() # Your image size
+y = torch.ones(32, dtype=torch.int64).cuda()
+y_pred = model(x)
+loss = criterion(y_pred, y)
+loss.backward()
+optimizer.step()
+log.info(f"Ran dummy batch to initialize torch.")
+
+volume = ...
+carm = ...
+with Projector(volume, carm=carm):
+  image = projector()
+  image = image.unsqueeze(0) # add batch dim
+  y_pred = model(image)
+  ...
+```
+
 ## Reference
 
 We hope this proves useful for medical imaging research. If you use our work, we would kindly ask you to reference our work.

deepdrr/annotations/line_annotation.py

@@ -1,13 +1,13 @@
 from __future__ import annotations
 
 import logging
-from typing import Optional
+from typing import List, Optional
 from pathlib import Path
 import numpy as np
 import json
 import pyvista as pv
 
-from .. import geo
+from .. import geo, utils
 from ..vol import Volume, AnyVolume
 
 logger = logging.getLogger(__name__)
@@ -65,6 +65,110 @@ def from_markup(cls, path: str, volume: AnyVolume) -> LineAnnotation:
 
         return cls(*points, volume)
 
+    def save(self, path: str, color: List[float] = [0.5, 0.5, 0.5]):
+        """Save the Line annotation to a mrk.json file, which can be opened by 3D Slicer.
+
+        Args:
+            path (str): Output path to the file.
+            color (List[int], optional): The color of the saved annotation.
+        """
+        path = Path(path).expanduser()
+
+        markup = {
+            "@schema": "https://raw.githubusercontent.com/slicer/slicer/master/Modules/Loadable/Markups/Resources/Schema/markups-schema-v1.0.0.json#",
+            "markups": [
+                {
+                    "type": "Line",
+                    "coordinateSystem": self.volume.anatomical_coordinate_system,
+                    "locked": True,
+                    "labelFormat": "%N-%d",
+                    "controlPoints": [
+                        {
+                            "id": "1",
+                            "label": "entry",
+                            "description": "",
+                            "associatedNodeID": "",
+                            "position": utils.jsonable(self.startpoint),
+                            "orientation": [
+                                -1.0,
+                                -0.0,
+                                -0.0,
+                                -0.0,
+                                -1.0,
+                                -0.0,
+                                0.0,
+                                0.0,
+                                1.0,
+                            ],
+                            "selected": True,
+                            "locked": False,
+                            "visibility": True,
+                            "positionStatus": "defined",
+                        },
+                        {
+                            "id": "2",
+                            "label": "exit",
+                            "description": "",
+                            "associatedNodeID": "",
+                            "position": utils.jsonable(self.endpoint),
+                            "orientation": [
+                                -1.0,
+                                -0.0,
+                                -0.0,
+                                -0.0,
+                                -1.0,
+                                -0.0,
+                                0.0,
+                                0.0,
+                                1.0,
+                            ],
+                            "selected": True,
+                            "locked": False,
+                            "visibility": True,
+                            "positionStatus": "defined",
+                        },
+                    ],
+                    "measurements": [
+                        {
+                            "name": "length",
+                            "enabled": True,
+                            "value": 124.90054351814699,
+                            "printFormat": "%-#4.4gmm",
+                        }
+                    ],
+                    "display": {
+                        "visibility": True,
+                        "opacity": 1.0,
+                        "color": color,
+                        "selectedColor": [1.0, 0.5000076295109484, 0.5000076295109484],
+                        "activeColor": [0.4, 1.0, 0.0],
+                        "propertiesLabelVisibility": True,
+                        "pointLabelsVisibility": True,
+                        "textScale": 3.0,
+                        "glyphType": "Sphere3D",
+                        "glyphScale": 5.800000000000001,
+                        "glyphSize": 5.0,
+                        "useGlyphScale": True,
+                        "sliceProjection": False,
+                        "sliceProjectionUseFiducialColor": True,
+                        "sliceProjectionOutlinedBehindSlicePlane": False,
+                        "sliceProjectionColor": [1.0, 1.0, 1.0],
+                        "sliceProjectionOpacity": 0.6,
+                        "lineThickness": 0.2,
+                        "lineColorFadingStart": 1.0,
+                        "lineColorFadingEnd": 10.0,
+                        "lineColorFadingSaturation": 1.0,
+                        "lineColorFadingHueOffset": 0.0,
+                        "handlesInteractive": False,
+                        "snapMode": "toVisibleSurface",
+                    },
+                }
+            ],
+        }
+
+        with open(path, "w") as file:
+            json.dump(markup, file)
+
     @property
     def startpoint_in_world(self) -> geo.Point:
         return self.volume.world_from_anatomical @ self.startpoint

deepdrr/device.py

@@ -1,4 +1,4 @@
-from typing import Optional, Tuple, Union, List
+from typing import Any, Dict, Optional, Tuple, Union, List
 
 import logging
 import numpy as np
@@ -96,14 +96,14 @@ def __init__(
         isocenter: geo.Point3D = [0, 0, 0],
         alpha: float = 0,
         beta: float = 0,
+        degrees: bool = True,
         horizontal_movement: float = 200,  # width of window in X and Y planes.
         vertical_travel: float = 430,  # width of window in Z plane.
         min_alpha: float = -40,
         max_alpha: float = 110,
         # note that this would collide with the patient. Suggested to limit to +/- 45
         min_beta: float = -225,
         max_beta: float = 225,
-        degrees: bool = True,
         source_to_detector_distance: float = 1020,
         # vertical component of the source point offset from the isocenter of rotation, in -Z. Previously called `isocenter_distance`
         source_to_isocenter_vertical_distance: float = 530,
@@ -117,7 +117,7 @@ def __init__(
         sensor_width: int = 1536,
         pixel_size: float = 0.194,
         rotate_camera_left: bool = True,  # make it so that down in the image corresponds to -x, so that patient images appear as expected.
-        enforce_isocenter_bounds: bool = True,
+        enforce_isocenter_bounds: bool = False,  # Allow the isocenter to travel arbitrarily far from the device origin
     ) -> None:
         """A simulated C-arm imaging device with orbital movement (alpha), angulation (beta) and 3D translation.
 
@@ -201,6 +201,34 @@ def __str__(self):
             f"alpha={np.degrees(self.alpha)}, beta={np.degrees(self.beta)}, degrees=True)"
         )
 
+    def to_config(self) -> Dict[str, Any]:
+        """Get a json-safe dictionary that can be used to initialize the C-arm in its current pose."""
+        return utils.jsonable(
+            dict(
+                world_from_device=self.world_from_device,
+                isocenter=self.isocenter,
+                alpha=self.alpha,
+                beta=self.beta,
+                degrees=False,
+                horizontal_movement=self.horizontal_movement,
+                vertical_travel=self.vertical_travel,
+                min_alpha=self.min_alpha,
+                max_alpha=self.max_alpha,
+                min_beta=self.min_beta,
+                max_beta=self.max_beta,
+                source_to_detector_distance=self.source_to_detector_distance,
+                source_to_isocenter_vertical_distance=self.source_to_isocenter_vertical_distance,
+                source_to_isocenter_horizontal_offset=self.source_to_isocenter_horizontal_offset,
+                immersion_depth=self.immersion_depth,
+                free_space=self.free_space,
+                sensor_height=self.sensor_height,
+                sensor_width=self.sensor_width,
+                pixel_size=self.pixel_size,
+                rotate_camera_left=self.rotate_camera_left,
+                enforce_isocenter_bounds=self.enforce_isocenter_bounds,
+            )
+        )
+
     @property
     def max_isocenter(self) -> np.ndarray:
         return (
@@ -257,10 +285,10 @@ def camera3d_from_device(self) -> geo.FrameTransform:
     @property
     def camera3d_from_world(self) -> geo.FrameTransform:
         """Rigid transformation of the C-arm camera pose."""
-        return self.camera3d_from_device @ self.device_from_world
+        return self.get_camera3d_from_world()
 
     def get_camera3d_from_world(self) -> geo.FrameTransform:
-        return self.camera3d_from_world
+        return self.camera3d_from_device @ self.device_from_world
 
     def get_camera_projection(self) -> geo.CameraProjection:
         return geo.CameraProjection(

deepdrr/geo/camera_projection.py

@@ -1,7 +1,7 @@
 from typing import Union, Optional, Any, TYPE_CHECKING
 import numpy as np
 
-from .core import Transform, FrameTransform, point, Point3D
+from .core import Transform, FrameTransform, point, Point3D, get_data
 from .camera_intrinsic_transform import CameraIntrinsicTransform
 from ..vol import AnyVolume
 
@@ -66,17 +66,25 @@ def extrinsic(self) -> FrameTransform:
         return self.camera3d_from_world
 
     @property
-    def index_from_world(self) -> FrameTransform:
+    def index_from_world(self) -> Transform:
         proj = np.concatenate([np.eye(3), np.zeros((3, 1))], axis=1)
         camera2d_from_camera3d = Transform(proj, _inv=proj.T)
         return (
             self.index_from_camera2d @ camera2d_from_camera3d @ self.camera3d_from_world
         )
 
     @property
-    def world_from_index(self) -> FrameTransform:
+    def world_from_index(self) -> Transform:
         return self.index_from_world.inv
 
+    @property
+    def world_from_index_on_image_plane(self) -> FrameTransform:
+        """Get the transform to points in world on the image (detector) plane from image indices."""
+        proj = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 0], [0, 0, 1]])
+        proj = Transform(proj, _inv=proj.T)
+        index_from_world_3d = proj @ self.index_from_world
+        return FrameTransform(data=get_data(index_from_world_3d.inv))
+
     @property
     def sensor_width(self) -> int:
         return self.intrinsic.sensor_width

deepdrr/geo/core.py

@@ -92,6 +92,10 @@ def to_array(self, is_point):
         """
         pass
 
+    def tolist(self) -> List:
+        """Get a json-save list with the data in this object."""
+        return self.data.tolist()
+
     def __array__(self, dtype=None):
         return self.to_array()