diff --git a/core/opengate_core/opengate_core.cpp b/core/opengate_core/opengate_core.cpp index 669857fa75..c6eaa47534 100644 --- a/core/opengate_core/opengate_core.cpp +++ b/core/opengate_core/opengate_core.cpp @@ -550,6 +550,10 @@ void init_GatePhaseSpaceSource(py::module &); void init_GateGANPairSource(py::module &); +void init_GateWindowTurboSource(py::module &); + +void init_GateVoxelWTSource(py::module &); + // Gate misc void init_GateExceptionHandler(py::module &); @@ -794,6 +798,8 @@ PYBIND11_MODULE(opengate_core, m) { init_GateGANSource(m); init_GatePhaseSpaceSource(m); init_GateGANPairSource(m); + init_GateWindowTurboSource(m); + init_GateVoxelWTSource(m); init_GateSPSPosDistribution(m); init_GateSPSVoxelsPosDistribution(m); init_GateRunAction(m); diff --git a/core/opengate_core/opengate_lib/GateSingleParticleSourceWindowTurbo.cpp b/core/opengate_core/opengate_lib/GateSingleParticleSourceWindowTurbo.cpp new file mode 100644 index 0000000000..02fe56c34b --- /dev/null +++ b/core/opengate_core/opengate_lib/GateSingleParticleSourceWindowTurbo.cpp @@ -0,0 +1,287 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include "GateSingleParticleSourceWindowTurbo.h" +#include "G4Threading.hh" +#include "GateHelpersDict.h" +#include "Randomize.hh" +#include +#include +#include +#include +#include + +GateSingleParticleSourceWindowTurbo::GateSingleParticleSourceWindowTurbo( + std::string mother_volume) + : GateSingleParticleSource(mother_volume) {} + +void GateSingleParticleSourceWindowTurbo::SetParameters( + G4double a1, G4double a2, G4double b1, G4double b2, G4double plane_distance, + G4double plane_phi) { + fA1 = a1; + fA2 = a2; + fB1 = b1; + fB2 = b2; + fPlaneDistance = plane_distance; + fPlanePhi = plane_phi; + fSinPlanePhi = sin(plane_phi); + fCosPlanePhi = cos(plane_phi); +} + +G4double solid_angle_pyramid(G4double a, G4double b, G4double d) { + return 4 * atan(a * b / (2 * d * sqrt(a * a + b * b + 4 * d * d))); +} + +G4double GateSingleParticleSourceWindowTurbo::GetSolidAngle( + const G4ThreeVector &pos) const { + + // rotate with -plane_phi + if (pos.x() * fCosPlanePhi + pos.y() * fSinPlanePhi >= fPlaneDistance) { + G4String error_msg = fmt::format( + "position ({}, {}, {}) is outside the plane distance {} for source: {}", + pos.x(), pos.y(), pos.z(), fPlaneDistance, fSourceName); + G4Exception("GateSingleParticleSourceWindowTurbo::GetSolidAngle", + "GetSolidAngleError", FatalException, error_msg); + } + + G4double x0 = pos.x() * fCosPlanePhi + pos.y() * fSinPlanePhi; + G4double y0 = -pos.x() * fSinPlanePhi + pos.y() * fCosPlanePhi; + G4double a1_rel = fA1 - y0; + G4double a2_rel = fA2 - y0; + G4double b1_rel = fB1 - pos.z(); + G4double b2_rel = fB2 - pos.z(); + G4double d_rel = fPlaneDistance - x0; + G4double sa11 = solid_angle_pyramid(2 * a1_rel, 2 * b1_rel, d_rel); + G4double sa12 = solid_angle_pyramid(2 * a1_rel, 2 * b2_rel, d_rel); + G4double sa21 = solid_angle_pyramid(2 * a2_rel, 2 * b1_rel, d_rel); + G4double sa22 = solid_angle_pyramid(2 * a2_rel, 2 * b2_rel, d_rel); + G4double sa = sa11 + sa22 - sa12 - sa21; + return fabs(sa * 0.25); +} + +void GateSingleParticleSourceWindowTurbo::ThreadFunc( + size_t count, G4double *act_ratio_all_thread, + G4double *max_solid_angle_thread) { + + G4ThreeVector pos; + for (size_t i = 0; i < count; i++) { + pos = fPositionGenerator->VGenerateOne(); + G4double solid_angle = GetSolidAngle(pos); + if (solid_angle > *max_solid_angle_thread) + *max_solid_angle_thread = solid_angle; + *act_ratio_all_thread += solid_angle / 4 / M_PI; + } +} + +G4double GateSingleParticleSourceWindowTurbo::InitializeBeforeRun( + G4double &act_ratio, G4double &max_solid_angle) { + auto start_time = std::chrono::high_resolution_clock::now(); + G4double act_ratio_all = 0; + fMaxSolidAngle = 0; + const size_t sampling_count_per_thread = + fSamplingCountInit / fThreadCountInit; + std::vector threads(fThreadCountInit); + std::vector act_ratio_all_thread(fThreadCountInit, 0); + std::vector max_solid_angle_thread(fThreadCountInit, 0); + for (G4int i = 0; i < fThreadCountInit; i++) { + threads[i] = + std::thread(&GateSingleParticleSourceWindowTurbo::ThreadFunc, this, + sampling_count_per_thread, &act_ratio_all_thread[i], + &max_solid_angle_thread[i]); + } + for (G4int i = 0; i < fThreadCountInit; i++) { + threads[i].join(); + act_ratio_all += act_ratio_all_thread[i]; + if (max_solid_angle_thread[i] > fMaxSolidAngle) + fMaxSolidAngle = max_solid_angle_thread[i]; + } + + act_ratio = act_ratio_all / sampling_count_per_thread / fThreadCountInit; + max_solid_angle = fMaxSolidAngle; + auto end_time = std::chrono::high_resolution_clock::now(); + auto duration = std::chrono::duration_cast( + end_time - start_time); + return duration.count() / 1e6; +} + +void GateSingleParticleSourceWindowTurbo::InitializeUserInfo( + py::dict &user_info) { + // TODO: make this run in master thread before each run + // and make the worker thread get info properly before each run + // fSourceName = source_name; + fSamplingCountInit = DictGetInt(user_info, "init_sampling_count"); + fThreadCountInit = DictGetInt(user_info, "init_number_of_threads"); + fSkip = DictGetBool(user_info, "skip_mode"); + + // fActRatio = DictGetDouble(user_info, "act_ratio"); + // fMaxSolidAngle = DictGetDouble(user_info, "max_solid_angle"); + // if (not isnan(fActRatio) && not isnan(fMaxSolidAngle) and fActRatio >= 0 + // and + // fActRatio <= 1 and fMaxSolidAngle >= 0 and fMaxSolidAngle <= 4 * M_PI) + // { + // return; + // } else { + // fMaxSolidAngle = 0; + // fActRatio = 0; + // } + + // if (not G4Threading::IsMasterThread()) { + // TBD: should I check validity of act_ratio and max_solid_angle here or in + // python side? + // if (isnan(fActRatio) || isnan(fMaxSolidAngle)) { + // G4String error_msg = + // "activity ratio or max solid angle not set for source: "; + // error_msg += fSourceName; + // G4Exception("GateSingleParticleSourceWindowTurbo::Initialize", + // "InitializeError", FatalException, error_msg); + // } + // return; + // } + + // TODO: check paramenters in python + // if (a1 != a1 || a2 != a2 || b1 != b1 || b2 != b2 || + // plane_distance != plane_distance || plane_phi != plane_phi) { + // G4Exception("GateWindowTurboSource::SetActRatio", "SetActRatioError", + // FatalException, "Not all parameters needed points are set"); + // } + + // if (a1 >= a2 || b1 >= b2) { + // G4Exception("GateWindowTurboSource::SetActRatio", "SetActRatioError", + // FatalException, "a1 >= a2 or b1 >= b2"); + // } + + // VerifyPhiTheta(samplingCount, 0.01); +} + +void GateSingleParticleSourceWindowTurbo::SetPhiTheta( + const G4ThreeVector &pos) const { + // compare theta with cos2, to avoid complex calculation + // G4double cot2theta = std::copysign(1.0, dir.z()) * dir.z() * dir.z() / + // (dir.x() * dir.x() + dir.y() * dir.y()); + // relationship between angular vector and theta and phi in Geant4 + // px = -sintheta * cosphi; + // py = -sintheta * sinphi; + // pz = -costheta; + + G4double x0 = pos.x() * fCosPlanePhi + pos.y() * fSinPlanePhi; + G4double y0 = -pos.x() * fSinPlanePhi + pos.y() * fCosPlanePhi; + G4double a1_rel = fA1 - y0; + G4double a2_rel = fA2 - y0; + G4double b1_rel = fB1 - pos.z(); + G4double b2_rel = fB2 - pos.z(); + G4double d_rel = fPlaneDistance - x0; + + G4double aamax = std::max(a1_rel * a1_rel, a2_rel * a2_rel); + G4double aamin = std::min(a1_rel * a1_rel, a2_rel * a2_rel); + G4double thetamax, thetamin; + + if (a1_rel < 0 and a2_rel > 0 and b2_rel > 0) + thetamax = M_PI - atan2(d_rel, b2_rel); + else + thetamax = M_PI - atan2(sqrt((b2_rel > 0 ? aamin : aamax) + d_rel * d_rel), + b2_rel); + + if (a1_rel < 0 and a2_rel > 0 and b1_rel < 0) + // in this case, need to check minmum/maxmum of the hyperbola + thetamin = M_PI - atan2(d_rel, b1_rel); + else + thetamin = M_PI - atan2(sqrt((b1_rel > 0 ? aamax : aamin) + d_rel * d_rel), + b1_rel); + + fDirectionGenerator->SetMinTheta(thetamin); + fDirectionGenerator->SetMaxTheta(thetamax); + G4double phimin = atan2(a1_rel, d_rel) + fPlanePhi; + G4double phimax = atan2(a2_rel, d_rel) + fPlanePhi; + + fDirectionGenerator->SetMinPhi(phimin + M_PI); + fDirectionGenerator->SetMaxPhi(phimax + M_PI); +} + +G4bool GateSingleParticleSourceWindowTurbo::CheckPosDirValid( + const G4ThreeVector &pos, const G4ThreeVector &dir) const { + // compare theta with cos2, to avoid complex calculation + // G4double cot2theta = std::copysign(1.0, dir.z()) * dir.z() * dir.z() / + // (dir.x() * dir.x() + dir.y() * dir.y()); + + G4double x0 = pos.x() * fCosPlanePhi + pos.y() * fSinPlanePhi; + G4double y0 = -pos.x() * fSinPlanePhi + pos.y() * fCosPlanePhi; + G4double a1_rel = fA1 - y0; + G4double a2_rel = fA2 - y0; + G4double b1_rel = fB1 - pos.z(); + G4double b2_rel = fB2 - pos.z(); + G4double d_rel = fPlaneDistance - x0; + G4double dir_x_rotated = dir.x() * fCosPlanePhi + dir.y() * fSinPlanePhi; + G4double dir_y_rotated = -dir.x() * fSinPlanePhi + dir.y() * fCosPlanePhi; + + G4double intersect_b = d_rel / dir_x_rotated * dir.z() + pos.z(); + G4double intersect_a = d_rel / dir_x_rotated * dir_y_rotated + y0; + return intersect_a <= fA2 && intersect_a >= fA1 && intersect_b <= fB2 && + intersect_b >= fB1; +} + +void GateSingleParticleSourceWindowTurbo::GeneratePos() { + fCurrentPos = fPositionGenerator->VGenerateOne(); + + // probability of the position is valid should be proportional to the solid + // angle + while (true) { + G4double solid_angle = GetSolidAngle(fCurrentPos); + if (solid_angle > fMaxSolidAngle * 1.1) { + G4String error_msg = "solid angle of position"; + error_msg += fmt::format(" ({}, {}, {}): {}", fCurrentPos.x(), + fCurrentPos.y(), fCurrentPos.z(), solid_angle); + error_msg += " is larger than max solid angle "; + error_msg += std::to_string(fMaxSolidAngle); + error_msg += " for source: "; + error_msg += fSourceName; + error_msg += "\nyou may increase max solid angle and try again"; + G4Exception("GateWindowTurboSource::GeneratePrimaryVertex", + "GeneratePrimaryVertexError", FatalException, error_msg); + } + if (G4UniformRand() < solid_angle / fMaxSolidAngle / 1.1) { + fCurrentSolidAngle = solid_angle; + break; + } + fCurrentPos = fPositionGenerator->VGenerateOne(); + } + fPosGenerated = true; +} + +void GateSingleParticleSourceWindowTurbo::GeneratePrimaryVertex( + G4Event *event) { + if (not fSkip) + GeneratePos(); + fPosGenerated = false; + SetPhiTheta(fCurrentPos); + G4ThreeVector direction; + while (true) { + direction = fDirectionGenerator->VGenerateOne(); + if (CheckPosDirValid(fCurrentPos, direction)) { + break; + } + } + G4PrimaryVertex *vertex = new G4PrimaryVertex(fCurrentPos, particle_time); + + // Set placement relative to attached volume + // DD(particle_momentum_direction); + + G4double energy = fEnergyGenerator->VGenerateOne(fParticleDefinition); + + // one single particle + auto *particle = new G4PrimaryParticle(fParticleDefinition); + particle->SetKineticEnergy(energy); + particle->SetMass(fMass); + particle->SetMomentumDirection(direction); + particle->SetCharge(fCharge); + particle->SetWeight(1.0); + if (fPolarizationFlag) + particle->SetPolarization(fPolarization); + + // set vertex + vertex->SetPrimary(particle); + event->AddPrimaryVertex(vertex); +} diff --git a/core/opengate_core/opengate_lib/GateSingleParticleSourceWindowTurbo.h b/core/opengate_core/opengate_lib/GateSingleParticleSourceWindowTurbo.h new file mode 100644 index 0000000000..eb2bbff4ae --- /dev/null +++ b/core/opengate_core/opengate_lib/GateSingleParticleSourceWindowTurbo.h @@ -0,0 +1,61 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#ifndef GateSingleParticleSourceWindowTurbo_h +#define GateSingleParticleSourceWindowTurbo_h + +#include "GateSingleParticleSource.h" +#include +#include +#include + +namespace py = pybind11; + +class GateSingleParticleSourceWindowTurbo : public GateSingleParticleSource { +public: + explicit GateSingleParticleSourceWindowTurbo(std::string mother_volume); + ~GateSingleParticleSourceWindowTurbo() override = default; + void GeneratePrimaryVertex(G4Event *event) override; + void InitializeUserInfo(py::dict &user_info); + G4double GetCurrentSolidAngle() const { return fCurrentSolidAngle; } + void GeneratePos(); + void SetSkipMode(G4bool skip) { fSkip = skip; } + G4bool PosGenerated() const { return fPosGenerated; } + void SetParameters(G4double a1, G4double a2, G4double b1, G4double b2, + G4double plane_distance, G4double plane_phi); + void SetMaxSolidAngle(G4double max_solid_angle) { + fMaxSolidAngle = max_solid_angle; + } + G4double GetMaxSolidAngle() const { return fMaxSolidAngle; } + G4double InitializeBeforeRun(G4double &act_ratio, G4double &max_solid_angle); + +private: + G4double GetSolidAngle( + const G4ThreeVector &pos) const; // get solid angle for the window + G4bool CheckPosDirValid(const G4ThreeVector &pos, + const G4ThreeVector &dir) + const; // check if the ray can pass through the window + void SetPhiTheta( + const G4ThreeVector &pos) const; // set the phi and theta of the direction + // distribution according to the position + G4double fPlaneDistance{NAN}; + G4double fPlanePhi{NAN}; + G4double fSinPlanePhi{NAN}, fCosPlanePhi{NAN}; + G4double fA1{NAN}, fA2{NAN}, fB1{NAN}, fB2{NAN}; + G4double fMaxSolidAngle = 0; + G4String fSourceName; + G4double fCurrentSolidAngle; + G4ThreeVector fCurrentPos; + G4bool fSkip; + G4int fThreadCountInit; + G4int fSamplingCountInit; + G4bool fPosGenerated = false; + void ThreadFunc(size_t count, G4double *act_ratio_all_thread, + G4double *max_solid_angle_thread); +}; + +#endif // GateSingleParticleSourceWindowTurbo_h diff --git a/core/opengate_core/opengate_lib/GateVoxelWTSource.cpp b/core/opengate_core/opengate_lib/GateVoxelWTSource.cpp new file mode 100644 index 0000000000..83b982e037 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateVoxelWTSource.cpp @@ -0,0 +1,19 @@ +#include "GateVoxelWTSource.h" +#include "GateSPSVoxelsPosDistribution.h" +#include "GateWindowTurboSource.h" + +GateVoxelWTSource::GateVoxelWTSource() : GateWindowTurboSource() { + fVoxelPositionGenerator = new GateSPSVoxelsPosDistribution(); +} + +void GateVoxelWTSource::PrepareNextRun() { + GateWindowTurboSource::PrepareNextRun(); + + fVoxelPositionGenerator->fGlobalRotation = fGlobalRotation; + fVoxelPositionGenerator->fGlobalTranslation = fGlobalTranslation; +} + +void GateVoxelWTSource::InitializePosition(py::dict) { + fSPS->SetPosGenerator(fVoxelPositionGenerator); + fVoxelPositionGenerator->SetPosDisType("Point"); +} diff --git a/core/opengate_core/opengate_lib/GateVoxelWTSource.h b/core/opengate_core/opengate_lib/GateVoxelWTSource.h new file mode 100644 index 0000000000..04827c39a1 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateVoxelWTSource.h @@ -0,0 +1,16 @@ +#include "GateWindowTurboSource.h" +class GateSPSVoxelsPosDistribution; + +class GateVoxelWTSource : public GateWindowTurboSource { +public: + GateVoxelWTSource(); + ~GateVoxelWTSource() = default; + void PrepareNextRun() override; + GateSPSVoxelsPosDistribution *GetSPSVoxelPosDistribution() { + return fVoxelPositionGenerator; + } + +protected: + void InitializePosition(py::dict user_info) override; + GateSPSVoxelsPosDistribution *fVoxelPositionGenerator; +}; diff --git a/core/opengate_core/opengate_lib/GateWindowTurboSource.cpp b/core/opengate_core/opengate_lib/GateWindowTurboSource.cpp new file mode 100644 index 0000000000..d2ceeb8201 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateWindowTurboSource.cpp @@ -0,0 +1,305 @@ +#include "GateWindowTurboSource.h" +#include "G4CallbackModel.hh" +#include "GateGenericSource.h" +#include "GateHelpersDict.h" +#include "GateSingleParticleSourceWindowTurbo.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +void GateWindowTurboSource::SetSharedCache( + std::shared_ptr cache) { + fSharedCache = std::move(cache); +} + +void GateWindowTurboSource::CreateSPS() { + fSPS = new GateSingleParticleSourceWindowTurbo(fAttachedToVolumeName); +} + +void GateWindowTurboSource::InitializeUserInfo(py::dict &user_info) { + GateGenericSource::InitializeUserInfo(user_info); + fWeight = -1; + fWeightSigma = -1; + fDirectionRelativeToAttachedVolume = false; + fUserInfo = user_info; +} + +double GateWindowTurboSource::CalcNextTime(double current_simulation_time) { + GateSingleParticleSourceWindowTurbo *spswt = + reinterpret_cast(fSPS); + G4double act_ratio; + if (fSkip) { + if (not spswt->PosGenerated()) { + spswt->GeneratePos(); + } + act_ratio = spswt->GetCurrentSolidAngle() / (4 * M_PI); + } else + act_ratio = fCurrentActRatio; + + double next_time = current_simulation_time; + if ((fMaxN <= 0)) { + next_time = current_simulation_time - + log(G4UniformRand()) * (1.0 / fActivity / act_ratio); + } + return next_time; +} + +void GateWindowTurboSource::InitializeSharedCache(py::dict &user_info) { + if (!fSharedCache->fActRatio.empty()) + return; + std::lock_guard lock(fSharedCache->fMutex); + if (py::isinstance(user_info["act_ratio"])) { + fSharedCache->fActRatio = {DictGetDouble(user_info, "act_ratio")}; + fSharedCache->fMaxSolidAngle = { + DictGetDouble(user_info, "max_solid_angle")}; + } else { + fSharedCache->fActRatio = DictGetVecDouble(user_info, "act_ratio"); + fSharedCache->fMaxSolidAngle = + DictGetVecDouble(user_info, "max_solid_angle"); + } +} + +void GateWindowTurboSource::WriteBackUserInfo() { + if (!fUserInfo) + return; + py::gil_scoped_acquire acquire; // write back needs to acquire gil + auto direction = py::dict(fUserInfo["direction"]); + direction["act_ratio"] = fSharedCache->fActRatio; + direction["max_solid_angle"] = fSharedCache->fMaxSolidAngle; + direction["init_duration"] = fSharedCache->fInitDuration; +} + +void GateWindowTurboSource::PrepareSharedBeforeRun() { + fCurrentRunId = G4RunManager::GetRunManager()->GetCurrentRun()->GetRunID(); + GateSingleParticleSourceWindowTurbo *spswt = + reinterpret_cast(fSPS); + fCurrentActRatio = GetValueThisRun(fSharedCache->fActRatio); + G4double max_solid_angle = GetValueThisRun(fSharedCache->fMaxSolidAngle); + if (max_solid_angle >= 0 and fCurrentActRatio >= 0) + return; + if (fSkip) + return; + + std::lock_guard lock(fSharedCache->fMutex); + + G4double a1 = GetValueThisRun(fA1); + G4double a2 = GetValueThisRun(fA2); + G4double b1 = GetValueThisRun(fB1); + G4double b2 = GetValueThisRun(fB2); + G4double plane_distance = GetValueThisRun(fPlaneDistance); + G4double plane_phi = GetValueThisRun(fPlanePhi); + spswt->SetParameters(a1, a2, b1, b2, plane_distance, plane_phi); + G4double duration_sec = + spswt->InitializeBeforeRun(fCurrentActRatio, max_solid_angle); + + fSharedCache->fInitDuration.push_back(duration_sec); + SetValueThisRun(fSharedCache->fActRatio, fCurrentActRatio); + SetValueThisRun(fSharedCache->fMaxSolidAngle, max_solid_angle); + WriteBackUserInfo(); +} + +void GateWindowTurboSource::PrepareNextRun() { + GateGenericSource::PrepareNextRun(); + // TBD: voxelized source prepare next run here + auto *ang = fSPS->GetAngDist(); + auto *spswt = reinterpret_cast(fSPS); + ang->fGlobalRotation = G4RotationMatrix(); + + // setup act ratio and max solid angle once per shared cache/run + PrepareSharedBeforeRun(); + + G4double a1 = GetValueThisRun(fA1); + G4double a2 = GetValueThisRun(fA2); + G4double b1 = GetValueThisRun(fB1); + G4double b2 = GetValueThisRun(fB2); + G4double plane_distance = GetValueThisRun(fPlaneDistance); + G4double plane_phi = GetValueThisRun(fPlanePhi); + spswt->SetParameters(a1, a2, b1, b2, plane_distance, plane_phi); + const G4double max_solid_angle = + GetValueThisRun(fSharedCache->fMaxSolidAngle); + spswt->SetMaxSolidAngle(max_solid_angle); +} + +void GateWindowTurboSource::Visualize() const { + if (G4Threading::GetNumberOfRunningWorkerThreads() > 0 and fVisCount > 0) { + G4Exception("GateWindowTurboSource::Visualize", "VisualizeWTSourceInMTMode", + JustWarning, + "Visualize for GateWindowTurboSource is not supported in MT " + "mode. The origin of the source will be wrong."); + } + GateGenericSource::Visualize(); + if (visualization_window_color.size() > 0) { + for (size_t i = 0; i < visualization_window_color.size(); i++) { + VisualizeOneWindow(visualization_window_color[i], + visualization_window_width[i], + visualization_window_run_id[i]); + } + } +} + +void GateWindowTurboSource::InitializeDirection(py::dict puser_info) { + + auto *ang = fSPS->GetAngDist(); + ang->SetAngDistType("iso"); + auto user_info = py::dict(puser_info["direction"]); + if (py::isinstance(user_info["a1"])) { + fA1 = {DictGetDouble(user_info, "a1")}; + fA2 = {DictGetDouble(user_info, "a2")}; + fB1 = {DictGetDouble(user_info, "b1")}; + fB2 = {DictGetDouble(user_info, "b2")}; + fPlaneDistance = {DictGetDouble(user_info, "plane_distance")}; + fPlanePhi = {DictGetDouble(user_info, "plane_phi")}; + } else { + fA1 = DictGetVecDouble(user_info, "a1"); + fA2 = DictGetVecDouble(user_info, "a2"); + fB1 = DictGetVecDouble(user_info, "b1"); + fB2 = DictGetVecDouble(user_info, "b2"); + fPlaneDistance = DictGetVecDouble(user_info, "plane_distance"); + fPlanePhi = DictGetVecDouble(user_info, "plane_phi"); + } + InitializeSharedCache(user_info); + + fSkip = DictGetBool(user_info, "skip_mode"); + + GateSingleParticleSourceWindowTurbo *spswt = + reinterpret_cast(fSPS); + spswt->InitializeUserInfo(user_info); + if (fAAManager == nullptr) { + fAAManager = new GateAcceptanceAngleManager; + fSPS->SetAAManager(fAAManager); + } + if (fFDManager == nullptr) { + fFDManager = new GateForcedDirectionManager; + fSPS->SetFDManager(fFDManager); + } +} +void GateWindowTurboSource::GetWindowVertex(G4ThreeVector &pos1, + G4ThreeVector &pos2, + G4ThreeVector &pos3, + G4ThreeVector &pos4, + G4int run_id) const { + G4double a1 = GetValueThisRun(fA1, run_id); + G4double a2 = GetValueThisRun(fA2, run_id); + G4double b1 = GetValueThisRun(fB1, run_id); + G4double b2 = GetValueThisRun(fB2, run_id); + G4double plane_distance = GetValueThisRun(fPlaneDistance, run_id); + G4double plane_phi = GetValueThisRun(fPlanePhi, run_id); + pos1 = {plane_distance, a1, b1}; + pos2 = {plane_distance, a1, b2}; + pos3 = {plane_distance, a2, b1}; + pos4 = {plane_distance, a2, b2}; + // rotate with plane_phi + G4double s = sin(plane_phi); + G4double c = cos(plane_phi); + G4RotationMatrix rot({{c, s, 0}, {-s, c, 0}, {0, 0, 1}}); + pos1 = rot * pos1; + pos2 = rot * pos2; + pos3 = rot * pos3; + pos4 = rot * pos4; +} + +namespace { +G4Color GetColor(const py::handle &color_py) { + if (py::isinstance(color_py)) { + G4Color color; + const std::string color_str = color_py.cast(); + G4Color::GetColor(color_str, color); + return color; + } + + const auto rgba = color_py.cast>(); + + if (rgba.size() == 3) + return {rgba[0], rgba[1], rgba[2], 1.0}; + return {rgba[0], rgba[1], rgba[2], rgba[3]}; +} +std::vector DictGetVecColor(py::dict &user_info, + const std::string &key) { + std::vector l; + auto color_list = py::list(user_info[key.c_str()]); + for (const auto color : color_list) { + l.push_back(GetColor(color)); + } + return l; +} + +} // namespace + +void GateWindowTurboSource::InitializeVisualization(py::dict puser_info) { + GateGenericSource::InitializeVisualization(puser_info); + auto user_info = py::dict(puser_info["visualization"]); + visualization_window_run_id = DictGetVecInt(user_info, "window_run_id"); + visualization_window_width = DictGetVecDouble(user_info, "window_width"); + visualization_window_color = DictGetVecColor(user_info, "window_color"); +} + +void GateWindowTurboSource::VisualizeOneWindow(G4Color color, G4double width, + int run_id) const { + G4ThreeVector pos1, pos2, pos3, pos4; + GetWindowVertex(pos1, pos2, pos3, pos4, run_id); + + VisWindow *window = new VisWindow(pos1, pos2, pos3, pos4, color, width); + G4VModel *model = + new G4CallbackModel(window); + model->SetType("Turbo Window"); + model->SetGlobalTag("Turbo Window"); + G4String description = "Turbo Window: "; + description += "(" + std::to_string(pos1.x()) + " " + + std::to_string(pos1.y()) + " " + std::to_string(pos1.z()) + + "), "; + description += "(" + std::to_string(pos2.x()) + " " + + std::to_string(pos2.y()) + " " + std::to_string(pos2.z()) + + "), "; + description += "(" + std::to_string(pos3.x()) + " " + + std::to_string(pos3.y()) + " " + std::to_string(pos3.z()) + + "), "; + description += "(" + std::to_string(pos4.x()) + " " + + std::to_string(pos4.y()) + " " + std::to_string(pos4.z()) + + ")"; + model->SetGlobalDescription(description); + G4cout << "Visualizing window for run " << run_id << ": " << description + << G4endl; + + G4VisManager *fpVisManager = G4VisManager::GetInstance(); + G4Scene *pScene = fpVisManager->GetCurrentScene(); + const G4String ¤tSceneName = pScene->GetName(); + G4bool successful = pScene->AddRunDurationModel(model, true); + G4UImanager::GetUIpointer()->ApplyCommand("/vis/scene/notifyHandlers"); +} + +GateWindowTurboSource::VisWindow::VisWindow(const G4Vector3D &pos1, + const G4Vector3D &pos2, + const G4Vector3D &pos3, + const G4Vector3D &pos4, + G4Color color, G4double width) { + + fPolyline.push_back(pos1); + fPolyline.push_back(pos3); + fPolyline.push_back(pos4); + fPolyline.push_back(pos2); + fPolyline.push_back(pos1); + G4VisAttributes va; + va.SetLineWidth(width); + va.SetColor(color); + fPolyline.SetVisAttributes(va); +} + +void GateWindowTurboSource::VisWindow::operator()( + G4VGraphicsScene &sceneHandler, const G4ModelingParameters *) { + sceneHandler.BeginPrimitives(); + sceneHandler.AddPrimitive(fPolyline); + sceneHandler.EndPrimitives(); +} diff --git a/core/opengate_core/opengate_lib/GateWindowTurboSource.h b/core/opengate_core/opengate_lib/GateWindowTurboSource.h new file mode 100644 index 0000000000..0940abb9eb --- /dev/null +++ b/core/opengate_core/opengate_lib/GateWindowTurboSource.h @@ -0,0 +1,99 @@ +#pragma once +#include "GateGenericSource.h" +#include "GateSingleParticleSource.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/* +author: LiKun (likun@dotuai.com/tontyoutoure@gmail.com) +source only generated particles that fulfill the following conditions: + +Given four points in space, pth1, pth2, pphi1, pphi2 +then elevation angle of the particle, theta from the source point, should +between the elevation angles of source point seeing of pth1 and pth2. Also the +azimuthal angle of the particle, phi from the source point, should between the +azimuthal angles of source point seeing pphi1 and pphi2. + +*/ + +class GateSingleParticleSourceWindowTurbo; + +struct GateWindowTurboSharedCache { + std::mutex fMutex; + std::vector fActRatio; + std::vector fMaxSolidAngle; + std::vector fInitDuration; +}; + +class GateWindowTurboSource : public GateGenericSource { +public: + GateWindowTurboSource() = default; + ~GateWindowTurboSource() override = default; + + virtual void PrepareNextRun() override; + void SetSharedCache(std::shared_ptr cache); + + void VisualizeOneWindow(G4Color color, G4double width, int run_id) const; + void InitializeUserInfo(py::dict &user_info) override; + virtual double CalcNextTime(double current_simulation_time) override; + virtual void Visualize() const override; + +protected: + virtual void CreateSPS() override; + virtual void InitializeDirection(py::dict puser_info) override; + virtual void InitializeVisualization(py::dict user_info) override; + +private: + std::vector fA1, fA2, fB1, fB2, fPlaneDistance, fPlanePhi; + G4int fCurrentRunId; + G4double fCurrentActRatio; + G4double GetValueThisRun(const std::vector &vec) const { + return GetValueThisRun(vec, fCurrentRunId); + } + G4double GetValueThisRun(const std::vector &vec, + G4int run_id) const { + if (vec.size() == 1) + return vec[0]; + else + return vec[run_id]; + } + void SetValueThisRun(std::vector &vec, G4double value) { + if (vec.size() == 1) + vec[0] = value; + else + vec[fCurrentRunId] = value; + } + std::vector visualization_window_color; + std::vector visualization_window_width; + std::vector visualization_window_run_id; + + struct VisWindow { + VisWindow(const G4Vector3D &pos1, const G4Vector3D &pos2, + const G4Vector3D &pos3, const G4Vector3D &pos4, G4Color color, + G4double width); + void operator()(G4VGraphicsScene &, const G4ModelingParameters *); + G4Polyline fPolyline; + G4Color fColor; + G4double fWidth; + }; + void GetWindowVertex(G4ThreeVector &pos1, G4ThreeVector &pos2, + G4ThreeVector &pos3, G4ThreeVector &pos4, + int run_id) const; + + void PrepareSharedBeforeRun(); + void InitializeSharedCache(py::dict &user_info); + void WriteBackUserInfo(); + py::dict fUserInfo; // to write back act ratio and max solid angle + + std::shared_ptr fSharedCache; + bool fSkip; // true for act ratio, false for event skipping base on + // solid angle +}; diff --git a/core/opengate_core/opengate_lib/pyGateVoxelWTSource.cpp b/core/opengate_core/opengate_lib/pyGateVoxelWTSource.cpp new file mode 100644 index 0000000000..6f82412052 --- /dev/null +++ b/core/opengate_core/opengate_lib/pyGateVoxelWTSource.cpp @@ -0,0 +1,23 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include + +namespace py = pybind11; + +#include "GateSPSVoxelsPosDistribution.h" +#include "GateVoxelWTSource.h" + +void init_GateVoxelWTSource(py::module &m) { + + py::class_(m, "GateVoxelWTSource") + .def(py::init()) + .def("GetSPSVoxelPosDistribution", + &GateVoxelWTSource::GetSPSVoxelPosDistribution, + py::return_value_policy::reference_internal) + .def("InitializeUserInfo", &GateVoxelWTSource::InitializeUserInfo); +} diff --git a/core/opengate_core/opengate_lib/pyGateWindowTurboSource.cpp b/core/opengate_core/opengate_lib/pyGateWindowTurboSource.cpp new file mode 100644 index 0000000000..fb2cae326b --- /dev/null +++ b/core/opengate_core/opengate_lib/pyGateWindowTurboSource.cpp @@ -0,0 +1,26 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include + +namespace py = pybind11; + +#include "GateWindowTurboSource.h" + +void init_GateWindowTurboSource(py::module &m) { + + py::class_>( + m, "GateWindowTurboSharedCache") + .def(py::init<>()); + + py::class_(m, + "GateWindowTurboSource") + .def(py::init()) + .def("SetSharedCache", &GateWindowTurboSource::SetSharedCache) + .def("InitializeUserInfo", &GateWindowTurboSource::InitializeUserInfo); +} diff --git a/docs/source/figures/visualize_window_turbo_source.png b/docs/source/figures/visualize_window_turbo_source.png new file mode 100755 index 0000000000..3d920c5138 Binary files /dev/null and b/docs/source/figures/visualize_window_turbo_source.png differ diff --git a/docs/source/figures/window_turbo_source.png b/docs/source/figures/window_turbo_source.png new file mode 100644 index 0000000000..c4b44fafca Binary files /dev/null and b/docs/source/figures/window_turbo_source.png differ diff --git a/docs/source/figures/window_turbo_source_definition.png b/docs/source/figures/window_turbo_source_definition.png new file mode 100644 index 0000000000..aa0517cfdf Binary files /dev/null and b/docs/source/figures/window_turbo_source_definition.png differ diff --git a/docs/source/user_guide/user_guide_reference_sources.rst b/docs/source/user_guide/user_guide_reference_sources.rst index 32b81670e2..ec02ce98a3 100644 --- a/docs/source/user_guide/user_guide_reference_sources.rst +++ b/docs/source/user_guide/user_guide_reference_sources.rst @@ -16,4 +16,5 @@ Details: Sources user_guide_reference_sources_gan_source.rst user_guide_reference_sources_phid_source.rst user_guide_reference_sources_phase_space_source.rst + .. user_guide_reference_sources_window_turbo_source.rst diff --git a/docs/source/user_guide/user_guide_reference_sources_window_turbo_source.rst b/docs/source/user_guide/user_guide_reference_sources_window_turbo_source.rst new file mode 100644 index 0000000000..a9bf43a487 --- /dev/null +++ b/docs/source/user_guide/user_guide_reference_sources_window_turbo_source.rst @@ -0,0 +1,261 @@ +.. _source-window-turbo-source: + +Window Turbo Source +=================== + +Description +----------- + +``WindowTurboSource`` is used to make a source emit photons only toward a rectangular window, +thereby reducing the generation of useless particles and the tracking time. +In scenarios such as multi-pinhole collimator SPECT, where only a small number of directions can reach the detector, +it can significantly shorten the simulation time; +the actual speedup depends on the window size, the distance from the source to the window, the source position distribution, +and the downstream geometry. + +The purpose of this source is to replace ``GenericSource`` emission with ``direction.type = "iso"`` +while preserving the spatial and temporal statistical properties of the resulting counts. +See ``test100`` for related comparisons. Each ``WindowTurboSource`` can define only one rectangular window. +The window must lie on a reference plane parallel to the z axis, and its two pairs of edges must be parallel to the z axis +and to the x-y plane, respectively. + +.. image:: ../figures/window_turbo_source.png + +``WindowTurboSource`` inherits the particle, energy, and position settings from ``GenericSource``, +but ``direction`` is replaced by the window parameters. +Usually, ``activity`` should be used to control the emission rate; the ``n`` parameter is not compatible with +``WindowTurboSource``. +This source is mainly intended for ``gamma`` particles, and the ``back_to_back`` particle type is not available. + +Basic usage +----------- + +The following example creates a cylindrical volume source and makes the photons emit only toward a rectangular window +in the y direction. +The window size is given by ``a1``, ``a2``, ``b1``, and ``b2``; the plane on which the window lies is given by +``plane_distance`` and ``plane_phi``. + +.. code:: python + + import numpy as np + import opengate as gate + + Bq = gate.g4_units.Bq + keV = gate.g4_units.keV + mm = gate.g4_units.mm + + source = sim.add_source("WindowTurboSource", "wts") + source.particle = "gamma" + source.activity = 1e6 * Bq + source.energy.mono = 141 * keV + + source.position.type = "cylinder" + source.position.translation = [0, -100 * mm, 0] + source.position.radius = 100 * mm + source.position.dz = 100 * mm + + radius = 13.6 * mm + source.direction.a1 = -radius + source.direction.a2 = radius + source.direction.b1 = -radius + source.direction.b2 = radius + source.direction.plane_distance = 86 * mm + source.direction.plane_phi = np.pi / 2 + +.. note:: + + Although ``WindowTurboSource`` inherits from ``GenericSource``, the ``direction`` parameters of ``GenericSource`` + (for example, ``type``, ``theta``, ``phi``, ``momentum``, ``focus_point``, and so on) are not used. + +Window parameters +----------------- + +The additional parameters of ``WindowTurboSource`` are all located in ``source.direction``: + +.. list-table:: + :header-rows: 1 + + * - **Parameter** + - **Description** + * - ``a1`` + - Left boundary of the window, relative to the center point of the reference plane; must be smaller than ``a2``. + * - ``a2`` + - Right boundary of the window, relative to the center point of the reference plane. + * - ``b1`` + - Lower boundary of the window, relative to the center point of the reference plane; must be smaller than ``b2``. + * - ``b2`` + - Upper boundary of the window, relative to the center point of the reference plane. + * - ``plane_distance`` + - Distance from the reference plane to the system center; must be positive. + * - ``plane_phi`` + - Angle between the reference plane normal vector and the positive x axis, in radians, in the range ``[0, 2*pi)``. + * - ``init_sampling_count`` + - Number of samples used during initialization to estimate the window acceptance ratio. The default value is ``1000000``. + * - ``init_number_of_threads`` + - Number of threads used during initialization. The default value ``0`` means that the number of simulation threads is used. + * - ``act_ratio`` + - Window acceptance ratio. Usually left empty at its default value, so it is estimated automatically during initialization. + * - ``max_solid_angle`` + - Maximum solid angle of the window over the source position distribution. Usually left empty at its default value, so it is estimated automatically during initialization. + * - ``skip_mode`` + - Advanced option, with a default value of ``False``. Its usage is not described in this document. + +``a1``, ``a2``, ``b1``, ``b2``, ``plane_distance``, and ``plane_phi`` jointly define the position and size of the window. +``a1`` and ``a2`` define the window extent within the reference plane in the direction perpendicular to the z axis; +``b1`` and ``b2`` define the window extent in the z direction. +``plane_distance`` defines the distance from the reference plane to the system center, and ``plane_phi`` defines the +orientation of this plane around the z axis. + +.. image:: ../figures/window_turbo_source_definition.png + +If a simulation contains multiple timing intervals, different window parameters can be used for each time interval. +The window parameters listed above, as well as ``act_ratio`` and ``max_solid_angle``, can be written as a single value +or as a list. +The list length can be ``1``, which means that the same value is used for all time intervals, or it can be equal to the +number of ``sim.run_timing_intervals``, which means that values are set time interval by time interval. + +.. code:: python + + sec = gate.g4_units.second + sim.run_timing_intervals = [[0, 1 * sec], [1 * sec, 2 * sec]] + + source.direction.a1 = [-10 * mm, -15 * mm] + source.direction.a2 = [10 * mm, 15 * mm] + source.direction.b1 = -10 * mm + source.direction.b2 = 10 * mm + source.direction.plane_distance = [80 * mm, 90 * mm] + source.direction.plane_phi = np.pi / 2 + +Initialization +-------------- + +``WindowTurboSource`` needs to estimate two quantities: + +* ``act_ratio``: on average, the fraction of ``GenericSource`` isotropic events that would pass through the window. +* ``max_solid_angle``: the maximum solid angle subtended by the window over the source position distribution. + +If the user does not provide these two parameters, the simulation estimates them automatically before the start of each +timing interval. +The estimation process samples from the source position distribution; the number of samples is controlled by +``init_sampling_count`` and the number of initialization threads is controlled by ``init_number_of_threads``. +After the simulation ends, the estimated values are written back into ``source.direction.act_ratio`` and +``source.direction.max_solid_angle``. +The initialization duration is also written back into ``source.direction.init_duration``, in seconds. This is an output +value for diagnostics, not an input parameter to configure the source. +For repeated simulations with the same configuration, the written-back ``act_ratio`` and ``max_solid_angle`` values can +be saved and set directly in the next run to reduce initialization time. + +.. note:: + + If the simulation is run with ``sim.run(start_new_process=True)``, the values written back in the child process are + not retained in the source object of the current Python process. + +Voxelization +------------ + +``VoxelWTSource`` is the voxelized version of ``WindowTurboSource``. +Like ``VoxelSource``, it samples the initial position from a 3D activity image, while the emission direction is +controlled by the window parameters of ``WindowTurboSource``. +Therefore, when using ``VoxelWTSource``, the ``image``, particle, and energy parameters must be set, as well as the same +set of ``direction`` window parameters. + +.. code:: python + + source = sim.add_source("VoxelWTSource", "voxel_wts") + source.image = "activity.mhd" + source.particle = "gamma" + source.activity = 1e6 * Bq + source.energy.mono = 141 * keV + + source.direction.a1 = -13.6 * mm + source.direction.a2 = 13.6 * mm + source.direction.b1 = -13.6 * mm + source.direction.b2 = 13.6 * mm + source.direction.plane_distance = 86 * mm + source.direction.plane_phi = np.pi / 2 + +Except for ``direction``, the meaning, normalization, and position settings of the voxel source image are the same as +for ``VoxelSource``. +If the activity image needs to be aligned with a CT image or another voxelized volume, see :ref:`source-voxel-source`. + +Dynamic activity image +~~~~~~~~~~~~~~~~~~~~~~ + +``VoxelWTSource`` can use a dynamic activity image in the same way as ``VoxelSource``. +After configuring the source as above, set ``sim.run_timing_intervals`` and provide one image per timing interval with +``add_dynamic_parametrisation(image=[...])``. + +.. code:: python + + source.add_dynamic_parametrisation( + image=[ + "activity_0.mhd", + "activity_1.mhd", + ] + ) + +The number of images must match the number of ``sim.run_timing_intervals``. +See :ref:`source-voxel-source` and :doc:`user_guide_dynamic_parametrisations` for the corresponding ``VoxelSource`` +usage. + +Visualization +------------- + +``WindowTurboSource`` and ``VoxelWTSource`` use the ``visualization`` attribute +to display source-position sampling points, like ``GenericSource``. +In addition, the ``visualization.window_run_id``, ``visualization.window_color``, +and ``visualization.window_width`` fields can draw one or more rectangular window +outlines: + +.. code:: python + + sim.visu = True + sim.visu_type = "qt" + + source.visualization.count = 1000 + source.visualization.color = "red" + source.visualization.size = 2 + + source.visualization.window_run_id = 0 + source.visualization.window_color = "red" + source.visualization.window_width = 2 + +The following visualization result is generated by ``test100_window_turbo_source_visu_wip.py``. + +.. image:: ../figures/visualize_window_turbo_source.png + +The visualization parameters specific to ``WindowTurboSource`` and +``VoxelWTSource`` are: + +.. list-table:: + :header-rows: 1 + + * - **Parameter** + - **Description** + * - ``window_color`` + - Window line color. It must be set when ``window_run_id`` is not empty. Use a color name such as ``"red"``, ``"green"``, or ``"blue"`` to reuse one color for all displayed windows, or use a list with one color per window. RGB or RGBA colors must be provided as color entries, for example ``[[1, 0, 0]]`` for one red window. + * - ``window_width`` + - Window line width. It must be set when ``window_run_id`` is not empty. Use a single number to reuse one width for all displayed windows, or use a list with one number per window. Values must be in the range ``(0, 10]``. + * - ``window_run_id`` + - Index, or list of indexes, of the timing intervals whose windows should be displayed. Indexes start from ``0`` and must be smaller than the number of ``sim.run_timing_intervals``. The default value is an empty list, so no window outline is drawn unless this field is set. + +When displaying several timing intervals, either give ``window_color`` and ``window_width`` as scalar values to reuse +the same style for every window, or give lists with the same length as ``window_run_id``. +For example: + +.. code:: python + + source.visualization.window_run_id = [0, 1] + source.visualization.window_color = ["red", "blue"] + source.visualization.window_width = 2 + +.. note:: + + Visualization of ``WindowTurboSource`` and ``VoxelWTSource`` currently does not support multithreaded mode. + Use ``sim.number_of_threads = 1`` when ``sim.visu`` is enabled; otherwise, the displayed source origin may be + incorrect. + +Implementation details +---------------------- + +FIXME (details to be added) diff --git a/opengate/managers.py b/opengate/managers.py index 48b6128c80..05f846e478 100644 --- a/opengate/managers.py +++ b/opengate/managers.py @@ -41,6 +41,8 @@ from .sources.generic import GenericSource, SourceBase from .sources.lastvertexsources import LastVertexSource from .sources.phidsources import PhotonFromIonDecaySource +from .sources.windowturbosource import WindowTurboSource +from .sources.voxelwtsource import VoxelWTSource from .sources.phspsources import PhaseSpaceSource from .sources.voxelsources import VoxelizedPromptGammaTLESource, VoxelSource from .utility import ( @@ -64,6 +66,8 @@ "PhotonFromIonDecaySource": PhotonFromIonDecaySource, "TreatmentPlanPBSource": TreatmentPlanPBSource, "VoxelizedPromptGammaTLESource": VoxelizedPromptGammaTLESource, + "WindowTurboSource": WindowTurboSource, + "VoxelWTSource": VoxelWTSource, } from .actors.chemistryactors import ChemistryActorBase, ChemicalCountingActor diff --git a/opengate/sources/voxelwtsource.py b/opengate/sources/voxelwtsource.py new file mode 100644 index 0000000000..2f8f2b4c2f --- /dev/null +++ b/opengate/sources/voxelwtsource.py @@ -0,0 +1,70 @@ +import itk + +import opengate_core as g4 + +from .windowturbosource import WindowTurboSource +from .voxelsources import VoxelSource +from ..utility import ensure_filename_is_str +from ..base import process_cls +from ..actors.dynamicactors import SourceActivityImageChanger + + +class VoxelWTSource(WindowTurboSource): + # basically the same as VoxelSource, just avoiding diamond inheritance + + # hints for IDE + image: str + + user_info_defaults = VoxelSource.user_info_defaults + + def __init__(self, *args, **kwargs): + WindowTurboSource.__init__(self, *args, **kwargs) + # the loaded image + self._current_itk_image = None + # cached CDFs + self._cdf_x = None + self._cdf_y = None + self._cdf_z = None + + def create_g4_source(self): + g4_source = g4.GateVoxelWTSource() + g4_source.SetSharedCache(self._g4_shared_cache) + return g4_source + + def create_changers(self): + changers = super().create_changers() + for dp in self.dynamic_params.values(): + if dp["extra_params"]["auto_changer"] is True: + if "image" in dp: + new_changer = SourceActivityImageChanger( + name=f"{self.name}_source_activity_changer_{len(changers)}", + activity_images=dp["image"], + attached_to=self, + simulation=self.simulation, + ) + changers.append(new_changer) + else: + self.warning( + f"You need to manually create a changer for dynamic parametrisation {dp} " + f"of source '{self.name}'." + ) + return changers + + def set_transform_from_user_info(self, g4_source): + VoxelSource.set_transform_from_user_info(self, g4_source) + + def cumulative_distribution_functions(self, g4_source): + VoxelSource.cumulative_distribution_functions(self, g4_source) + + def update_activity_image(self, filename): + VoxelSource.update_activity_image(self, filename) + + def initialize_g4_source(self, g4_source, run_timing_intervals): + if self._current_itk_image is None: + self._current_itk_image = itk.imread(ensure_filename_is_str(self.image)) + self.set_transform_from_user_info(g4_source) + self.cumulative_distribution_functions(g4_source) + WindowTurboSource.initialize_g4_source(self, g4_source, run_timing_intervals) + + +process_cls(VoxelWTSource) diff --git a/opengate/sources/windowturbosource.py b/opengate/sources/windowturbosource.py new file mode 100644 index 0000000000..4e791a4b22 --- /dev/null +++ b/opengate/sources/windowturbosource.py @@ -0,0 +1,316 @@ +import opengate_core as g4 +from .generic import GenericSource, VisualizationValidator +from box import Box +from ..base import UserInfoValidatorBase +from ..base import process_cls +from ..exception import fatal, warning +import numpy as np +from ..logger import logger +import os + + +def _wts_direction_parameters(): + return Box( + { + "a1": [], + "a2": [], + "b1": [], + "b2": [], + "plane_distance": [], + "plane_phi": [], + "init_sampling_count": 1000000, + "init_number_of_threads": 0, + "act_ratio": [], + "max_solid_angle": [], + "skip_mode": False, + } + ) + + +def _wts_visualization_parameters(): + return Box( + { + "window_color": [], + "window_width": [], + "window_run_id": [], + "count": 2000, + "color": "yellow", + "size": 2, + } + ) + + +class WTSDirectionValidator(UserInfoValidatorBase): + """Validates the 'direction' Box.""" + + __schema__ = set(_wts_direction_parameters().keys()) + + def set_simulation(self, simulation): + self.simulation = simulation + self.act_ratio_inited = True + self.max_solid_angle_inited = True + + def is_integer(self, val): + return isinstance(val, (int, np.integer)) + + def is_numeric(self, val): + return isinstance(val, (int, float, np.number)) + + def validate_attr_against_nti(self, b, attr_name): + attr = b[attr_name] + if isinstance(attr, list): + if len(attr) != self.num_intervals and len(attr) != 1: + fatal( + f"'{self.context_name}.{attr_name}' must be a list of length 1 or {self.num_intervals} (number of timing intervals)." + ) + for i, v in enumerate(attr): + if not self.is_numeric(v): + fatal( + f"All elements of '{self.context_name}.{attr_name}' must be numbers." + ) + b[attr_name][i] = float(v) + elif self.is_numeric(attr): + logger.debug( + f"'{self.context_name}.{attr_name}' is converted to a list of length 1." + ) + b[attr_name] = [float(attr)] + else: + fatal( + f"'{self.context_name}.{attr_name}' must be a number or a list of numbers." + ) + + def get_attr_interval(self, b, attr_name, interval_index): + attr = b[attr_name] + if len(attr) == 1: + return attr[0] + else: + return attr[interval_index] + + def validate(self, parent_obj, attr_name: str, parent_context: str = None): + self.context_name = super().validate(parent_obj, attr_name, parent_context) + self.num_intervals = len(self.simulation.run_timing_intervals) + b = getattr(parent_obj, attr_name) + + if isinstance(b.act_ratio, list) and len(b.act_ratio) == 0: + b.act_ratio = [-1] * self.num_intervals # default value indicating not set + self.act_ratio_inited = False + if isinstance(b.max_solid_angle, list) and len(b.max_solid_angle) == 0: + b.max_solid_angle = [ + -1 + ] * self.num_intervals # default value indicating not set + self.max_solid_angle_inited = False + + nti_sized_attrs = [ + "a1", + "a2", + "b1", + "b2", + "act_ratio", + "max_solid_angle", + "plane_distance", + "plane_phi", + ] + for attr in nti_sized_attrs: + self.validate_attr_against_nti(b, attr) + + for int_idx in range(self.num_intervals): + a1 = self.get_attr_interval(b, "a1", int_idx) + a2 = self.get_attr_interval(b, "a2", int_idx) + if a1 >= a2: + fatal( + f"'a1' must be less than 'a2' for timing interval {int_idx} in '{self.context_name}'." + ) + b1 = self.get_attr_interval(b, "b1", int_idx) + b2 = self.get_attr_interval(b, "b2", int_idx) + if b1 >= b2: + fatal( + f"'b1' must be less than 'b2' for timing interval {int_idx} in '{self.context_name}'." + ) + max_solid_angle = self.get_attr_interval(b, "max_solid_angle", int_idx) + if ( + max_solid_angle <= 0 or max_solid_angle > 2 * np.pi + ) and self.max_solid_angle_inited: + fatal( + f"'max_solid_angle' must be in the range (0, 2*pi] for timing interval {int_idx} in '{self.context_name}'." + ) + act_ratio = self.get_attr_interval(b, "act_ratio", int_idx) + if (act_ratio < 0 or act_ratio > 1) and self.act_ratio_inited: + fatal( + f"'act_ratio' must be in the range [0, 1] for timing interval {int_idx} in '{self.context_name}'." + ) + plane_distance = self.get_attr_interval(b, "plane_distance", int_idx) + if plane_distance <= 0: + fatal(f"'plane_distance' must be positive in '{self.context_name}'.") + plane_phi = self.get_attr_interval(b, "plane_phi", int_idx) + if plane_phi < 0 or plane_phi >= np.pi * 2: + fatal( + f"'plane_phi' must be in the range [0, np.pi * 2) in '{self.context_name}'." + ) + if not self.is_integer(b.init_sampling_count) or b.init_sampling_count <= 0: + fatal(f"'init_sampling_count' must be positive in '{self.context_name}'.") + + if b.init_number_of_threads == 0: + b.init_number_of_threads = self.simulation.number_of_threads + logger.debug( + f"'init_number_of_threads' is set to the number of CPU cores: {b.init_number_of_threads}." + ) + + if ( + not self.is_integer(b.init_number_of_threads) + or b.init_number_of_threads < 0 + or b.init_number_of_threads > os.cpu_count() + ): + warning( + f"'init_number_of_threads' must be a positive integer less than or equal to the number of CPU cores. Setting it to {os.cpu_count()}." + ) + + if not isinstance(b.skip_mode, bool): + fatal(f"'skip_mode' must be a boolean in '{self.context_name}'.") + + +class WTSVisualizationValidator(VisualizationValidator): + """Validates the visualization parameters for WindowTurboSource.""" + + __schema__ = set(_wts_visualization_parameters().keys()) + + def set_simulation(self, simulation): + self.simulation = simulation + + def validate_width(self, width, context): + if not isinstance(width, (int, float)): + fatal(f"'window_width' must be a number in '{context}'.") + if width <= 0 or width > 10: + fatal(f"'window_width' must be in the range (0, 10] in '{context}'.") + + def validate_run_id(self, run_id, context): + if not isinstance(run_id, int): + fatal(f"'window_run_id' must be an integer in '{context}'.") + + if run_id < 0 or run_id >= len(self.simulation.run_timing_intervals): + fatal( + f"'window_run_id' must be between 0 and {len(self.simulation.run_timing_intervals) - 1} in '{context}'." + ) + + def validate_list_length(self, lst, context): + if len(lst) == 1: + lst = lst * self.max_list_length + elif len(lst) != self.max_list_length: + fatal(f"Length of list must be 1 or {self.max_list_length} in '{context}'.") + + def validate(self, parent_obj, attr_name: str, parent_context: str = None): + self.context_name = super().validate(parent_obj, attr_name, parent_context) + b = getattr(parent_obj, attr_name) + self.max_list_length = ( + 1 if not isinstance(b.window_run_id, list) else len(b.window_run_id) + ) + + if isinstance(b.window_run_id, list): + for run_id in b.window_run_id: + self.validate_run_id(run_id, f"{self.context_name}.window_run_id") + else: + self.validate_run_id(b.window_run_id, f"{self.context_name}.window_run_id") + b.window_run_id = [b.window_run_id] + + if not isinstance(b.window_color, list): + self.validate_color(b.window_color, f"{self.context_name}.window_color[0]") + b.window_color = [b.window_color] * self.max_list_length + else: + for i, color in enumerate(b.window_color): + self.validate_color(color, f"{self.context_name}.window_color[{i}]") + self.validate_list_length( + b.window_color, f"{self.context_name}.window_color" + ) + + if not isinstance(b.window_width, list): + self.validate_width(b.window_width, f"{self.context_name}") + b.window_width = [b.window_width] * self.max_list_length + else: + for i, width in enumerate(b.window_width): + self.validate_width(width, f"{self.context_name}") + self.validate_list_length( + b.window_width, f"{self.context_name}.window_width" + ) + + +class WindowTurboSource(GenericSource): + direction: Box + + user_info_defaults = { + "direction": ( + _wts_direction_parameters(), + {"doc": "Define the direction of the primary particles.", "override": True}, + ), + "visualization": ( + _wts_visualization_parameters(), + { + "doc": "Define the visualization parameters for the source.", + "override": True, + }, + ), + } + + def __init__(self, *args, **kwargs): + GenericSource.__init__(self, *args, **kwargs) + self._g4_shared_cache = g4.GateWindowTurboSharedCache() + self._dir_validator = WTSDirectionValidator() + self._dir_validator.set_simulation(self.simulation) + self._visu_validator = WTSVisualizationValidator() + self._visu_validator.set_simulation(self.simulation) + + def create_g4_source(self): + g4_source = g4.GateWindowTurboSource() + g4_source.SetSharedCache(self._g4_shared_cache) + return g4_source + + def initialize_g4_source(self, g4_source, run_timing_intervals): + + if self.particle == "back_to_back": + fatal( + "The 'back_to_back' particle type is not compatible with WindowTurboSource." + ) + if isinstance(self.n, list): + if any(n_i != 0 for n_i in self.n): + fatal( + "The 'n' parameter must be 0 for all timing intervals in WindowTurboSource." + ) + elif isinstance(self.n, (int, float)): + if self.n != 0: + fatal("The 'n' parameter is not compatible with WindowTurboSource.") + elif isinstance(self.n, np.ndarray): + if np.any(self.n != 0): + fatal( + "The 'n' parameter must be 0 for all timing intervals in WindowTurboSource." + ) + else: + fatal("Invalid type for 'n' parameter in WindowTurboSource.") + + if self.particle != "gamma": + warning( + f"Particle type '{self.particle}' is not 'gamma'. WindowTurboSource is designed for gamma primary purpose only. Proceed ONLY if you know what you are doing." + ) + + GenericSource.initialize_g4_source(self, g4_source, run_timing_intervals) + + def can_predict_number_of_events(self): + # Actually, can, but initialization is needed. However act ratio is dependent on mother volume movement, therefore cannot be predicted before the simulation starts. + return False + + def _visualize_window( + self, color, width: float = 2.0, timing_interval_index: int = 0 + ): + if ( + timing_interval_index >= len(self.simulation.run_timing_intervals) + or timing_interval_index < 0 + ): + fatal( + f"Invalid timing interval index {timing_interval_index}. Must be between 0 and {len(self.simulation.run_timing_intervals) - 1}." + ) + self._visu_validator.validate_color(color) + if isinstance(color, str): + color = color.lower() + self.visualization.window_color.append(color) + self.visualization.window_width.append(width) + self.visualization.window_run_id.append(timing_interval_index) + + +process_cls(WindowTurboSource) diff --git a/opengate/tests/src/source/test100_window_turbo_source_base.py b/opengate/tests/src/source/test100_window_turbo_source_base.py new file mode 100644 index 0000000000..6297e73343 --- /dev/null +++ b/opengate/tests/src/source/test100_window_turbo_source_base.py @@ -0,0 +1,158 @@ +import opengate as gate +import numpy as np +from opengate.tests import utility +import matplotlib.pyplot as plt + + +def build_collimator(sim, head, pin_radius_up=3.6, pin_radius_down=13.6): + world = sim.world + gcm3 = gate.g4_units.g_cm3 + mm = gate.g4_units.mm + sim.volume_manager.material_database.add_material_weights( + "Tungsten", + ["W"], + [1], + 19.3 * gcm3, + ) + + pinboard_inner = sim.add_volume("Box", "pinboard_inner") + pinboard_inner.mother = head + pinboard_inner.material = "Tungsten" + pinboard_inner.translation = [0, -59.5 * mm, 0] + pinboard_inner.size = [500 * mm, 2 * mm, 500 * mm] + pinboard_inner.color = [0.5, 0.5, 0.5, 0.5] + + # kill_actor_inner = sim.add_actor("KillActor", "kill_inner") + # kill_actor_inner.attached_to = "pinboard_inner" + + pinboard_cylinder = sim.add_volume("Tubs", "pin_cylinder") + pinboard_cylinder.mother = pinboard_inner + pinboard_cylinder.rmin = 0 + pinboard_cylinder.rmax = pin_radius_up * mm + pinboard_cylinder.dz = 1 * mm + pinboard_cylinder.material = "G4_AIR" + pinboard_cylinder.rotation = np.array([[1, 0, 0], [0, 0, -1], [0, 1, 0]]) + pinboard_cylinder.color = [1, 1, 1, 1] + + pinboard_outer = sim.add_volume("Box", "pinboard_outer") + pinboard_outer.mother = head + pinboard_outer.material = "Tungsten" + pinboard_outer.translation = [0, -67 * mm, 0] + pinboard_outer.size = [500 * mm, 12.99999 * mm, 500 * mm] + pinboard_outer.color = [0.5, 0.5, 0.5, 0.5] + + pin_cone = sim.add_volume("Cons", "pin_cone") + pin_cone.mother = pinboard_outer + pin_cone.rmax1 = pin_radius_up * mm + pin_cone.rmax2 = pin_radius_down * mm + pin_cone.rmin1 = 0 + pin_cone.rmin2 = 0 + pin_cone.dz = (12.99999 / 2) * mm + pin_cone.material = "G4_AIR" + pin_cone.rotation = np.array([[1, 0, 0], [0, 0, -1], [0, 1, 0]]) + pin_cone.color = [1, 1, 1, 1] + pin_cone.dphi = 2 * np.pi + + # kill_actor_outer = sim.add_actor("KillActor", "kill_outer") + # kill_actor_outer.attached_to = "pinboard_outer" + + +def build_crystal(sim, head, prj_name): + world = sim.world + gcm3 = gate.g4_units.g_cm3 + sim.volume_manager.material_database.add_material_weights( + "CsI", + ["Cs", "I"], + [1, 1], + 4.51 * gcm3, + ) + + mm = gate.g4_units.mm + + head_crystal = sim.add_volume("Box", "head_crystal") + head_crystal.mother = head + head_crystal.size = [160 * mm, 8 * mm, 160 * mm] + head_crystal.material = "CsI" + head_crystal.translation = [0, 69.5 * mm, 0] + head_crystal.color = [0, 0, 1, 1] + hc = sim.add_actor("DigitizerHitsCollectionActor", "Hits") + hc.attributes = [ + "TotalEnergyDeposit", + "KineticEnergy", + "PostPosition", + "TrackCreatorProcess", + "GlobalTime", + "TrackVolumeName", + "RunID", + "ThreadID", + "TrackID", + "PreStepUniqueVolumeID", + ] + hc.attached_to = ["head_crystal"] + sc = sim.add_actor("DigitizerAdderActor", "Singles") + sc.input_digi_collection = "Hits" + sc.policy = "EnergyWeightedCentroidPosition" + sc.group_volume = "head_crystal" + proj = sim.add_actor("DigitizerProjectionActor", "Projection") + proj.attached_to = "head_crystal" + proj.input_digi_collections = ["Singles"] + proj.spacing = [1.5 * mm, 1.5 * mm] + proj.size = [100, 100] + proj.origin_as_image_center = False + proj.output_filename = f"{prj_name}.mhd" + proj.detector_orientation_matrix = np.array([[1, 0, 0], [0, 0, -1], [0, 1, 0]]) + + +def build_geometry( + sim, prj_name, pin_radius_up=3.6, pin_radius_down=13.6, head_y_pos=100 +): + # head_y_pos is the distance between center of the collimator and the system center in y direction + mm = gate.g4_units.mm + world = sim.world + world.size = [500 * mm, 1000 * mm, 500 * mm] + world.color = [1, 1, 1, 0.05] + + head = sim.add_volume("Box", "head") + head.mother = world + head.size = [500 * mm, 147 * mm, 500 * mm] + head.material = "G4_AIR" + head.translation = [0, (73.5 + head_y_pos - 14) * mm, 0] + head.color = [0, 1, 0, 0.05] + + build_collimator(sim, head, pin_radius_up, pin_radius_down) + build_crystal(sim, head, prj_name) + + +def calculate_profile(image_path): + import SimpleITK as sitk + + image = sitk.ReadImage(image_path) + array = sitk.GetArrayFromImage(image) + profile = np.sum(array, axis=(0, 1)) + return profile + + +def compare_profiles(ref, test, tolerance=8.0, fig_name=None): + ref = np.asarray(ref, dtype=float) + test = np.asarray(test, dtype=float) + + if ref.shape != test.shape: + utility.print_test(False, f"Profile shapes differ: {ref.shape} vs {test.shape}") + return False + + sad = np.abs(ref - test).sum() / (ref.sum() + test.sum()) * 100 + is_ok = sad < tolerance + utility.print_test( + is_ok, f"Profile relative SAD = {sad:.2f}% (tol {tolerance:.2f}%)" + ) + + if fig_name is not None: + plt.figure(figsize=(10, 5)) + plt.plot(ref / ref.sum(), label="reference") + plt.plot(test / test.sum(), label="test") + plt.legend() + plt.xlabel("Pixel") + plt.ylabel("Normalized counts") + plt.savefig(fig_name) + + return is_ok diff --git a/opengate/tests/src/source/test100_window_turbo_source_skip_wip.py b/opengate/tests/src/source/test100_window_turbo_source_skip_wip.py new file mode 100644 index 0000000000..dbd1f16306 --- /dev/null +++ b/opengate/tests/src/source/test100_window_turbo_source_skip_wip.py @@ -0,0 +1,193 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +from opengate.tests import utility +import pathlib +import numpy as np +from box import Box +import matplotlib.pyplot as plt +from opengate.tests import utility +from test100_window_turbo_source_base import build_geometry + +paths = utility.get_default_test_paths(__file__, output_folder="test100") + + +def change_source_parameters(source_back, source_1, source_2, source_3): + keV = gate.g4_units.keV + mm = gate.g4_units.mm + source_back.particle = "gamma" + source_back.energy.mono = 141 * keV + source_back.position.type = "cylinder" + source_back.position.translation = [0, -100, 0] + source_back.position.radius = 100 * mm + source_back.position.dz = 100 * mm + + source_1.particle = "gamma" + source_1.energy.mono = 141 * keV + source_1.position.type = "cylinder" + source_1.position.translation = [-50, -100, 0] + source_1.position.radius = 10 * mm + source_1.position.dz = 100 * mm + + source_2.particle = "gamma" + source_2.energy.mono = 141 * keV + source_2.position.type = "cylinder" + source_2.position.translation = [0, -100, 0] + source_2.position.radius = 15 * mm + source_2.position.dz = 100 * mm + + source_3.particle = "gamma" + source_3.energy.mono = 141 * keV + source_3.position.type = "cylinder" + source_3.position.translation = [50, -100, 0] + source_3.position.radius = 20 * mm + source_3.position.dz = 100 * mm + + +def initialize(duration=10): + sim = gate.Simulation() + sec = gate.g4_units.second + sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3" + sim.physics_manager.global_production_cuts.all = 1 * gate.g4_units.mm + # main options + sim.g4_verbose = False + sim.g4_verbose_level = 1 + sim.visu = False + sim.visu_type = "qt" + sim.number_of_threads = 32 + sim.progress_bar = True + sim.run_timing_intervals = [[0, duration * sec]] + sim.add_actor("SimulationStatisticsActor", "Stats") + return sim + + +def calculate_profile(image_path): + import SimpleITK as sitk + + image = sitk.ReadImage(image_path) + array = sitk.GetArrayFromImage(image) + profile = np.sum(array, axis=(0, 1)) + return profile + + +def compare_profiles(ref, test, tolerance=8.0, fig_name=None): + ref = np.asarray(ref, dtype=float) + if fig_name is not None: + plt.figure(figsize=(10, 5)) + plt.plot(ref / ref.sum(), label="reference") + is_ok = True + for i, t in enumerate(test): + test_1 = np.asarray(t, dtype=float) + + if ref.shape != test_1.shape: + utility.print_test( + False, f"Profile shapes differ: {ref.shape} vs {test_1.shape}" + ) + return False + + sad = np.abs(ref - test_1).sum() / (ref.sum() + test_1.sum()) * 100 + is_ok = is_ok and (sad < tolerance) + + if fig_name is not None: + plt.plot(test_1 / test_1.sum(), label=f"test {i}") + if fig_name is not None: + plt.legend() + plt.xlabel("Pixel") + plt.ylabel("Normalized counts") + plt.savefig(fig_name) + + utility.print_test( + is_ok, f"Profile relative SAD = {sad:.2f}% (tol {tolerance:.2f}%)" + ) + return is_ok + + +def run_window_turbo_source(activity=1000000, skip_mode=False): + Bq = gate.g4_units.Bq + mm = gate.g4_units.mm + NoT = 4 + duration = 320 / NoT + sim = initialize(duration) + sim.g4_verbose = False + sim.number_of_threads = NoT + sim.progress_bar = False + sim.random_seed = 1 + radius_down = 13.6 + build_geometry( + sim, paths.output / f"window_turbo_{skip_mode}", pin_radius_down=radius_down + ) + source_back = sim.add_source("WindowTurboSource", "source_back") + source_1 = sim.add_source("WindowTurboSource", "source_1") + source_2 = sim.add_source("WindowTurboSource", "source_2") + source_3 = sim.add_source("WindowTurboSource", "source_3") + source_1.activity = activity * Bq + source_2.activity = activity * Bq + source_3.activity = activity * Bq + source_back.activity = activity * Bq + source_back.direction.a1 = -radius_down * mm + source_back.direction.a2 = radius_down * mm + source_back.direction.b1 = -radius_down * mm + source_back.direction.b2 = radius_down * mm + source_back.direction.plane_distance = 86 * mm + source_back.direction.plane_phi = np.pi / 2 + source_back.direction.skip_mode = skip_mode + source_2.direction = source_back.direction.copy() + source_3.direction = source_back.direction.copy() + source_1.direction = source_back.direction.copy() + change_source_parameters(source_back, source_1, source_2, source_3) + if skip_mode: + source_back.direction.max_solid_angle = [0.09743848102142992] + source_1.direction.max_solid_angle = [0.021343358734360683] + source_2.direction.max_solid_angle = [0.025128504107295623] + source_3.direction.max_solid_angle = [0.023893613241897496] + sim.run(start_new_process=True) + + stats = sim.get_actor("Stats") + print(stats) + print("-" * 80) + + +def run_generic_source(activity=1000000): + Bq = gate.g4_units.Bq + + sim = initialize(10) + sim.number_of_threads = 32 + build_geometry(sim, paths.output / "generic") + + # physic list + # print('Phys lists :', sim.get_available_physicLists()) + + source_back = sim.add_source("GenericSource", "source_back") + source_1 = sim.add_source("GenericSource", "source_1") + source_2 = sim.add_source("GenericSource", "source_2") + source_3 = sim.add_source("GenericSource", "source_3") + source_1.activity = activity * Bq + source_2.activity = activity * Bq + source_3.activity = activity * Bq + source_back.activity = activity * Bq + change_source_parameters(source_back, source_1, source_2, source_3) + + sim.run() + + stats = sim.get_actor("Stats") + print(stats) + print("-" * 80) + + +if __name__ == "__main__": + pathFile = pathlib.Path(__file__).parent.resolve() + # run_generic_source() + run_window_turbo_source(skip_mode=False) + run_window_turbo_source(skip_mode=True) + profile_wt_true = calculate_profile(paths.output / "window_turbo_True_counts.mhd") + profile_wt_false = calculate_profile(paths.output / "window_turbo_False_counts.mhd") + profile_generic = calculate_profile(paths.output_ref / "generic.mhd") + compare_result = compare_profiles( + profile_generic, + [profile_wt_true, profile_wt_false], + tolerance=4.0, + fig_name=paths.output / "profile_comparison_skip.png", + ) + + utility.test_ok(compare_result) diff --git a/opengate/tests/src/source/test100_window_turbo_source_visu_wip.py b/opengate/tests/src/source/test100_window_turbo_source_visu_wip.py new file mode 100644 index 0000000000..bde664b42a --- /dev/null +++ b/opengate/tests/src/source/test100_window_turbo_source_visu_wip.py @@ -0,0 +1,179 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +from opengate.tests import utility +import pathlib +import numpy as np +from box import Box +import matplotlib.pyplot as plt +from opengate.tests import utility +from test100_window_turbo_source_base import build_geometry + +paths = utility.get_default_test_paths(__file__, output_folder="test100") + + +def change_source_parameters(source_back, source_1, source_2, source_3): + keV = gate.g4_units.keV + mm = gate.g4_units.mm + source_back.particle = "gamma" + source_back.energy.mono = 141 * keV + source_back.position.type = "cylinder" + source_back.position.translation = [0, -100, 0] + source_back.position.radius = 100 * mm + source_back.position.dz = 100 * mm + + source_1.particle = "gamma" + source_1.energy.mono = 141 * keV + source_1.position.type = "cylinder" + source_1.position.translation = [-50, -100, 0] + source_1.position.radius = 10 * mm + source_1.position.dz = 100 * mm + + source_2.particle = "gamma" + source_2.energy.mono = 141 * keV + source_2.position.type = "cylinder" + source_2.position.translation = [0, -100, 0] + source_2.position.radius = 15 * mm + source_2.position.dz = 100 * mm + + source_3.particle = "gamma" + source_3.energy.mono = 141 * keV + source_3.position.type = "cylinder" + source_3.position.translation = [50, -100, 0] + source_3.position.radius = 20 * mm + source_3.position.dz = 100 * mm + + +def initialize(duration=10): + sim = gate.Simulation() + sec = gate.g4_units.second + sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3" + sim.physics_manager.global_production_cuts.all = 1 * gate.g4_units.mm + # main options + sim.g4_verbose = True + sim.g4_verbose_level = 1 + sim.visu = True + sim.visu_type = "qt" + sim.number_of_threads = 1 + sim.progress_bar = True + sim.run_timing_intervals = [[0, duration * sec]] + sim.add_actor("SimulationStatisticsActor", "Stats") + return sim + + +def calculate_profile(image_path): + import SimpleITK as sitk + + image = sitk.ReadImage(image_path) + array = sitk.GetArrayFromImage(image) + profile = np.sum(array, axis=(0, 1)) + return profile + + +def compare_profiles(ref, test, tolerance=8.0, fig_name=None): + ref = np.asarray(ref, dtype=float) + test = np.asarray(test, dtype=float) + + if ref.shape != test.shape: + utility.print_test(False, f"Profile shapes differ: {ref.shape} vs {test.shape}") + return False + + sad = np.abs(ref - test).sum() / (ref.sum() + test.sum()) * 100 + is_ok = sad < tolerance + utility.print_test( + is_ok, f"Profile relative SAD = {sad:.2f}% (tol {tolerance:.2f}%)" + ) + + if fig_name is not None: + plt.figure(figsize=(10, 5)) + plt.plot(ref / ref.sum(), label="reference") + plt.plot(test / test.sum(), label="test") + plt.legend() + plt.xlabel("Pixel") + plt.ylabel("Normalized counts") + plt.savefig(fig_name) + + return is_ok + + +def run_window_turbo_source(activity=1): + Bq = gate.g4_units.Bq + mm = gate.g4_units.mm + sim = initialize(0) + sim.g4_verbose = True + sim.number_of_threads = 1 + sim.random_seed = 1 + radius_down = 13.6 + build_geometry(sim, paths.output / "window_turbo", pin_radius_down=radius_down) + source_back = sim.add_source("WindowTurboSource", "source_back") + source_1 = sim.add_source("WindowTurboSource", "source_1") + source_2 = sim.add_source("WindowTurboSource", "source_2") + source_3 = sim.add_source("WindowTurboSource", "source_3") + source_1.activity = activity * Bq + source_2.activity = activity * Bq + source_3.activity = activity * Bq + source_back.activity = activity * Bq + source_back.direction.a1 = -radius_down * mm + source_back.direction.a2 = radius_down * mm + source_back.direction.b1 = -radius_down * mm + source_back.direction.b2 = radius_down * mm + source_back.direction.plane_distance = 86 * mm + source_back.direction.plane_phi = np.pi / 2 + source_2.direction = source_back.direction.copy() + source_3.direction = source_back.direction.copy() + source_1.direction = source_back.direction.copy() + change_source_parameters(source_back, source_1, source_2, source_3) + source_back.visualization.window_run_id = 0 + source_back.visualization.window_width = 2.0 + source_back.visualization.window_color = "red" + source_back.visualization.color = "red" + source_1.visualization.color = "green" + source_2.visualization.color = "blue" + source_3.visualization.color = "cyan" + sim.run() + stats = sim.get_actor("Stats") + print(stats) + print("-" * 80) + + +def run_generic_source(activity=1000000): + Bq = gate.g4_units.Bq + + sim = initialize() + build_geometry(sim, "generic") + + # physic list + # print('Phys lists :', sim.get_available_physicLists()) + + source_back = sim.add_source("GenericSource", "source_back") + source_1 = sim.add_source("GenericSource", "source_1") + source_2 = sim.add_source("GenericSource", "source_2") + source_3 = sim.add_source("GenericSource", "source_3") + source_1.activity = activity * Bq + source_2.activity = activity * Bq + source_3.activity = activity * Bq + source_back.activity = activity * Bq + change_source_parameters(source_back, source_1, source_2, source_3) + + sim.run() + + stats = sim.get_actor("Stats") + print(stats) + print("-" * 80) + + +if __name__ == "__main__": + pathFile = pathlib.Path(__file__).parent.resolve() + # run_generic_source() + run_window_turbo_source() + # profile_wt = calculate_profile(paths.output / "window_turbo_counts.mhd") + # profile_generic = calculate_profile(paths.output_ref / "generic.mhd") + # compare_result = compare_profiles( + # profile_generic, + # profile_wt, + # tolerance=4.0, + # fig_name=paths.output / "profile_comparison.png", + # ) + + # utility.test_ok(compare_result) diff --git a/opengate/tests/src/source/test100_window_turbo_source_voxel_dynamic_wip.py b/opengate/tests/src/source/test100_window_turbo_source_voxel_dynamic_wip.py new file mode 100644 index 0000000000..11ad479d33 --- /dev/null +++ b/opengate/tests/src/source/test100_window_turbo_source_voxel_dynamic_wip.py @@ -0,0 +1,191 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import numpy as np +import SimpleITK as sitk + +import opengate as gate +from opengate.tests import utility +from test100_window_turbo_source_base import ( + build_geometry, + calculate_profile, + compare_profiles, +) + +paths = utility.get_default_test_paths(__file__, output_folder="test100") + + +def make_itk_source(output_path, zero_first_voxel=False): + data = np.ones((13, 2)) + if zero_first_voxel: + data[:, 0] = 0 + else: + data[:, 1] = 0 + data = data.flatten()[:-1].reshape((5, 1, 5)) + itk_image = sitk.GetImageFromArray(data) + itk_image.SetSpacing([50, 50, 50]) + itk_image.SetOrigin([-125 + 25, -125 + 25, -125 + 25]) + sitk.WriteImage(itk_image, output_path) + return output_path + + +def initialize(run_timing_intervals): + sim = gate.Simulation() + sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3" + sim.physics_manager.global_production_cuts.all = 1 * gate.g4_units.mm + sim.g4_verbose = False + sim.g4_verbose_level = 1 + sim.visu = False + sim.visu_type = "qt" + sim.number_of_threads = 4 + sim.progress_bar = False + sim.run_timing_intervals = run_timing_intervals + sim.add_actor("SimulationStatisticsActor", "Stats") + return sim + + +def change_source_parameters(source): + keV = gate.g4_units.keV + mm = gate.g4_units.mm + source.particle = "gamma" + source.energy.mono = 141 * keV + source.position.type = "cylinder" + source.position.translation = [0, -100, 0] + source.position.radius = 100 * mm + source.position.dz = 100 * mm + + +def configure_voxel_wt_source(sim, image_path, activity=1000000): + Bq = gate.g4_units.Bq + mm = gate.g4_units.mm + radius_down = 13.6 + head_y_pos = 100 + + voxel_wt_source = sim.add_source("VoxelWTSource", "source_back") + voxel_wt_source.image = str(image_path) + voxel_wt_source.activity = activity * Bq + voxel_wt_source.direction.a1 = -radius_down * mm + voxel_wt_source.direction.a2 = radius_down * mm + voxel_wt_source.direction.b1 = -radius_down * mm + voxel_wt_source.direction.b2 = radius_down * mm + voxel_wt_source.direction.plane_distance = (head_y_pos - 14) * mm + voxel_wt_source.direction.plane_phi = np.pi / 2 + change_source_parameters(voxel_wt_source) + return voxel_wt_source + + +def projection_exists(image_path): + if not image_path.exists(): + return False + try: + sitk.ReadImage(image_path) + except RuntimeError: + return False + return True + + +def run_static_voxel_wt_source(image_path, output_name): + output_path = paths.output / f"{output_name}_counts.mhd" + if projection_exists(output_path): + print(f"Reuse existing static projection: {output_path}") + return output_path + + sec = gate.g4_units.second + interval_duration = 1920 // 4 + sim = initialize([[0, interval_duration * sec]]) + build_geometry( + sim, + paths.output / output_name, + pin_radius_down=13.6, + head_y_pos=100, + ) + configure_voxel_wt_source(sim, image_path) + sim.run(start_new_process=True) + stats = sim.get_actor("Stats") + print(stats) + print("-" * 80) + return output_path + + +def run_dynamic_voxel_wt_source(image_paths, output_name): + sec = gate.g4_units.second + interval_duration = 1920 // 4 + sim = initialize( + [ + [0, interval_duration * sec], + [interval_duration * sec, 2 * interval_duration * sec], + ] + ) + build_geometry( + sim, + paths.output / output_name, + pin_radius_down=13.6, + head_y_pos=100, + ) + voxel_wt_source = configure_voxel_wt_source(sim, image_paths[0]) + voxel_wt_source.add_dynamic_parametrisation(image=image_paths) + sim.run(start_new_process=True) + stats = sim.get_actor("Stats") + print(stats) + print("-" * 80) + + expected_runs = sim.number_of_threads * len(sim.run_timing_intervals) + stats_ok = stats.counts.runs == expected_runs + utility.print_test( + stats_ok, + f"Stats runs count {stats.counts.runs} matches threads x timing intervals {expected_runs}", + ) + return paths.output / f"{output_name}_counts.mhd", stats_ok + + +def calculate_dynamic_profiles(image_path): + image = sitk.ReadImage(image_path) + array = sitk.GetArrayFromImage(image) + shape_ok = array.shape == (2, 100, 100) + utility.print_test( + shape_ok, + f"Dynamic projection shape is {array.shape}, expected (2, 100, 100)", + ) + if not shape_ok: + return [], False + return [np.sum(array[i], axis=0) for i in range(array.shape[0])], True + + +if __name__ == "__main__": + source_image_0 = make_itk_source(paths.output / "voxel_source_dynamic_0.mhd") + source_image_1 = make_itk_source( + paths.output / "voxel_source_dynamic_1.mhd", + zero_first_voxel=True, + ) + + static_projection_0 = run_static_voxel_wt_source( + source_image_0, + "voxel_wt_dynamic_static_0", + ) + static_projection_1 = run_static_voxel_wt_source( + source_image_1, + "voxel_wt_dynamic_static_1", + ) + dynamic_projection, stats_ok = run_dynamic_voxel_wt_source( + [source_image_0, source_image_1], + "voxel_wt_dynamic", + ) + + static_profile_0 = calculate_profile(static_projection_0) + static_profile_1 = calculate_profile(static_projection_1) + dynamic_profiles, shape_ok = calculate_dynamic_profiles(dynamic_projection) + + compare_result_0 = shape_ok and compare_profiles( + static_profile_0, + dynamic_profiles[0], + tolerance=4.0, + fig_name=paths.output / "profile_comparison_voxel_dynamic_0.png", + ) + compare_result_1 = shape_ok and compare_profiles( + static_profile_1, + dynamic_profiles[1], + tolerance=4.0, + fig_name=paths.output / "profile_comparison_voxel_dynamic_1.png", + ) + + utility.test_ok(stats_ok and shape_ok and compare_result_0 and compare_result_1) diff --git a/opengate/tests/src/source/test100_window_turbo_source_voxel_wip.py b/opengate/tests/src/source/test100_window_turbo_source_voxel_wip.py new file mode 100644 index 0000000000..8017da41ac --- /dev/null +++ b/opengate/tests/src/source/test100_window_turbo_source_voxel_wip.py @@ -0,0 +1,134 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +from opengate.tests import utility +import pathlib +import numpy as np +from box import Box +import matplotlib.pyplot as plt +import SimpleITK as sitk +from test100_window_turbo_source_base import ( + build_geometry, + compare_profiles, + calculate_profile, +) + +paths = utility.get_default_test_paths(__file__, output_folder="test100") + + +def change_source_parameters(source_back): + keV = gate.g4_units.keV + mm = gate.g4_units.mm + source_back.particle = "gamma" + source_back.energy.mono = 141 * keV + source_back.position.type = "cylinder" + source_back.position.translation = [0, -100, 0] + source_back.position.radius = 100 * mm + source_back.position.dz = 100 * mm + + +def make_itk_source(): + data = np.ones((13, 2)) + data[:, 1] = 0 + data = data.flatten()[:-1].reshape((5, 1, 5)) + itk_image = sitk.GetImageFromArray(data) + itk_image.SetSpacing([50, 50, 50]) + itk_image.SetOrigin([-125 + 25, -125 + 25, -125 + 25]) + sitk.WriteImage(itk_image, paths.output / "voxel_source.mhd") + # print(data) + + +def initialize(duration=10): + sim = gate.Simulation() + sec = gate.g4_units.second + sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3" + sim.physics_manager.global_production_cuts.all = 1 * gate.g4_units.mm + # main options + sim.g4_verbose = False + sim.g4_verbose_level = 1 + sim.visu = False + sim.visu_type = "qt" + sim.number_of_threads = 32 + sim.progress_bar = True + sim.run_timing_intervals = [[0, duration * sec]] + sim.add_actor("SimulationStatisticsActor", "Stats") + return sim + + +def run_window_turbo_source(activity=1000000): + Bq = gate.g4_units.Bq + mm = gate.g4_units.mm + NoT = 4 + sim = initialize(1920 // NoT) + sim.g4_verbose = False + + sim.progress_bar = False + sim.number_of_threads = NoT + # sim.random_seed = 1 + radius_down = 13.6 + head_y_pos = 100 + build_geometry( + sim, + paths.output / "voxel_wt", + pin_radius_down=radius_down, + head_y_pos=head_y_pos, + ) + voxel_wt_source = sim.add_source("VoxelWTSource", "source_back") + voxel_wt_source.image = str(paths.output / "voxel_source.mhd") + voxel_wt_source.activity = activity * Bq + voxel_wt_source.direction.a1 = -radius_down * mm + voxel_wt_source.direction.a2 = radius_down * mm + voxel_wt_source.direction.b1 = -radius_down * mm + voxel_wt_source.direction.b2 = radius_down * mm + voxel_wt_source.direction.plane_distance = (head_y_pos - 14) * mm + voxel_wt_source.direction.plane_phi = np.pi / 2 + # voxel_wt_source.visualize(2000,"red",5) + # voxel_wt_source.visualize_window("red",2) + change_source_parameters(voxel_wt_source) + sim.run() + stats = sim.get_actor("Stats") + print(stats) + print("-" * 80) + + +def run_generic_source(activity=1000000): + Bq = gate.g4_units.Bq + total_num = 1920 + thread_count = 10 + duration = 1920 // thread_count + sim = initialize(duration) + sim.number_of_threads = thread_count + build_geometry(sim, paths.output / f"voxel") + + # physic list + # print('Phys lists :', sim.get_available_physicLists()) + + voxel_source = sim.add_source("VoxelSource", "source_back") + voxel_source.image = str(paths.output / "voxel_source.mhd") + voxel_source.activity = activity * Bq + change_source_parameters(voxel_source) + # voxel_source.visualize(2000,"red",5) + + sim.run() + + stats = sim.get_actor("Stats") + print(stats) + print("-" * 80) + + +if __name__ == "__main__": + pathFile = pathlib.Path(__file__).parent.resolve() + make_itk_source() + # run_generic_source(1000000) + run_window_turbo_source() + profile_voxel_wt = calculate_profile(paths.output / "voxel_wt_counts.mhd") + profile_voxel = calculate_profile(paths.output_ref / "voxel.mhd") + compare_result = compare_profiles( + profile_voxel, + profile_voxel_wt, + tolerance=4.0, + fig_name=paths.output / "profile_comparison_voxel.png", + ) + + utility.test_ok(compare_result) diff --git a/opengate/tests/src/source/test100_window_turbo_source_wip.py b/opengate/tests/src/source/test100_window_turbo_source_wip.py new file mode 100644 index 0000000000..db22a22221 --- /dev/null +++ b/opengate/tests/src/source/test100_window_turbo_source_wip.py @@ -0,0 +1,142 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +from opengate.tests import utility +import pathlib +import numpy as np +from box import Box +import matplotlib.pyplot as plt +from opengate.tests import utility +from test100_window_turbo_source_base import ( + build_geometry, + compare_profiles, + calculate_profile, +) + +paths = utility.get_default_test_paths(__file__, output_folder="test100") + + +def change_source_parameters(source_back, source_1, source_2, source_3): + keV = gate.g4_units.keV + mm = gate.g4_units.mm + source_back.particle = "gamma" + source_back.energy.mono = 141 * keV + source_back.position.type = "cylinder" + source_back.position.translation = [0, -100, 0] + source_back.position.radius = 100 * mm + source_back.position.dz = 100 * mm + + source_1.particle = "gamma" + source_1.energy.mono = 141 * keV + source_1.position.type = "cylinder" + source_1.position.translation = [-50, -100, 0] + source_1.position.radius = 10 * mm + source_1.position.dz = 100 * mm + + source_2.particle = "gamma" + source_2.energy.mono = 141 * keV + source_2.position.type = "cylinder" + source_2.position.translation = [0, -100, 0] + source_2.position.radius = 15 * mm + source_2.position.dz = 100 * mm + + source_3.particle = "gamma" + source_3.energy.mono = 141 * keV + source_3.position.type = "cylinder" + source_3.position.translation = [50, -100, 0] + source_3.position.radius = 20 * mm + source_3.position.dz = 100 * mm + + +def initialize(duration=10): + sim = gate.Simulation() + sec = gate.g4_units.second + sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3" + sim.physics_manager.global_production_cuts.all = 1 * gate.g4_units.mm + # main options + sim.g4_verbose = False + sim.g4_verbose_level = 1 + sim.visu = False + sim.visu_type = "qt" + sim.number_of_threads = 32 + sim.progress_bar = True + sim.run_timing_intervals = [[0, duration * sec]] + sim.add_actor("SimulationStatisticsActor", "Stats") + return sim + + +def run_window_turbo_source(activity=1000000): + Bq = gate.g4_units.Bq + mm = gate.g4_units.mm + sim = initialize(80) + sim.g4_verbose = False + sim.number_of_threads = 4 + sim.random_seed = 1 + sim.progress_bar = False + radius_down = 13.6 + build_geometry(sim, paths.output / "window_turbo", pin_radius_down=radius_down) + source_back = sim.add_source("WindowTurboSource", "source_back") + source_1 = sim.add_source("WindowTurboSource", "source_1") + source_2 = sim.add_source("WindowTurboSource", "source_2") + source_3 = sim.add_source("WindowTurboSource", "source_3") + source_1.activity = activity * Bq + source_2.activity = activity * Bq + source_3.activity = activity * Bq + source_back.activity = activity * Bq + source_back.direction.a1 = -radius_down * mm + source_back.direction.a2 = radius_down * mm + source_back.direction.b1 = -radius_down * mm + source_back.direction.b2 = radius_down * mm + source_back.direction.plane_distance = 86 * mm + source_back.direction.plane_phi = np.pi / 2 + source_2.direction = source_back.direction.copy() + source_3.direction = source_back.direction.copy() + source_1.direction = source_back.direction.copy() + change_source_parameters(source_back, source_1, source_2, source_3) + sim.run() + stats = sim.get_actor("Stats") + print(stats) + print("-" * 80) + + +def run_generic_source(activity=1000000): + Bq = gate.g4_units.Bq + + sim = initialize() + build_geometry(sim, "generic") + + # physic list + # print('Phys lists :', sim.get_available_physicLists()) + + source_back = sim.add_source("GenericSource", "source_back") + source_1 = sim.add_source("GenericSource", "source_1") + source_2 = sim.add_source("GenericSource", "source_2") + source_3 = sim.add_source("GenericSource", "source_3") + source_1.activity = activity * Bq + source_2.activity = activity * Bq + source_3.activity = activity * Bq + source_back.activity = activity * Bq + change_source_parameters(source_back, source_1, source_2, source_3) + + sim.run() + + stats = sim.get_actor("Stats") + print(stats) + print("-" * 80) + + +if __name__ == "__main__": + pathFile = pathlib.Path(__file__).parent.resolve() + # run_generic_source() + run_window_turbo_source() + profile_wt = calculate_profile(paths.output / "window_turbo_counts.mhd") + profile_generic = calculate_profile(paths.output_ref / "generic.mhd") + compare_result = compare_profiles( + profile_generic, + profile_wt, + tolerance=4.0, + fig_name=paths.output / "profile_comparison.png", + ) + + utility.test_ok(compare_result)