ENH: Add parallel_beam_geometry, closes #711 by adler-j · Pull Request #775 · odlgroup/odl

adler-j · 2016-12-20T11:21:32Z

Adds a simple utility to create parallel beam geometries. Useful for testing and such cases where you don't really care about what the geometry is but simply "want it to work".

kohr-h

Some minor points, and having a reference would be nice.

kohr-h · 2016-12-20T12:35:14Z

odl/tomo/geometry/parallel.py

+
+    Returns
+    -------
+    geometry : ParallelGeometry


`ParallelGeometry`

kohr-h · 2016-12-20T12:35:26Z

odl/tomo/geometry/parallel.py

+    -------
+    geometry : ParallelGeometry
+        If ``space`` is 2d, returns a ``Parallel2dGeometry``.
+        If ``space`` is 3d, returns a ``Parallel3dAxisGeometry``.


Single ticks on classes

kohr-h · 2016-12-20T12:36:19Z

odl/tomo/geometry/parallel.py

+    """Create default parallel beam geometry from space.
+
+    This default geometry gives a fully sampled sinogram according to the
+    nyquist criterium. In general this gives a very large number of samples.


Nyquist
~~criterium~~ criterion

Do we have a reference?

Kak A C and Slaney M 2001 Principles of Computerized Tomographic Imaging (Philadelphia, PA: SIAM)

Thank you for the reference @chongchenmath!

Kak A C and Slaney M 2001 Principles of Computerized Tomographic Imaging (Philadelphia, PA: SIAM)

In general it is not a good idea to give an entire book as a reference for a specific result. The sampling theory based on the Nyquist criteria, as outlined in section 4.2and pp. 73-78 in "Mathematical Methods in Image Reconstruction" by F. Natterer and F. Wübbeling, SIAM, 2001 is only developed for the 2D case. Of course, 3D parallel beam geometry reduces to a sequence of 2D problems, so it also makes sense in that setting.

kohr-h · 2016-12-20T12:37:44Z

odl/tomo/geometry/parallel.py

+    Parameters
+    ----------
+    space : `DiscreteLp`
+        Space in which the image should live. Needs to be 2d or 3d.


Perhaps
Reconstruction or "image" space.

kohr-h · 2016-12-20T12:43:37Z

odl/tomo/geometry/parallel.py

+    space : `DiscreteLp`
+        Space in which the image should live. Needs to be 2d or 3d.
+    angles : int, optional
+        Number of angles. Default: Enough to fully sample the data.


"fully sample the data" - I get what you mean, but it could be clearer. Why not refer to the Nyquist criterion as you did before and perhaps include a link to a source above? Like this one.

adler-j · 2016-12-20T13:26:12Z

Will fix comments and I agree a reference is a must.

chongchenmath · 2016-12-20T13:27:58Z

odl/tomo/geometry/parallel.py

+    nyquist criterium. In general this gives a very large number of samples.
+
+    This is intended for simple test cases where users do not need the full
+    flexibility of the geometries, but simply wants a geomoetry that works.


adler-j · 2016-12-25T12:07:53Z

Fixed comments, updated doc, added tests. Next review!

kohr-h

Some more remarks, fix them and merge.

kohr-h · 2017-01-01T20:59:37Z

examples/tomo/filtered_backprojection_parallel_2d.py

 """
-Example using a filtered back-projection in 2d using the ray transform and a
+Example creating a filtered back-projection in 2d using the ray transform and a
 ramp filter. The ramp filter is implemented in fourier space.


kohr-h · 2017-01-01T21:00:43Z

examples/tomo/filtered_backprojection_parallel_2d.py

 Also note that ODL has a utility function, `fbp_op` that can be used to perform
-filtered back-projection.
+filtered back-projection, and that this example is intended to show how this
+could be implemented in ODL.


perhaps "by hand using ODL"

kohr-h · 2017-01-01T21:06:25Z

odl/tomo/geometry/parallel.py



+def parallel_beam_geometry(space, angles=None, det_shape=None):
+    """Create default parallel beam geometry from space.


kohr-h · 2017-01-01T21:06:49Z

odl/tomo/geometry/parallel.py

+    """Create default parallel beam geometry from space.
+
+    This default geometry gives a fully sampled sinogram according to the
+    nyquist criterion. In general this gives a very large number of samples.


Nyquist
"criterion, which in general results in a very ..." (reason: too many "this")

kohr-h · 2017-01-01T21:08:14Z

odl/tomo/geometry/parallel.py

+    nyquist criterion. In general this gives a very large number of samples.
+
+    This is intended for simple test cases where users do not need the full
+    flexibility of the geometries, but simply want a geomoetry that works.


