Make reproject_to_healpix Compatible with Multiple images with the same WCS

reproject_to_healpix can be made faster if it allows for the same input format as retroject_interp, where states:

> If the data array contains more dimensions than are described by the input header or WCS, the extra dimensions (assumed to be the first dimensions) are taken to represent multiple images with the same coordinate information. The coordinate transformation will be computed once and then each image will be reprojected, offering a speedup over reprojecting each image individually.

I have multiple images, formatted in a numpy array where each index of the first array corresponds to different energy ranges, in which each image has the same WCS. I currently have to iterate over the first axis to get the HEALPIX map for each energy range and then combine all the maps appropriately. Since I know that these images all correspond to the same WCS, it would be ideal to save the computation of the coordinate transformation.

To be explicit my code looks something like this:
```
new_array = np.zeros((energy_nbins, healpix_npixels))
for e in range(energy_nbins):
array, footprint = reproject_to_healpix(
(image[e, :, :], wcs), coordsys,
nside=healpix_nside)
new_array[e, :] = array
```

And it would be great for the reprojection call to be:
```
new_array= reproject_to_healpix((image, wcs), coordsys, nside=healpix_nside, footprint=False)
```

Currently, my version of reproject is 0.14.0.

Read the reproject_to_healpix entry point and compare its input handling with retroject_interp; check how leading dimensions and the WCS are currently treated. Done means the shown (image, wcs) call accepts multiple leading energy-image dimensions, reuses the coordinate transformation, and preserves the requested footprint=False output.

Written by the indexing model from the issue text.