Skip to content

Commit

Permalink
Merge pull request #1059 from kohr-h/vectorize_geometries
Browse files Browse the repository at this point in the history 
Vectorize geometries
  • Loading branch information
Holger Kohr authored Sep 26, 2017
2 parents b14b962 + 4e42a1b commit ee4ad6d
Show file tree
Hide file tree
Showing 7 changed files with 1,939 additions and 676 deletions.
113 changes: 57 additions & 56 deletions odl/tomo/backends/astra_setup.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,8 @@

from odl.discr import DiscreteLp, DiscreteLpElement
from odl.tomo.geometry import (
Geometry, DivergentBeamGeometry, ParallelBeamGeometry, FlatDetector)
Geometry, DivergentBeamGeometry, ParallelBeamGeometry,
Flat1dDetector, Flat2dDetector)
from odl.tomo.util.utility import euler_matrix


Expand Down Expand Up @@ -245,27 +246,26 @@ def astra_conebeam_3d_geom_to_vec(geometry):
angles = geometry.angles
vectors = np.zeros((angles.size, 12))

for ang_idx, angle in enumerate(angles):
# Source position
vectors[ang_idx, 0:3] = geometry.src_position(angle)

# Center of detector in 3D space
mid_pt = geometry.det_params.mid_pt
vectors[ang_idx, 3:6] = geometry.det_point_position(angle, mid_pt)

# Vectors from detector pixel (0, 0) to (1, 0) and (0, 0) to (0, 1)
det_axes = geometry.det_axes(angle)
px_sizes = geometry.det_partition.cell_sides

# Swap detector axes to have better memory layout in projection data.
# ASTRA produces `(v, theta, u)` layout, and to map to ODL layout
# `(theta, u, v)` a complete roll must be performed, which is the
# worst case (compeltely discontiguous).
# Instead we swap `u` and `v`, resulting in the effective ASTRA result
# `(u, theta, v)`. Here we only need to swap axes 0 and 1, which
# keeps at least contiguous blocks in `v`.
vectors[ang_idx, 9:12] = det_axes[0] * px_sizes[0]
vectors[ang_idx, 6:9] = det_axes[1] * px_sizes[1]
# Source position
vectors[:, 0:3] = geometry.src_position(angles)

# Center of detector in 3D space
mid_pt = geometry.det_params.mid_pt
vectors[:, 3:6] = geometry.det_point_position(angles, mid_pt)

# Vectors from detector pixel (0, 0) to (1, 0) and (0, 0) to (0, 1)
# `det_axes` gives shape (N, 2, 3), swap to get (2, N, 3)
det_axes = np.moveaxis(geometry.det_axes(angles), -2, 0)
px_sizes = geometry.det_partition.cell_sides
# Swap detector axes to have better memory layout in projection data.
# ASTRA produces `(v, theta, u)` layout, and to map to ODL layout
# `(theta, u, v)` a complete roll must be performed, which is the
# worst case (compeltely discontiguous).
# Instead we swap `u` and `v`, resulting in the effective ASTRA result
# `(u, theta, v)`. Here we only need to swap axes 0 and 1, which
# keeps at least contiguous blocks in `v`.
vectors[:, 9:12] = det_axes[0] * px_sizes[0]
vectors[:, 6:9] = det_axes[1] * px_sizes[1]

# ASTRA has (z, y, x) axis convention, in contrast to (x, y, z) in ODL,
# so we need to adapt to this by changing the order.
Expand Down Expand Up @@ -317,19 +317,21 @@ def astra_conebeam_2d_geom_to_vec(geometry):
angles = geometry.angles
vectors = np.zeros((angles.size, 6))

for ang_idx, angle in enumerate(angles):
# Source position
vectors[ang_idx, 0:2] = rot_minus_90.dot(geometry.src_position(angle))
# Source position
src_pos = geometry.src_position(angles)
vectors[:, 0:2] = rot_minus_90.dot(src_pos.T).T # dot along 2nd axis

# Center of detector
mid_pt = geometry.det_params.mid_pt
vectors[ang_idx, 2:4] = rot_minus_90.dot(
geometry.det_point_position(angle, mid_pt))
# Center of detector
mid_pt = geometry.det_params.mid_pt
# Need to cast `mid_pt` to float since otherwise the empty axis is
# not removed
centers = geometry.det_point_position(angles, float(mid_pt))
vectors[:, 2:4] = rot_minus_90.dot(centers.T).T

# Vector from detector pixel 0 to 1
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
# Vector from detector pixel 0 to 1
det_axis = rot_minus_90.dot(geometry.det_axis(angles).T).T
px_size = geometry.det_partition.cell_sides[0]
vectors[:, 4:6] = det_axis * px_size

