diff --git a/diffdrr/__init__.py b/diffdrr/__init__.py index 98a433b3..3dd3d2d5 100644 --- a/diffdrr/__init__.py +++ b/diffdrr/__init__.py @@ -1 +1 @@ -__version__ = "0.4.5" +__version__ = "0.4.6" diff --git a/diffdrr/data.py b/diffdrr/data.py index c7c5e481..b9464291 100644 --- a/diffdrr/data.py +++ b/diffdrr/data.py @@ -85,35 +85,37 @@ def read( # Frame-of-reference change if orientation == "AP": # Rotates the C-arm about the x-axis by 90 degrees - # Rotates the C-arm about the z-axis by -90 degrees reorient = torch.tensor( [ - [0.0, 1.0, 0.0, 0.0], - [0.0, 0.0, -1.0, 0.0], - [-1.0, 0.0, 0.0, 0.0], - [0.0, 0.0, 0.0, 1.0], - ] + [1, 0, 0, 0], + [0, 0, -1, 0], + [0, 1, 0, 0], + [0, 0, 0, 1], + ], + dtype=torch.float32, ) elif orientation == "PA": # Rotates the C-arm about the x-axis by 90 degrees - # Rotates the C-arm about the z-axis by 90 degrees + # Reverses the direction of the y-axis reorient = torch.tensor( [ - [0.0, 1.0, 0.0, 0.0], - [0.0, 0.0, 1.0, 0.0], - [-1.0, 0.0, 0.0, 0.0], - [0.0, 0.0, 0.0, 1.0], - ] + [1, 0, 0, 0], + [0, 0, 1, 0], + [0, 1, 0, 0], + [0, 0, 0, 1], + ], + dtype=torch.float32, ) elif orientation is None: # Identity transform reorient = torch.tensor( [ - [1.0, 0.0, 0.0, 0.0], - [0.0, 1.0, 0.0, 0.0], - [0.0, 0.0, 1.0, 0.0], - [0.0, 0.0, 0.0, 1.0], - ] + [1, 0, 0, 0], + [0, 1, 0, 0], + [0, 0, 1, 0], + [0, 0, 0, 1], + ], + dtype=torch.float32, ) else: raise ValueError(f"Unrecognized orientation {orientation}") @@ -122,6 +124,7 @@ def read( subject = Subject( volume=volume, mask=mask, + orientation=orientation, reorient=reorient, density=density, fiducials=fiducials, @@ -161,9 +164,13 @@ def read( dim=0, ) - subject.volume.data = subject.volume.data * mask - subject.mask.data = subject.mask.data * mask - subject.density.data = subject.density.data * mask + # Mask all volumes, unless error, then just mask the density + try: + subject.volume.data = subject.volume.data * mask + subject.mask.data = subject.mask.data * mask + subject.density.data = subject.density.data * mask + except: + subject.density.data = subject.density.data * mask return subject diff --git a/diffdrr/detector.py b/diffdrr/detector.py index e9afeb47..1bed607d 100644 --- a/diffdrr/detector.py +++ b/diffdrr/detector.py @@ -51,8 +51,8 @@ def __init__( "_calibration", torch.tensor( [ - [dely, 0, 0, -y0], - [0, delx, 0, -x0], + [delx, 0, 0, x0], + [0, dely, 0, y0], [0, 0, sdd, 0], [0, 0, 0, 1], ] @@ -65,19 +65,19 @@ def sdd(self): @property def delx(self): - return self._calibration[1, 1].item() + return self._calibration[0, 0].item() @property def dely(self): - return self._calibration[0, 0].item() + return self._calibration[1, 1].item() @property def x0(self): - return -self._calibration[1, -1].item() + return -self._calibration[0, -1].item() @property def y0(self): - return -self._calibration[0, -1].item() + return -self._calibration[1, -1].item() @property def reorient(self): @@ -107,7 +107,7 @@ def _initialize_carm(self: Detector): center = torch.tensor([[0.0, 0.0, 1.0]], device=device) # Use the standard basis for the detector plane - basis = torch.tensor([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], device=device) + basis = torch.tensor([[0.0, 1.0, 0.0], [1.0, 0.0, 0.0]], device=device) # Construct the detector plane with different offsets for even or odd heights # These ensure that the detector plane is centered around (0, 0, 1) @@ -117,8 +117,12 @@ def _initialize_carm(self: Detector): # Construct equally spaced points along the basis vectors t = torch.arange(-self.height // 2, self.height // 2, device=device) + h_off s = torch.arange(-self.width // 2, self.width // 2, device=device) + w_off - if self.reverse_x_axis: + + t = -t + s = -s + if not self.reverse_x_axis: s = -s + coefs = torch.cartesian_prod(t, s).reshape(-1, 2) target = torch.einsum("cd,nc->nd", basis, coefs) target += center diff --git a/diffdrr/drr.py b/diffdrr/drr.py index 1eced0f9..0ef37d36 100644 --- a/diffdrr/drr.py +++ b/diffdrr/drr.py @@ -230,8 +230,8 @@ def set_intrinsics_( width if width is not None else self.detector.width, delx if delx is not None else self.detector.delx, dely if dely is not None else self.detector.dely, - x0 if x0 is not None else self.detector.x0, - y0 if y0 is not None else self.detector.y0, + x0 if x0 is not None else -self.detector.x0, + y0 if y0 is not None else -self.detector.y0, self.subject.reorient, n_subsample if n_subsample is not None else self.detector.n_subsample, reverse_x_axis if reverse_x_axis is not None else self.detector.reverse_x_axis, @@ -256,14 +256,21 @@ def perspective_projection( pts: torch.Tensor, ): """Project points in world coordinates (3D) onto the pixel plane (2D).""" + # Poses in DiffDRR are world2camera, but perspective transforms use camera2world, so invert extrinsic = (self.detector.reorient.compose(pose)).inverse() x = extrinsic(pts) + + # Project onto the detector plane x = torch.einsum("ij, bnj -> bni", self.detector.intrinsic, x) z = x[..., -1].unsqueeze(-1).clone() x = x / z + + # Move origin to upper-left corner + x[..., 1] = self.detector.height - x[..., 1] if self.detector.reverse_x_axis: - x[..., 1] = self.detector.width - x[..., 1] - return x[..., :2].flip(-1) + x[..., 0] = self.detector.width - x[..., 0] + + return x[..., :2] # %% ../notebooks/api/00_drr.ipynb 14 from torch.nn.functional import pad @@ -276,7 +283,7 @@ def inverse_projection( pts: torch.Tensor, ): """Backproject points in pixel plane (2D) onto the image plane in world coordinates (3D).""" - pts = pts.flip(-1) + # pts = pts.flip(-1) if self.detector.reverse_x_axis: pts[..., 1] = self.detector.width - pts[..., 1] x = self.detector.sdd * torch.einsum( diff --git a/diffdrr/utils.py b/diffdrr/utils.py index 9caa4c87..0481accc 100644 --- a/diffdrr/utils.py +++ b/diffdrr/utils.py @@ -54,6 +54,7 @@ def resample( # %% ../notebooks/api/07_utils.ipynb 6 from kornia.geometry.camera.pinhole import PinholeCamera as KorniaPinholeCamera +from torchio import Subject from diffdrr.detector import Detector @@ -66,9 +67,11 @@ def __init__( height: torch.Tensor, width: torch.Tensor, detector: Detector, + subject: Subject, ): super().__init__(intrinsics, extrinsics, height, width) - self.f = detector.sdd + multiplier = -1 if subject.orientation == "PA" else 1 + self.sdd = multiplier * detector.sdd self.delx = detector.delx self.dely = detector.dely self.x0 = detector.x0 @@ -94,9 +97,9 @@ def pose(self): from kornia.geometry.calibration import solve_pnp_dlt +from .detector import make_intrinsic_matrix from .drr import DRR from .pose import RigidTransform -from .detector import make_intrinsic_matrix def get_pinhole_camera( @@ -107,14 +110,12 @@ def get_pinhole_camera( pose = deepcopy(pose).to(device="cpu", dtype=dtype) # Make the intrinsic matrix (in pixels) - fx = drr.detector.sdd / drr.detector.delx - fy = drr.detector.sdd / drr.detector.dely + multiplier = -1 if drr.subject.orientation == "PA" else 1 + fx = multiplier * drr.detector.sdd / drr.detector.delx + fy = multiplier * drr.detector.sdd / drr.detector.dely u0 = drr.detector.x0 / drr.detector.delx + drr.detector.width / 2 v0 = drr.detector.y0 / drr.detector.dely + drr.detector.height / 2 - intrinsics = torch.eye(4)[None] - intrinsics[0, :3, :3] = make_intrinsic_matrix(drr.detector) - - torch.tensor( + intrinsics = torch.tensor( [ [ [fx, 0.0, u0, 0.0], @@ -156,6 +157,7 @@ def get_pinhole_camera( torch.tensor([drr.detector.height]), torch.tensor([drr.detector.width]), drr.detector, + drr.subject, ) return camera diff --git a/notebooks/api/00_drr.ipynb b/notebooks/api/00_drr.ipynb index 0eb2d1c4..d43dcb91 100644 --- a/notebooks/api/00_drr.ipynb +++ b/notebooks/api/00_drr.ipynb @@ -357,8 +357,8 @@ " width if width is not None else self.detector.width,\n", " delx if delx is not None else self.detector.delx,\n", " dely if dely is not None else self.detector.dely,\n", - " x0 if x0 is not None else self.detector.x0,\n", - " y0 if y0 is not None else self.detector.y0,\n", + " x0 if x0 is not None else -self.detector.x0,\n", + " y0 if y0 is not None else -self.detector.y0,\n", " self.subject.reorient,\n", " n_subsample if n_subsample is not None else self.detector.n_subsample,\n", " reverse_x_axis if reverse_x_axis is not None else self.detector.reverse_x_axis,\n", @@ -399,14 +399,21 @@ " pts: torch.Tensor,\n", "):\n", " \"\"\"Project points in world coordinates (3D) onto the pixel plane (2D).\"\"\"\n", + " # Poses in DiffDRR are world2camera, but perspective transforms use camera2world, so invert\n", " extrinsic = (self.detector.reorient.compose(pose)).inverse()\n", " x = extrinsic(pts)\n", + "\n", + " # Project onto the detector plane\n", " x = torch.einsum(\"ij, bnj -> bni\", self.detector.intrinsic, x)\n", " z = x[..., -1].unsqueeze(-1).clone()\n", " x = x / z\n", + "\n", + " # Move origin to upper-left corner\n", + " x[..., 1] = self.detector.height - x[..., 1]\n", " if self.detector.reverse_x_axis:\n", - " x[..., 1] = self.detector.width - x[..., 1]\n", - " return x[..., :2].flip(-1)" + " x[..., 0] = self.detector.width - x[..., 0]\n", + " \n", + " return x[..., :2]" ] }, { @@ -427,7 +434,7 @@ " pts: torch.Tensor,\n", "):\n", " \"\"\"Backproject points in pixel plane (2D) onto the image plane in world coordinates (3D).\"\"\"\n", - " pts = pts.flip(-1)\n", + " # pts = pts.flip(-1)\n", " if self.detector.reverse_x_axis:\n", " pts[..., 1] = self.detector.width - pts[..., 1]\n", " x = self.detector.sdd * torch.einsum(\n", diff --git a/notebooks/api/02_detector.ipynb b/notebooks/api/02_detector.ipynb index 86417b51..732bcb7a 100644 --- a/notebooks/api/02_detector.ipynb +++ b/notebooks/api/02_detector.ipynb @@ -106,8 +106,8 @@ " \"_calibration\",\n", " torch.tensor(\n", " [\n", - " [dely, 0, 0, -y0],\n", - " [0, delx, 0, -x0],\n", + " [delx, 0, 0, x0],\n", + " [0, dely, 0, y0],\n", " [0, 0, sdd, 0],\n", " [0, 0, 0, 1],\n", " ]\n", @@ -120,19 +120,19 @@ "\n", " @property\n", " def delx(self):\n", - " return self._calibration[1, 1].item()\n", + " return self._calibration[0, 0].item()\n", "\n", " @property\n", " def dely(self):\n", - " return self._calibration[0, 0].item()\n", + " return self._calibration[1, 1].item()\n", "\n", " @property\n", " def x0(self):\n", - " return -self._calibration[1, -1].item()\n", + " return -self._calibration[0, -1].item()\n", "\n", " @property\n", " def y0(self):\n", - " return -self._calibration[0, -1].item()\n", + " return -self._calibration[1, -1].item()\n", "\n", " @property\n", " def reorient(self):\n", @@ -170,7 +170,7 @@ " center = torch.tensor([[0.0, 0.0, 1.0]], device=device)\n", "\n", " # Use the standard basis for the detector plane\n", - " basis = torch.tensor([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], device=device)\n", + " basis = torch.tensor([[0.0, 1.0, 0.0], [1.0, 0.0, 0.0]], device=device)\n", "\n", " # Construct the detector plane with different offsets for even or odd heights\n", " # These ensure that the detector plane is centered around (0, 0, 1)\n", @@ -180,8 +180,12 @@ " # Construct equally spaced points along the basis vectors\n", " t = torch.arange(-self.height // 2, self.height // 2, device=device) + h_off\n", " s = torch.arange(-self.width // 2, self.width // 2, device=device) + w_off\n", - " if self.reverse_x_axis:\n", + "\n", + " t = -t\n", + " s = -s\n", + " if not self.reverse_x_axis:\n", " s = -s\n", + "\n", " coefs = torch.cartesian_prod(t, s).reshape(-1, 2)\n", " target = torch.einsum(\"cd,nc->nd\", basis, coefs)\n", " target += center\n", diff --git a/notebooks/api/03_data.ipynb b/notebooks/api/03_data.ipynb index c1b3b943..b97e67d2 100644 --- a/notebooks/api/03_data.ipynb +++ b/notebooks/api/03_data.ipynb @@ -156,35 +156,37 @@ " # Frame-of-reference change\n", " if orientation == \"AP\":\n", " # Rotates the C-arm about the x-axis by 90 degrees\n", - " # Rotates the C-arm about the z-axis by -90 degrees\n", " reorient = torch.tensor(\n", " [\n", - " [0.0, 1.0, 0.0, 0.0],\n", - " [0.0, 0.0, -1.0, 0.0],\n", - " [-1.0, 0.0, 0.0, 0.0],\n", - " [0.0, 0.0, 0.0, 1.0],\n", - " ]\n", + " [1, 0, 0, 0],\n", + " [0, 0, -1, 0],\n", + " [0, 1, 0, 0],\n", + " [0, 0, 0, 1],\n", + " ],\n", + " dtype=torch.float32,\n", " )\n", " elif orientation == \"PA\":\n", - " # Rotates the C-arm about the x-axis by 90 degrees\n", - " # Rotates the C-arm about the z-axis by 90 degrees\n", + " # Rotates the C-arm about the x-axis by 90 degrees \n", + " # Reverses the direction of the y-axis\n", " reorient = torch.tensor(\n", " [\n", - " [0.0, 1.0, 0.0, 0.0],\n", - " [0.0, 0.0, 1.0, 0.0],\n", - " [-1.0, 0.0, 0.0, 0.0],\n", - " [0.0, 0.0, 0.0, 1.0],\n", - " ]\n", + " [1, 0, 0, 0],\n", + " [0, 0, 1, 0],\n", + " [0, 1, 0, 0],\n", + " [0, 0, 0, 1],\n", + " ],\n", + " dtype=torch.float32,\n", " )\n", " elif orientation is None:\n", " # Identity transform\n", " reorient = torch.tensor(\n", " [\n", - " [1.0, 0.0, 0.0, 0.0],\n", - " [0.0, 1.0, 0.0, 0.0],\n", - " [0.0, 0.0, 1.0, 0.0],\n", - " [0.0, 0.0, 0.0, 1.0],\n", - " ]\n", + " [1, 0, 0, 0],\n", + " [0, 1, 0, 0],\n", + " [0, 0, 1, 0],\n", + " [0, 0, 0, 1],\n", + " ],\n", + " dtype=torch.float32,\n", " )\n", " else:\n", " raise ValueError(f\"Unrecognized orientation {orientation}\")\n", @@ -193,6 +195,7 @@ " subject = Subject(\n", " volume=volume,\n", " mask=mask,\n", + " orientation=orientation,\n", " reorient=reorient,\n", " density=density,\n", " fiducials=fiducials,\n", @@ -232,9 +235,13 @@ " dim=0,\n", " )\n", "\n", - " subject.volume.data = subject.volume.data * mask\n", - " subject.mask.data = subject.mask.data * mask\n", - " subject.density.data = subject.density.data * mask\n", + " # Mask all volumes, unless error, then just mask the density\n", + " try:\n", + " subject.volume.data = subject.volume.data * mask\n", + " subject.mask.data = subject.mask.data * mask\n", + " subject.density.data = subject.density.data * mask\n", + " except:\n", + " subject.density.data = subject.density.data * mask\n", "\n", " return subject" ] diff --git a/notebooks/api/07_utils.ipynb b/notebooks/api/07_utils.ipynb index 4b6d5178..2bdd0599 100644 --- a/notebooks/api/07_utils.ipynb +++ b/notebooks/api/07_utils.ipynb @@ -119,6 +119,7 @@ "source": [ "#| exporti\n", "from kornia.geometry.camera.pinhole import PinholeCamera as KorniaPinholeCamera\n", + "from torchio import Subject\n", "\n", "from diffdrr.detector import Detector\n", "\n", @@ -131,9 +132,12 @@ " height: torch.Tensor,\n", " width: torch.Tensor,\n", " detector: Detector,\n", + " subject: Subject,\n", + " \n", " ):\n", " super().__init__(intrinsics, extrinsics, height, width)\n", - " self.f = detector.sdd\n", + " multiplier = -1 if subject.orientation == \"PA\" else 1\n", + " self.sdd = multiplier * detector.sdd\n", " self.delx = detector.delx\n", " self.dely = detector.dely\n", " self.x0 = detector.x0\n", @@ -167,9 +171,9 @@ "\n", "from kornia.geometry.calibration import solve_pnp_dlt\n", "\n", + "from diffdrr.detector import make_intrinsic_matrix\n", "from diffdrr.drr import DRR\n", "from diffdrr.pose import RigidTransform\n", - "from diffdrr.detector import make_intrinsic_matrix\n", "\n", "\n", "def get_pinhole_camera(\n", @@ -180,14 +184,12 @@ " pose = deepcopy(pose).to(device=\"cpu\", dtype=dtype)\n", "\n", " # Make the intrinsic matrix (in pixels)\n", - " fx = drr.detector.sdd / drr.detector.delx\n", - " fy = drr.detector.sdd / drr.detector.dely\n", + " multiplier = -1 if drr.subject.orientation == \"PA\" else 1\n", + " fx = multiplier * drr.detector.sdd / drr.detector.delx\n", + " fy = multiplier * drr.detector.sdd / drr.detector.dely\n", " u0 = drr.detector.x0 / drr.detector.delx + drr.detector.width / 2\n", " v0 = drr.detector.y0 / drr.detector.dely + drr.detector.height / 2\n", - " intrinsics = torch.eye(4)[None]\n", - " intrinsics[0, :3, :3] = make_intrinsic_matrix(drr.detector)\n", - " \n", - " torch.tensor(\n", + " intrinsics = torch.tensor(\n", " [\n", " [\n", " [fx, 0.0, u0, 0.0],\n", @@ -229,6 +231,7 @@ " torch.tensor([drr.detector.height]),\n", " torch.tensor([drr.detector.width]),\n", " drr.detector,\n", + " drr.subject,\n", " )\n", "\n", " return camera" diff --git a/settings.ini b/settings.ini index eb0107be..867a785b 100644 --- a/settings.ini +++ b/settings.ini @@ -1,7 +1,7 @@ [DEFAULT] repo = DiffDRR lib_name = diffdrr -version = 0.4.5 +version = 0.4.6 min_python = 3.8 license = mit black_formatting = True