added nrrd capabilities

.gitignore

@@ -9,3 +9,4 @@ docs/build
 .vscode
 **flycheck*.py
 /.DS_Store
+output

deepdrr/device.py

@@ -83,36 +83,7 @@ def pose_vector_angles(pose: geo.Vector3D) -> Tuple[float, float]:
 
 
 class MobileCArm(object):
-    """A simulated C-arm imaging device with orbital movement (alpha), angulation (beta) and 3D translation.
-
-    A MobileCArm has its own `device` frame, which moves independent of the `arm` frame and the
-    `camera` frame, which actually defines the projection.
-
-    The geometry and default values, are based on the Cios Fusion device, specifically this_
-    document, with some exceptions. Rather than incorporating a swivel, this device allows for
-    simple horizontal translations. Additionally, by default the principle ray is not offset from
-    the isocenter.
-
-    .. _this: https://www.lomisa.com/app/download/10978681598/ARCO_EN_C_SIEMENS_CIOS_FUSION.pdf?t=1490962065
-
-    Additionall, the orbital movement and angulation are reversed from the convention used in `Kausch
-    et al`_, Figure 2.
-
-    .. _Kausch et al: https://pubmed.ncbi.nlm.nih.gov/32533315/
-
-    All length units are in millimeters.
-
-    Args:
-        isocenter (geo.Point3D): the initial isocenter of in the device frame. This is the point
-            about which rotations are performed.
-        isocenter_distance (float): the distance from the X-ray source to the isocenter of the CAarm. (The center of rotation).
-        alpha (float): initial LAO/RAO angle of the C-Arm. alpha > 0 is in the RAO direction. This is the angle along arm of the C-arm.
-        beta (float): initial CRA/CAU angulation of the C-arm. beta > 0 is in the CAU direction.
-        degrees (bool, optional): Whether given angles are in degrees. Defaults to False.
-        world_from_device: (Optional[geo.FrameTransform], optional): Transform that defines the device coordinate space in world coordinates. None is the identity transform. Defaults to None.
-        camera_intrinsics: (Optional[Union[geo.CameraIntrinsicTransform, dict]], optional): either a CameraIntrinsicTransform instance or kwargs for CameraIntrinsicTransform.from_sizes
-
-    """
+
 
     # basic parameters which can be safely set by user, but move_by() and reposition() are recommended.
     isocenter: geo.Point3D  # the isocenter point in the device frame
@@ -142,6 +113,36 @@ def __init__(
         sensor_width: int = 1536,
         pixel_size: float = 0.194,
     ) -> None:
+        """A simulated C-arm imaging device with orbital movement (alpha), angulation (beta) and 3D translation.
+
+        A MobileCArm has its own `device` frame, which moves independent of the `arm` frame and the
+        `camera` frame, which actually defines the projection.
+
+        The geometry and default values, are based on the Cios Fusion device, specifically this_
+        document, with some exceptions. Rather than incorporating a swivel, this device allows for
+        simple horizontal translations. Additionally, by default the principle ray is not offset from
+        the isocenter.
+
+        .. _this: https://www.lomisa.com/app/download/10978681598/ARCO_EN_C_SIEMENS_CIOS_FUSION.pdf?t=1490962065
+
+        Additionall, the orbital movement and angulation are reversed from the convention used in `Kausch
+        et al`_, Figure 2.
+
+        .. _Kausch et al: https://pubmed.ncbi.nlm.nih.gov/32533315/
+
+        All length units are in millimeters.
+
+        Args:
+            isocenter (geo.Point3D): the initial isocenter of in the device frame. This is the point
+                about which rotations are performed.
+            isocenter_distance (float): the distance from the X-ray source to the isocenter of the CAarm. (The center of rotation).
+            alpha (float): initial LAO/RAO angle of the C-Arm. alpha > 0 is in the RAO direction. This is the angle along arm of the C-arm.
+            beta (float): initial CRA/CAU angulation of the C-arm. beta > 0 is in the CAU direction.
+            degrees (bool, optional): Whether given angles are in degrees. Defaults to False.
+            world_from_device: (Optional[geo.FrameTransform], optional): Transform that defines the device coordinate space in world coordinates. None is the identity transform. Defaults to None.
+            camera_intrinsics: (Optional[Union[geo.CameraIntrinsicTransform, dict]], optional): either a CameraIntrinsicTransform instance or kwargs for CameraIntrinsicTransform.from_sizes
+
+        """
         self.isocenter = geo.point(isocenter)
         self.alpha = utils.radians(alpha, degrees=degrees)
         self.beta = utils.radians(beta, degrees=degrees)