return vectors

Expand Down Expand Up @@ -371,28 +373,27 @@ def astra_parallel_3d_geom_to_vec(geometry):
angles = geometry.angles
vectors = np.zeros((angles.shape[0], 12))

for ang_idx, angle in enumerate(angles):
mid_pt = geometry.det_params.mid_pt

# Ray direction = -(detector-to-source normal vector)
vectors[ang_idx, 0:3] = -geometry.det_to_src(angle, mid_pt)
mid_pt = geometry.det_params.mid_pt

# Center of the detector in 3D space
vectors[ang_idx, 3:6] = geometry.det_point_position(angle, mid_pt)
# Ray direction = -(detector-to-source normal vector)
vectors[:, 0:3] = -geometry.det_to_src(angles, mid_pt)

# Vectors from detector pixel (0, 0) to (1, 0) and (0, 0) to (0, 1)
det_axes = geometry.det_axes(angle)
px_sizes = geometry.det_partition.cell_sides
# Center of the detector in 3D space
vectors[:, 3:6] = geometry.det_point_position(angles, mid_pt)

# Swap detector axes to have better memory layout in projection data.
# ASTRA produces `(v, theta, u)` layout, and to map to ODL layout
# `(theta, u, v)` a complete roll must be performed, which is the
# worst case (compeltely discontiguous).
# Instead we swap `u` and `v`, resulting in the effective ASTRA result
# `(u, theta, v)`. Here we only need to swap axes 0 and 1, which
# keeps at least contiguous blocks in `v`.
vectors[ang_idx, 9:12] = det_axes[0] * px_sizes[0]
vectors[ang_idx, 6:9] = det_axes[1] * px_sizes[1]
# Vectors from detector pixel (0, 0) to (1, 0) and (0, 0) to (0, 1)
# `det_axes` gives shape (N, 2, 3), swap to get (2, N, 3)
det_axes = np.moveaxis(geometry.det_axes(angles), -2, 0)
px_sizes = geometry.det_partition.cell_sides
# Swap detector axes to have better memory layout in projection data.
# ASTRA produces `(v, theta, u)` layout, and to map to ODL layout
# `(theta, u, v)` a complete roll must be performed, which is the
# worst case (compeltely discontiguous).
# Instead we swap `u` and `v`, resulting in the effective ASTRA result
# `(u, theta, v)`. Here we only need to swap axes 0 and 1, which
# keeps at least contiguous blocks in `v`.
vectors[:, 9:12] = det_axes[0] * px_sizes[0]
vectors[:, 6:9] = det_axes[1] * px_sizes[1]

# ASTRA has (z, y, x) axis convention, in contrast to (x, y, z) in ODL,
# so we need to adapt to this by changing the order.
Expand Down Expand Up @@ -431,7 +432,7 @@ def astra_projection_geometry(geometry):
raise ValueError('non-uniform detector sampling is not supported')

if (isinstance(geometry, ParallelBeamGeometry) and
isinstance(geometry.detector, FlatDetector) and
isinstance(geometry.detector, (Flat1dDetector, Flat2dDetector)) and
geometry.ndim == 2):
# TODO: change to parallel_vec when available
det_width = geometry.det_partition.cell_sides[0]
Expand All @@ -444,14 +445,14 @@ def astra_projection_geometry(geometry):
angles)

elif (isinstance(geometry, DivergentBeamGeometry) and
isinstance(geometry.detector, FlatDetector) and
isinstance(geometry.detector, (Flat1dDetector, Flat2dDetector)) and
geometry.ndim == 2):
det_count = geometry.detector.size
vec = astra_conebeam_2d_geom_to_vec(geometry)
proj_geom = astra.create_proj_geom('fanflat_vec', det_count, vec)

elif (isinstance(geometry, ParallelBeamGeometry) and
isinstance(geometry.detector, FlatDetector) and
isinstance(geometry.detector, (Flat1dDetector, Flat2dDetector)) and
geometry.ndim == 3):
# Swap detector axes (see astra_*_3d_to_vec)
det_row_count = geometry.det_partition.shape[0]
Expand All @@ -461,7 +462,7 @@ def astra_projection_geometry(geometry):
det_col_count, vec)

elif (isinstance(geometry, DivergentBeamGeometry) and
isinstance(geometry.detector, FlatDetector) and
isinstance(geometry.detector, (Flat1dDetector, Flat2dDetector)) and
geometry.ndim == 3):
# Swap detector axes (see astra_*_3d_to_vec)
det_row_count = geometry.det_partition.shape[0]
Expand Down
Loading

0 comments on commit ee4ad6d

Please sign in to comment.