diff --git a/src/pyronn_torch/conebeam.py b/src/pyronn_torch/conebeam.py
index 816845d6a45675a21d37f94ed5c5371010b7501c..95791748653e23fdce21546c223f78dcaa353277 100644
--- a/src/pyronn_torch/conebeam.py
+++ b/src/pyronn_torch/conebeam.py
@@ -2,20 +2,20 @@
# Copyright © 2020 Stephan Seitz <stephan.seitz@fau.de>
#
# Distributed under terms of the GPLv3 license.
-
"""
"""
import numpy as np
-import sympy as sp
import torch
import pyronn_torch
+import sympy as sp
class _ForwardProjection(torch.autograd.Function):
-
- def __init__(self, projection_shape,
+ def __init__(self,
+ projection_shape,
+ volume_shape,
source_points,
inverse_matrices,
projection_matrices,
@@ -25,6 +25,7 @@ class _ForwardProjection(torch.autograd.Function):
step_size=1.,
with_texture=True):
self.projection_shape = projection_shape
+ self.volume_shape = volume_shape
self.source_points = source_points
self.inverse_matrices = inverse_matrices
self.projection_matrices = projection_matrices
@@ -36,41 +37,31 @@ class _ForwardProjection(torch.autograd.Function):
def forward(self, volume):
volume = volume.cuda().contiguous()
- projection = torch.zeros(self.projection_shape, device='cuda', requires_grad=volume.requires_grad)
+ projection = torch.zeros(self.projection_shape,
+ device='cuda',
+ requires_grad=volume.requires_grad)
assert pyronn_torch.cpp_extension
if self.with_texture:
pyronn_torch.cpp_extension.call_Cone_Projection_Kernel_Tex_Interp_Launcher(
- self.inverse_matrices,
- projection,
- self.source_points,
- self.step_size,
- volume,
- *self.volume_spacing)
+ self.inverse_matrices, projection, self.source_points,
+ self.step_size, volume, *self.volume_spacing)
else:
pyronn_torch.cpp_extension.call_Cone_Projection_Kernel_Launcher(
- self.inverse_matrices,
- projection,
- self.source_points,
- self.step_size,
- volume,
- *self.volume_spacing)
+ self.inverse_matrices, projection, self.source_points,
+ self.step_size, volume, *self.volume_spacing)
return projection,
def backward(self, *projection_grad):
projection_grad = projection_grad[0]
- self.projection_matrices
- volume_grad = torch.zeros(self.volume_shape, device='cuda')
+ volume_grad = torch.zeros(self.volume_shape, device='cuda', requires_grad=projection_grad.requires_grad)
assert pyronn_torch.cpp_extension
pyronn_torch.cpp_extension.call_Cone_Backprojection3D_Kernel_Launcher(
- self.projection_matrices,
- projection_grad,
- self.projection_multiplier,
- volume_grad,
- self.volume_origin,
- self.volume_spacing)
+ self.projection_matrices, projection_grad,
+ self.projection_multiplier, volume_grad, *self.volume_origin,
+ *self.volume_spacing)
return volume_grad,
@@ -82,14 +73,8 @@ class _BackwardProjection(torch.autograd.Function):
class ConeBeamProjector:
-
- def __init__(self,
- volume_shape,
- volume_spacing,
- volume_origin,
- projection_shape,
- projection_spacing,
- projection_origin,
+ def __init__(self, volume_shape, volume_spacing, volume_origin,
+ projection_shape, projection_spacing, projection_origin,
projection_matrices):
self._volume_shape = volume_shape
self._volume_origin = volume_origin
@@ -108,74 +93,82 @@ class ConeBeamProjector:
volume_spacing = pyconrad.config.get_reco_spacing()
volume_origin = pyconrad.config.get_reco_origin()
projection_shape = pyconrad.config.get_sino_shape()
- projection_spacing = [pyconrad.config.get_geometry().getPixelDimensionY(),
- pyconrad.config.get_geometry().getPixelDimensionX()]
- projection_origin = [pyconrad.config.get_geometry().getDetectorOffsetV(),
- pyconrad.config.get_geometry().getDetectorOffsetU()]
+ projection_spacing = [
+ pyconrad.config.get_geometry().getPixelDimensionY(),
+ pyconrad.config.get_geometry().getPixelDimensionX()
+ ]
+ projection_origin = [
+ pyconrad.config.get_geometry().getDetectorOffsetV(),
+ pyconrad.config.get_geometry().getDetectorOffsetU()
+ ]
projection_matrices = pyconrad.config.get_projection_matrices()
- obj = cls(volume_shape,
- volume_spacing,
- volume_origin,
- projection_shape,
- projection_spacing,
- projection_origin,
+ obj = cls(volume_shape, volume_spacing, volume_origin,
+ projection_shape, projection_spacing, projection_origin,
projection_matrices)
return obj
def new_volume_tensor(self, requires_grad=False):
- return torch.zeros(self._volume_shape, requires_grad=requires_grad).cuda()
+ return torch.zeros(self._volume_shape,
+ requires_grad=requires_grad).cuda()
def new_projection_tensor(self, requires_grad=False):
- return torch.zeros(self._projection_shape, requires_grad=requires_grad).cuda()
+ return torch.zeros(self._projection_shape,
+ requires_grad=requires_grad).cuda()
def project_forward(self, volume, step_size=1., use_texture=True):
- return _ForwardProjection(self._projection_shape,
- self._source_points,
- self._inverse_matrices,
+ return _ForwardProjection(self._projection_shape, self._volume_shape,
+ self._source_points, self._inverse_matrices,
self._projection_matrices,
- self._volume_origin,
- self._volume_shape,
- self._projection_multiplier,
- step_size,
+ self._volume_origin, self._volume_spacing,
+ self._projection_multiplier, step_size,
use_texture).forward(volume)[0]
- def project_backward(self, projection_stack, step_size=1., use_texture=True):
- return _BackwardProjection(self._projection_shape,
- self._source_points,
- self._inverse_matrices,
+ def project_backward(self,
+ projection_stack,
+ step_size=1.,
+ use_texture=True):
+ return _BackwardProjection(self._projection_shape, self._volume_shape,
+ self._source_points, self._inverse_matrices,
self._projection_matrices,
- self._volume_origin,
- self._volume_shape,
- self._projection_multiplier,
- step_size,
- use_texture).backward(projection_stack)[0]
+ self._volume_origin, self._volume_spacing,
+ self._projection_multiplier, step_size,
+ use_texture).forward(projection_stack)[0]
def _calc_inverse_matrices(self):
if self._projection_matrices_numpy is None:
return
- self._projection_matrices = torch.stack(tuple(
- map(torch.from_numpy, self._projection_matrices_numpy))).cuda().contiguous()
-
- inv_spacing = np.array([1/s for s in reversed(self._volume_spacing)], np.float32)
-
- camera_centers = map(lambda x: np.array(sp.Matrix(x).nullspace(), np.float32), self._projection_matrices_numpy)
-
- source_points = map(lambda x: (x[0, :3] / x[0, 3] * inv_spacing - np.array(list(reversed(self._volume_origin)))
- * inv_spacing).astype(np.float32), camera_centers)
-
- inv_matrices = map(lambda x: (np.linalg.inv(x[:3, :3]) *
- inv_spacing).astype(np.float32), self._projection_matrices_numpy)
-
- self._inverse_matrices = torch.stack(tuple(map(torch.from_numpy, inv_matrices))).cuda().contiguous()
- self._source_points = torch.stack(tuple(map(torch.from_numpy, source_points))).cuda().contiguous()
+ self._projection_matrices = torch.stack(
+ tuple(torch.from_numpy(p.astype(np.float32)) for p in self._projection_matrices_numpy)).cuda().contiguous()
+
+ inv_spacing = np.array([1 / s for s in reversed(self._volume_spacing)],
+ np.float32)
+
+ camera_centers = map(
+ lambda x: np.array(sp.Matrix(x).nullspace(), np.float32),
+ self._projection_matrices_numpy)
+
+ source_points = map(
+ lambda x: (x[0, :3] / x[0, 3] * inv_spacing - np.array(
+ list(reversed(self._volume_origin))) * inv_spacing).astype(
+ np.float32), camera_centers)
+
+ inv_matrices = map(
+ lambda x:
+ (np.linalg.inv(x[:3, :3]) * inv_spacing).astype(np.float32),
+ self._projection_matrices_numpy)
+
+ self._inverse_matrices = torch.stack(
+ tuple(map(torch.from_numpy, inv_matrices))).cuda().contiguous()
+ self._source_points = torch.stack(
+ tuple(map(torch.from_numpy, source_points))).cuda().contiguous()
self._projection_multiplier = 1.
@property
def projection_matrices(self):
- return self._projection_matrices_numpy
+ return self._projection_matrices_numpy
@projection_matrices.setter
def projection_matrices(self, numpy_matrices):
- self._projection_matrices_numpy = numpy_matrices
- self._calc_inverse_matrices()
+ self._projection_matrices_numpy = numpy_matrices
+ self._calc_inverse_matrices()
diff --git a/tests/test_projection.py b/tests/test_projection.py
index 96f815920c1d2d0c186b14c7b653b4cf28c5cda7..35ff398a5842ab16e5617289284ef5b6afa87305 100644
--- a/tests/test_projection.py
+++ b/tests/test_projection.py
@@ -6,7 +6,6 @@
"""
"""
-import os
import numpy as np
import pytest
@@ -73,11 +72,14 @@ def test_projection_backward(with_texture, with_backward):
[-9.99847710e-01, -1.74524058e-2, 0.00000000e+0,
6.00000000e+2]]])
)
-
projection = projector.new_projection_tensor(requires_grad=True if with_backward else False)
projection += 1.
+
result = projector.project_backward(projection, use_texture=with_texture)
+ import pyconrad.autoinit
+ pyconrad.imshow(result)
+ assert result.shape == projector._volume_shape
assert result is not None
if with_backward:
@@ -99,6 +101,8 @@ def test_conrad_config(with_backward, with_texture=True):
volume += 1.
result = projector.project_forward(volume, use_texture=with_texture)
+ import pyconrad.autoinit
+ pyconrad.imshow(result)
assert result is not None
if with_backward:
@@ -107,3 +111,49 @@ def test_conrad_config(with_backward, with_texture=True):
loss = result.mean()
loss.backward()
+
+
+def test_projection_backward_conrad(with_texture=True, with_backward=True):
+ import pytest
+ pytest.importorskip("pyconrad")
+
+ projector = pyronn_torch.ConeBeamProjector.from_conrad_config()
+
+ projection = projector.new_projection_tensor(requires_grad=True if with_backward else False)
+
+ projection += 1000.
+
+ result = projector.project_backward(projection, use_texture=with_texture)
+ import pyconrad.autoinit
+ pyconrad.imshow(result)
+ assert result.shape == projector._volume_shape
+
+ assert result is not None
+ if with_backward:
+ assert projection.requires_grad
+ assert result.requires_grad
+
+ loss = result.mean()
+ loss.backward()
+
+
+def test_conrad_forward_backward():
+ import pytest
+ pytest.importorskip("pyconrad")
+
+ projector = pyronn_torch.ConeBeamProjector.from_conrad_config()
+
+ volume = projector.new_volume_tensor()
+
+ volume += 1.
+ result = projector.project_forward(volume)
+
+ reco = projector.project_backward(result)
+ import pyconrad.autoinit
+ pyconrad.imshow(reco)
+
+ assert result is not None
+ assert reco is not None
+
+ while True:
+ pass