Vectorize geometries #1059

kohr-h · 2017-07-04T22:50:37Z

~~Incomplete~~ PR for vectorization of geometry methods.

TODO:

Propagate check_bounds to detectors
Update doc for detectors
Improve repr of detectors
Decide on what to do with nonsense parameters (e.g. normalized when only one choice is valid, but the parameter is there for API consistency)
Solution: I removed normalized from ParallelGeometry.det_to_src, can't seem to remember what other cases I had in mind.
Replace super() calls
Solution: I did like suggested n Use super(cls, self) instead of super() or hard-coded class? #1061
Check if the custom all_contained() is really needed, seems fishy
Solution: Turned out a simple transpose did the job
Fix methods depending on two variables (as in 56393e6)
~~Add unit tests~~ Have already a bunch of doctests, and ASTRA vec geometry generation already uses vectorization, should be good enough?

kohr-h · 2017-07-05T21:59:15Z

Ready for review. As often before when working on geometries, there's a lot of new documentation, including doctests, which accounts for the biggest part of the PR. Documentation is way improved now IMO.

kohr-h · 2017-07-05T22:21:02Z

odl/tomo/geometry/detector.py

-        """Parametrization of the detector reference surface.
+        """Return the detector surface point corresponding to ``param``.
+
+        The surface point lies on a circle around ``radius * center_dir``


kohr-h · 2017-07-05T22:30:18Z

odl/tomo/geometry/geometry.py

+            infinity).
+            The shape of the returned array is as follows:
+
+            - ``mparam`` and ``dparam`` single: ``(ndim,)``


kohr-h · 2017-07-05T22:34:56Z

odl/tomo/geometry/parallel.py

-            contained in `motion_params`.
+        angle : float or `array-like`
+            Angle(s) in radians describing the counter-clockwise
+            rotation of the detector.


Needs more general formulation that includes multiple angles.

kohr-h · 2017-07-05T22:38:03Z

odl/tomo/geometry/parallel.py

+
+        # Need custom check here until #861 is in because arrays are
+        # assumed to be flat in the "point enumeration" axis
+        # TODO: replace when #861 is in


Missed this old check, needs to be replaced.

adler-j

One big overall comment about how we should handle arrays of input that needs to be discussed.

Other than that only minor stuff, good code quality as usual.

adler-j · 2017-07-06T08:24:42Z

odl/tomo/backends/astra_setup.py

@@ -431,7 +430,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


Wait wat. The next line guarantees that its a 1d detector, no?

Well not really. It just says that the whole geometry is embedded in a 2D space, the detector can be anything. This was probably supposed to be the duck-typed version of isinstance(geometry, Parallel2dGeometry).
But honestly if somebody decides to write a new geometry that doesn't subclass Parallel2dGeometry but still fulfills these conditions (whatever that may be), it would require a new conversion implementation anyway. So we might as well just do isinstance(geometry, Parallel2dGeometry).

adler-j · 2017-07-06T08:27:03Z

odl/tomo/geometry/conebeam.py

-        if angle not in self.motion_params:
-            raise ValueError('`angle` {} is not in the valid range {}'
-                             ''.format(angle, self.motion_params))
+        squeeze_out = np.isscalar(angle)


>>> np.isscalar(np.array(1)) False

is this intended?

Good point. I guess I would like to change this to np.shape(angle) == () which covers both cases.

adler-j · 2017-07-06T08:28:27Z

odl/tomo/geometry/conebeam.py

-        if np.linalg.norm(src_to_det_init) <= 1e-10:
-            raise ValueError('`src_to_det_init` norm {} too close to 0'
-                             ''.format(np.linalg.norm(src_to_det_init)))
+        if np.linalg.norm(src_to_det_init) == 0:


we removed the close-ness check for what reason?

I expect issues regarding scales. If somebody works on electron tomography where lengths are in the range of nanometers, and they decide to not (or forget to) normalize the vector, the old check will produce an unreasonable error.

With all these closeness checks on user input my attitude is to let the user figure out what to provide and how to deal with numerical accuracy issues, instead of us imposing some arbitrary bounds.

Agree. Good change.

adler-j · 2017-07-08T18:08:54Z

odl/tomo/geometry/conebeam.py

+
+        Returns
+        -------
+        axes : tuple of `numpy.ndarray`'s


Why return a tuple of arrays instead of simply one large array? Seems much more reasonable and makes the note much more simple

Makes sense, yes. I mainly did it in the original version to resemble the input for axes parameters, i.e., a sequence of vectors.
I'll change it to a big array with shape (num_axes, num_params, ndim).

Do we want the same for the axes property of detectors? Single array instead of tuple?

For consistency yeah i think so

OK I'll go with it

adler-j · 2017-07-09T13:07:46Z

odl/tomo/geometry/geometry.py

+        Other Parameters
+        ----------------
+        check_bounds : bool, optional
+            If ``True``, methods perform sanity checks on provided input


Somewhat ambiguous what parameters are refered to here, I'd read it as the parameters to this method.

