From 1483e8aab671826acf45d6a520223812cecb6f44 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Aur=C3=A9lien=20Coussat?=
 <aurelien.coussat@creatis.insa-lyon.fr>
Date: Thu, 2 Apr 2026 13:07:04 +0200
Subject: [PATCH] DOC: Add full Python example

Add a new Python example that illustrate the full reconstruction workflow, from the GATE simulation until the final FDK reconstruction.
---
 documentation/docs/getting_started.md     |  2 +
 examples/README.md                        |  9 +++
 examples/Reconstruction/README.md         | 12 ++++
 examples/Reconstruction/Reconstruction.py | 67 +++++++++++++++++++++++
 index.md                                  |  1 +
 5 files changed, 91 insertions(+)
 create mode 100644 examples/README.md
 create mode 100644 examples/Reconstruction/README.md
 create mode 100644 examples/Reconstruction/Reconstruction.py

diff --git a/documentation/docs/getting_started.md b/documentation/docs/getting_started.md
index 9c3aedd..1db1dca 100644
--- a/documentation/docs/getting_started.md
+++ b/documentation/docs/getting_started.md
@@ -136,6 +136,8 @@ The resulting image can be visualized by any MHD image viewer, such as [vv](http
 
 Note that to keep the running times in this guide fast enough, not enough statistics are generated to yield a nice output image. This is just an illustration of the PCT reconstruction workflow.
 
+PCT also provides Python functions that behave exactly like the applications described above. For illustration, the exact same simulation as the one described in this page is reimplemented in Python [`here`](../../examples/Reconstruction/README.md).
+
 ## Conclusion
 
 This guide gives an overview of a complete workflow that uses GATE to generate proton CT data, and PCT to process the data all the way to image reconstruction. To further familiarize yourself with PCT, you can explore the other applications offered by PCT, customize the current workflow with additional parameters using the `--help` flag, or implement your own simulation in GATE and try to produce a reconstruction.
diff --git a/examples/README.md b/examples/README.md
new file mode 100644
index 0000000..4d6fab2
--- /dev/null
+++ b/examples/README.md
@@ -0,0 +1,9 @@
+# Python examples
+
+This section provides a collection of Python code examples to demonstrate how to effectively use PCT in various applications.
+
+```{toctree}
+:maxdepth: 1
+
+./Reconstruction/README.md
+```
diff --git a/examples/Reconstruction/README.md b/examples/Reconstruction/README.md
new file mode 100644
index 0000000..eb45ac9
--- /dev/null
+++ b/examples/Reconstruction/README.md
@@ -0,0 +1,12 @@
+# Reconstruction
+
+A simple but complete pCT reconstruction workflow, starting from GATE Monte Carlo data, illustrating how to:
+- pair protons between the two detectors;
+- create distance-driven projections from the paired protons;
+- reconstruct an image using distance-driven FDK.
+
+## Code
+
+```{literalinclude} ./Reconstruction.py
+:language: python
+```
diff --git a/examples/Reconstruction/Reconstruction.py b/examples/Reconstruction/Reconstruction.py
new file mode 100644
index 0000000..abfb139
--- /dev/null
+++ b/examples/Reconstruction/Reconstruction.py
@@ -0,0 +1,67 @@
+#!/usr/bin/env python
+
+import os
+import sys
+import warnings
+from itk import PCT as pct
+from itk import RTK as rtk
+from opengate.contrib.protonct.protonct import protonct
+
+if len(sys.argv) < 1:
+    print("Usage: python Reconstruction.py <outputfolder>")
+    sys.exit(1)
+
+output_folder = sys.argv[1]
+
+number_of_projections = 90
+
+# Generate some data
+gate_folder = os.path.join(output_folder, "gate")
+protonct(gate_folder, projections=number_of_projections, verbose=False)
+
+# TODO example on how to make data noisy
+
+# Convert GATE data to PCT list-mode
+pairs_folder = os.path.join(output_folder, "pairs")
+os.makedirs(pairs_folder, exist_ok=True)
+pct.pctpairprotons(
+    input_in=os.path.join(gate_folder, "PhaseSpaceIn.root"),
+    input_out=os.path.join(gate_folder, "PhaseSpaceOut.root"),
+    output=os.path.join(pairs_folder, "pairs.mhd"),
+    psin="PhaseSpaceIn",
+    psout="PhaseSpaceOut",
+    plane_in=-110.0,
+    plane_out=110.0,
+    verbose=True,
+)
+
+# TODO cut the pairs (pctpaircuts is not converted to Python yet)
+
+# Bin the pairs into projections
+projections_folder = os.path.join(output_folder, "projections")
+os.makedirs(projections_folder, exist_ok=True)
+for p in range(number_of_projections):
+    pct.pctbinning(
+        input=os.path.join(pairs_folder, f"pairs{p:04d}.mhd"),
+        output=os.path.join(projections_folder, f"projections{p}.mhd"),
+        source=-1000.0,
+        size=[200, 1, 200],
+        spacing=[2.0, 1.0, 1.0],
+        verbose=True,
+    )
+
+# Build geometry
+geometry = os.path.join(output_folder, "geometry.xml")
+rtk.rtksimulatedgeometry(
+    nproj=number_of_projections, output=geometry, sdd=1000.0 + 110.0, sid=1000.0
+)
+
+# Reconstruct
+pct.pctfdk(
+    geometry=geometry,
+    path=projections_folder,
+    regexp=r"projections.*\.mhd",
+    output=os.path.join(output_folder, "recon.mhd"),
+    size=[210, 1, 210],
+    verbose=True,
+)
diff --git a/index.md b/index.md
index 84a76cc..7fbb9a8 100644
--- a/index.md
+++ b/index.md
@@ -22,6 +22,7 @@ documentation/docs/installation.md
 documentation/docs/getting_started.md
 documentation/docs/pct_format.md
 documentation/docs/wepl_calibration.md
+examples/README.md
 ```
 
 ```{toctree}