def surface(self, param):
        """Return the detector surface point corresponding to ``param``.

        For parameter value ``p``, the surface point is given by ::

            surf = p[0] * axes[0] + p[1] * axes[1]

        Parameters
        ----------
        param : `array-like` or sequence
            Parameter value(s) at which to evaluate. A sequence of
            parameters must have length 2.

        Returns
        -------
        point : `numpy.ndarray`
            Vector(s) pointing from the origin to the detector surface
            point at ``param``.
            If ``param`` is a single parameter, the returned array has
            shape ``(3,)``, otherwise ``broadcast(*param).shape + (3,)``.

        Examples
        --------
        The method works with a single parameter, resulting in a single
        vector:

        >>> part = odl.uniform_partition([0, 0], [1, 1], (10, 10))
        >>> det = Flat2dDetector(part, axes=[(1, 0, 0), (0, 0, 1)])
        >>> det.surface([0, 0])
        array([ 0.,  0.,  0.])
        >>> det.surface([0, 1])
        array([ 0.,  0.,  1.])
        >>> det.surface([1, 1])
        array([ 1.,  0.,  1.])

        It is also vectorized, i.e., it can be called with multiple
        parameters at once (or n-dimensional arrays of parameters):

        >>> # 3 pairs of parameters, resulting in 3 vectors
        >>> det.surface([[0, 0, 1],
        ...              [0, 1, 1]])
        array([[ 0.,  0.,  0.],
               [ 0.,  0.,  1.],
               [ 1.,  0.,  1.]])
        >>> # Pairs of parameters in a (4, 5) array each
        >>> param = (np.zeros((4, 5)), np.zeros((4, 5)))
        >>> det.surface(param).shape
        (4, 5, 3)
        >>> # Using broadcasting for "outer product" type result
        >>> param = (np.zeros((4, 1)), np.zeros((1, 5)))
        >>> det.surface(param).shape
        (4, 5, 3)
        """
        squeeze_out = (np.broadcast(*param).shape == ())
        param_in = param
        param = tuple(np.array(p, dtype=float, copy=False, ndmin=1)
                      for p in param)
        if self.check_bounds and not is_inside_bounds(param, self.params):
            raise ValueError('`param` {} not in the valid range '
                             '{}'.format(param_in, self.params))

        # Compute outer product of the i-th spatial component of the
        # parameter and sum up the contributions
        surf = sum(np.multiply.outer(p, ax) for p, ax in zip(param, self.axes))
        if squeeze_out:
            surf = surf.squeeze()

        return surf