docs: def template_ellipsoid(shape):
      r"""
      Returns an ellipsoid binary structure of a of the supplied radius that can be used as
      template input to the generalized hough transform.
      
      Parameters
      ----------
      shape : tuple of integers
          The main axes of the ellipsoid in voxel units.
      
      Returns
      -------
      template_sphere : ndarray
          A boolean array containing an ellipsoid.
      """
      # prepare template array
      template = numpy.zeros([int(x // 2 + (x % 2)) for x in shape], dtype=numpy.bool) # in odd shape cases, this will include the ellipses middle line, otherwise not
  
      # get real world offset to compute the ellipsoid membership
      rw_offset = []
      for s in shape:
          if int(s) % 2 == 0: rw_offset.append(0.5 - (s % 2) / 2.) # number before point is even 
          else: rw_offset.append(-1 * (s % int(s)) / 2.) # number before point is odd 
  
      # prepare an array containing the squares of the half axes to avoid computing inside the loop
      shape_pow = numpy.power(numpy.asarray(shape) / 2., 2)
  
      # we use the ellipse normal form to find all point in its surface as well as volume
      # e.g. for 2D, all voxels inside the ellipse (or on its surface) with half-axes a and b
      #      follow x^2/a^2 + y^2/b^2 <= 1; for higher dimensions accordingly
      # to not have to iterate over each voxel, we make use of the ellipsoids symmetry
      # and construct just a part of the whole ellipse here
      for idx in numpy.ndindex(template.shape):
          distance = sum((math.pow(coordinate + rwo, 2) / axes_pow for axes_pow, coordinate, rwo in zip(shape_pow, idx, rw_offset))) # plus once since ndarray is zero based, but real-world coordinates not
          if distance <= 1: template[idx] = True
          
      # we take now our ellipse part and flip it once along each dimension, concatenating it in each step
      # the slicers are constructed to flip in each step the current dimension i.e. to behave like arr[...,::-1,...]
      for i in range(template.ndim):
          slicers = [(slice(None, None, -1) if i == j else slice(None)) for j in range(template.ndim)]
          if 0 == int(shape[i]) % 2: # even case
              template = numpy.concatenate((template[slicers], template), i)
          else: # odd case, in which an overlap has to be created
              slicers_truncate = [(slice(None, -1) if i == j else slice(None)) for j in range(template.ndim)]
              template = numpy.concatenate((template[slicers][slicers_truncate], template), i)
  
      return template