Add full Python example

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.

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
+```

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
+```

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,
+)

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}