Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WIP: Add vector geometries #1418

Open
wants to merge 24 commits into
base: master
Choose a base branch
from

Conversation

kohr-h
Copy link
Member

@kohr-h kohr-h commented Aug 22, 2018

I picked up an old WIP branch that I currently need. The code adds support for ASTRA vec geometries in a fairly direct way.

Since the coordinate systems of ODL and ASTRA differ, I've taken the following approach:

  • The geometries ParallelVecGeometry and ConeVecGeometry take a matrix of shape (N, 6) for 2D or (N, 12) for 3D and determine ndim accordingly.
  • The vectors are stored as arrays without further modification.
  • The axes X,Y,Z in the definition of the vectors refer to ODL coordinates, e.g., (rayX, rayY) is defined with respect to the ODL coordinate system.
    As a consequence, a conversion step is necessary before handing these vectors over to ASTRA (or, inversely, when we get vectors from ASTRA or from somewhere with ASTRA coordinates).
    The reasoning behind this is that ODL should use consistent coordinates internally and translate to other systems when interfacing with other libraries. The "list of vectors" format is just treated as a representation for geometries, without any prescribed coordinate system.

I've added a couple of "check" examples to make sure that the transformations are defined correctly. It's looking good, and the existing geometries are still correct.

TODO:

  • Move axis swapping in 3D to a separate function so that *3d_geom_to_vecs returns the ASTRA vectors "as expected", and downstream code does the swapping instead.
  • Add adjointness test for vec geometries
  • Clean up commits

Closes #1023

Sorry, something went wrong.

angle_partition, detector_partition, src_radius=1000, det_radius=100,
axis=[1, 0, 0])

circle_vecs = odl.tomo.astra_conebeam_3d_geom_to_vec(circle_geometry)
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've renamed this function, TODO change here

@@ -54,6 +54,7 @@ def sparse_meshgrid(*x):
xi = np.asarray(xi)
slc = [None] * n
slc[ax] = slice(None)
slc = tuple(slc)
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Code like this without conversion from list to tuple raises a DeprecationWarning in NumPy 1.15, we should apply the fix throughout.

'DetectorCount': 10, 'DetectorWidth': 0.2}

assert is_subdict(correct_subdict, proj_geom)
assert 'ProjectionAngles' in proj_geom
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This test was kind of pointless since it really tests ASTRA rather than ODL, and some implementation detail of ODL. Not sure if we really need a replacement.

elif geometry.ndim == 3:
det_px_area = geometry.det_partition.cell_volume
scaling_factor *= (det_px_area ** 2 /
reco_space.cell_volume ** 2)
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Some testing needed here

# ODL x = ASTRA -y, ODL y = ASTRA x
# Thus: ODL->ASTRA conversion is a rotation by +90 degrees
rot_90 = euler_matrix(np.pi / 2)
# Bulk matrix-vector product
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

with the transposed matrix

if vecs.shape[1] == 6:
# 2D geometry
# ODL x = ASTRA -y, ODL y = ASTRA x
# Thus: ODL->ASTRA conversion is a rotation by +90 degrees
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe these comments are redundant

# TODO: shapes are not quite right here. Scalar works, but
# vectorized goes wrong
ray_dir = np.empty(bcast_shape + (self.ndim,), dtype=float)
print('result shape:', ray_dir.shape)
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Debug statement, TODO remove

Regarding vectorization, I'd suggest putting it in a separate PR.

@adler-j
Copy link
Member

adler-j commented Aug 23, 2018

Should this be reviewed?

@kohr-h
Copy link
Member Author

kohr-h commented Aug 24, 2018

Should this be reviewed?

Hm, I think some things are still a bit in the flow, so I'll make it a WIP in the title. I'll give a heads-up when it's ready.

@kohr-h kohr-h changed the title Add vector geometries WIP: Add vector geometries Aug 24, 2018
@kohr-h kohr-h force-pushed the issue-1023__vec_geometry branch from 621e03a to 8d211d6 Compare September 10, 2018 08:56
@adler-j adler-j mentioned this pull request Sep 11, 2018
20 tasks
@pep8speaks
Copy link

pep8speaks commented Feb 13, 2019

Checking updated PR...

Line 508:5: E303 too many blank lines (2)

Line 1650:9: E303 too many blank lines (2)

Comment last updated at 2020-11-08 00:20:59 UTC

@kohr-h kohr-h force-pushed the issue-1023__vec_geometry branch from f6bf263 to 4052454 Compare March 13, 2019 21:15
@kohr-h kohr-h force-pushed the issue-1023__vec_geometry branch from 4052454 to 521d7cf Compare November 5, 2019 13:19
@kohr-h
Copy link
Member Author

kohr-h commented Apr 13, 2020

Deferring until #1540 has landed. Too many overlapping changes.

@kohr-h
Copy link
Member Author

kohr-h commented Sep 6, 2020

Okay, time to pick this one up again.

@adriaangraas
Copy link
Contributor

adriaangraas commented Sep 6, 2020

I did a merge with a few minor fixes, may or may not be useful. See your repo.

@kohr-h
Copy link
Member Author

kohr-h commented Sep 6, 2020

I did a merge with a few minor fixes, may or may not be useful. See your repo.

Thanks, I saw it. I'll take it as reference. 👍

@kohr-h kohr-h force-pushed the issue-1023__vec_geometry branch from 1ca6fed to 252b7d0 Compare September 6, 2020 10:59
Comment on lines 443 to 452
elif isinstance(geometry, VecGeometry):
scaling_factor = (geometry.det_partition.cell_volume /
vol_space.cell_volume)
elif isinstance(geometry, ParallelVecGeometry):
if geometry.ndim == 2:
# Scales with 1 / cell_volume
scaling_factor *= float(vol_space.cell_volume)
elif geometry.ndim == 3:
# Scales with cell volume
scaling_factor /= vol_space.cell_volume
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Some part of this logic makes no sense...

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Vector geometries
4 participants