diff --git a/CHANGELOG.md b/CHANGELOG.md
index dfbf7aef..3807f9ba 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -21,6 +21,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 * Creation of a MOC from a zone (defined by min/max ra and dec)`MOC.from_zone`
 * Creation of a single MOC from a lot of small cones is faster with the new option in
 `MOC.from_cones`: the keyword 'union_strategy' can now take the value 'small_cones'.
+* `MOC.from_fits_images` can now loop through the HDUList to only keep images with the
+parameter `hdu_index` set to -1.
 
 ## [0.16.2]
 
diff --git a/python/mocpy/moc/moc.py b/python/mocpy/moc/moc.py
index a3cdce07..895c4aea 100644
--- a/python/mocpy/moc/moc.py
+++ b/python/mocpy/moc/moc.py
@@ -499,7 +499,7 @@ def from_fits_image(cls, hdu, max_norder, mask=None):
 
         Info
         ----
-        When giving a mask the MOC computed will only take into account the center
+        When giving a mask, the MOC computed will only take into account the center
         of the image pixels and not the whole pixel borders.
         This leads to an approximate resulting MOC.
 
@@ -510,12 +510,12 @@ def from_fits_image(cls, hdu, max_norder, mask=None):
         max_norder : int
             The moc resolution.
         mask : `numpy.ndarray`, optional
-            A boolean array of the same shape of the image where True valued pixels are part of
+            A boolean array of the same shape than the image where True valued pixels are part of
             the final MOC and False valued pixels are not.
 
         Returns
         -------
-        moc : `~mocpy.moc.MOC`
+        `~mocpy.moc.MOC`
             The resulting MOC.
         """
         # Only take the first HDU
@@ -573,7 +573,7 @@ def from_fits_image(cls, hdu, max_norder, mask=None):
             )  # in steradians
 
             # Division by 0 case
-            if px_sky_area == 0.0:
+            if px_sky_area == 0:
                 healpix_order_computed = False
             else:
                 depth_px = np.uint8(
@@ -585,7 +585,8 @@ def from_fits_image(cls, hdu, max_norder, mask=None):
 
         if not healpix_order_computed:
             warnings.warn(
-                "MOC precision HEALPix order could not be determined because sky coordinates from the corners of the image has not have been correctly retrieved. "
+                "MOC precision HEALPix order could not be determined because sky coordinates "
+                "from the corners of the image has not have been correctly retrieved. "
                 "Therefore MOC precision will be set to max_norder",
                 UserWarning,
                 stacklevel=2,
@@ -603,33 +604,46 @@ def from_fits_images(cls, path_l, max_norder, hdu_index=0):
         """
         Load a MOC from a set of FITS file images.
 
-        Assumes the data of the image is stored in the first HDU of the FITS file.
-        Please call `~mocpy.moc.MOC.from_fits_image` for passing another hdu than the first one.
-
         Parameters
         ----------
         path_l : [str]
-            A list of path where the fits image are located.
+            A list of path where the fits images are located.
         max_norder : int
             The MOC resolution.
-        hdu_index : int
-            Index of the the HDU containing the image in each FITS file (default = 0)
+        hdu_index : int, optional
+            Index of the the HDUs containing the image in each FITS file (default = 0)
+            If set to -1, all the HUD will be taken in account, and only the ones
+            corresponding to images will be kept.
 
         Returns
         -------
         moc : `~mocpy.moc.MOC`
             The union of all the MOCs created from the paths found in ``path_l``.
         """
-        moc = MOC.new_empty(max_norder)
-        for filename in path_l:
-            with fits.open(filename) as hdul:
-                current_moc = MOC.from_fits_image(
-                    hdu=hdul[hdu_index],
-                    max_norder=max_norder,
-                )
-                moc = moc.union(current_moc)
+        if not isinstance(path_l, list):
+            path_l = [path_l]
+
+        mocs = []
+        if hdu_index == -1:
+            for filename in path_l:
+                with fits.open(filename) as hdul:
+                    for hdu in hdul:
+                        if (
+                            isinstance(hdu, (fits.ImageHDU, fits.PrimaryHDU))
+                            and len(hdu.data.shape) == 2
+                        ):
+                            mocs.append(MOC.from_fits_image(hdu, max_norder))
+
+        else:
+            for filename in path_l:
+                with fits.open(filename) as hdul:
+                    mocs.append(MOC.from_fits_image(hdu=hdul[hdu_index],
+                                                    max_norder=max_norder,))
+
+        if len(mocs) == 1:
+            return mocs[0]
+        return mocs[0].union(*mocs[1:]) # this is the fastest way to do multi union
 
-        return moc
 
     @classmethod
     def from_vizier_table(cls, table_id, max_depth=None, nside=None):
diff --git a/python/mocpy/tests/test_moc.py b/python/mocpy/tests/test_moc.py
index f6538274..862425b9 100644
--- a/python/mocpy/tests/test_moc.py
+++ b/python/mocpy/tests/test_moc.py
@@ -344,8 +344,7 @@ def test_moc_consistent_with_aladin():
 
 def test_moc_from_fits_images():
     image_path = "resources/image_with_mask.fits.gz"
-
-    MOC.from_fits_images([image_path], max_norder=15)
+    MOC.from_fits_images([image_path], max_norder=15, hdu_index=-1)
 
 
 def test_from_fits_images_2():
@@ -353,7 +352,7 @@ def test_from_fits_images_2():
 
 
 def test_from_fits_image_without_cdelt():
-    MOC.from_fits_images(["resources/horsehead.fits"], max_norder=15)
+    MOC.from_fits_images(["resources/horsehead.fits"], max_norder=5, hdu_index=-1)
 
 
 @pytest.fixture()
@@ -853,8 +852,7 @@ def test_moc_complement_consistency():
 
 
 def test_from_fits_old():
-    moc = MOC.from_fits("resources/V_147_sdss12.moc.fits")
-    assert moc.complement().complement() == moc
+    MOC.from_fits("resources/V_147_sdss12.moc.fits")
 
 
 @pytest.mark.parametrize(