Skip to content

Commit

Permalink
fixed a bug where multi voxel grids in a comp domain would cause pybi…
Browse files Browse the repository at this point in the history 
…nd errror
  • Loading branch information
@jmeneghini
jmeneghini committed Feb 3, 2024
1 parent ed4727f commit 16df802
Show file tree
Hide file tree
Showing 5 changed files with 79 additions and 184 deletions.
158 changes: 24 additions & 134 deletions cpp_simulations/lead_attenuation/lead_attenuation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,19 +8,20 @@
namespace py = pybind11;

struct DIM {
static constexpr const double X = 107.0;
static constexpr const double Y = 62.0;
static constexpr const double Z = 82.0;
static constexpr const double X = 10.7;
static constexpr const double Y = 9.2;
static constexpr const double Z = 8.2;
};

const double N_SLICES = 10;

struct SCIN_DIM {
static constexpr const double X = 10;
static constexpr const double Y = 50.0;
static constexpr const double Z = 50.0;
static constexpr const double X = 1.0;
static constexpr const double Y = 5.0;
static constexpr const double Z = 5.0;
};

Eigen::Vector3d scin_origin(9.65, 0.2, 1.6);



std::vector<std::unique_ptr<SurfaceTally>> initializeSurfaceTallies() {
Expand All @@ -29,71 +30,6 @@ std::vector<std::unique_ptr<SurfaceTally>> initializeSurfaceTallies() {
auto surface_container = SurfaceQuantityContainer();
surface_container.addVectorQuantity(VectorSurfaceQuantity(VectorSurfaceQuantityType::IncidentEnergy));

// TEST
// tallies.emplace_back(std::make_unique<RectangularSurfaceTally>(
// Eigen::Vector3d(39.0/2 - 1.5, 39.0/2 - 1.5, 0),
// Eigen::Vector3d(3, 0, 0),
// Eigen::Vector3d(0, 3, 0),
// std::move(SurfaceQuantityContainerFactory::AllQuantities())));
//

// FULL FIELD ROIS:
Eigen::Vector3d slice_y(0, SCIN_DIM::Y/N_SLICES, 0);
Eigen::Vector3d slice_z(0, 0, SCIN_DIM::Z);
Eigen::Vector3d scin_origin(12.5 + 84, 6, 16);

for (int i = 0; i < N_SLICES; i++) {
tallies.emplace_back(std::make_unique<RectangularSurfaceTally>(
scin_origin + i*slice_y,
slice_y,
slice_z,
surface_container));
}


// PENCIL BEAM ROIS:

// tallies.emplace_back(std::make_unique<RectangularSurfaceTally>(
// Eigen::Vector3d(39.0/2 - 6 - 1.5, 39.0/2 - 6 - 1.5, 180),
// Eigen::Vector3d(3, 0, 0),
// Eigen::Vector3d(0, 3, 0),
// surface_container));
//
// tallies.emplace_back(std::make_unique<RectangularSurfaceTally>(
// Eigen::Vector3d(39.0/2 - 1.5, 39.0/2 - 6 - 1.5, 180),
// Eigen::Vector3d(3, 0, 0),
// Eigen::Vector3d(0, 3, 0),
// surface_container));
//
// tallies.emplace_back(std::make_unique<RectangularSurfaceTally>(
// Eigen::Vector3d(39.0/2 - 3 - 1.5, 39.0/2 - 3 - 1.5, 180),
// Eigen::Vector3d(3, 0, 0),
// Eigen::Vector3d(0, 3, 0),
// surface_container));
//
// tallies.emplace_back(std::make_unique<RectangularSurfaceTally>(
// Eigen::Vector3d(39.0/2 - 1.5, 39.0/2 - 3 - 1.5, 180),
// Eigen::Vector3d(3, 0, 0),
// Eigen::Vector3d(0, 3, 0),
// surface_container));
//
// tallies.emplace_back(std::make_unique<RectangularSurfaceTally>(
// Eigen::Vector3d(39.0/2 - 1.5, 39.0/2 - 1.5, 180),
// Eigen::Vector3d(3, 0, 0),
// Eigen::Vector3d(0, 3, 0),
// surface_container));
//
// tallies.emplace_back(std::make_unique<RectangularSurfaceTally>(
// Eigen::Vector3d(39.0/2 + 3 - 1.5, 39.0/2 + 3 - 1.5, 180),
// Eigen::Vector3d(3, 0, 0),
// Eigen::Vector3d(0, 3, 0),
// surface_container));
//
// tallies.emplace_back(std::make_unique<RectangularSurfaceTally>(
// Eigen::Vector3d(39.0/2 + 6 - 1.5, 39.0/2 + 6 - 1.5, 180),
// Eigen::Vector3d(3, 0, 0),
// Eigen::Vector3d(0, 3, 0),
// surface_container));

return tallies;
}
Expand All @@ -103,60 +39,14 @@ std::vector<std::unique_ptr<VolumeTally>> initializeVolumeTallies() {

auto volume_container = VolumeQuantityContainer();
volume_container.addVectorQuantity(VectorVolumeQuantity(VectorVolumeQuantityType::EnergyDeposition));
volume_container.addCountQuantity(CountVolumeQuantity(CountVolumeQuantityType::NumberOfPhotons));
volume_container.addCountQuantity(CountVolumeQuantity(CountVolumeQuantityType::NumberOfInteractions));

// WHOLE BODY VOI:

tallies.emplace_back(std::make_unique<AACuboidVolumeTally>(
Eigen::Vector3d(0, 0, 155),
Eigen::Vector3d(39.0, 39.0, 175),
volume_container));


// BODY VOIS:

tallies.emplace_back(std::make_unique<AACuboidVolumeTally>(
Eigen::Vector3d(39.0/2 - 1.5, 39.0/2 - 15 - 1.5, 155 + 10 - 1.5),
Eigen::Vector3d(39.0/2 + 1.5, 39.0/2 - 15 + 1.5, 155 + 10 + 1.5),
volume_container));

tallies.emplace_back(std::make_unique<AACuboidVolumeTally>(
Eigen::Vector3d(39.0/2 - 15 - 1.5, 39.0/2 - 1.5, 155 + 10 - 1.5),
Eigen::Vector3d(39.0/2 - 15 + 1.5, 39.0/2 + 1.5, 155 + 10 + 1.5),
volume_container));

tallies.emplace_back(std::make_unique<AACuboidVolumeTally>(
Eigen::Vector3d(39.0/2 - 1.5, 39.0/2 - 1.5, 155 + 10 - 1.5),
Eigen::Vector3d(39.0/2 + 1.5, 39.0/2 + 1.5, 155 + 10 + 1.5),
volume_container));

tallies.emplace_back(std::make_unique<AACuboidVolumeTally>(
Eigen::Vector3d(39.0/2 + 15 - 1.5, 39.0/2 - 1.5, 155 + 10 - 1.5),
Eigen::Vector3d(39.0/2 + 15 + 1.5, 39.0/2 + 1.5, 155 + 10 + 1.5),
volume_container));
// WHOLE SCIN VOI:

tallies.emplace_back(std::make_unique<AACuboidVolumeTally>(
Eigen::Vector3d(39.0/2 - 1.5, 39.0/2 + 15 - 1.5, 155 + 10 - 1.5),
Eigen::Vector3d(39.0/2 + 1.5, 39.0/2 + 15 + 1.5, 155 + 10 + 1.5),
volume_container));

tallies.emplace_back(std::make_unique<AACuboidVolumeTally>(
Eigen::Vector3d(39.0/2 - 1.5, 39.0/2 - 1.5, 155 + 10 - 6 - 1.5),
Eigen::Vector3d(39.0/2 + 1.5, 39.0/2 + 1.5, 155 + 10 - 6 + 1.5),
volume_container));

tallies.emplace_back(std::make_unique<AACuboidVolumeTally>(
Eigen::Vector3d(39.0/2 - 1.5, 39.0/2 - 1.5, 155 + 10 - 3 - 1.5),
Eigen::Vector3d(39.0/2 + 1.5, 39.0/2 + 1.5, 155 + 10 - 3 + 1.5),
volume_container));

tallies.emplace_back(std::make_unique<AACuboidVolumeTally>(
Eigen::Vector3d(39.0/2 - 1.5, 39.0/2 - 1.5, 155 + 10 + 3 - 1.5),
Eigen::Vector3d(39.0/2 + 1.5, 39.0/2 + 1.5, 155 + 10 + 3 + 1.5),
volume_container));

tallies.emplace_back(std::make_unique<AACuboidVolumeTally>(
Eigen::Vector3d(39.0/2 - 1.5, 39.0/2 - 1.5, 155 + 10 + 6 - 1.5),
Eigen::Vector3d(39.0/2 + 1.5, 39.0/2 + 1.5, 155 + 10 + 6 + 1.5),
scin_origin,
scin_origin + Eigen::Vector3d(SCIN_DIM::X, SCIN_DIM::Y, SCIN_DIM::Z),
volume_container));

return tallies;
Expand All @@ -171,17 +61,15 @@ Eigen::MatrixXd processEnergySpectrum() {
PhotonSource initializeSource() {
auto energy_spectrum = processEnergySpectrum();

PolyenergeticSpectrum poly_spectrum(energy_spectrum);
std::unique_ptr<EnergySpectrum> spectrum = std::make_unique<PolyenergeticSpectrum>(poly_spectrum);
// PolyenergeticSpectrum poly_spectrum(energy_spectrum);
// std::unique_ptr<EnergySpectrum> spectrum = std::make_unique<PolyenergeticSpectrum>(poly_spectrum);
MonoenergeticSpectrum mono_spectrum(661.7E3);
std::unique_ptr<EnergySpectrum> spectrum = std::make_unique<MonoenergeticSpectrum>(mono_spectrum);

// std::unique_ptr<Directionality> directionality = std::make_unique<BeamDirectionality>(BeamDirectionality(Eigen::Vector3d(
// 39.0/2, 39.0/2, 180)));
// essentially isotropic over the whole field of view
// std::unique_ptr<Directionality> directionality = std::make_unique<BeamDirectionality>(BeamDirectionality(Eigen::Vector3d(1, DIM::Y/2, DIM::Z/2)));
std::unique_ptr<Directionality> directionality = std::make_unique<RectangularIsotropicDirectionality>(
RectangularIsotropicDirectionality(Eigen::Vector3d(1, 0, 0)
Eigen::Vector3d(0, DIM::Y, 0),
Eigen::Vector3d(0, 0, DIM::Z)));
RectangularIsotropicDirectionality(Eigen::Vector3d(10, 0, 0), Eigen::Vector3d(0, DIM::Y, 0), Eigen::Vector3d(0, 0, DIM::Z)));

std::unique_ptr<SourceGeometry> geometry = std::make_unique<PointGeometry>(PointGeometry(Eigen::Vector3d(0, DIM::Y/2, DIM::Z/2)));

PhotonSource source(std::move(spectrum), std::move(directionality), std::move(geometry));
Expand Down Expand Up @@ -285,22 +173,24 @@ void displayVoxelData(ComputationalDomain& comp_domain, int N_photons) {
std::cout << "Voxel Grid " << i << ":" << std::endl;
std::cout << "----------------" << std::endl;
VoxelGrid& voxel_grid = comp_domain.getVoxelGridN(i);
Eigen::Vector3d dim_space = voxel_grid.getDimSpace();
auto material_deposition = voxel_grid.getEnergyDepositedInMaterials();
for (auto &material: material_deposition) {
std::cout << "Material " << material.first << ": " << material.second.getSum()/N_photons << " +- "
<< material.second.getSumSTD()/N_photons << std::endl << "\n";
std::cout << "Dim Space: " << dim_space.x() << " " << dim_space.y() << " " << dim_space.z() << std::endl;
}
}
}

int main() {
ComputationalDomain comp_domain("radiography_0_degrees.json");
ComputationalDomain comp_domain("/home/john/Documents/MIDSX/cpp_simulations/lead_attenuation/lead_attenuation.json");
InteractionData interaction_data = comp_domain.getInteractionData();
PhysicsEngine physics_engine(comp_domain, interaction_data);

PhotonSource source = initializeSource();

const int NUM_OF_PHOTONS = 1000000;
const int NUM_OF_PHOTONS = 10000000;

std::cout << std::fixed << std::setprecision(15);

Expand All @@ -313,7 +203,7 @@ int main() {
auto surface_tally_results = physics_engine.getSurfaceQuantityContainers();
auto volume_tally_results = physics_engine.getVolumeQuantityContainers();

displaySurfaceTallyResults(surface_tally_results, NUM_OF_PHOTONS, comp_domain);
// displaySurfaceTallyResults(surface_tally_results, NUM_OF_PHOTONS, comp_domain);
displayVolumeTallyResults(volume_tally_results, NUM_OF_PHOTONS, comp_domain);
displayVoxelData(comp_domain, NUM_OF_PHOTONS);

Expand Down
18 changes: 13 additions & 5 deletions cpp_simulations/lead_attenuation/lead_attenuation.json
View file
Original file line number Diff line number Diff line change
@@ -1,14 +1,22 @@
{
"dim_space": [
107,
62,
82
10.7,
9.2,
8.2
],
"background_material_name": "Air, Dry (near sea level)",
"voxel_grids": [
{
"file_path": "../../data/voxels/lead_attenuation.nii.gz",
"origin": [0, 0, 0]
"file_path": "../../data/voxels/lead_block.nii.gz",
"origin": [0.2, 1.5, 0]
},
{
"file_path": "../../data/voxels/lead_block.nii.gz",
"origin": [0.2, 5.2, 0]
},
{
"file_path": "../../data/voxels/scintillator.nii.gz",
"origin": [9.65, 0.2, 1.6]
}
]
}
4 changes: 4 additions & 0 deletions src/MIDSX/Core/computational_domain.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,6 +97,10 @@ void ComputationalDomain::setVoxelGrids(const json &json_object, const std::stri
getNIFTIFilePaths(voxel_grid_json, json_directory_path, nifti_file_paths);
getOrigins(voxel_grid_json, origins);
}
std::unique_ptr<py::scoped_interpreter> guard;
if (!is_python_environment_) {
guard = std::make_unique<py::scoped_interpreter>(); // start python interpreter only if not already started
}
for (int i = 0; i < nifti_file_paths.size(); i++) {
voxel_grids_.emplace_back(VoxelGrid(nifti_file_paths[i], is_python_environment_), origins[i]);
}
Expand Down
6 changes: 1 addition & 5 deletions src/MIDSX/Core/voxel_grid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -107,11 +107,7 @@ std::unordered_map<int, VectorValue> VoxelGrid::getEnergyDepositedInMaterials()


void VoxelGrid::initializeVoxels() {
std::unique_ptr<py::scoped_interpreter> guard;
if (!is_python_environment_) {
guard = std::make_unique<py::scoped_interpreter>(); // start python interpreter only if not already started
}
py::module nibabel = py::module::import("nibabel"); // import nibabel
py::module nibabel = py::module::import("nibabel");
py::object img = nibabel.attr("load")(filename_); // load the nifti file into img
py::object data = img.attr("get_fdata")(); // get the data from img
py::object header = img.attr("header");
Expand Down
Loading

0 comments on commit 16df802

Please sign in to comment.