deepdrr/utils/testing.py

@@ -1,2 +1,43 @@
 """Util functions for automated tests."""
 
+import logging
+from pathlib import Path
+from torchvision.datasets.utils import download_url
+
+logger = logging.getLogger(__name__)
+
+_sampledata = {
+    'CT-chest': 'http://www.slicer.org/w/img_auth.php/3/31/CT-chest.nrrd',
+}
+
+def download_sampledata(name: str = 'CT-chest', root: str = "~/datasets/") -> Path:
+    """Download the sample volumes for testing, if they are not already present.
+
+    Args:
+        name (str, optional): The name of the volume, used as a key to the public domain downloadable data. Defaults to 'CT-chest'.
+        root (str, optional): The root where datasets are kept. A new dataset directory will be created in `root`, called `DeepDRR_SampleData`. Defaults to "~/datasets/".
+
+    Returns:
+        Path: The path to the downloaded file.
+    """
+    if name not in _sampledata:
+        logger.error(f"unrecognized sample data name: {name}. Options are:\n{list(_sampledata.keys())}")
+
+    root = Path(root).expanduser()
+    if not root.exists():
+        root.mkdir()
+
+    dataset_dir = root / 'DeepDRR_SampleData'
+    if not dataset_dir.exists():
+        dataset_dir.mkdir()
+
+    fpath = dataset_dir / f"{name}.nrrd"
+    if fpath.exists():
+        return fpath
+
+    logger.info(f"Downloading sample volume to {fpath}")
+
+    download_url(_sampledata[name], root=dataset_dir, filename=fpath.name)
+
+    logger.info("Done.")
+    return fpath

deepdrr/vol.py

@@ -10,6 +10,7 @@
 from pathlib import Path
 import nibabel as nib
 from pydicom.filereader import dcmread
+import nrrd
 
 from . import load_dicom
 from . import geo
