diff --git a/applications/pctfilterprotons/pctfilterprotons.py b/applications/pctfilterprotons/pctfilterprotons.py new file mode 100644 index 0000000..548bb45 --- /dev/null +++ b/applications/pctfilterprotons/pctfilterprotons.py @@ -0,0 +1,102 @@ +#!/usr/bin/env python +import argparse +import itk +import numpy as np + + +def build_parser(): + + parser = argparse.ArgumentParser( + description="Filter protons according to various criteria" + ) + parser.add_argument("-i", "--input", help="Input list-mode file", required=True) + parser.add_argument("-o", "--output", help="Output list-mode file", required=True) + parser.add_argument("--roi-radius", help="Radius of the ROI in mm", type=float) + parser.add_argument( + "--fluence", help="Fluence level as fraction of full fluence", type=float + ) + parser.add_argument("--min-wepl", help="Minimum WEPL", type=float) + parser.add_argument("--max-wepl", help="Maximum WEPL", type=float) + parser.add_argument( + "--verbose", "-v", help="Verbose execution", default=False, action="store_true" + ) + + return parser + + +def process(args_info: argparse.Namespace): + + if args_info.verbose: + + def verbose(message): + print(message) + + else: + + def verbose(message): + pass + + pairs_itk = itk.imread(args_info.input) + pairs = itk.GetArrayFromImage(pairs_itk) + + u_in = np.s_[:, 0, 0] + w_in = np.s_[:, 0, 2] + u_out = np.s_[:, 1, 0] + w_out = np.s_[:, 1, 2] + e_in = np.s_[:, 4, 0] + wepl = np.s_[:, 4, 1] + + if np.any(pairs[e_in] != 0.0) and ( + args_info.min_wepl is not None or args_info.max_wepl is not None + ): + raise ValueError( + "Cannot filter WEPLs because input file does not contain WEPL (e_in != 0)! Aborting." + ) + + interception = np.full_like(pairs[wepl], True) + + if args_info.roi_radius is not None: + # Circular ROI intersection + verbose("Filtering ions outside of the ROI…") + radius = args_info.roi_radius + du = pairs[u_out] - pairs[u_in] + dt = pairs[w_out] - pairs[w_in] + dr_sq = (du * du) + (dt * dt) + d = (pairs[u_in] * pairs[w_out]) - (pairs[u_out] * pairs[w_in]) + delta = (radius * radius * dr_sq) - (d * d) + interception = np.logical_and(interception, delta >= 0.0) + + if args_info.fluence is not None: + # Fluence filtering + verbose("Applying fluence level…") + rng = np.random.default_rng() + fluence_filter = rng.random(interception.shape) + interception = np.logical_and(interception, fluence_filter < args_info.fluence) + + if args_info.min_wepl is not None: + verbose("Filtering low WEPLs…") + interception = np.logical_and(interception, pairs[wepl] >= args_info.min_wepl) + + if args_info.max_wepl is not None: + verbose("Filtering high WEPLs…") + interception = np.logical_and(interception, pairs[wepl] <= args_info.max_wepl) + + verbose("Applying interception filter…") + + ComponentType = itk.ctype("float") + PixelType = itk.Vector[ComponentType, 3] + ImageType = itk.Image[PixelType, 2] + + output_itk = itk.GetImageFromArray(pairs[interception], ttype=ImageType) + itk.imwrite(output_itk, args_info.output) + verbose(f"Wrote file {args_info.output}.") + + +def main(argv=None): + parser = build_parser() + args_info = parser.parse_args(argv) + process(args_info) + + +if __name__ == "__main__": + main() diff --git a/applications/pctlomalinda/pctlomalinda.py b/applications/pctlomalinda/pctlomalinda.py new file mode 100644 index 0000000..cb7a687 --- /dev/null +++ b/applications/pctlomalinda/pctlomalinda.py @@ -0,0 +1,128 @@ +#!/usr/bin/env python +import argparse +import uproot +import itk + +import numpy as np + + +def build_parser(): + + parser = argparse.ArgumentParser(description="Convert Loma Linda data to PCT data") + parser.add_argument( + "-i", "--input", help="Root phase space file of particles", required=True + ) + parser.add_argument("-o", "--output", help="Output file name", required=True) + parser.add_argument( + "--plane-in", + help="Plane position of incoming protons", + required=True, + type=float, + ) + parser.add_argument( + "--plane-out", + help="Plane position of outgoing protons", + required=True, + type=float, + ) + parser.add_argument( + "--min-run", help="Minimum run (inclusive)", default=0, type=int + ) + parser.add_argument( + "--max-run", help="Maximum run (exclusive)", default=1e6, type=int + ) + parser.add_argument( + "--verbose", "-v", help="Verbose execution", default=False, action="store_true" + ) + parser.add_argument( + "--ps", help="Name of tree in input phase space", default="PhaseSpace" + ) + + return parser + + +def process(args_info: argparse.Namespace): + + if args_info.verbose: + + def verbose(message): + print(message) + + else: + + def verbose(message): + pass + + tree = uproot.open(args_info.input)[args_info.ps] + pairs = tree.arrays(library="np") + verbose("Read input phase space:\n" + str(pairs)) + + # should match what it in the file, no checks is performed + pairs["u_hit1"] = np.full_like(pairs["u_hit1"], args_info.plane_in) + pairs["u_hit2"] = np.full_like(pairs["u_hit2"], args_info.plane_out) + + verbose("Filter RunIDs…") + # "projection_angle" acts as the run ID here + interception = np.logical_and( + pairs["projection_angle"] < args_info.max_run, + pairs["projection_angle"] >= args_info.min_run, + ) + for k, v in pairs.items(): + pairs[k] = v[interception] + + verbose("Calculating directions…") + dir_in_t = pairs["t_hit1"] - pairs["t_hit0"] + dir_in_v = pairs["v_hit1"] - pairs["v_hit0"] + dir_in_u = pairs["u_hit1"] - pairs["u_hit0"] + norm_in = np.sqrt(dir_in_t**2 + dir_in_v**2 + dir_in_u**2) + dir_out_t = pairs["t_hit3"] - pairs["t_hit2"] + dir_out_v = pairs["v_hit3"] - pairs["v_hit2"] + dir_out_u = pairs["u_hit3"] - pairs["u_hit2"] + norm_out = np.sqrt(dir_out_t**2 + dir_out_v**2 + dir_out_u**2) + + ComponentType = itk.ctype("float") + PixelType = itk.Vector[ComponentType, 3] + ImageType = itk.Image[PixelType, 2] + + runs = np.unique(pairs["projection_angle"]) + number_of_runs = len(runs) + verbose("Identified number of runs: " + str(number_of_runs)) + run_range = range(args_info.min_run, min(number_of_runs, args_info.max_run)) + for r in run_range: + ps_run = pairs["projection_angle"] == runs[r] + len_ps_run = ps_run.sum() + if len_ps_run == 0: + continue + + ps_np = np.empty(shape=(len_ps_run, 5, 3), dtype=np.float32) + ps_np[:, 0, 0] = pairs["t_hit1"][ps_run] + ps_np[:, 0, 1] = pairs["v_hit1"][ps_run] + ps_np[:, 0, 2] = pairs["u_hit1"][ps_run] + ps_np[:, 1, 0] = pairs["t_hit2"][ps_run] + ps_np[:, 1, 1] = pairs["v_hit2"][ps_run] + ps_np[:, 1, 2] = pairs["u_hit2"][ps_run] + ps_np[:, 2, 0] = dir_in_t[ps_run] / norm_in[ps_run] + ps_np[:, 2, 1] = dir_in_v[ps_run] / norm_in[ps_run] + ps_np[:, 2, 2] = dir_in_u[ps_run] / norm_in[ps_run] + ps_np[:, 3, 0] = dir_in_t[ps_run] / norm_out[ps_run] + ps_np[:, 3, 1] = dir_in_v[ps_run] / norm_out[ps_run] + ps_np[:, 3, 2] = dir_in_u[ps_run] / norm_out[ps_run] + ps_np[:, 4, 0] = 0.0 # WEPL is already present in the ROOT file + ps_np[:, 4, 1] = pairs["calculated_WEPL"][ps_run] + ps_np[:, 4, 2] = pairs["timestamp"][ps_run] + + df_itk = itk.GetImageFromArray(ps_np, ttype=ImageType) + + output_file = args_info.output.replace(".", f"{r:04d}.") + itk.imwrite(df_itk, output_file) + verbose(f"Wrote file {output_file}.") + + +def main(argv=None): + parser = build_parser() + args_info = parser.parse_args(argv) + process(args_info) + + +if __name__ == "__main__": + main() diff --git a/pctpairprotonsLomaLinda.cxx b/pctpairprotonsLomaLinda.cxx deleted file mode 100644 index 6841292..0000000 --- a/pctpairprotonsLomaLinda.cxx +++ /dev/null @@ -1,256 +0,0 @@ -#include "pctpairprotonsLomaLinda_ggo.h" - -#include -#include - -#include -#include -#include -#include -#include - -// Root includes -#include -#include -#include - -#define MAX_RUNS 4096 - -struct ParticleData -{ - float wepl, time; - itk::Vector position0; - itk::Vector position1; - itk::Vector position2; - itk::Vector position3; -}; - - -struct ParticleDataFinal -{ - float wepl; - itk::Vector position; - itk::Vector direction; - float time; -}; - -struct ParticleInfo -{ - int runID; -}; - - -bool -SetTreeBranch(TChain * tree, std::string branchName, void * add, bool mandatory = true) -{ - unsigned int found = 0; - tree->SetBranchStatus(branchName.c_str(), 1, &found); - if (!found) - { - if (mandatory) - { - std::cerr << "Could not load branch " << branchName << std::endl; - exit(EXIT_FAILURE); - } - } - else - tree->SetBranchAddress(branchName.c_str(), add); - return found; -} - -void -BranchParticleToPhaseSpace(struct ParticleInfo & piInput, struct ParticleData & pdInput, TChain * tree) -{ - std::cout << "In BranchParticleToPhaseSpace" << std::endl; - tree->GetListOfBranches(); // force reading of chain - SetTreeBranch(tree, "projection_angle", &piInput.runID); - ; - SetTreeBranch(tree, "calculated_WEPL", &pdInput.wepl); - SetTreeBranch(tree, "timestamp", &pdInput.time); - - // WARNING: X and Z are purposely swap... - SetTreeBranch(tree, "v_hit0", pdInput.position0.GetDataPointer() + 1); - SetTreeBranch(tree, "t_hit0", pdInput.position0.GetDataPointer()); - SetTreeBranch(tree, "u_hit0", pdInput.position0.GetDataPointer() + 2); - SetTreeBranch(tree, "v_hit1", pdInput.position1.GetDataPointer() + 1); - SetTreeBranch(tree, "t_hit1", pdInput.position1.GetDataPointer()); - SetTreeBranch(tree, "u_hit1", pdInput.position1.GetDataPointer() + 2); - SetTreeBranch(tree, "v_hit2", pdInput.position2.GetDataPointer() + 1); - SetTreeBranch(tree, "t_hit2", pdInput.position2.GetDataPointer()); - SetTreeBranch(tree, "u_hit2", pdInput.position2.GetDataPointer() + 2); - SetTreeBranch(tree, "v_hit3", pdInput.position3.GetDataPointer() + 1); - SetTreeBranch(tree, "t_hit3", pdInput.position3.GetDataPointer()); - SetTreeBranch(tree, "u_hit3", pdInput.position3.GetDataPointer() + 2); -} - - -void -WritePairs(const std::vector> & pairs, std::string fileName) -{ - itk::ImageRegion<2> region; - region.SetSize(itk::MakeSize(5, pairs.size())); - - using PixelType = itk::Vector; - using ImageType = itk::Image; - ImageType::Pointer img = ImageType::New(); - img->SetRegions(region); - img->Allocate(); - - itk::ImageRegionIterator it(img, region); - PixelType eet; - for (size_t i = 0; i < pairs.size(); i++) - { - eet[0] = pairs[i].first.wepl; - eet[1] = pairs[i].second.wepl; - eet[2] = pairs[i].second.time - pairs[i].first.time; - - it.Set(pairs[i].first.position); - ++it; - it.Set(pairs[i].second.position); - ++it; - it.Set(pairs[i].first.direction); - ++it; - it.Set(pairs[i].second.direction); - ++it; - it.Set(eet); - ++it; - } - - // Write - using WriterType = itk::ImageFileWriter; - WriterType::Pointer writer = WriterType::New(); - writer->SetFileName(fileName); - writer->SetInput(img); - TRY_AND_EXIT_ON_ITK_EXCEPTION(writer->Update()); -} - -int -main(int argc, char * argv[]) -{ - GGO(pctpairprotonsLomaLinda, args_info); // RTK macro parsing options from .ggo file (rtkMacro.h) - - // Create root trees - TChain * treeIn = new TChain("recoENTRY"); - treeIn->AddFile(args_info.input_arg); - std::cout << "Reading in file:" << args_info.input_arg << std::endl; - if (args_info.fmpct_flag) - std::cout << "In fmpct mode with ROI radius = " << args_info.roiR_arg << " mm" << std::endl; - - // Branch particles - struct ParticleInfo pi; - struct ParticleData pd; - struct ParticleDataFinal pdIn; - struct ParticleDataFinal pdOut; - - BranchParticleToPhaseSpace(pi, pd, treeIn); - - // Init - std::vector>> pairs(MAX_RUNS); - size_t nparticulesIn = treeIn->GetEntries(); - std::cout << "Number of entries = " << nparticulesIn << std::endl; - size_t iIn = 0; - size_t counterPairs = 0; - - std::cout << iIn << " particles of input phase space processed (" << 100 * iIn / nparticulesIn << "%)" << std::flush; - // Go over root files - while (iIn < nparticulesIn) - { - if (iIn % 10000 == 0) - std::cout << '\r' << iIn << " particles of input phase space processed (" << 100 * iIn / nparticulesIn << "%)" - << std::flush; - - treeIn->GetEntry(iIn); - -#if 0 - std::cout << "Input Entry = " << iIn << " with data : " - << pd.position0[0] << " " << pd.position0[1] << " " << pd.position0[2] << " " - << pd.position1[0] << " " << pd.position1[1] << " " << pd.position1[2] << " " - << pd.position2[0] << " " << pd.position2[1] << " " << pd.position2[2] << " " - << pd.position3[0] << " " << pd.position3[1] << " " << pd.position3[2] << " " - << " " << pd.wepl << std::endl; -#endif - - - // Part for proton track segment - corcular ROI intersection - // http://mathworld.wolfram.com/Circle-LineIntersection.html - bool interceptionFlag = true; - if (args_info.fmpct_flag) - { - double radius = args_info.roiR_arg; // ROI radius in mm - double dx = pd.position2[2] - pd.position1[2]; - double dy = pd.position2[0] - pd.position1[0]; - double dr_sq = (dx * dx) + (dy * dy); - double D = (pd.position1[2] * pd.position2[0]) - (pd.position2[2] * pd.position1[0]); - double Delta = (radius * radius * dr_sq) - (D * D); - double randomNumber = ((double)rand() / (RAND_MAX)); - // std::cout << "RAND_MAX = " << RAND_MAX << std::endl; - // std::cout << "randomNumber = " << randomNumber << std::endl; - // no interception, then discard proton with randomNumber > modF_arg - if (Delta < 0) - { - if (randomNumber > args_info.modF_arg) - interceptionFlag = false; - } - // std::cout << "Delta, randomNumber, interceptionFlag = " << Delta << " " << randomNumber << " " << - // interceptionFlag <= -20. && pdOut.wepl <= 900. && interceptionFlag) - { - pairs[args_info.runID_arg].push_back(std::pair(pdIn, pdOut)); - /* std::cout << "Test = " << iIn << " with data : " - << pdIn.position[0] << " " << pdIn.position[1] << " " << pdIn.position[2] << " " - << pdOut.position[0] << " " << pdOut.position[1] << " " << pdOut.position[2] << " " - << pdIn.direction[0] << " " << pdIn.direction[1] << " " << pdIn.direction[2] << " " - << pdOut.direction[0] << " " << pdOut.direction[1] << " " << pdOut.direction[2] << " " - << " " << pdIn.wepl << " " << pdOut.wepl << std::endl; */ - /*if (pdIn.position[0]>-45. && pdIn.position[0]<-10. && pdIn.position[1]>-20. && pdIn.position[1]<30.)*/ - counterPairs++; - } - - iIn++; - } - - - std::cout << "\r" << nparticulesIn << " particles of input phase space processed (" << 100 << "%)" << std::endl - << "Writing..." << std::endl; - std::cout << "number of accepted pairs = " << counterPairs << std::endl; - - for (unsigned int i = 0; i < MAX_RUNS; i++) - { - if (pairs[i].size()) - { - std::ostringstream os; - os << itksys::SystemTools::GetFilenamePath(args_info.output_arg) << "/" - << itksys::SystemTools::GetFilenameWithoutLastExtension(args_info.output_arg) std::cout - << "Writing into file:" << os.str() << std::endl; - WritePairs(pairs[i], os.str()); - } - } - - - return EXIT_SUCCESS; -} diff --git a/pctpairprotonsLomaLinda.ggo b/pctpairprotonsLomaLinda.ggo deleted file mode 100644 index 514f2a4..0000000 --- a/pctpairprotonsLomaLinda.ggo +++ /dev/null @@ -1,14 +0,0 @@ -package "pct" -version "Pair corresponding protons from Gate root files" - -option "verbose" v "Verbose execution" flag off -option "config" - "Config file" string no -option "input" i "Root phase space file of particles (entrance/exit)" string yes -option "output" o "Output file name" string yes -option "runID" - "Run number, indicating a projection in the sim" int yes -option "minRun" - "Minimum run (inclusive)" int no default="0" -option "maxRun" - "Maximum run (exclusive)" int no default="1000000" -option "nonuclear" - "Remove inelastic nuclear collisions" flag off -option "fmpct" - "Checks intersection with a cylindrical ROI" flag off -option "roiR" - "Radius of the above ROI in mm" double no default="0" -option "modF" - "Fluence level as fraction of full fluence" double no default="0" diff --git a/test/pctapplicationtest.py b/test/pctapplicationtest.py new file mode 100644 index 0000000..bd9acef --- /dev/null +++ b/test/pctapplicationtest.py @@ -0,0 +1,33 @@ +import pytest +from itk import PCT as pct +import urllib.request +import numpy as np + + +def download_file_fixture(file_key, filename): + + @pytest.fixture(scope="session") + def fixture(tmp_path_factory): + path = tmp_path_factory.getbasetemp() / filename + url = f"https://data.kitware.com/api/v1/file/{file_key}/download" + with urllib.request.urlopen(url) as response, open(path, "wb") as out_file: + out_file.write(response.read()) + return path + + return fixture + + +lomalinda_data = download_file_fixture( + "69e21803ed08a1c077afd077", "projection_045.root" +) + + +def test_lomalinda_application(tmp_path, lomalinda_data): + output = tmp_path / "baseline_lomalinda.mhd" + pct.pctlomalinda( + f"-i {lomalinda_data} -o {str(output)} --plane-in -167.2 --plane-out 167.2 --ps recoENTRY -v" + ) + # TODO finish after rebase onto main + # test_lomalinda = itk.array_from_image(itk.imread(output)) + # reference_lomalinda = itk.array_from_image(itk.imread(baseline_lomalinda)) + # assert np.array_equal(test_lomalinda, reference_lomalinda) diff --git a/wrapping/CMakeLists.txt b/wrapping/CMakeLists.txt index f245cfd..ec7e80b 100644 --- a/wrapping/CMakeLists.txt +++ b/wrapping/CMakeLists.txt @@ -8,6 +8,8 @@ list( ${PCT_SOURCE_DIR}/applications/pctfdktwodweights/pctfdktwodweights.py ${PCT_SOURCE_DIR}/applications/pctpairprotons/pctpairprotons.py ${PCT_SOURCE_DIR}/applications/pctweplfit/pctweplfit.py + ${PCT_SOURCE_DIR}/applications/pctlomalinda/pctlomalinda.py + ${PCT_SOURCE_DIR}/applications/pctfilterprotons/pctfilterprotons.py ${PCT_SOURCE_DIR}/gate/protonct.py )