docs: def construct_zernike_polynomials(x, y, zernike_indexes, mask=None, weight=None):
      """Return the zerike polynomials for all objects in an image
      
      x - the X distance of a point from the center of its object
      y - the Y distance of a point from the center of its object
      zernike_indexes - an Nx2 array of the Zernike polynomials to be computed.
      mask - a mask with same shape as X and Y of the points to consider
      weight - weightings of points with the same shape as X and Y (default
               weight on each point is 1).
      
      returns a height x width x N array of complex numbers which are the
      e^i portion of the sine and cosine of the Zernikes
      """
      if x.shape != y.shape:
          raise ValueError("X and Y must have the same shape")
      if mask is None:
          pass
      elif mask.shape != x.shape:
          raise ValueError("The mask must have the same shape as X and Y")
      else:
          x = x[mask]
          y = y[mask]
          if weight is not None:
              weight = weight[mask]
      lut = construct_zernike_lookuptable(zernike_indexes) # precompute poly. coeffs.
      nzernikes = zernike_indexes.shape[0]
      # compute radii
      r_square = np.square(x) # r_square = x**2
      np.add(r_square, np.square(y), out=r_square) # r_square = x**2 + y**2
      # z = y + 1j*x
      # each Zernike polynomial is poly(r)*(r**m * np.exp(1j*m*phi)) ==
      #                            poly(r)*(y + 1j*x)**m
      z = np.empty(x.shape, np.complex)
      np.copyto(z.real, y)
      np.copyto(z.imag, x)
      # preallocate buffers
      s = np.empty_like(x)
      zf = np.zeros((nzernikes,) + x.shape, np.complex)
      z_pows = {}
      for idx, (n, m) in enumerate(zernike_indexes):
          s[:]=0
          if not m in z_pows:
              if m == 0:
                  z_pows[m] = np.complex(1.0)
              else:
                  z_pows[m] = z if m == 1 else (z ** m)
          z_pow = z_pows[m]
          # use Horner scheme
          for k in range((n-m)//2+1):
              s *= r_square
              s += lut[idx, k]
          s[r_square>1]=0
          if weight is not None:
              s *= weight.astype(s.dtype)
          if m == 0:
              np.copyto(zf[idx], s) # zf[idx] = s
          else:
              np.multiply(s, z_pow, out=zf[idx]) # zf[idx] = s*exp_term
      
      if mask is None:
          result = zf.transpose( tuple(range(1, 1+x.ndim)) + (0, ))
      else:
          result = np.zeros( mask.shape + (nzernikes,), np.complex)
          result[mask] = zf.transpose( tuple(range(1, 1 + x.ndim)) + (0, ))
      return result