def det_point_position(self, mparam, dparam):
        """Return the detector point at ``(mparam, dparam)``.

        The position is computed as follows::

            pos = refpoint(mparam) +
                  rotation_matrix(mparam).dot(detector.surface(dparam))

        In other words, the motion parameter ``mparam`` is used to move the
        detector reference point, and the detector parameter ``dparam``
        defines an intrinsic shift that is added to the reference point.

        Parameters
        ----------
        mparam : `array-like` or sequence
            Motion parameter(s) at which to evaluate. If
            ``motion_params.ndim >= 2``, a sequence of that length must be
            provided.
        dparam : `array-like` or sequence
            Detector parameter(s) at which to evaluate. If
            ``det_params.ndim >= 2``, a sequence of that length must be
            provided.

        Returns
        -------
        pos : `numpy.ndarray`
            Vector(s) pointing from the origin to the detector point.
            The shape of the returned array is obtained from the
            (broadcast) shapes of ``mparam`` and ``dparam``, and
            broadcasting is supported within both parameters and between
            them. The precise definition of the shape is
            ``broadcast(bcast_mparam, bcast_dparam).shape + (ndim,)``,
            where ``bcast_mparam`` is

            - ``mparam`` if `motion_params` is 1D,
            - ``broadcast(*mparam)`` otherwise,

            and ``bcast_dparam`` defined analogously.

        Examples
        --------
        The method works with single parameter values, in which case
        a single vector is returned:

        >>> apart = odl.uniform_partition(0, np.pi, 10)
        >>> dpart = odl.uniform_partition(-1, 1, 20)
        >>> geom = odl.tomo.Parallel2dGeometry(apart, dpart)
        >>> geom.det_point_position(0, 0)  # (0, 1) + 0 * (1, 0)
        array([ 0.,  1.])
        >>> geom.det_point_position(0, 1)  # (0, 1) + 1 * (1, 0)
        array([ 1.,  1.])
        >>> pt = geom.det_point_position(np.pi / 2, 0)  # (-1, 0) + 0 * (0, 1)
        >>> np.allclose(pt, [-1, 0])
        True
        >>> pt = geom.det_point_position(np.pi / 2, 1)  # (-1, 0) + 1 * (0, 1)
        >>> np.allclose(pt, [-1, 1])
        True

        Both variables support vectorized calls, i.e., stacks of
        parameters can be provided. The order of axes in the output (left
        of the ``ndim`` axis for the vector dimension) corresponds to the
        order of arguments:

        >>> geom.det_point_position(0, [-1, 0, 0.5, 1])
        array([[-1. ,  1. ],
               [ 0. ,  1. ],
               [ 0.5,  1. ],
               [ 1. ,  1. ]])
        >>> pts = geom.det_point_position([0, np.pi / 2, np.pi], 0)
        >>> np.allclose(pts, [[0, 1],
        ...                   [-1, 0],
        ...                   [0, -1]])
        True
        >>> # Providing 3 pairs of parameters, resulting in 3 vectors
        >>> pts = geom.det_point_position([0, np.pi / 2, np.pi],
        ...                               [-1, 0, 1])
        >>> pts[0]  # Corresponds to angle = 0, dparam = -1
        array([-1.,  1.])
        >>> pts.shape
        (3, 2)
        >>> # Pairs of parameters arranged in arrays of same size
        >>> geom.det_point_position(np.zeros((4, 5)), np.zeros((4, 5))).shape
        (4, 5, 2)
        >>> # "Outer product" type evaluation using broadcasting
        >>> geom.det_point_position(np.zeros((4, 1)), np.zeros((1, 5))).shape
        (4, 5, 2)

        More complicated 3D geometry with 2 angle variables and 2 detector
        variables:

        >>> apart = odl.uniform_partition([0, 0], [np.pi, 2 * np.pi],
        ...                               (10, 20))
        >>> dpart = odl.uniform_partition([-1, -1], [1, 1], (20, 20))
        >>> geom = odl.tomo.Parallel3dEulerGeometry(apart, dpart)
        >>> # 2 values for each variable, resulting in 2 vectors
        >>> angles = ([0, np.pi / 2], [0, np.pi])
        >>> dparams = ([-1, 0], [-1, 0])
        >>> pts = geom.det_point_position(angles, dparams)
        >>> pts[0]  # Corresponds to angle = (0, 0), dparam = (-1, -1)
        array([-1.,  1., -1.])
        >>> pts.shape
        (2, 3)
        >>> # 4 x 5 parameters for both
        >>> angles = dparams = (np.zeros((4, 5)), np.zeros((4, 5)))
        >>> geom.det_point_position(angles, dparams).shape
        (4, 5, 3)
        >>> # Broadcasting angles to shape (4, 5, 1, 1)
        >>> angles = (np.zeros((4, 1, 1, 1)), np.zeros((1, 5, 1, 1)))
        >>> # Broadcasting dparams to shape (1, 1, 6, 7)
        >>> dparams = (np.zeros((1, 1, 6, 1)), np.zeros((1, 1, 1, 7)))
        >>> # Total broadcast parameter shape is (4, 5, 6, 7)
        >>> geom.det_point_position(angles, dparams).shape
        (4, 5, 6, 7, 3)
        """
        # Always call the downstream methods with vectorized arguments
        # to be able to reliably manipulate the final axes of the result
        if self.motion_params.ndim == 1:
            squeeze_mparam = (np.shape(mparam) == ())
            mparam = np.array(mparam, dtype=float, copy=False, ndmin=1)
            matrix = self.rotation_matrix(mparam)  # shape (m, ndim, ndim)
        else:
            squeeze_mparam = (np.broadcast(*mparam).shape == ())
            mparam = tuple(np.array(a, dtype=float, copy=False, ndmin=1)
                           for a in mparam)
            matrix = self.rotation_matrix(mparam)  # shape (m, ndim, ndim)

        if self.det_params.ndim == 1:
            squeeze_dparam = (np.shape(dparam) == ())
            dparam = np.array(dparam, dtype=float, copy=False, ndmin=1)
        else:
            squeeze_dparam = (np.broadcast(*dparam).shape == ())
            dparam = tuple(np.array(p, dtype=float, copy=False, ndmin=1)
                           for p in dparam)

        surf = self.detector.surface(dparam)  # shape (d, ndim)

        # Perform matrix-vector multiplication along the last axis of both
        # `matrix` and `surf` while "zipping" all axes that do not
        # participate in the matrix-vector product. In other words, the axes
        # are labelled
        # [0, 1, ..., r-1, r, r+1] for `matrix` and
        # [0, 1, ..., r-1, r+1] for `surf`, and the output axes are set to
        # [0, 1, ..., r-1, r]. This automatically supports broadcasting
        # along the axes 0, ..., r-1.
        matrix_axes = list(range(matrix.ndim))
        surf_axes = list(range(matrix.ndim - 2)) + [matrix_axes[-1]]
        out_axes = list(range(matrix.ndim - 1))
        det_part = np.einsum(matrix, matrix_axes, surf, surf_axes, out_axes)

        refpt = self.det_refpoint(mparam)
        det_pt_pos = refpt + det_part
        if squeeze_mparam and squeeze_dparam:
            det_pt_pos = det_pt_pos.squeeze()

        return det_pt_pos