Tried to come up with a better yet short explanation.

adler-j · 2017-07-09T13:26:33Z

odl/tomo/geometry/parallel.py

@@ -1014,9 +1196,9 @@ def frommatrix(cls, apart, dpart, init_matrix):
        # Use the standard constructor with these vectors
        axis, det_pos, det_axis_0, det_axis_1 = transformed_vecs
        if translation.size == 0:
-            kwargs = {}
+            pass


just make an if translation.size != 0:

Sure, this was the work-in-progress version just in case I want to switch back. Forgot to change to final.

adler-j · 2017-07-09T13:28:09Z

odl/tomo/util/utility.py

+    if np.isscalar(angles[0]):
+        return mat.squeeze()
+    else:
+        # Move third axis to the beginning


This is obvious from the code. But why do we do it?

We want mat[0] to be the first matrix etc. But you're right, the comment should say that rather than state the obvious.

adler-j · 2017-07-09T13:29:01Z

odl/tomo/util/utility.py

        raise ValueError('zero vector')

    result = np.zeros(vec.shape)
+    cond = np.any(vec[:, :2] != 0, axis=1)
+    if np.any(cond):


does this "any" make a difference in any way? I guess performance but I'd not over-optimize this right now.

Well this any goes along the first axis whose size is arbitrary (and can be large) so I think using np.any makes sense there.

My point was if the if statement is needed, or if we can simply go to the next line, where we get an empty assignment. I guess performance wise it shouldn't matter?

~~Ah, no that actually doesn't work. If the condition is empty, the array cannot be used for indexing.~~
Hm, now it works. Must have been some issue during the implementation. Will remove the checks.

adler-j · 2017-07-09T13:29:17Z

odl/tomo/util/utility.py

+        result[cond, 1] = vec[cond, 0]
+    if np.any(~cond):
+        result[~cond, 0] = 1
+    result / np.linalg.norm(result, axis=1, keepdims=True)


We're not afraid of zero division here, right?

This cannot happen. The rows in which there are nonzero entries in the first two coordinates we switch them, so the norm cannot be zero. In the other rows we insert a 1 so no way to get zero norm either (including rows that are all zeros).

adler-j · 2017-07-09T13:31:05Z

odl/trafos/backends/pyfftw_bindings.py

 except ImportError:
    PYFFTW_AVAILABLE = False
+else:
+    import pip


I don't like adding pip as a dependency to ODL, and further you mention this taking some performance.

I'd rather remove this, or as mentioned somewhere else explcitly check against something like

if parse_version(pyfftw.__version__) < parse_version('0.10.3.devSOMETHING'): warnings.warn('PyFFTW < 0.10.4 is known to cause problems with some ' 'ODL functionality, see issue #1002.', RuntimeWarning)

Yeah this was just a fix for myself to reduce annoyance. I'll just remove the commit.

kohr-h · 2017-07-13T22:27:41Z

odl/tomo/geometry/parallel.py

+            Unit vector(s) along which the detector is aligned.
+            If ``angles`` is a single pair (or triplet) of Euler angles,
+            the returned array has shape ``(2, 3)``, otherwise
+            is ``(2,) + broadcast(*angles).shape + (3,)``.


Not quite sure about this ordering. Opinion? Is it better to just add the extra dimensions to the left?

I'd say that they should go to the left, more consistent with other stuff

kohr-h · 2017-07-14T08:30:41Z

Rebased, ready for the final round.

adler-j

Looking mighty fine!

adler-j · 2017-07-11T11:31:10Z

odl/tomo/geometry/detector.py

+        This default (``space_ndim = ndim + 1``) can be overridden by
+        subclasses.
+        """
+        return self.ndim + 1


I would prefer if __init__ took an argument space_ndim, that would make overriding this easier and consistent with e.g. ndim

adler-j · 2017-07-11T11:31:40Z

odl/tomo/geometry/detector.py

-            raise NotImplementedError('normal not defined for ndim >= 3')
+            raise NotImplementedError(
+                'no default implementation of `surface_normal` available '
+                'for spaces with 4 or more dimensions')


Also 0d and 1d space ^.^

adler-j · 2017-07-11T11:40:00Z

odl/tomo/geometry/detector.py

-            zeros_shape = [m if n == 1 and m != 1 else 1
-                           for n, m in zip(axis_arr.shape, bcast_shape)]
-            return axis_arr + np.zeros(zeros_shape)
+            # TODO: use broadcast_to from Numpy when v1.10 is required


we have a utility for that

adler-j · 2017-07-11T11:41:42Z

odl/tomo/geometry/detector.py

+        # to broadcast along all axes
+        center_part = np.multiply.outer(1 - np.cos(param), self.center_dir)
+        tangent_part = np.multiply.outer(np.sin(param), self.tangent_at_0)
+        surf = self.radius * (center_part + tangent_part)
        if squeeze_out:


will this ever be needed given the current format?

adler-j · 2017-08-31T12:54:01Z

odl/tomo/geometry/detector.py

+            - ``param.shape[:-1] + (space_ndim,)`` otherwise.
+        """
+        if self.ndim == 1 and self.space_ndim == 2:
+            return -perpendicular_vector(self.surface_deriv(param))


