OpenGATE / opengate

Gate 10 (beta)
http://www.opengatecollaboration.org
GNU Lesser General Public License v3.0
42 stars 38 forks source link

Projection Exception #339

Closed jizhang02 closed 7 months ago

jizhang02 commented 7 months ago

Hello,

I'm using the Gate beta07 version to simulate a pet image. Everything works fine except for adding a DigitizerProjectionActor. The output is like this: image

Does anyone have a clue about it? Thanks a lot.

Here is my code:

import opengate as gate
import opengate.contrib.pet.siemensbiograph as pet_biograph # pet system
from opengate.userhooks import check_production_cuts
import numpy as np

'''
General
'''
sim = gate.Simulation()
sim.number_of_threads = 1
ui = sim.user_info
ui.g4_verbose = False
ui.g4_verbose_level = 1
ui.visu = False
ui.visu_type = "vrml_file_only"
ui.visu_filename = "VisuFile.wrl"
ui.visu_verbose = False
sim.random_engine = "MersenneTwister"
sim.random_seed = 123456
m = gate.g4_units.m
mm = gate.g4_units.mm
Bq = gate.g4_units.Bq
kBq = Bq * 1e3
sec = gate.g4_units.second
keV = gate.g4_units.keV
MeV = gate.g4_units.MeV

'''
Geometry
'''
# world
sim.world.size = [2 * m, 2 * m, 2 * m]
sim.world.material = "G4_AIR"
sim.volume_manager.add_material_database("MC/GateMaterials_pet.db") # add specified materials
# add a PET Biograph system
pet = pet_biograph.add_pet(sim, "pet")

# CT image
ct = sim.add_volume('Image', 'ct')
ct.mother = 'world'
ct.image = 'data/ct_cropped.mhd'
ct_info = gate.image.read_image_info(ct.image)
ct.voxel_materials = [
    [0,   4,  "G4_AIR"],
    [5,   6,  "G4_EYE_LENS_ICRP"],
    [7,   7,  "G4_SKIN_ICRP"],
    [8,  43,  "G4_BRAIN_ICRP"],
    [44, 51,  "G4_BLOOD_ICRP"],
    [52, 54,  "G4_LUNG_ICRP"],
    [55, 61,  "G4_MUSCLE_SKELETAL_ICRP"],
    [62, 63,  "G4_B-100_BONE"],
    [64, 66,  "G4_BONE_COMPACT_ICRU"],
    [67, 111, "G4_TISSUE_SOFT_ICRP"]
]

'''
Physics
'''
sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option4"
sim.physics_manager.global_production_cuts.all = 1 * m
sim.physics_manager.set_production_cut("world", "all", 10 * mm)
sim.physics_manager.set_production_cut(f"{pet.name}_crystal", "all", 0.1 * mm)

'''
Sources
'''
source = sim.add_source("VoxelsSource", "source")
source.mother = ct.name
source.image = 'data/activity_ga_lesion_norm_cropped.mhd'
total_yield = gate.sources.generic.get_rad_yield("Ga68")
source.activity = 28.118 * kBq * total_yield / sim.number_of_threads
source.half_life = 4062.6 * sec
source.energy.type = "Ga68"

'''
Actors
'''

# get crystal volume
crystal = sim.volume_manager.volumes[f"{pet.name}_crystal"]
block = sim.volume_manager.volumes[f"{pet.name}_block"]
ring = sim.volume_manager.volumes[f"{pet.name}_ring"]

# hits collection
hc = sim.add_actor("DigitizerHitsCollectionActor", "Hits")
hc.mother = crystal.name
hc.output = f"pet.root"
hc.attributes = ["PostPosition", "TotalEnergyDeposit", "PreStepUniqueVolumeID", "GlobalTime",]

# singles collection
sc = sim.add_actor("DigitizerReadoutActor", "Singles")
sc.output = hc.output
sc.input_digi_collection = hc.name
sc.group_volume = block.name
# sc.group_volume = ring.name # (I checked that is it different from block)
# sc.group_volume = crystal.name # (I checked that is it different from block)
sc.discretize_volume = crystal.name
sc.policy = "EnergyWeightedCentroidPosition"

# EnergyWindows
ew = sim.add_actor("DigitizerEnergyWindowsActor", "EnergyWindows")
ew.mother = crystal.name
ew.input_digi_collection = "Singles"
ew.channels = [
    {"name": "scatter", "min": 114 * keV, "max": 126 * keV},
    {"name": "peak140", "min": 126 * keV, "max": 154.55 * keV},
    # {'name': 'spectrum', 'min': 0 * keV, 'max': 5000 * keV}  # should be strictly equal to 'Singles'
]
ew.output = hc.output

# projection-not working
proj = sim.add_actor("DigitizerProjectionActor", "Projection")
proj.mother = hc.mother
proj.input_digi_collections = ["Singles", "scatter", "peak140"]
proj.spacing = [ct_info.spacing[0] * mm, ct_info.spacing[2] * mm]
proj.size = [ct_info.size[0], ct_info.size[2]]
proj.origin_as_image_center = False
# proj.plane = 'XY' # not implemented yet
proj.output = "projection.mhd"

# add stat actor
stats = sim.add_actor("SimulationStatisticsActor", "Stats")
stats.track_types_flag = True

'''
Simulation
'''
sim.run_timing_intervals = [[0, 1 * sec]]
sim.user_hook_after_init = check_production_cuts # set user hook to dump production cuts from G4
sim.run()
stats = sim.output.get_actor("Stats")
print(stats)
nkrah commented 7 months ago

I need to check in more detail, but I think your DigitizerProjectionActor is not hooked up correctly. Shouldn't its mother be the DigitizerEnergyWindowsActor?

proj.mother = ew.name

You are setting the mother to hc.mother, which points to the crystal volume.

Again, I have not tried to reproduce your error yet, just glanced over your simulation script.

jizhang02 commented 7 months ago

I need to check in more detail, but I think your DigitizerProjectionActor is not hooked up correctly. Shouldn't its mother be the DigitizerEnergyWindowsActor?

proj.mother = ew.name

You are setting the mother to hc.mother, which points to the crystal volume.

Again, I have not tried to reproduce your error yet, just glanced over your simulation script.

Hi, I tried to change the mother, but it still didn't work.

Exception: Cannot find volume EnergyWindows. Volumes included in this simulation are: dict_keys(['world', 'pet', 'pet_ring', 'pet_block', 'pet_crystal', 'ct'])
dsarrut commented 7 months ago

I think all "mother" must be "crystal.name" here. However, the error says that pet_block is not a volume, while it is defined in the add_pet, so something is weird.

jizhang02 commented 7 months ago

Hello,

I commented out:

#proj.mother = crystal.name

Then it works, thank you!