diff --git a/src/Cxx/ImageData/ImageIteratorDemo.cxx b/src/Cxx/ImageData/ImageIteratorDemo.cxx
index f43f75fd8bce273ac5dc0fe8202b2a2add4530e3..df5fda89aea14f67a7b3365efab63b4b7dc5a7da 100644
--- a/src/Cxx/ImageData/ImageIteratorDemo.cxx
+++ b/src/Cxx/ImageData/ImageIteratorDemo.cxx
@@ -67,9 +67,9 @@ int main(int, char*[])
     it.NextSpan();
     ++counter;
   }
-  std::cout << "# of spans: " << counter << std::endl;
+  std::cout << "Number of spans:      " << counter << std::endl;
 
-  std::cout << "Increments: " << imageData->GetIncrements()[0] << ", "
+  std::cout << "Increments:           " << imageData->GetIncrements()[0] << ", "
             << imageData->GetIncrements()[1] << ", "
             << imageData->GetIncrements()[2] << std::endl;
   vtkIdType incX, incY, incZ;
diff --git a/src/PythonicAPI.md b/src/PythonicAPI.md
index 23d53de6a355f6d65f10d5d38d655c5de4ba9c19..0bacdcfe561b2dcb2853db215ff1c6fc0bf72267 100644
--- a/src/PythonicAPI.md
+++ b/src/PythonicAPI.md
@@ -290,6 +290,8 @@ This section includes examples of manipulating meshes.
 | Example Name | Description | Image |
 | -------------- | ------------- | ------- |
 [ClipVolume](/PythonicAPI/ImageData/ClipVolume) | Clip a volume and produce a vtkUnhstructuredGrid.
+[ImageIterator](/PythonicAPI/ImageData/ImageIterator) |
+[ImageIteratorDemo](/PythonicAPI/ImageData/ImageIteratorDemo) | Demonstrate using indexing to access pixels in a region.
 [ImageWeightedSum](/PythonicAPI/ImageData/ImageWeightedSum) | Add two or more images.
 
 #### ?vtkExplicitStructuredGrid?
