diff --git a/applications/pctdoublelut/pctdoublelut.py b/applications/pctdoublelut/pctdoublelut.py new file mode 100644 index 00000000..1ed7189c --- /dev/null +++ b/applications/pctdoublelut/pctdoublelut.py @@ -0,0 +1,307 @@ +#!/usr/bin/env python +import argparse +import json +import sys +from multiprocessing import Pool, Manager +import os +import warnings + +import numpy as np +import itk +from itk import PCT as pct + + +def build_parser(): + parser = pct.PCTArgumentParser( + description="Convert TOF to WEPL using a fit on Monte Carlo data", + ) + + parser.add_argument( + "-o", "--output", help="Path of outputs", default="pctdoublelut" + ) + parser.add_argument( + "-n", + "--number-of-particles", + help="Number of generated particles", + default=1000, + type=int, + ) + parser.add_argument("-m", "--material", help="Material", default="G4_WATER") + parser.add_argument( + "--wepl-samples", + help="Number of WEPL samples", + default=(2 * 260) + 1, + type=int, + ) + parser.add_argument( + "--max-wepl", + help="Maximum WEPL", + default=260.0, + type=float, + ) + parser.add_argument( + "-w", + "--phantom-width", + help="Phantom width", + default=40.0, + type=float, + ) + parser.add_argument( + "-e", + "--initial-energy", + help="Initial energy of the protons (in MeV)", + default=200.0, + type=float, + ) + parser.add_argument("--seed", help="Seed for random number generator", type=int) + parser.add_argument( + "--verbose", + "-v", + help="Verbose execution", + action="store_true", + ) + + return parser + + +epsilon_mm = 1e-5 + + +def pv(verbose, *args, **kwargs): + if verbose: + print(*args, **kwargs) + + +def get_gate_output(output): + return os.path.join(output, "gate") + + +def tof_fit_mc( + wepl, + output, + material, + phantom_width_cm, + number_of_particles, + initial_energy, + seed, + verbose, +): + import opengate as gate + + u = gate.g4_units + nm, mm, cm, m, sec, MeV = u.nm, u.mm, u.cm, u.m, u.second, u.MeV + + # Simulation + sim = gate.Simulation() + + sim.random_engine = "MersenneTwister" + sim.random_seed = seed + sim.run_timing_intervals = [[0.0 * sec, 1.0 * sec]] + sim.check_volumes_overlap = False + sim.g4_verbose = True + sim.progress_bar = verbose + + # Misc + yellow = [1, 1, 0, 1] + blue = [0, 0, 1, 1] + + # Geometry + sim.world.material = "G4_AIR" + sim.world.size = [4 * m, 4 * m, 4 * m] + + # Phantom + if wepl > 0.0: + phantom = sim.add_volume("Box", name="Phantom") + phantom.size = [ + phantom_width_cm * cm, + phantom_width_cm * cm, + wepl * mm, + ] + phantom.material = material + phantom.color = blue + phantom.set_max_step_size(1.0 * mm) + + # Beam + source = sim.add_source("GenericSource", "mybeam") + source.particle = "proton" + source.energy.mono = initial_energy * MeV + source.energy.type = "mono" + source.position.type = "box" + source.position.size = [1 * nm, 1 * nm, 1 * nm] + source.position.translation = [ + 0 * mm, + 0 * mm, + (-wepl / 2 - 0.1) * mm, + ] + source.direction.type = "momentum" + source.direction.momentum = [0, 0, 1] + source.n = number_of_particles + + # Physics list + sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option4" + sim.physics_manager.set_user_limits_particles(["proton"]) + + # Phase spaces + def add_detector(name, translation, attach_to_phantom=False): + plane = sim.add_volume("Box", "PlanePhaseSpace" + name) + plane.size = [phantom_width_cm * cm, phantom_width_cm * cm, 1 * nm] + plane.translation = translation + plane.material = "G4_AIR" + plane.color = yellow + if attach_to_phantom: + plane.mother = phantom.name + + phase_space = sim.add_actor("PhaseSpaceActor", "PhaseSpace" + name) + phase_space.attached_to = plane.name + phase_space.output_filename = os.path.join( + get_gate_output(output), f"wepl{wepl:.3f}_ps{name}.root" + ) + phase_space.attributes = [ + "EventID", + "TrackID", + "Position", + "LocalTime", + "PostVelocity", + ] + if int(gate.utility.version("opengate").split(".")[1]) > 0: + F = gate.actors.filters.GateFilterBuilder() + phase_space.filter = F.ParticleName == "proton" + else: + particle_filter = sim.add_filter("ParticleFilter", "Filter" + name) + particle_filter.particle = "proton" + + add_detector("In", [0.0 * mm, 0.0 * mm, (-wepl / 2) * mm]) + add_detector("Out", [0.0 * mm, 0.0 * mm, (wepl / 2) * mm]) + + sim.run() + + +def process_wepl(wepl, output, verbose): + import uproot + + wepls = [] + tofs_wepl = [] + velocities_wepl = [] + + data = uproot.concatenate( + os.path.join(get_gate_output(output), f"wepl{wepl:.3f}_*.root"), library="np" + ) + pv( + verbose, + "Loaded", + len(data["EventID"]), + "events for WEPL", + wepl, + ) + + # Sort data if needed + ws = data["Position_Z"] + if not np.all(ws[1:] > ws[:-1]): + pv(verbose, "Sorting input data…") + index_sorted = np.argsort(ws) + for key in data.keys(): + data[key] = data[key][index_sorted] + + for n in np.unique(data["EventID"]): + event_mask = data["EventID"] == n + + number_of_hits = np.sum(event_mask) + if number_of_hits == 0: + continue + + times = data["LocalTime"][event_mask] + tof = times[-1] - times[0] + + velocity = data["PostVelocity"][event_mask][-1] + + wepls.append(wepl) + tofs_wepl.append(tof) + velocities_wepl.append(velocity) + + return np.mean(wepls), np.mean(tofs_wepl), np.mean(velocities_wepl) + + +def tof_fit( + wepls, + output, + verbose, +): + + results = [] + points = [] + + with Pool() as pool: + for wepl in wepls: + result = pool.apply_async( + process_wepl, + ( + wepl, + output, + verbose, + ), + ) + results.append(result) + pool.close() + pool.join() + for result in results: + point = result.get() + points.append(point) + + np.savetxt(os.path.join(output, "data.txt"), points) + + data = np.array(points) + wepls = data[:, 0] + tofs = data[:, 1] + velocities = data[:, 2] + + degs = range(3, 10) + + for deg in degs: + tof_coeffs = np.polyfit(wepls, tofs, deg) + vel_coeffs = np.polyfit(wepls, velocities, deg) + np.savetxt(os.path.join(output, f"tof_coeffs_{deg}.txt"), tof_coeffs) + np.savetxt(os.path.join(output, f"vel_coeffs_{deg}.txt"), vel_coeffs) + + +def process(args_info: argparse.Namespace): + wepls = np.linspace(1.0, args_info.max_wepl, args_info.wepl_samples) + + results = [] + + # Catch DeprecationWarnings coming from opengate and making PCT tests fail + # To be removed once the warnings are gone from opengate + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", category=DeprecationWarning) + + with Pool(maxtasksperchild=1) as pool: + for wepl in wepls: + result = pool.apply_async( + tof_fit_mc, + ( + wepl, + args_info.output, + args_info.material, + args_info.phantom_width, + args_info.number_of_particles, + args_info.initial_energy, + args_info.seed, + args_info.verbose, + ), + ) + results.append(result) + pool.close() + pool.join() + for result in results: + result.get() + + tof_fit(wepls, args_info.output, args_info.verbose) + + +def main(argv=None): + parser = build_parser() + args_info = parser.parse_args(argv) + process(args_info) + + +if __name__ == "__main__": + main() diff --git a/applications/pctpairprotons/pctpairprotons.py b/applications/pctpairprotons/pctpairprotons.py index 07d4c9a4..d9259168 100644 --- a/applications/pctpairprotons/pctpairprotons.py +++ b/applications/pctpairprotons/pctpairprotons.py @@ -48,14 +48,6 @@ def build_parser(): default=False, action="store_true", ) - parser.add_argument( - "--fit", help="Fit file used to convert from energy loss or TOF to WEPL" - ) - parser.add_argument( - "--fit-kind", - help="Whether to convert to WEPL using energy loss or TOF", - choices=["tof", "energy"], - ) parser.add_argument( "--store-time", help="Store time instead of energy in the output list-mode", @@ -72,9 +64,132 @@ def build_parser(): "--psout", help="Name of tree in output phase space", default="PhaseSpace" ) + single_lut = parser.add_argument_group("Single LUT for conversion to WEPL") + single_lut.add_argument( + "--fit", help="Fit file used to convert from energy loss or TOF to WEPL" + ) + single_lut.add_argument( + "--fit-kind", + help="Whether to convert to WEPL using energy loss or TOF", + choices=["tof", "energy"], + ) + + double_lut = parser.add_argument_group("Double LUT for conversion to WEPL") + double_lut.add_argument( + "--lut-tof", help="LUT corresponding to TOF in double LUT method" + ) + double_lut.add_argument( + "--lut-vel", help="LUT corresponding to velocity in double LUT method" + ) + double_lut.add_argument( + "--quadric", + help="Quadric representing the hull of the object", + type=float, + nargs=10, + required=False, + ) + double_lut.add_argument( + "--angle", help="Angle of the current projection", type=float, required=False + ) + return parser +m0 = 938.27208943 # MeV/c2 +c = 299.792458 # mm/ns + + +def make_hull(quadric, angle): + from itk import RTK as rtk + + hull = rtk.QuadricShape.New() + hull.SetA(quadric[0]) + hull.SetB(quadric[1]) + hull.SetC(quadric[2]) + hull.SetD(quadric[3]) + hull.SetE(quadric[4]) + hull.SetF(quadric[5]) + hull.SetG(quadric[6]) + hull.SetH(quadric[7]) + hull.SetI(quadric[8]) + hull.SetJ(quadric[9]) + + theta = np.deg2rad(angle) + hull.Rotate( + itk.matrix_from_array( + [ + [np.cos(theta), 0.0, np.sin(theta)], + [0.0, 1.0, 0.0], + [-np.sin(theta), 0.0, np.cos(theta)], + ] + ) + ) + + return hull + + +def tof_in_vacuum(d, e): + v = c * np.sqrt(1 - (m0**2) / ((m0 + e) ** 2)) + return d / v + + +def get_tof_and_distance(hull, p_u, d_u, p_d, d_d, energy, tof): + + d_d_r = [-d for d in d_d] + + intersect_u, d_uo, _ = hull.IsIntersectedByRay(p_u, d_u) # assuming d_u is unitary + intersect_d, d2, _ = hull.IsIntersectedByRay(p_d, d_d_r) + if not (intersect_u and intersect_d): + raise ValueError + + int1 = [p_u[i] + d_u[i] * d_uo for i in range(3)] + int2 = [p_d[i] + d_d_r[i] * d2 for i in range(3)] + d1 = np.sqrt(np.sum([(int1[i] - int2[i]) ** 2 for i in range(3)])) + + tof_od = tof - tof_in_vacuum(d_uo, energy) + + return tof_od, d1, d2 + + +def convert_tof_to_wepl(fit_tof, fit_vel, pairs, quadric, angle): + tof_coeffs = np.loadtxt(fit_tof) + vel_coeffs = np.loadtxt(fit_vel) + + def tof_to_wepl(tof, d1, d2): + t1 = d1 * np.polymul(tof_coeffs, vel_coeffs) + t2 = [d2, 0.0] + t3 = tof * np.polymul(vel_coeffs, [1.0, 0.0]) + p = np.polysub(np.polyadd(t1, t2), t3) + roots = [ + np.real(r) + for r in np.roots(p) + if np.imag(r) == 0.0 + and np.real(r) >= 0.0 + and (d1 / r) * np.polyval(tof_coeffs, r) < tof + ] + return roots[-1] + + names = pairs.dtype.names + hull = make_hull(quadric, angle) + + for pi, pair in enumerate(pairs): + p = dict(zip(names, pair.tolist())) + p_u = [p["u_in"], p["v_in"], p["w_in"]] + p_d = [p["u_out"], p["v_out"], p["w_out"]] + d_u = [p["du_in"], p["dv_in"], p["dw_in"]] + d_d = [p["du_out"], p["dv_out"], p["dw_out"]] + energy = p["KineticEnergy_in"] + tof = p["LocalTime_out"] - p["LocalTime_in"] + try: + tof, d1, d2 = get_tof_and_distance(hull, p_u, d_u, p_d, d_d, energy, tof) + wepl = tof_to_wepl(tof, d1, d2) + except ValueError: + wepl = 0.0 + + pairs[pi]["KineticEnergy_out"] = wepl + pairs[pi]["KineticEnergy_in"] = 0.0 + + def process(args_info: argparse.Namespace): import uproot @@ -194,8 +309,8 @@ def load_tree_as_df(root_file, tree_name): np.recarray.sort(pairs, order=["RunID", "EventID", "TrackID_in", "TrackID_out"]) verbose("Merged input and output phase spaces.") - if args_info.fit is not None: - verbose("Converting energy loss or TOF to WEPL…") + if args_info.fit is not None: # Single LUT + verbose("Converting energy loss or TOF to WEPL using single LUT technique…") with open(args_info.fit, encoding="utf-8") as f: p = json.load(f) if args_info.fit_kind == "tof": @@ -207,6 +322,15 @@ def load_tree_as_df(root_file, tree_name): wepls = np.polyval(p, xs) pairs["KineticEnergy_in"] = 0.0 pairs["KineticEnergy_out"] = wepls + elif args_info.lut_tof is not None and args_info.lut_vel is not None: # Double LUT + verbose("Converting energy loss or TOF to WEPL using double LUT technique…") + convert_tof_to_wepl( + args_info.lut_tof, + args_info.lut_vel, + pairs, + args_info.quadric, + args_info.angle, + ) number_of_runs = pairs["RunID"].max() + 1 verbose("Identified number of runs: " + str(number_of_runs)) diff --git a/documentation/docs/wepl_calibration.md b/documentation/docs/wepl_calibration.md index 512ad271..b29b4f81 100644 --- a/documentation/docs/wepl_calibration.md +++ b/documentation/docs/wepl_calibration.md @@ -1,6 +1,12 @@ # WEPL calibration -PCT provides the application `pctweplfit` that generates calibration curves from energy loss and TOF to WEPL. The resulting polynomial fit can then be used when pairing protons with `pctpairprotons` to convert the energy loss or the TOF directly to WEPL. Please note that `pctweplfit` requires GATE, which can be installed following the instructions in the [installation guide](installation.md). +PCT supports several calibration curve methods to convert from a given measurement (energy loss or TOF) to the corresponding WEPL. +First, the calibration curve(s) must be generated using one of the applications provided by PCT. Generating the calibration curves typically requires GATE, which can be installed following the instructions in the [installation guide](installation.md). +Then, the resulting curves can be passed to `pctpairprotons` in order to properly convert the measurements to WEPL when forming proton pairs. The resulting pairs are directly associated to the corresponding WEPL following the convention explained in the [PCT data format](pct_format.md), that is that $e_\text{in}=0$ and $e_\text{out}=\text{WEPL}$. + +## Single 1D calibration curve + +PCT provides the application `pctweplfit` that generates 1D calibration curves from energy loss or TOF to WEPL, following the technique described in [Coussat et al, Fully 3D, 2025](https://hal.science/hal-05375602). The resulting polynomial fit can then be used when pairing protons with `pctpairprotons` to convert the energy loss or the TOF directly to WEPL. Below is an example of how to run `pctweplfit`: ```bash @@ -10,7 +16,7 @@ pctweplfit \ -v ``` -In the resuling `output` folder, along with all intermediate results (ROOT files from the GATE simulations) will be two files `eloss_to_wepl_fit_deg3.json` and `tof_to_wepl_fit_deg3.json` that contain the coefficients for the polynomials. These files can then be passed to `pctpairprotons` with the corresponding `--fit-kind` parameters. +In the resulting `output` folder, along with all intermediate results (ROOT files from the GATE simulations) will be two files `eloss_to_wepl_fit_deg3.json` and `tof_to_wepl_fit_deg3.json` that contain the coefficients of the polynomials. These files can then be passed to `pctpairprotons` with the corresponding `--fit-kind` parameters. Here is an example of a `pctpairprotons` invocations with energy-loss fit: ```bash @@ -24,4 +30,34 @@ pctpairprotons \ --fit-kind energy ``` -The resulting pairs are directly associated to the corresponding WEPL following the convention explained in the [PCT data format](pct_format.md), that is that $e_\text{in}=0$ and $e_\text{out}=\text{WEPL}$. +## Double 1D calibration curve + +PCT provides the application `pctdoublelut` that generates two 1D calibration curves from TOF to WEPL. The resulting polynomial fits can then be used when pairing protons with `pctpairprotons` to convert the TOF directly to WEPL. This conversion is not available for energy loss measurements. + +Below is an example of how to run `pctdoublelut` to generate calibration curves for a 200 MeV beam: +```bash +pctdoublelut \ + -o output \ + -e 200 \ + -v +``` + +In the resulting `output` folder, along with all intermediate results (ROOT files from the GATE simulations) will be files named `tof_coeffs_X.json` and `vel_coeffs_X.json` that contain the coefficients of the polynomials, where `X` is the polynomial degree. These files can then be passed to `pctpairprotons` with the appropriate parameters. + +Here is an example of a `pctpairprotons` invocations that uses the double LUT fits with polynomials of degree 9: +```bash +pctpairprotons \ + -i PhaseSpaceIn_0.root \ + -j PhaseSpaceOut_0.root \ + -o pairs.mhd \ + --plane-in -110 \ + --plane-out 110 \ + --lut-tof output/tof_coeffs_9.json \ + --lut-vel output/vel_coeffs_9.json \ + --quadric 1 0 1 0 0 0 0 0 0 -10000 \ + --angle 0 +``` + +Additional notes: +- This method requires the knowledge of the (potentially approximate) support of the object, which is passed to the script as a quadric shape. For additional details about how to define quadric shapes, you can refer to [the documentation of `rtkQuadricShape`](https://www.openrtk.org/Doxygen/classrtk_1_1QuadricShape.html). +- This method can only process list-mode files that contain data for a single acquisition angle, which is specified using the `--angle` parameter. It cannot process data that contain more than one angle. diff --git a/pyproject.toml b/pyproject.toml index cbc5503f..23be7277 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -47,6 +47,7 @@ pctpairprotons = "itk.pctpairprotons:main" pctweplfit = "itk.pctweplfit:main" pctlomalinda = "itk.pctlomalinda:main" pctaddnoise = "itk.pctaddnoise:main" +pctdoublelut = "itk.pctdoublelut:main" [project.urls] Download = "https://github.com/RTKConsortium/PCT" diff --git a/test/pct_application_test.py b/test/pct_application_test.py index a7f0918c..9043d9cf 100644 --- a/test/pct_application_test.py +++ b/test/pct_application_test.py @@ -77,6 +77,105 @@ def test_weplfit_application(tmp_path): assert np.allclose(eloss_to_wepl_fit, reference) +def test_doublelut_application(tmp_path): + output = tmp_path / "doublelut" + pct.pctdoublelut(f"-o {output} -n 1000 --wepl-samples 10 -e 200 --seed 1234") + + tof_coeffs = np.loadtxt(output / "tof_coeffs_9.txt") + tof_reference = np.array( + [ + -1.290921042580240798e-19, + 1.351207981731195971e-16, + -5.963202845355236675e-14, + 1.443362133502381465e-11, + -2.084980429067098285e-09, + 1.829330424219974866e-07, + -9.443542359443581318e-06, + 2.631908622195730887e-04, + 2.921848578064085208e-03, + 2.715972855380448421e-03, + ] + ) + assert np.allclose(tof_coeffs, tof_reference) + + vel_coeffs = np.loadtxt(output / "vel_coeffs_9.txt") + vel_reference = np.array( + [ + 5.023900900539231649e-18, + -5.304589281207117926e-15, + 2.358303897263661076e-12, + -5.748047473922940780e-10, + 8.360830241840195663e-08, + -7.391334678087265226e-06, + 3.842573655671269075e-04, + -1.096847900269657186e-02, + -1.773574286549965684e-02, + 1.696161111400466552e02, + ] + ) + assert np.allclose(vel_reference, vel_coeffs) + + +baseline_pairs_doublelut_mhd = download_file_fixture( + "6a8f008e92f283f838800623", "baseline_pairs_doublelut.mhd" +) +baseline_pairs_doublelut_raw = download_file_fixture( + "6a8f009092f283f838800626", "baseline_pairs_doublelut.raw" +) + + +def test_pairprotons_doublelut_application( + tmp_path, + phasespacein_root, + phasespaceout_root, + baseline_pairs_doublelut_mhd, + baseline_pairs_doublelut_raw, +): + output = tmp_path / "pairs_doublelut.mhd" + + tof_coeffs = tmp_path / "tof_coeffs.txt" + np.savetxt( + tof_coeffs, + [ + 2.796830534907338879e-22, + -1.788065254786836053e-19, + 4.591994205164617155e-17, + -5.662289174905150623e-15, + 3.899965359778268460e-13, + -4.224885502649042288e-12, + 4.745768765394230579e-09, + 2.432353352721748280e-06, + 5.889704508365986406e-03, + -9.553979783485672557e-07, + ], + ) + vel_coeffs = tmp_path / "vel_coeffs.txt" + np.savetxt( + vel_coeffs, + [ + -1.811548969775848679e-19, + 1.487181098408996225e-16, + -5.245527890416944471e-14, + 1.021392695091731880e-11, + -1.214073609154529679e-09, + 8.787192875077365930e-08, + -4.584117446224705820e-06, + -1.447480766127036390e-04, + -1.423014493400832636e-01, + 1.697332540903398126e02, + ], + ) + + pct.pctpairprotons( + f"-i {phasespacein_root} -j {phasespaceout_root} -o {output} --plane-in -110 --plane-out 110 --psin PhaseSpaceIn --psout PhaseSpaceOut --lut-tof {tof_coeffs} --lut-vel {vel_coeffs} --quadric 1 0 1 0 0 0 0 0 0 -10000 --angle 0" + ) + + output0000 = tmp_path / "pairs_doublelut0000.mhd" + pairs_test = itk.array_from_image(itk.imread(output0000)) + pairs_baseline = itk.array_from_image(itk.imread(baseline_pairs_doublelut_mhd)) + assert np.array_equal(pairs_test, pairs_baseline) + + lomalinda_data = download_file_fixture( "69e21803ed08a1c077afd077", "projection_045.root" ) diff --git a/wrapping/__init_pct__.py b/wrapping/__init_pct__.py index 191c94ed..3611fb4d 100644 --- a/wrapping/__init_pct__.py +++ b/wrapping/__init_pct__.py @@ -27,6 +27,7 @@ "pctweplfit", "pctlomalinda", "pctaddnoise", + "pctdoublelut", ] # Dynamically access make_application_func from pctExtras