code: 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.0
              )  # number before point is even
          else:
              rw_offset.append(
                  -1 * (s % int(s)) / 2.0
              )  # 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.0, 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
  
imports: import math
  import numpy
  from .utilities import pad
label_desc: Use image node_id 4 for calling a global function with example usage: template_ellipsoid(shape) and returns: template