diff --git a/src/PythonicAPI/ImageData/ImageIterator.md b/src/PythonicAPI/ImageData/ImageIterator.md
new file mode 100644
index 0000000000000000000000000000000000000000..d383815ac28a96a31a4af16c8f9c6b300a5dfe0e
--- /dev/null
+++ b/src/PythonicAPI/ImageData/ImageIterator.md
@@ -0,0 +1,7 @@
+### Description
+
+Extracts an extent from an image.
+
+The differences from the C++ example are:
+- We cannot use pointers in Python so we use `SetScalarComponentFromDoubl(...)` and `GetScalarComponentFromDouble(...)` instead of `GetScalarPointer(...)`.
+- We cannot use vtkImageIterator as it is not wrapped in Python. So we use indexing instead.
diff --git a/src/PythonicAPI/ImageData/ImageIterator.py b/src/PythonicAPI/ImageData/ImageIterator.py
new file mode 100755
index 0000000000000000000000000000000000000000..d0cf9e47887bbf92e31d3d536552fbbe7c8dfc6f
--- /dev/null
+++ b/src/PythonicAPI/ImageData/ImageIterator.py
@@ -0,0 +1,62 @@
+#!/usr/bin/env python3
+
+# noinspection PyUnresolvedReferences
+import vtkmodules.vtkInteractionStyle
+# noinspection PyUnresolvedReferences
+import vtkmodules.vtkRenderingOpenGL2
+from vtkmodules.vtkCommonCore import VTK_DOUBLE
+from vtkmodules.vtkCommonDataModel import vtkImageData
+
+
+def main():
+    # We cannot use vtkImageIterator as it is not wrapped in the Python API
+    #  so we use indexing to produce the same result as in the C++ example.
+
+    # Create an image data and specify the size and type of the image data.
+    image_data = vtkImageData(dimensions=(10, 20, 30))
+    image_data.AllocateScalars(VTK_DOUBLE, 3)
+
+    # Fill every entry of the image data with x,y,z.
+    dims = image_data.GetDimensions()
+    for z in range(dims[2]):
+        for y in range(dims[1]):
+            for x in range(dims[0]):
+                image_data.SetScalarComponentFromDouble(x, y, z, 0, x)
+                image_data.SetScalarComponentFromDouble(x, y, z, 1, y)
+                image_data.SetScalarComponentFromDouble(x, y, z, 2, z)
+
+    # Define the extent to be extracted.
+    extent = [2, 5, 2, 5, 15, 15]
+
+    # Retrieve the entries from the image data and print them to the screen.
+    # We cannot use vtkImageIterator as it is not wrapped in the Python API so we use indexing.
+    # Note that we add 1 to the upper index in extent so that range() works correctly.
+    for z in range(extent[4], extent[5] + 1):
+        for y in range(extent[2], extent[3] + 1):
+            res = list()
+            for x in range(extent[0], extent[1] + 1):
+                zyx = list()
+                for i in reversed(range(0, 3)):
+                    zyx.append(image_data.GetScalarComponentAsDouble(x, y, z, i))
+                res.append(f'({fmt_floats(zyx)})')
+            print(' '.join(res))
+
+
+def fmt_floats(v, w=0, d=6, pt='g'):
+    """
+    Pretty print a list or tuple of floats.
+
+    :param v: The list or tuple of floats.
+    :param w: Total width of the field.
+    :param d: The number of decimal places.
+    :param pt: The presentation type, 'f', 'g' or 'e'.
+    :return: A string.
+    """
+    pt = pt.lower()
+    if pt not in ['f', 'g', 'e']:
+        pt = 'f'
+    return ', '.join([f'{element:{w}.{d}{pt}}' for element in v])
+
+
+if __name__ == '__main__':
+    main()
diff --git a/src/PythonicAPI/ImageData/ImageIteratorDemo.md b/src/PythonicAPI/ImageData/ImageIteratorDemo.md
new file mode 100644
index 0000000000000000000000000000000000000000..487aaedb6b3b11bc68b78b35d55cd182f1fb0457
--- /dev/null
+++ b/src/PythonicAPI/ImageData/ImageIteratorDemo.md
@@ -0,0 +1,8 @@
+### Description
+
+Unlike in C++ where vtkImageIterator is an efficient way to access the regions of a vtkImageData. We have to use indexing instead.
+
+The differences from the C++ example are:
+- We cannot use pointers in Python so we use `SetScalarComponentFromDoubl(...)` and `GetScalarComponentFromDouble(...)` instead of `GetScalarPointer(...)`.
+- We cannot use vtkImageIterator as it is not wrapped in Python.
+So we use indexing instead.
diff --git a/src/PythonicAPI/ImageData/ImageIteratorDemo.py b/src/PythonicAPI/ImageData/ImageIteratorDemo.py
new file mode 100755
index 0000000000000000000000000000000000000000..eab92065d9a5effa7d4d97a9e3dd8a34a0b28083
--- /dev/null
+++ b/src/PythonicAPI/ImageData/ImageIteratorDemo.py
@@ -0,0 +1,110 @@
+#!/usr/bin/env python3
+
+# noinspection PyUnresolvedReferences
+import vtkmodules.vtkInteractionStyle
+# noinspection PyUnresolvedReferences
+import vtkmodules.vtkRenderingOpenGL2
+from vtkmodules.vtkCommonColor import vtkNamedColors
+from vtkmodules.vtkCommonCore import VTK_UNSIGNED_CHAR
+from vtkmodules.vtkCommonDataModel import vtkImageData
+from vtkmodules.vtkCommonCore import reference
+from vtkmodules.vtkInteractionImage import vtkImageViewer2
+from vtkmodules.vtkRenderingCore import vtkRenderWindowInteractor
+from vtkmodules.vtkInteractionStyle import vtkInteractorStyleImage
+
+
+def main():
+    # We cannot use vtkImageIterator as it is not wrapped in the Python API
+    #  so we use indexing to produce the same result as in the C++ example.
+
+    # Create an image data and specify the size and type of the image data.
+    image_data = vtkImageData(dimensions=(100, 200, 30))
+    image_data.AllocateScalars(VTK_UNSIGNED_CHAR, 3)
+
+    colors = vtkNamedColors()
+
+    rgba = tuple(colors.GetColor4ub('Banana'))
+
+    # Fill every entry of the image data with r, g, b.
+    dims = image_data.GetDimensions()
+    for z in range(dims[2]):
+        for y in range(dims[1]):
+            for x in range(dims[0]):
+                for i in range(0, 3):
+                    image_data.SetScalarComponentFromDouble(x, y, z, i, rgba[i])
+
+    # whole_extent = image_data.GetExtent()
+    # Define the extent to be extracted.
+    extent = [20, 50, 30, 60, 10, 20]
+
+    # We cannot use vtkImageIterator as it is not wrapped in the Python API so we use indexing.
+    # Note that we add 1 to the upper index in extent so that range() works correctly.
+    rgba = tuple(colors.GetColor4ub('Tomato'))
+    counter=0
+    for z in range(extent[4], extent[5] + 1):
+        for y in range(extent[2], extent[3] + 1):
+            counter += 1
+            for x in range(extent[0], extent[1] + 1):
+                for i in range(0, 3):
+                    image_data.SetScalarComponentFromDouble(x, y, z, i, rgba[i])
+
+    print(f'Number of spans:      {counter}')
+    print(f'Increments:           {fmt_ints(image_data.increments)}')
+    inc_x = reference(0)
+    inc_y = reference(0)
+    inc_z = reference(0)
+    image_data.GetContinuousIncrements(extent, inc_x, inc_y, inc_z)
+    print(f'ContinuousIncrements: {inc_x}, {inc_y}, {inc_z}')
+
+    # Visualize
+    image_viewer = vtkImageViewer2(input_data=image_data)
+
+    style = vtkInteractorStyleImage()
+    style.SetInteractionModeToImageSlicing()
+
+    render_window_interactor = vtkRenderWindowInteractor()
+    render_window_interactor.interactor_style = style
+    image_viewer.SetupInteractor(render_window_interactor)
+    slice = (extent[5] - extent[4]) // 2 + extent[4]
+    image_viewer.SetSlice(slice)
+
+    image_viewer.renderer.background = colors.GetColor3d('Slate_grey')
+    image_viewer.image_actor.interpolate = False
+
+    image_viewer.Render()
+    image_viewer.renderer.ResetCamera()
+    image_viewer.render_window.window_name = 'ImageIteratorDemo'
+
+    image_viewer.Render()
+
+    render_window_interactor.Start()
+
+
+def fmt_ints(v, w=0):
+    """
+    Pretty print a list or tuple of ints.
+
+    :param v: The list or tuple of ints.
+    :param w: Total width of the field.
+    :return: A string.
+    """
+    return ', '.join([f'{element:{w}}' for element in v])
+
+def fmt_floats(v, w=0, d=6, pt='g'):
+    """
+    Pretty print a list or tuple of floats.
+
+    :param v: The list or tuple of floats.
+    :param w: Total width of the field.
+    :param d: The number of decimal places.
+    :param pt: The presentation type, 'f', 'g' or 'e'.
+    :return: A string.
+    """
+    pt = pt.lower()
+    if pt not in ['f', 'g', 'e']:
+        pt = 'f'
+    return ', '.join([f'{element:{w}.{d}{pt}}' for element in v])
+
+
+if __name__ == '__main__':
+    main()
diff --git a/src/Testing/Baseline/PythonicAPI/ImageData/TestImageIteratorDemo.png b/src/Testing/Baseline/PythonicAPI/ImageData/TestImageIteratorDemo.png
new file mode 100644
index 0000000000000000000000000000000000000000..a2f4bddc0dc28b108218967c2c1d4d7635da1413
--- /dev/null
+++ b/src/Testing/Baseline/PythonicAPI/ImageData/TestImageIteratorDemo.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:2bca8453952b026aaadb5a48c7695a1875ed1d872109946319a292b881893708
+size 604