diff --git a/.gitignore b/.gitignore
index 047e028..8fa32e7 100644
--- a/.gitignore
+++ b/.gitignore
@@ -26,3 +26,5 @@ plots/*
 *.html
 analysis_config.json
 .Rhistory
+
+.vscode/*
diff --git a/README.md b/README.md
index b494859..3e9bb80 100644
--- a/README.md
+++ b/README.md
@@ -82,6 +82,21 @@ sq.pl.spatial_scatter(
 
 ![Image](docs/_static/img/example_img_sq.png)
 
+#### Convert AnnData objects to ion image arrays
+
+```python
+from metaspace_converter import metaspace_to_anndata, anndata_to_image_array
+
+# Download data
+adata2 = metaspace_to_anndata(dataset_id="2023-11-14_21h58m39s", fdr=0.1)
+
+ion_images = anndata_to_image_array(adata2)
+
+# 6 ion images of shape 61x78
+print(ion_images.shape)
+# > (6, 61, 78)
+```
+
 ### SpatialData
 
 Here using a reversed colormap which better represents intense values on bright background.
diff --git a/docs/api.rst b/docs/api.rst
index 07b55b6..96da798 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -5,6 +5,6 @@ metaspace_converter
 -------------------
 
 .. automodule:: metaspace_converter
-   :members: metaspace_to_anndata, metaspace_to_spatialdata
+   :members:
    :undoc-members:
    :show-inheritance:
diff --git a/metaspace2anndata/__pycache__/__init__.cpython-310.pyc b/metaspace2anndata/__pycache__/__init__.cpython-310.pyc
deleted file mode 100644
index 7ad2e89..0000000
Binary files a/metaspace2anndata/__pycache__/__init__.cpython-310.pyc and /dev/null differ
diff --git a/metaspace2anndata/__pycache__/converter.cpython-310.pyc b/metaspace2anndata/__pycache__/converter.cpython-310.pyc
deleted file mode 100644
index 1c2c834..0000000
Binary files a/metaspace2anndata/__pycache__/converter.cpython-310.pyc and /dev/null differ
diff --git a/metaspace2anndata/__pycache__/utils.cpython-310.pyc b/metaspace2anndata/__pycache__/utils.cpython-310.pyc
deleted file mode 100644
index 368da7f..0000000
Binary files a/metaspace2anndata/__pycache__/utils.cpython-310.pyc and /dev/null differ
diff --git a/metaspace_converter/__init__.py b/metaspace_converter/__init__.py
index cb4cdb8..87c9805 100644
--- a/metaspace_converter/__init__.py
+++ b/metaspace_converter/__init__.py
@@ -1,4 +1,5 @@
+from metaspace_converter.anndata_to_array import anndata_to_image_array
 from metaspace_converter.to_anndata import metaspace_to_anndata
 from metaspace_converter.to_spatialdata import metaspace_to_spatialdata
 
-__all__ = ["metaspace_to_anndata", "metaspace_to_spatialdata"]
+__all__ = ["metaspace_to_anndata", "metaspace_to_spatialdata", "anndata_to_image_array"]
diff --git a/metaspace_converter/anndata_to_array.py b/metaspace_converter/anndata_to_array.py
new file mode 100644
index 0000000..61d98a2
--- /dev/null
+++ b/metaspace_converter/anndata_to_array.py
@@ -0,0 +1,65 @@
+import numpy as np
+from anndata import AnnData
+
+from metaspace_converter.constants import COL, METASPACE_KEY, X, Y
+from metaspace_converter.to_anndata import all_image_pixel_coordinates
+
+
+def _check_pixel_coordinates(adata: AnnData) -> bool:
+    sorted_obs = adata.obs.sort_values([COL.ion_image_pixel_y, COL.ion_image_pixel_x])
+
+    pixel_list = sorted_obs[[COL.ion_image_pixel_y, COL.ion_image_pixel_x]].values
+
+    img_size = adata.uns[METASPACE_KEY]["image_size"]
+    required_pixels = all_image_pixel_coordinates((img_size[Y], img_size[X]))
+
+    return np.all(np.equal(pixel_list, required_pixels))
+
+
+def anndata_to_image_array(adata: AnnData) -> np.ndarray:
+    """
+    Extracts an array of ion images from an AnnData object
+    (that has been generated through the ``metaspace_to_anndata`` function).
+
+    Args:
+        adata: An AnnData object.
+
+    Returns:
+        A three-dimensional Numpy array in the following shape
+
+            * Dimension 0: Number of ion images in the order of ``adata.var_names``
+            * Dimension 1: Image height ``adata.uns["metaspace"]["image_size"]["y"]``
+            * Dimension 2: Image width ``adata.uns["metaspace"]["image_size"]["x"]``
+
+    Raises:
+        ValueError: If the AnnData object has been modified.
+            E.g. Pixel have been removed/added and the number of pixels
+            and their coordinates do not match the original image dimensions.
+    """
+
+    pixel_array = adata.X.transpose().copy()
+    img_size = adata.uns[METASPACE_KEY]["image_size"]
+
+    # Check if image dimensions are okay
+    if img_size[X] * img_size[Y] != pixel_array.shape[1]:
+        raise ValueError("Number of observations does not match the original image dimensions")
+
+    # Check if all pixels are available
+    if not _check_pixel_coordinates(adata):
+        raise ValueError("Not all pixels for ion images are available")
+
+    # Sort indices, in case of modified order of pixels (obs)
+    image_sorting = adata.obs.sort_values(
+        [COL.ion_image_pixel_y, COL.ion_image_pixel_x]
+    ).index.values.astype(int)
+    pixel_array = pixel_array[:, image_sorting]
+
+    image_array = pixel_array.reshape(
+        (
+            pixel_array.shape[0],
+            adata.obs[COL.ion_image_pixel_y].max() + 1,
+            adata.obs[COL.ion_image_pixel_x].max() + 1,
+        )
+    )
+
+    return image_array
diff --git a/metaspace_converter/constants.py b/metaspace_converter/constants.py
new file mode 100644
index 0000000..99bf75e
--- /dev/null
+++ b/metaspace_converter/constants.py
@@ -0,0 +1,30 @@
+from dataclasses import dataclass
+from typing import Final
+
+SPATIAL_KEY: Final = "spatial"
+METASPACE_KEY: Final = "metaspace"
+OBS_INDEX_NAME: Final = "ion_image_pixel"
+VAR_INDEX_NAME: Final = "formula_adduct"
+Shape2d = tuple[int, int]
+
+
+@dataclass
+class COL:
+    ion_image_shape_y: Final = "ion_image_shape_y"
+    ion_image_shape_x: Final = "ion_image_shape_x"
+    ion_image_pixel_y: Final = "ion_image_pixel_y"
+    ion_image_pixel_x: Final = "ion_image_pixel_x"
+
+
+COORD_SYS_GLOBAL: Final = "global"
+MICROMETER: Final = "micrometer"
+OPTICAL_IMAGE_KEY: Final = "optical_image"
+POINTS_KEY: Final = "maldi_points"
+REGION_KEY: Final = "region"
+INSTANCE_KEY: Final = "instance_id"
+X: Final = "x"
+Y: Final = "y"
+C: Final = "c"
+XY: Final = (X, Y)
+YX: Final = (Y, X)
+YXC: Final = (Y, X, C)
diff --git a/metaspace_converter/tests/anndata_to_array_test.py b/metaspace_converter/tests/anndata_to_array_test.py
new file mode 100644
index 0000000..2956193
--- /dev/null
+++ b/metaspace_converter/tests/anndata_to_array_test.py
@@ -0,0 +1,75 @@
+from typing import TYPE_CHECKING
+
+import numpy as np
+import pandas as pd
+import pytest
+from anndata import AnnData
+
+if TYPE_CHECKING:
+    import _pytest.fixtures
+
+from metaspace_converter import anndata_to_image_array
+from metaspace_converter.constants import COL, METASPACE_KEY, X, Y
+from metaspace_converter.to_anndata import all_image_pixel_coordinates
+
+
+@pytest.fixture
+def adata_with_coordinates_and_image_size(
+    request: "_pytest.fixtures.SubRequest",
+) -> tuple[AnnData, np.ndarray]:
+    """
+    Create an AnnData object filled with parametrized shape
+
+    Args:
+        request: Request passed by Pytest, can be parametrized with a dictionary containing
+            ``num_ions``, ``height``, ``width``
+
+    Returns:
+        A tuple of an AnnData object and the expected array of stacked ion images
+    """
+    num_ions = request.param.get("num_ions", 0)
+    height = request.param.get("height", 0)
+    width = request.param.get("width", 0)
+    shape = (num_ions, height, width)
+    # Create a stack of ion images where every pixel has a different value
+    ion_images_stack = np.arange(np.prod(shape)).reshape(shape)
+    coordinates = all_image_pixel_coordinates(shape[1:])
+    # Create an AnnData with pixel coordinates
+    obs = pd.DataFrame(
+        {COL.ion_image_pixel_y: coordinates[:, 0], COL.ion_image_pixel_x: coordinates[:, 1]}
+    )
+    adata = AnnData(
+        X=ion_images_stack.reshape((height * width, num_ions)),
+        obs=obs,
+        uns={METASPACE_KEY: {"image_size": {X: width, Y: height}}},
+    )
+    expected = ion_images_stack
+    return adata, expected
+
+
+@pytest.mark.parametrize(
+    "adata_with_coordinates_and_image_size",
+    [
+        # Non-square ion images
+        dict(num_ions=4, height=2, width=3),
+        # Edge case: No annotations found
+        dict(num_ions=0, height=2, width=3),
+    ],
+    indirect=["adata_with_coordinates_and_image_size"],
+)
+def test_anndata_to_image_array(adata_with_coordinates_and_image_size: AnnData):
+    adata, expected = adata_with_coordinates_and_image_size
+
+    actual = anndata_to_image_array(adata)
+
+    assert actual.shape == (
+        adata.shape[1],
+        adata.obs[COL.ion_image_pixel_y].max() + 1,
+        adata.obs[COL.ion_image_pixel_x].max() + 1,
+    )
+
+    assert actual.shape == (
+        adata.shape[1],
+        adata.uns[METASPACE_KEY]["image_size"][Y],
+        adata.uns[METASPACE_KEY]["image_size"][X],
+    )
diff --git a/metaspace_converter/tests/to_anndata_test.py b/metaspace_converter/tests/to_anndata_test.py
index 2362900..cc6aac9 100644
--- a/metaspace_converter/tests/to_anndata_test.py
+++ b/metaspace_converter/tests/to_anndata_test.py
@@ -6,13 +6,8 @@
 import pytest
 from metaspace import SMInstance
 
-from metaspace_converter.to_anndata import (
-    COL,
-    METASPACE_KEY,
-    SPATIAL_KEY,
-    get_ion_image_shape,
-    metaspace_to_anndata,
-)
+from metaspace_converter.constants import COL, METASPACE_KEY, SPATIAL_KEY
+from metaspace_converter.to_anndata import get_ion_image_shape, metaspace_to_anndata
 
 METASPACE_DEFAULT_CONFIG_FILE = "~/.metaspace"
 METASPACE_EMAIL_ENV = "METASPACE_EMAIL"
diff --git a/metaspace_converter/tests/to_spatialdata_test.py b/metaspace_converter/tests/to_spatialdata_test.py
index 11dd47c..88e51ed 100644
--- a/metaspace_converter/tests/to_spatialdata_test.py
+++ b/metaspace_converter/tests/to_spatialdata_test.py
@@ -2,17 +2,12 @@
 import pytest
 from multiscale_spatial_image import MultiscaleSpatialImage
 
+from metaspace_converter.constants import COORD_SYS_GLOBAL, OPTICAL_IMAGE_KEY, POINTS_KEY, X, Y
+
 # The import metaspace_credentials is need by the fixture sm.
 from metaspace_converter.tests.to_anndata_test import metaspace_credentials, sm  # noqa
 from metaspace_converter.to_anndata import COL, METASPACE_KEY, get_ion_image_shape
-from metaspace_converter.to_spatialdata import (
-    COORD_SYS_GLOBAL,
-    OPTICAL_IMAGE_KEY,
-    POINTS_KEY,
-    X,
-    Y,
-    metaspace_to_spatialdata,
-)
+from metaspace_converter.to_spatialdata import metaspace_to_spatialdata
 
 
 @pytest.mark.parametrize(
diff --git a/metaspace_converter/to_anndata.py b/metaspace_converter/to_anndata.py
index d597450..e2fcb87 100644
--- a/metaspace_converter/to_anndata.py
+++ b/metaspace_converter/to_anndata.py
@@ -1,5 +1,4 @@
-from dataclasses import dataclass
-from typing import Container, Final, Optional
+from typing import Container, Optional
 
 import numpy as np
 import pandas as pd
@@ -12,6 +11,14 @@
     SMDataset,
 )
 
+from metaspace_converter.constants import (
+    COL,
+    METASPACE_KEY,
+    OBS_INDEX_NAME,
+    SPATIAL_KEY,
+    VAR_INDEX_NAME,
+    Shape2d,
+)
 from metaspace_converter.utils import (
     empty_columns_to_nan,
     has_optical_image,
@@ -20,20 +27,6 @@
     transform_coords,
 )
 
-SPATIAL_KEY: Final = "spatial"
-METASPACE_KEY: Final = "metaspace"
-OBS_INDEX_NAME: Final = "ion_image_pixel"
-VAR_INDEX_NAME: Final = "formula_adduct"
-Shape2d = tuple[int, int]
-
-
-@dataclass
-class COL:
-    ion_image_shape_y: Final = "ion_image_shape_y"
-    ion_image_shape_x: Final = "ion_image_shape_x"
-    ion_image_pixel_y: Final = "ion_image_pixel_y"
-    ion_image_pixel_x: Final = "ion_image_pixel_x"
-
 
 def metaspace_to_anndata(
     dataset: Optional[SMDataset] = None,
diff --git a/metaspace_converter/to_spatialdata.py b/metaspace_converter/to_spatialdata.py
index b508000..b478c2e 100644
--- a/metaspace_converter/to_spatialdata.py
+++ b/metaspace_converter/to_spatialdata.py
@@ -1,5 +1,5 @@
 import warnings
-from typing import Final, Optional
+from typing import Optional
 
 import numpy as np
 import pandas as pd
@@ -11,22 +11,23 @@
 from spatialdata.models import Image2DModel, PointsModel, TableModel
 from spatialdata.transformations import Affine, Scale, Sequence as SequenceTransform, Translation
 
-from metaspace_converter.to_anndata import COL, metaspace_to_anndata
+from metaspace_converter.constants import (
+    COL,
+    COORD_SYS_GLOBAL,
+    INSTANCE_KEY,
+    MICROMETER,
+    OPTICAL_IMAGE_KEY,
+    POINTS_KEY,
+    REGION_KEY,
+    XY,
+    YX,
+    YXC,
+    X,
+    Y,
+)
+from metaspace_converter.to_anndata import metaspace_to_anndata
 from metaspace_converter.utils import has_optical_image
 
-COORD_SYS_GLOBAL: Final = "global"
-MICROMETER: Final = "micrometer"
-OPTICAL_IMAGE_KEY: Final = "optical_image"
-POINTS_KEY: Final = "maldi_points"
-REGION_KEY: Final = "region"
-INSTANCE_KEY: Final = "instance_id"
-X: Final = "x"
-Y: Final = "y"
-C: Final = "c"
-XY: Final = (X, Y)
-YX: Final = (Y, X)
-YXC: Final = (Y, X, C)
-
 
 def metaspace_to_spatialdata(
     dataset: Optional[SMDataset] = None,