We could extend perpendicular_vector to take several arguments as input and cover both of these cases

adler-j · 2017-08-31T14:41:08Z

odl/tomo/geometry/detector.py

+        else:
+            # Produce array of shape `param.shape + (ndim,)` by broadcasting
+            axis_slc = (None,) * param.ndim + (slice(None),)
+            # TODO: use broadcast_to from Numpy when v1.10 is required


Don't we have a util for this somewhere?

No, it's currently buried somewhere internally. It only exists on the tensor branch.

adler-j · 2017-08-31T15:23:38Z

odl/tomo/geometry/parallel.py

+            Unit vector(s) along which the detector is aligned.
+            If ``angles`` is a single pair (or triplet) of Euler angles,
+            the returned array has shape ``(2, 3)``, otherwise
+            is ``(2,) + broadcast(*angles).shape + (3,)``.


I'd say that they should go to the left, more consistent with other stuff

kohr-h · 2017-09-26T16:45:17Z

Added the space_ndim thing to the abstract Detector class and reordered the det_axes and surface_deriv things.

I'll rebase and if CI is happy, I'll merge.

The following changes are part of this commit: - All detector methods are now vectorized. - Detectors now take a `check_bounds` parameter to switch off bounds checking for cases where milliseconds count. - The majority of methods has doctests now. - `Geometry` has a new method `surface_normal` which also replaces `normal` in flat detectors. - Implementation of `repr` for detectors is improved. - `CircleSectionDetector` has gained some flexibility (can be defined with a custom circle radius now), plus doctests.

adler-j · 2017-09-26T16:57:28Z

Great work!

kohr-h force-pushed the vectorize_geometries branch from 1e65f17 to e4176b1 Compare July 5, 2017 14:59

kohr-h changed the title ~~WIP: Vectorize geometries~~ Vectorize geometries Jul 5, 2017

kohr-h added the status: review needed label Jul 5, 2017

kohr-h commented Jul 5, 2017

View reviewed changes

adler-j reviewed Jul 9, 2017

View reviewed changes

kohr-h commented Jul 13, 2017

View reviewed changes

kohr-h force-pushed the vectorize_geometries branch 3 times, most recently from 58b021d to 116a545 Compare July 14, 2017 08:30

kohr-h force-pushed the vectorize_geometries branch from 116a545 to c412be8 Compare July 14, 2017 08:55

adler-j approved these changes Aug 31, 2017

View reviewed changes

kohr-h force-pushed the vectorize_geometries branch from c412be8 to 05777bf Compare September 5, 2017 11:10

kohr-h force-pushed the vectorize_geometries branch from 1fa0054 to 6292fa7 Compare September 26, 2017 14:51

Holger Kohr added 11 commits September 26, 2017 18:47

ENH: vectorize geometry methods

858de5a

ENH: vectorize perpendicular_vector

92e0316

MAINT: remove FlatDetector class

d50c459

MAINT: fix test failures of geometries

79004fa

MAINT: propagate check_bounds from geometry to detector

ca6373d

BUG: fix det_point_position and add a doctest

948bbdd

MAINT: replace super() by explicit super(cls, self)

1e3c496

ENH: vectorize geometry functions with 2 variables

696a01e

MAINT: fix minor things in geometry and detector

9de30be

ENH: add broadcasting to geometry methods

a7c6dc2

Holger Kohr added 4 commits September 26, 2017 18:55

ENH: add utility to check bounds for geometry methods

8ed437d

MAINT: add space_ndim arg to abstract Detector class

b843c16

MAINT: change det_axes in 3D to return (N, 2, 3) shape

030da5b

MAINT: change also surface_deriv from (2, N, 3) to (N, 2, 3)

4e42a1b

kohr-h force-pushed the vectorize_geometries branch from b72b65a to 4e42a1b Compare September 26, 2017 16:55

kohr-h merged commit ee4ad6d into odlgroup:master Sep 26, 2017

kohr-h deleted the vectorize_geometries branch September 26, 2017 17:04

kohr-h removed the status: review needed label Oct 6, 2017

kohr-h mentioned this pull request Nov 13, 2017

Input parameter checks in geometry methods #296

Closed

Vectorize geometries #1059

Vectorize geometries #1059

Conversation

kohr-h commented Jul 4, 2017 • edited Loading

TODO:

kohr-h commented Jul 5, 2017

kohr-h Jul 5, 2017 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

adler-j left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

adler-j Jul 9, 2017 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kohr-h commented Jul 14, 2017

adler-j left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kohr-h commented Sep 26, 2017

adler-j commented Sep 26, 2017

kohr-h commented Jul 4, 2017 •

edited

Loading

kohr-h Jul 5, 2017 •

edited

Loading

adler-j Jul 9, 2017 •

edited

Loading