typo "geomoetry"

kohr-h · 2017-01-01T21:12:44Z

odl/tomo/geometry/parallel.py

+
+    The geometry returned by this function satisfies these conditions exactly.
+
+    If the domain is 3 dimensional, the geometry is "separable", in that each


3-dimensional

kohr-h · 2017-01-01T21:13:05Z

odl/tomo/geometry/parallel.py

+    The geometry returned by this function satisfies these conditions exactly.
+
+    If the domain is 3 dimensional, the geometry is "separable", in that each
+    slice along the z-dimension of the data is treated as independed 2d data.


Good description!

kohr-h · 2017-01-01T21:14:14Z

odl/tomo/geometry/parallel.py

+    corners = space.domain.corners()[:, :2]
+    rho = np.max(np.linalg.norm(corners, axis=1))
+
+    # Find default values according to nyquist criterion


kohr-h · 2017-01-01T21:16:50Z

odl/tomo/geometry/parallel.py

+    rho = np.max(np.linalg.norm(corners, axis=1))
+
+    # Find default values according to nyquist criterion
+    min_side = min(space.partition.cell_sides[:2])


~~.partition~~

kohr-h · 2017-01-01T21:19:01Z

odl/tomo/geometry/parallel.py

+        elif space.ndim == 3:
+            width = int(2 * rho * omega / np.pi) + 1
+            height = space.shape[2]
+            det_shape = [width, height]


No natural place to check for ndim in (2, 3), but an else case at some place would be nice.

Added afterwards

adler-j · 2017-01-02T09:11:20Z

Fixed the comments but realized that I may have an error in the computations. At one point I compute:

min_side = min(space.partition.cell_sides[:2])
omega = 2 * np.pi / min_side

Shouldn't this be

min_side = min(space.partition.cell_sides[:2])
omega = np.pi / min_side

?

Edit: I just checked with what Fourier transform gives, and it supports this idea. Anyway this would be great since then we'd only need 1/4 of the samples.

Edit2: With this change, the total number of samples also becomes about 4 times the samples in the input, which makes perfect sense given the Nyqvist criterion.

Edit3: Fixed this in a commit.

kohr-h

One spelling thing left, merge after that. Good catch with the definition of the sampling, I didn't go and check the reference.

kohr-h · 2017-01-02T22:53:39Z

odl/tomo/geometry/parallel.py

+    flexibility of the geometries, but simply want a geometry that works.
+
+    This default geometry gives a fully sampled sinogram according to the
+    nyquist criterion, which in general results in a very large number of


adler-j · 2017-01-03T10:07:08Z

Fixed comment, merge after CI.

ENH: Add parallel_beam_geometry, closes #711

00c462f

kohr-h requested changes Dec 20, 2016

View reviewed changes

chongchenmath reviewed Dec 20, 2016

View reviewed changes

adler-j added 2 commits December 25, 2016 03:10

ENH: Bug fixes to parallel_beam_geometry, add tests

16698cf

MAINT: Fix wrongly sized angle intervals in some examples

7ed4b4a

adler-j added status: review needed area: tomo labels Dec 28, 2016

kohr-h approved these changes Jan 1, 2017

View reviewed changes

MAINT: Style fixes to parallel_beam_geometry

93497c7

kohr-h approved these changes Jan 2, 2017

View reviewed changes

adler-j force-pushed the issue-711__default_parallel_geometry branch from 444e6e5 to d605c64 Compare January 3, 2017 10:07

BUG: fix error in parallel_beam_geometry overestimating number of points

d605c64

adler-j merged commit 9cc183d into master Jan 3, 2017

adler-j deleted the issue-711__default_parallel_geometry branch January 3, 2017 10:32

kohr-h removed the status: review needed label Jan 22, 2017



		def parallel_beam_geometry(space, angles=None, det_shape=None):
		"""Create default parallel beam geometry from space.


		The geometry returned by this function satisfies these conditions exactly.

		If the domain is 3 dimensional, the geometry is "separable", in that each

Conversation

adler-j commented Dec 20, 2016

Uh oh!

kohr-h left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

ozanoktem Dec 20, 2016 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

adler-j commented Dec 20, 2016

Uh oh!

Choose a reason for hiding this comment

Uh oh!

adler-j commented Dec 25, 2016

Uh oh!

kohr-h left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

adler-j commented Jan 2, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

kohr-h left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

adler-j commented Jan 3, 2017

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

ozanoktem Dec 20, 2016 •

edited

Loading

adler-j commented Jan 2, 2017 •

edited

Loading