Return the detector point at ``(mparam, dparam)``.

        The position is computed as follows::

            pos = refpoint(mparam) +
                  rotation_matrix(mparam).dot(detector.surface(dparam))

        In other words, the motion parameter ``mparam`` is used to move the
        detector reference point, and the detector parameter ``dparam``
        defines an intrinsic shift that is added to the reference point.

        Parameters
        ----------
        mparam : `array-like` or sequence
            Motion parameter(s) at which to evaluate. If
            ``motion_params.ndim >= 2``, a sequence of that length must be
            provided.
        dparam : `array-like` or sequence
            Detector parameter(s) at which to evaluate. If
            ``det_params.ndim >= 2``, a sequence of that length must be
            provided.

        Returns
        -------
        pos : `numpy.ndarray`
            Vector(s) pointing from the origin to the detector point.
            The shape of the returned array is obtained from the
            (broadcast) shapes of ``mparam`` and ``dparam``, and
            broadcasting is supported within both parameters and between
            them. The precise definition of the shape is
            ``broadcast(bcast_mparam, bcast_dparam).shape + (ndim,)``,
            where ``bcast_mparam`` is

            - ``mparam`` if `motion_params` is 1D,
            - ``broadcast(*mparam)`` otherwise,

            and ``bcast_dparam`` defined analogously.

        Examples
        --------
        The method works with single parameter values, in which case
        a single vector is returned:

        >>> apart = odl.uniform_partition(0, np.pi, 10)
        >>> dpart = odl.uniform_partition(-1, 1, 20)
        >>> geom = odl.tomo.Parallel2dGeometry(apart, dpart)
        >>> geom.det_point_position(0, 0)  # (0, 1) + 0 * (1, 0)
        array([ 0.,  1.])
        >>> geom.det_point_position(0, 1)  # (0, 1) + 1 * (1, 0)
        array([ 1.,  1.])
        >>> pt = geom.det_point_position(np.pi / 2, 0)  # (-1, 0) + 0 * (0, 1)
        >>> np.allclose(pt, [-1, 0])
        True
        >>> pt = geom.det_point_position(np.pi / 2, 1)  # (-1, 0) + 1 * (0, 1)
        >>> np.allclose(pt, [-1, 1])
        True

        Both variables support vectorized calls, i.e., stacks of
        parameters can be provided. The order of axes in the output (left
        of the ``ndim`` axis for the vector dimension) corresponds to the
        order of arguments:

        >>> geom.det_point_position(0, [-1, 0, 0.5, 1])
        array([[-1. ,  1. ],
               [ 0. ,  1. ],
               [ 0.5,  1. ],
               [ 1. ,  1. ]])
        >>> pts = geom.det_point_position([0, np.pi / 2, np.pi], 0)
        >>> np.allclose(pts, [[0, 1],
        ...                   [-1, 0],
        ...                   [0, -1]])
        True
        >>> # Providing 3 pairs of parameters, resulting in 3 vectors
        >>> pts = geom.det_point_position([0, np.pi / 2, np.pi],
        ...                               [-1, 0, 1])
        >>> pts[0]  # Corresponds to angle = 0, dparam = -1
        array([-1.,  1.])
        >>> pts.shape
        (3, 2)
        >>> # Pairs of parameters arranged in arrays of same size
        >>> geom.det_point_position(np.zeros((4, 5)), np.zeros((4, 5))).shape
        (4, 5, 2)
        >>> # "Outer product" type evaluation using broadcasting
        >>> geom.det_point_position(np.zeros((4, 1)), np.zeros((1, 5))).shape
        (4, 5, 2)

        More complicated 3D geometry with 2 angle variables and 2 detector
        variables:

        >>> apart = odl.uniform_partition([0, 0], [np.pi, 2 * np.pi],
        ...                               (10, 20))
        >>> dpart = odl.uniform_partition([-1, -1], [1, 1], (20, 20))
        >>> geom = odl.tomo.Parallel3dEulerGeometry(apart, dpart)
        >>> # 2 values for each variable, resulting in 2 vectors
        >>> angles = ([0, np.pi / 2], [0, np.pi])
        >>> dparams = ([-1, 0], [-1, 0])
        >>> pts = geom.det_point_position(angles, dparams)
        >>> pts[0]  # Corresponds to angle = (0, 0), dparam = (-1, -1)
        array([-1.,  1., -1.])
        >>> pts.shape
        (2, 3)
        >>> # 4 x 5 parameters for both
        >>> angles = dparams = (np.zeros((4, 5)), np.zeros((4, 5)))
        >>> geom.det_point_position(angles, dparams).shape
        (4, 5, 3)
        >>> # Broadcasting angles to shape (4, 5, 1, 1)
        >>> angles = (np.zeros((4, 1, 1, 1)), np.zeros((1, 5, 1, 1)))
        >>> # Broadcasting dparams to shape (1, 1, 6, 7)
        >>> dparams = (np.zeros((1, 1, 6, 1)), np.zeros((1, 1, 1, 7)))
        >>> # Total broadcast parameter shape is (4, 5, 6, 7)
        >>> geom.det_point_position(angles, dparams).shape
        (4, 5, 6, 7, 3)