@@ -47,11 +48,11 @@ def __init__(
         """
         self.data = np.array(data).astype(np.float32)
         self.materials = self._format_materials(materials)
-        self.anatomical_from_ijk = anatomical_from_ijk
+        self.anatomical_from_ijk = geo.frame_transform(anatomical_from_ijk)
         self.world_from_anatomical = (
             geo.FrameTransform.identity(3)
             if world_from_anatomical is None
-            else world_from_anatomical
+            else geo.frame_transform(world_from_anatomical)
         )
 
     @classmethod
@@ -410,7 +411,7 @@ def from_dicom(
         # constructing the volume
         return cls(data, materials, lps_from_ijk, world_from_anatomical,)
 
-    def to_dicom(self, path: Union[str, Path]):
+    def to_dicom(self, path: str):
         """Write the volume to a DICOM file.
 
         Args:
@@ -420,6 +421,47 @@ def to_dicom(self, path: Union[str, Path]):
 
         raise NotImplementedError("save volume to dicom file")
 
+    @classmethod
+    def from_nrrd(
+        cls, 
+        path: str, 
+        world_from_anatomical: Optional[geo.FrameTransform] = None, 
+        use_thresholding: bool = True,
+        use_cached: bool = True,
+        cache_dir: Optional[Path] = None,
+    ):
+        """Load a volume from a nrrd file.
+
+        Args:
+            path (str): path to the file.
+            use_thresholding (bool, optional): segment the materials using thresholding (faster but less accurate). Defaults to True.
+            world_from_anatomical (Optional[geo.FrameTransform], optional): position the volume in world space. If None, uses identity. Defaults to None.
+            use_cached (bool, optional): Use a cached segmentation if available. Defaults to True.
+            cache_dir (Optional[Path], optional): Where to load/save the cached segmentation. If None, use the parent dir of `path`. Defaults to None.
+
+        Returns:
+            Volume: A volume formed from the NRRD.
+        """
+        hu_values, header = nrrd.read(path)
+        ijk_from_anatomical = np.concatenate(
+            [
+                header['space directions'],
+                header['space origin'].reshape(-1, 1),
+            ],
+            axis=1,
+        )
+        ijk_from_anatomical = np.concatenate([ijk_from_anatomical, [[0, 0, 0, 1]]], axis=0)
+        anatomical_from_ijk = np.linalg.inv(ijk_from_anatomical)
+        data = cls._convert_hounsfield_to_density(hu_values)
+        materials = cls.segment_materials(
+            hu_values,
+            use_thresholding=use_thresholding,
+            use_cached=use_cached,
+            cache_dir=cache_dir,
+        )
+
+        return cls(data, materials, anatomical_from_ijk, world_from_anatomical)
+
     @property
     def origin(self) -> geo.Point3D:
         """The origin of the volume in anatomical space."""

docs/requirements.txt

@@ -6,6 +6,7 @@ scikit-image
 pydicom
 scipy
 pyvista 
+pynrrd
 # pycuda # docs build machines don't have GPUs
 
 flake8

requirements.txt

@@ -8,3 +8,4 @@ pydicom
 scipy
 pyvista
 pycuda
+pynrrd

setup.py

@@ -5,7 +5,7 @@
 
 setuptools.setup(
     name="deepdrr",
-    version="1.1.0a1",
+    version="1.1.0a2",
     author="Benjamin D. Killeen",
     author_email="[email protected]",
     description="A Catalyst for Machine Learning in Fluoroscopy-guided Procedures.",
@@ -23,6 +23,7 @@
         "pyvista",
         "scipy",
         "pyvista",
+        "pynrrd"
     ],
     extras_require={"gpu": ["pycuda"]},
     include_package_data=True,

tests/test.py

@@ -1,11 +1,37 @@
 #!/usr/bin/env python3
 
-from pydicom.data import get_testdata_file
+import numpy as np
+from pathlib import Path
 import deepdrr
 from deepdrr import geo
+from deepdrr.utils import testing
+from PIL import Image
 
 
-def test_simple():
-    file_path = get_testdata_file("CT_small.dcm")
-    volume = deepdrr.Volume.from_dicom(file_path)
-    print(volume)
+class TestSingleVolume:
+    output_dir = Path.cwd() / 'output'
+    output_dir.mkdir(exist_ok=True)
+
+    def test_simple(self):
+        file_path = testing.download_sampledata("CT-chest")
+        volume = deepdrr.Volume.from_nrrd(file_path)
+        carm = deepdrr.MobileCArm(isocenter=volume.center_in_world)
+        with deepdrr.Projector(
+            volume=volume,
+            carm=carm,
+            step=0.1,  # stepsize along projection ray, measured in voxels
+            mode="linear",
+            max_block_index=200,
+            spectrum="90KV_AL40",
+            photon_count=100000,
+            add_scatter=False,
+            threads=8,
+            neglog=True,
+        ) as projector:
+            image = projector.project()
+
+        image = (image * 255).astype(np.uint8)
+        Image.fromarray(image).save(self.output_dir / 'test_simple.png')
+
+if __name__ == "__main__":
+    TestSingleVolume().test_simple()