ENH: rotate geometry instead of data in 2D

Author: Holger Kohr

odl/tomo/backends/astra_cpu.py

@@ -97,10 +97,7 @@ def astra_cpu_forward_projector(vol_data, geometry, proj_space, out=None):
                               impl='cpu')
 
     # Create ASTRA data structures
-    # Since ASTRA uses (rows, cols) == (-y, x) as coordinate system, we need
-    # to rotate by 90 degrees counter-clockwise to match the axes (x, y) as
-    # used in `astra_volume_geometry`.
-    vol_data_arr = np.rot90(np.asarray(vol_data), 1)
+    vol_data_arr = np.asarray(vol_data)
     vol_id = astra_data(vol_geom, datatype='volume', data=vol_data_arr,
                         allow_copy=True)
 
@@ -195,8 +192,8 @@ def astra_cpu_back_projector(proj_data, geometry, reco_space, out=None):
 
     # Convert out to correct dtype and order if needed.
     # Since we transpose, we need to use F-contiguousness.
-    with writable_array(out, dtype='float32', order='F') as out_arr:
-        vol_id = astra_data(vol_geom, datatype='volume', data=out_arr.T,
+    with writable_array(out, dtype='float32', order='C') as out_arr:
+        vol_id = astra_data(vol_geom, datatype='volume', data=out_arr,
                             ndim=reco_space.ndim)
         # Create algorithm
         algo_id = astra_algorithm('backward', ndim, vol_id, sino_id, proj_id,
@@ -205,13 +202,6 @@ def astra_cpu_back_projector(proj_data, geometry, reco_space, out=None):
         # Run algorithm
         astra.algorithm.run(algo_id)
 
-        # If we hadn't transposed `out_arr` we would have to rotate
-        # clockwise by 90 degrees to invert the transition from (rows, cols)
-        # to (x, y). Since transposition happens automatically when exiting
-        # this context, we need to flip vertically, because
-        # rot90(a, -1) == fliplr(a.T)
-        out_arr[:] = np.fliplr(out_arr)
-
     # Weight the adjoint by appropriate weights
     scaling_factor = float(proj_data.space.weighting.const)
     scaling_factor /= float(reco_space.weighting.const)

odl/tomo/backends/astra_cuda.py

@@ -87,10 +87,7 @@ def call_forward(self, vol_data, out=None):
 
         # Copy data to GPU memory
         if self.geometry.ndim == 2:
-            # Rotate 90 degrees counter-clockwise to transfer from coordinate
-            # system (x, y) to (rows, cols)
-            # TODO: rotate the geometry instead
-            astra.data2d.store(self.vol_id, np.rot90(vol_data.asarray(), 1))
+            astra.data2d.store(self.vol_id, vol_data.asarray())
         elif self.geometry.ndim == 3:
             astra.data3d.store(self.vol_id, vol_data.asarray())
         else:
@@ -127,8 +124,7 @@ def create_ids(self):
 
         if proj_ndim == 2:
             astra_proj_shape = proj_shape
-            astra_vol_shape = (self.reco_space.shape[1],
-                               self.reco_space.shape[0])
+            astra_vol_shape = self.reco_space.shape
         elif proj_ndim == 3:
             astra_proj_shape = (proj_shape[2], proj_shape[0], proj_shape[1])
             astra_vol_shape = self.reco_space.shape
@@ -247,12 +243,7 @@ def call_backward(self, proj_data, out=None):
         astra.algorithm.run(self.algo_id)
 
         # Copy result to CPU memory
-        if self.geometry.ndim == 2:
-            # Rotate 90 degrees clockwise from coordinate system
-            # (rows, cols) to (x, y)
-            out[:] = np.rot90(self.out_array, -1)
-        elif self.geometry.ndim == 3:
-            out[:] = self.out_array
+        out[:] = self.out_array
 
         # Fix scaling to weight by pixel/voxel size
         out *= astra_cuda_bp_scaling_factor(
@@ -273,8 +264,7 @@ def create_ids(self):
 
         if proj_ndim == 2:
             astra_proj_shape = proj_shape
-            astra_vol_shape = (self.reco_space.shape[1],
-                               self.reco_space.shape[0])
+            astra_vol_shape = self.reco_space.shape
         elif proj_ndim == 3:
             astra_proj_shape = (proj_shape[2], proj_shape[0], proj_shape[1])
             astra_vol_shape = self.reco_space.shape

odl/tomo/backends/astra_setup.py

@@ -53,6 +53,7 @@
 from odl.discr import DiscreteLp, DiscreteLpElement
 from odl.tomo.geometry import (
     Geometry, DivergentBeamGeometry, ParallelBeamGeometry, FlatDetector)
+from odl.tomo.util.utility import euler_matrix
 from odl.util.utility import pkg_supports
 
 
@@ -165,17 +166,21 @@ def astra_volume_geometry(discr_reco):
                 'non-isotropic pixels in 2d volumes not supported by ASTRA '
                 'v{}'.format(ASTRA_VERSION))
         # Given a 2D array of shape (x, y), a volume geometry is created as:
-        #    astra.create_vol_geom(y, x, x_min, x_max, y_min, y_max)
+        #    astra.create_vol_geom(x, y, y_min, y_max, x_min, x_max)
         # yielding a dictionary:
-        #   {'GridRowCount': y,
-        #    'GridColCount': x,
-        #    'WindowMinX': x_min,
-        #    'WindowMaxX': x_max,
-        #    'WindowMinY': y_min,
-        #    'WindowMaxY': y_max}
-        vol_geom = astra.create_vol_geom(vol_shp[1], vol_shp[0],
-                                         vol_min[0], vol_max[0],
-                                         vol_min[1], vol_max[1])
+        #   {'GridRowCount': x,
+        #    'GridColCount': y,
+        #    'WindowMinX': y_min,
+        #    'WindowMaxX': y_max,
+        #    'WindowMinY': x_min,
+        #    'WindowMaxY': x_max}
+        #
+        # NOTE: this setting is flipped with respect to x and y. We do this
+        # as part of a global rotation of the geometry by -90 degrees, which
+        # avoids rotating the data.
+        vol_geom = astra.create_vol_geom(vol_shp[0], vol_shp[1],
+                                         vol_min[1], vol_max[1],
+                                         vol_min[0], vol_max[0])
     elif discr_reco.ndim == 3:
         # Not supported in all versions of ASTRA
         if (not discr_reco.partition.has_isotropic_cells and
@@ -300,19 +305,24 @@ def astra_conebeam_2d_geom_to_vec(geometry):
     .. _ASTRA projection geometry documentation:
        http://www.astra-toolbox.com/docs/geom3d.html#projection-geometries
     """
+    # Instead of rotating the data by 90 degrees counter-clockwise,
+    # we subtract pi/2 from the geometry angles, thereby rotating the
+    # geometry by 90 degrees clockwise
+    rot_minus_90 = euler_matrix(-np.pi / 2)
     angles = geometry.angles
     vectors = np.zeros((angles.size, 6))
 
     for ang_idx, angle in enumerate(angles):
         # Source position
-        vectors[ang_idx, 0:2] = geometry.src_position(angle)
+        vectors[ang_idx, 0:2] = rot_minus_90.dot(geometry.src_position(angle))
 
         # Center of detector
         mid_pt = geometry.det_params.mid_pt
-        vectors[ang_idx, 2:4] = geometry.det_point_position(angle, mid_pt)
+        vectors[ang_idx, 2:4] = rot_minus_90.dot(
+            geometry.det_point_position(angle, mid_pt))
 
         # Vector from detector pixel 0 to 1
-        det_axis = geometry.det_axis(angle)
+        det_axis = rot_minus_90.dot(geometry.det_axis(angle))
         px_size = geometry.det_partition.cell_sides[0]
         vectors[ang_idx, 4:6] = det_axis * px_size
 
@@ -413,7 +423,10 @@ def astra_projection_geometry(geometry):
         # TODO: change to parallel_vec when available
         det_width = geometry.det_partition.cell_sides[0]
         det_count = geometry.detector.size
-        angles = geometry.angles
+        # Instead of rotating the data by 90 degrees counter-clockwise,
+        # we subtract pi/2 from the geometry angles, thereby rotating the
+        # geometry by 90 degrees clockwise
+        angles = geometry.angles - np.pi / 2
         proj_geom = astra.create_proj_geom('parallel', det_width, det_count,
                                            angles)