diff --git a/python/mocpy/moc/moc.py b/python/mocpy/moc/moc.py
index a3cdce07..7cd968e5 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,
@@ -599,27 +600,44 @@ def from_fits_image(cls, hdu, max_norder, mask=None):
         return moc  # noqa: RET504
 
     @classmethod
-    def from_fits_images(cls, path_l, max_norder, hdu_index=0):
+    def from_fits_images(cls, path_l, max_norder, hdu_index=0, detect_hdu=False):
         """
         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)
+            Index of the the HDUs containing the image in each FITS file (default = 0)
+        detect_hdu : bool, optional
+            If set to True, all the HUD will be taken in account, and only the ones
+            corresponding to images will be kept. 'hdu_index' is ignored if this is set to
+            True. Default is False.
 
         Returns
         -------
         moc : `~mocpy.moc.MOC`
             The union of all the MOCs created from the paths found in ``path_l``.
         """
+        if not isinstance(path_l, list):
+            path_l = [path_l]
+        mocs = []
+        if detect_hdu:
+            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))
+            if len(mocs) == 1:
+                return mocs[0]
+            return mocs[0].union(*mocs[1:])
+
         moc = MOC.new_empty(max_norder)
         for filename in path_l:
             with fits.open(filename) as hdul:
diff --git a/python/mocpy/tests/test_moc.py b/python/mocpy/tests/test_moc.py
index f6538274..e68b05a2 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, detect_hdu=True)
 
 
 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, detect_hdu=True)
 
 
 @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(