From c70fe318a71d457e5ce1319b41829785b547b67f Mon Sep 17 00:00:00 2001 From: "Sung, Chi-Hsun" Date: Thu, 16 Apr 2026 11:35:56 +0200 Subject: [PATCH 1/4] code changes by mammo dicussed on 20260121 --- gecatsim/pyfiles/CommonTools.py | 123 ++++-- gecatsim/pyfiles/Detector_RayAngles_2D.py | 2 +- gecatsim/pyfiles/PhantomProjectorWrapper.py | 457 ++++++++++++++++++-- gecatsim/pyfiles/Phantom_Analytic.py | 10 +- gecatsim/pyfiles/Xray_Filter.py | 44 +- 5 files changed, 541 insertions(+), 95 deletions(-) diff --git a/gecatsim/pyfiles/CommonTools.py b/gecatsim/pyfiles/CommonTools.py index 3f74878..2b206e6 100644 --- a/gecatsim/pyfiles/CommonTools.py +++ b/gecatsim/pyfiles/CommonTools.py @@ -26,33 +26,19 @@ def make_col(a): return a -def feval(func_name, *args): - module_paths = [ - func_name, - f"gecatsim.pyfiles.{func_name}", - f"gecatsim.pyfiles.FlatPanel.{func_name}" - ] - - module = None - for path in module_paths: - try: - module = importlib.import_module(path) - break - except ModuleNotFoundError: - continue - - if module is None: - raise ImportError(f"Could not import module: {func_name}. Tried paths: {module_paths}") - - # Extract the function name (the part after the last dot) - func_name_only = func_name.split('.')[-1] - +def feval(funcName, *args): try: - func = getattr(module, func_name_only) - except AttributeError: - raise AttributeError(f"Function '{func_name_only}' not found in module '{module.__name__}'") - - return func(*args) + md = __import__(funcName) + except: + md = __import__("xcist.src.app.pyfiles."+funcName, fromlist=[funcName]) # equal to: from gecatsim.foo import foo + strip_leading_module = '.'.join(funcName.split('.')[1:]) + func_name_only = funcName.split('.')[-1] + + if len(strip_leading_module) > 0: + eval_name = f"md.{strip_leading_module}.{func_name_only}" + else: + eval_name = f"md.{func_name_only}" + return eval(eval_name)(*args) def load_C_lib(): @@ -74,16 +60,35 @@ def load_C_lib(): class emptyCFG: - pass + + def __str__(self): + attrs = ', '.join(f"{key}={value!r}" for key, value in self.__dict__.items()) + return f"{self.__class__.__name__}({attrs})" + + def __repr__(self): + return self.__str__() + + +import sys class PathHelper: def __init__(self): self._base_dir = os.getcwd() # Locate paths of lib and data. self.paths = {} - self.paths["main"] = os.path.dirname(os.path.abspath(__file__)) - self.paths["top"] = os.path.split(self.paths["main"])[0] + + # Handle PyInstaller bundled environment + if getattr(sys, 'frozen', False) and hasattr(sys, '_MEIPASS'): + # Running as compiled executable + bundle_dir = sys._MEIPASS + self.paths["main"] = os.path.join(bundle_dir, 'app', 'pyfiles') + self.paths["top"] = os.path.join(bundle_dir, 'app') + else: + # Running as normal Python script + self.paths["main"] = os.path.dirname(os.path.abspath(__file__)) + self.paths["top"] = os.path.split(self.paths["main"])[0] + self.paths["bowtie"] = os.path.join(self.paths["top"], 'bowtie') self.paths["cfg"] = os.path.join(self.paths["top"], 'cfg') self.paths["dose"] = os.path.join(self.paths["top"], 'dose') @@ -259,16 +264,11 @@ def source_cfg(*para): # initialize structs in cfg and structs attrList = ['sim', 'det', 'detNew', 'src', 'srcNew', 'spec', 'protocol', 'scanner', 'phantom', 'physics', 'recon', 'dose'] - for attr in attrList: - # Ensure cfg has the attribute if not hasattr(cfg, attr): setattr(cfg, attr, emptyCFG()) - - # Define the variable in the local scope if not already defined - if attr not in locals(): - globals()[attr] = emptyCFG() - # exec("%s = emptyCFG()" % attr) # globals()[attr] = ... is more explicit and safer for dynamic variable creation + if not attr in dir(): + exec("%s = emptyCFG()" % attr) # execute scripts in cfg file exec(open(cfg_file).read()) @@ -381,27 +381,60 @@ def overlap2d(oldimg, old_pos_x, old_pos_y, pos_x, pos_y): def rawread(fname, dataShape, dataType): # dataType is for numpy, ONLY allows: 'float'/'single', 'double', 'int'/'int32', 'uint'/'uint32', 'int8', 'int16' # they are single, double, int32, uin32, int8, int16 - with open(fname, 'rb') as fin: - data = fin.read() - # https://docs.python.org/3/library/struct.html - switcher = {'float': ['f', 4, np.single], - 'single': ['f', 4, np.single], - 'double': ['d', 8, np.double], + switcher = {'float': ['f', 4, np.float32], + 'float32': ['f', 4, np.float32], + 'single': ['f', 4, np.float32], + 'double': ['d', 8, np.float64], 'int': ['i', 4, np.int32], 'uint': ['I', 4, np.uint32], 'int32': ['i', 4, np.int32], 'uint32': ['I', 4, np.uint32], + 'uint8': ['B', 1, np.uint8], 'int8': ['b', 1, np.int8], 'int16': ['h', 2, np.int16]} - fmt = switcher[dataType] - data = struct.unpack("%d%s" % (len(data)/fmt[1], fmt[0]), data) - data = np.array(data, dtype=fmt[2]) + fmt_char, byte_size, np_dtype = switcher[dataType] + + # Use numpy.fromfile for direct binary reading - much faster + try: + data = np.fromfile(fname, dtype=np_dtype) + except: + # Fallback to manual reading if fromfile fails + with open(fname, 'rb') as fin: + raw_data = fin.read() + num_elements = len(raw_data) // byte_size + data = struct.unpack(f"{num_elements}{fmt_char}", raw_data) + data = np.array(data, dtype=np_dtype) + if dataShape: data = data.reshape(dataShape) return data + + # with open(fname, 'rb') as fin: + # data = fin.read() + + # # https://docs.python.org/3/library/struct.html + # switcher = {'float': ['f', 4, np.single], + # 'float32': ['f', 4, np.float32], + # 'single': ['f', 4, np.single], + # 'double': ['d', 8, np.double], + # 'int': ['i', 4, np.int32], + # 'uint': ['I', 4, np.uint32], + # 'int32': ['i', 4, np.int32], + # 'uint32': ['I', 4, np.uint32], + # 'uint8': ['B', 1, np.uint8], + # 'int8': ['b', 1, np.int8], + # 'int16': ['h', 2, np.int16]} + # fmt = switcher[dataType] + # data = struct.unpack("%d%s" % (len(data)/fmt[1], fmt[0]), data) + + # data = np.array(data, dtype=fmt[2]) + # if dataShape: + # data = data.reshape(dataShape) + + # return data def rawwrite(fname, data): with open(fname, 'wb') as fout: diff --git a/gecatsim/pyfiles/Detector_RayAngles_2D.py b/gecatsim/pyfiles/Detector_RayAngles_2D.py index 562024c..e3d15eb 100644 --- a/gecatsim/pyfiles/Detector_RayAngles_2D.py +++ b/gecatsim/pyfiles/Detector_RayAngles_2D.py @@ -30,7 +30,7 @@ def Detector_RayAngles_2D(cfg): xyzDet = nm.repmat(det.modCoords[m, :], nCells, 1) + \ det.cellCoords @ (np.c_[det.uvecs[m, :], det.vvecs[m, :]].T) - xyzR = xyzDet - nm.repmat(xyzSrc, nCells, 1); + xyzR = xyzDet - nm.repmat(xyzSrc, nCells, 1) cellInd = range(startInd, startInd+nCells) rayDistance[cellInd] = vectornorm(xyzR.T) diff --git a/gecatsim/pyfiles/PhantomProjectorWrapper.py b/gecatsim/pyfiles/PhantomProjectorWrapper.py index 1a912f8..088aeba 100644 --- a/gecatsim/pyfiles/PhantomProjectorWrapper.py +++ b/gecatsim/pyfiles/PhantomProjectorWrapper.py @@ -1,16 +1,143 @@ +""" +PhantomProjectorWrapper.py - Phantom and Projector Configuration Management + +============================================================================= +OVERVIEW - WHERE THIS FITS IN THE XCIST PIPELINE +============================================================================= + +This module is called during the phantom_scan phase of CatSim.run_all(): + + CatSim.run_all() + └── phantom_scan() + └── one_scan() + └── RunModels() with callback_map including: + - 'Phantom': PhantomWrapper() [initialization phase] + - 'Projector': ProjectorWrapper() [per-view projection] + +The module handles THREE phantom configuration scenarios: + +1. SINGLE PHANTOM (simplest case): + cfg.phantom.projectorCallback = "C_Projector_Voxelized" + → One phantom, same for all views + +2. MULTIPLE PHANTOMS PER VIEW (e.g., breast + lesion + calcifications): + cfg.phantom.projectorCallback = ["C_Projector_Voxelized", "Analytic_Projector", ...] + cfg.phantom.filename = ["breast.json", "lesion.json", ...] + → Multiple phantoms combined, same set for all views + +3. DYNAMIC PHANTOM (nested list - different phantom per view): + cfg.phantom.projectorCallback = [ + ["C_Projector_Voxelized"], # view 0 + ["C_Projector_Voxelized"], # view 1 + ["Analytic_Projector", "C_Proj..."], # view 2 (multiple phantoms) + ... + ] + → Each view can have different phantom configuration (e.g., motion simulation) + +============================================================================= +MEMORY CONSIDERATIONS +============================================================================= + +For multiple phantoms per view, partial_subview is allocated: + Shape: [nSamples, totalNumCells, nEbin] + Size: nSamples × totalNumCells × nEbin × 4 bytes + Example: 9 × 170,572,500 × 74 × 4 = ~450 GB for full mammo resolution! + +This is the LARGEST memory allocation in the entire pipeline and occurs +because each phantom's transmission must be computed independently before +combining (multiplicative combination for overlapping phantoms). + +For single phantom, no partial_subview is needed - projection writes +directly to cfg.thisSubView. + +============================================================================= +FUNCTION CALL SEQUENCE +============================================================================= + +Per-view execution flow: + + RunModels(cfg, viewId, subViewId) + │ + ├── [First view only] PhantomWrapper(cfg) + │ └── Handles dynamic phantom: extracts cfg.phantom for viewId + │ + └── ProjectorWrapper(cfg, viewId, subViewId) + │ + ├── [Dynamic phantom] Re-extract cfg.phantom for viewId + │ + ├── [Multiple phantoms] + │ ├── SplitCfgPhantom() → list of individual phantom cfgs + │ ├── Allocate partial_subview [nSamples, nCells, nEbin] + │ ├── Loop over each phantom: + │ │ ├── CopyCfgPhantom() → copy phantom config + │ │ ├── feval(phantom.callback) → Phantom_Voxelized() + │ │ └── feval(projectorCallback) → C_Projector_Voxelized() + │ └── Combine: thisSubView *= sum(partial_subview) + │ + └── [Single phantom] + ├── feval(phantom.callback) → Phantom_Voxelized() + └── feval(projectorCallback) → C_Projector_Voxelized() + +============================================================================= +""" + import sys from copy import copy, deepcopy from gecatsim.pyfiles.CommonTools import * def CopyCfgPhantom(cfgfrom, cfgto=None): + """ + Copy phantom configuration from one CFG object to another. + + Used during multiple-phantom iteration to update the main cfg with + each individual phantom's configuration while preserving all other + cfg attributes (detector, source, spectrum, etc.). + + Args: + cfgfrom: Source CFG object containing phantom to copy + cfgto: Destination CFG object (created if None) + + Returns: + cfgto with phantom attributes deep-copied from cfgfrom + + Note: Uses deepcopy to avoid reference issues when phantom attributes + are modified during projection. + """ if cfgto is None: cfgto = CFG() cfgto.phantom = deepcopy(cfgfrom.phantom) return cfgto -# split cfg.phantom from a list into + def SplitCfgPhantom(cfg): + """ + Split a multi-phantom configuration into a list of single-phantom CFGs. + + When cfg.phantom contains lists (multiple phantoms to combine), this + function creates separate CFG objects for each phantom so they can be + processed individually by their respective projectors. + + Example input: + cfg.phantom.filename = ["breast.raw", "lesion.raw", "calc.raw"] + cfg.phantom.projectorCallback = ["C_Projector", "C_Projector", "Analytic"] + cfg.phantom.materialList = [["adipose"], ["tumor"], ["calcium"]] + + Output: + cfg_list[0].phantom.filename = "breast.raw" + cfg_list[0].phantom.projectorCallback = "C_Projector" + cfg_list[0].phantom.materialList = ["adipose"] + ...and so on for each phantom + + Args: + cfg: Main CFG with phantom attributes as lists + + Returns: + List of CFG objects, each with single phantom configuration + + Note: Only phantom attributes are copied; other cfg attributes (det, src, etc.) + are accessed from the main cfg via CopyCfgPhantom during iteration. + """ cfg_list = [] cfgphantom_attrs = [x for x in dir(cfg.phantom) if not x.startswith('__')] for i in range(len(cfg.phantom.filename)): @@ -22,44 +149,318 @@ def SplitCfgPhantom(cfg): return cfg_list -def PhantomWrapper(cfg): - if isinstance(cfg.phantom.callback, list): - # if duplicate phantom callback, raise an error - if (len(cfg.phantom.callback) != len(set(cfg.phantom.callback))) and len(cfg.phantom.callback) > 1: - print("Error! XCIST does not support repeated phantom types.\n") - sys.exit(1) - cfg_list = SplitCfgPhantom(cfg) - origCfgPhantom = CopyCfgPhantom(cfg) - cfg_nom_list = [] - for thiscfg in cfg_list: - cfg = CopyCfgPhantom(thiscfg, cfg) - cfg = feval(cfg.phantom.callback, cfg) - if hasattr(cfg.phantom, 'numberOfMaterials'): - cfg_nom_list.append(cfg.phantom.numberOfMaterials) - delattr(cfg.phantom, 'numberOfMaterials') +def PhantomWrapper(cfg): + """ + Prepare phantom configuration for the current view (initialization phase). + + ========================================================================== + WHEN THIS IS CALLED + ========================================================================== + Called by RunModels() during the FIRST view/subview only (initialization), + via the 'Phantom' callback in callback_map. The recalc_map typically has: + 'Phantom': ['viewId', 'subViewId'] + meaning it only recalculates when view or subview changes. + + ========================================================================== + PURPOSE + ========================================================================== + Handles DYNAMIC PHANTOM configuration where each view may have different + phantom data (useful for motion simulation, contrast agent wash-in, etc.). + + Detection: Checks if projectorCallback is a NESTED LIST: + - Nested list (dynamic): [["proj_view0"], ["proj_view1"], ...] + - Simple list (multiple phantoms): ["proj1", "proj2", ...] + - String (single phantom): "C_Projector_Voxelized" + + ========================================================================== + WHAT IT DOES FOR DYNAMIC PHANTOMS + ========================================================================== + 1. Stores original phantom config in cfg._original_dynamic_phantom + (needed to restore full config for subsequent views) + + 2. Extracts phantom attributes at index [viewId] from all list attributes + Example: cfg.phantom.filename[viewId] → temp_cfg.phantom.filename + + 3. Copies all other cfg attributes to temp_cfg (detector, source, etc.) + + 4. Returns temp_cfg which now has single-view phantom configuration + + ========================================================================== + FOR NON-DYNAMIC PHANTOMS + ========================================================================== + Simply returns cfg unchanged - no processing needed. + + Args: + cfg: Main configuration object with cfg.viewId set + + Returns: + cfg (possibly modified temp_cfg for dynamic phantoms) + + Note: The actual phantom data loading (Phantom_Voxelized) and projection + happen later in ProjectorWrapper, not here. This only prepares the config. + """ + # ======================================================================== + # DYNAMIC PHANTOM DETECTION + # ======================================================================== + # Check if projectorCallback is a nested list (list of lists) + # This indicates different phantom configurations per view + # Example: [["C_Projector"], ["C_Projector", "Analytic"], ...] + # ↑ view 0 ↑ view 1 (multiple phantoms) + if (isinstance(cfg.phantom.projectorCallback, list) and + len(cfg.phantom.projectorCallback) > 0 and + isinstance(cfg.phantom.projectorCallback[0], list)): + + # ==================================================================== + # PRESERVE ORIGINAL CONFIG + # ==================================================================== + # Store the full dynamic phantom config so ProjectorWrapper can + # restore it for subsequent views (each view needs its own config) + cfg._original_dynamic_phantom = deepcopy(cfg.phantom) + + # ==================================================================== + # EXTRACT VIEW-SPECIFIC PHANTOM ATTRIBUTES + # ==================================================================== + # Get all phantom attributes (filename, projectorCallback, materialList, etc.) + cfgphantom_attrs = [x for x in dir(cfg.phantom) if not x.startswith('__')] + + # Create a temporary cfg to hold the extracted single-view config + temp_cfg = CFG() + + # For each phantom attribute, extract the value at index [viewId] + # if it's a list with enough elements, otherwise keep as-is + for thisattr in cfgphantom_attrs: + phantom_attr = getattr(cfg.phantom, thisattr) + if isinstance(phantom_attr, list) and len(phantom_attr) > cfg.viewId: + # Extract view-specific value: e.g., filename[viewId] + setattr(temp_cfg.phantom, thisattr, phantom_attr[cfg.viewId]) else: - # ncat does not have such property - cfg_nom_list.append(None) - - # restore original CFG - cfg = CopyCfgPhantom(origCfgPhantom, cfg) - cfg.phantom.numberOfMaterials = cfg_nom_list - else: - cfg = feval(cfg.phantom.callback, cfg) + # Keep non-list or short-list attributes unchanged + setattr(temp_cfg.phantom, thisattr, phantom_attr) + + # ==================================================================== + # COPY NON-PHANTOM CFG ATTRIBUTES + # ==================================================================== + # Preserve all other configuration (detector, source, spectrum, etc.) + for attr in dir(cfg): + if not attr.startswith('__') and attr != 'phantom': + setattr(temp_cfg, attr, getattr(cfg, attr)) + + # Return the view-specific configuration + cfg = temp_cfg + # For non-dynamic phantoms, return cfg unchanged return cfg def ProjectorWrapper(cfg, viewId, subViewId): + """ + Execute phantom projection for the current view and subview. + + ========================================================================== + WHEN THIS IS CALLED + ========================================================================== + Called by RunModels() for EVERY view/subview during phantom_scan, via the + 'Projector' callback in callback_map. The recalc_map typically has: + 'Projector': ['viewId', 'subViewId'] + + This is called AFTER all other RunModels callbacks (detector, source, + gantry, spectrum, filter, flux) have updated cfg for the current geometry. + + ========================================================================== + INPUT STATE + ========================================================================== + At entry, cfg contains: + - cfg.thisSubView: [totalNumCells, nEbin] float32 array + Already populated with detection flux (from Mammo_Detection_Flux) + This represents the X-ray intensity reaching each detector cell + before passing through the phantom. + + ========================================================================== + OUTPUT STATE + ========================================================================== + At exit, cfg.thisSubView has been MULTIPLIED by phantom transmission: + cfg.thisSubView *= transmission + + Where transmission = exp(-sum of path integrals through all materials) + + ========================================================================== + THREE MODES OF OPERATION + ========================================================================== + + MODE 1: DYNAMIC PHANTOM (nested list projectorCallback) + -------------------------------------------------------- + projectorCallback = [["proj_view0"], ["proj_view1"], ...] + + - Restores original phantom config from _original_dynamic_phantom + - Extracts view-specific phantom attributes at index [viewId] + - Then continues with Mode 2 or Mode 3 logic + + MODE 2: MULTIPLE PHANTOMS PER VIEW (simple list projectorCallback) + ------------------------------------------------------------------- + projectorCallback = ["C_Projector_Voxelized", "Analytic_Projector", ...] + + - Allocates partial_subview: [nSamples, totalNumCells, nEbin] + ⚠️ WARNING: This can be ~450 GB for full mammo resolution! + + - For CPU projector: pre-multiplies by source weights + - Iterates over each phantom: + 1. SplitCfgPhantom → get individual phantom configs + 2. CopyCfgPhantom → update cfg with this phantom + 3. feval(phantom.callback) → load phantom data (e.g., Phantom_Voxelized) + 4. feval(projectorCallback) → run projector (e.g., C_Projector_Voxelized) + The projector accumulates into partial_subview + + - Combines results: + CPU: thisSubView *= sum(partial_subview, axis=0) # sum over samples + GPU: thisSubView *= partial_subview[0,:,:] # GPU handles samples internally + + - Restores original phantom config for next view + + MODE 3: SINGLE PHANTOM (string projectorCallback) + ------------------------------------------------- + projectorCallback = "C_Projector_Voxelized" + + - Simplest case, no partial_subview allocation needed + - feval(phantom.callback) → load phantom data + - feval(projectorCallback) → run projector + - Projector directly modifies cfg.thisSubView + + ========================================================================== + PROJECTOR CALLBACK DETAILS (C_Projector_Voxelized) + ========================================================================== + The projector callback performs: + + 1. Double loop: for srcId in nSamples, for matId in nMaterials + 2. Calls C library: clib.voxelized_projector(phantom, rayOrigins, rayDirections) + 3. Computes path integrals through voxelized phantom + 4. Accumulates transmission: trans += weight * exp(-path_integral) + 5. Final: thisSubView *= trans (or partial_subview for multiple phantoms) + + Memory in projector (per call, float64): + - matPVS: [totalNumCells, nEbin] ~ 100 GB + - trans: [totalNumCells, nEbin] ~ 100 GB + - pValueSpectrum: [totalNumCells, nEbin] ~ 100 GB + These are freed when the projector function returns. + + ========================================================================== + MEMORY IMPACT + ========================================================================== + Mode 2 (multiple phantoms) is the most memory-intensive due to partial_subview. + Peak memory during multiple-phantom projection: + partial_subview (450 GB) + projector buffers (300 GB) = ~750 GB + + Mode 3 (single phantom) avoids partial_subview: + Only projector buffers (~300 GB) + + Args: + cfg: Main configuration with thisSubView already containing detection flux + viewId: Current view index (0 to protocol.viewCount-1) + subViewId: Current subview index (0 to subViewCount-1) + + Returns: + cfg with cfg.thisSubView multiplied by phantom transmission + """ + # ========================================================================== + # MODE 1: DYNAMIC PHANTOM HANDLING + # ========================================================================== + # Check if projectorCallback is a nested list (different config per view) + # This check runs EVERY view to extract the correct phantom for this viewId + if (isinstance(cfg.phantom.projectorCallback, list) and + len(cfg.phantom.projectorCallback) > 0 and + isinstance(cfg.phantom.projectorCallback[0], list)): + + # Restore the full dynamic phantom config (may have been modified by previous view) + if hasattr(cfg, '_original_dynamic_phantom'): + cfg.phantom = deepcopy(cfg._original_dynamic_phantom) + + # Extract phantom attributes at index [viewId] + cfgphantom_attrs = [x for x in dir(cfg.phantom) if not x.startswith('__')] + + temp_cfg = CFG() + + for thisattr in cfgphantom_attrs: + phantom_attr = getattr(cfg.phantom, thisattr) + if isinstance(phantom_attr, list) and len(phantom_attr) > viewId: + # Extract: filename[viewId], projectorCallback[viewId], etc. + setattr(temp_cfg.phantom, thisattr, phantom_attr[viewId]) + else: + setattr(temp_cfg.phantom, thisattr, phantom_attr) + + # Preserve non-phantom cfg attributes + for attr in dir(cfg): + if not attr.startswith('__') and attr != 'phantom': + setattr(temp_cfg, attr, getattr(cfg, attr)) + + # Continue with view-specific phantom config + cfg = temp_cfg + + # ========================================================================== + # MODE 2: MULTIPLE PHANTOMS PER VIEW + # ========================================================================== + # Check if projectorCallback is a simple list (not nested) = multiple phantoms if isinstance(cfg.phantom.projectorCallback, list): + # Split the multi-phantom config into individual phantom configs cfg_list = SplitCfgPhantom(cfg) - origCfgPhantom = CopyCfgPhantom(cfg) + + # ====================================================================== + # ALLOCATE PARTIAL_SUBVIEW - LARGEST MEMORY ALLOCATION IN PIPELINE + # ====================================================================== + # Shape: [nSamples, totalNumCells, nEbin] + # This stores transmission for each source sample separately because + # multiple phantoms may have different spatial relationships to each sample + # ⚠️ WARNING: ~450 GB for full mammo (9 samples × 170M cells × 74 bins × 4B) + cfg.partial_subview = np.ones([cfg.src.nSamples, cfg.det.totalNumCells, cfg.spec.nEbin], dtype=np.single) + + if not cfg.protocol.gpu_optimization: + # ================================================================== + # CPU PROJECTOR PATH + # ================================================================== + # Pre-multiply by source sample weights for proper flux weighting + # Shape broadcast: [nSamples,1,1] * [nSamples, nCells, nEbin] + cfg.partial_subview *= cfg.src.weights[:, np.newaxis, np.newaxis] + + # ====================================================================== + # ITERATE OVER EACH PHANTOM + # ====================================================================== for thiscfg in cfg_list: + # Copy this phantom's config into main cfg cfg = CopyCfgPhantom(thiscfg, cfg) + + # Load phantom data into C library (e.g., Phantom_Voxelized) + # This reads the .raw file and sets up material attenuation tables + cfg = feval(cfg.phantom.callback, cfg) + + # Execute projection (e.g., C_Projector_Voxelized) + # This multiplies partial_subview by this phantom's transmission cfg = feval(cfg.phantom.projectorCallback, cfg, viewId, subViewId) - # restore original CFG - cfg = CopyCfgPhantom(origCfgPhantom, cfg) + + # ====================================================================== + # COMBINE PARTIAL RESULTS + # ====================================================================== + # After all phantoms processed, combine into thisSubView + if not cfg.protocol.gpu_optimization: + # CPU: Sum over source samples (weighted sum already applied above) + # Physics: Total transmission = weighted sum of sample transmissions + cfg.thisSubView *= np.sum(cfg.partial_subview, axis=0) + else: + # GPU: GPU projector handles source sample weighting internally + # Only need the first "sample" which contains combined result + cfg.thisSubView *= cfg.partial_subview[0,:,:] + + # Restore original phantom config for next view (important for dynamic phantom) + cfg.phantom = deepcopy(cfg._original_dynamic_phantom) + + # ========================================================================== + # MODE 3: SINGLE PHANTOM + # ========================================================================== else: + # Simplest case: one phantom, one projector + # No partial_subview needed - projector writes directly to thisSubView + + # Load phantom data into C library + cfg = feval(cfg.phantom.callback, cfg) + + # Execute projection - directly modifies cfg.thisSubView cfg = feval(cfg.phantom.projectorCallback, cfg, viewId, subViewId) - return cfg + return cfg \ No newline at end of file diff --git a/gecatsim/pyfiles/Phantom_Analytic.py b/gecatsim/pyfiles/Phantom_Analytic.py index fc7aeb4..f5350b2 100644 --- a/gecatsim/pyfiles/Phantom_Analytic.py +++ b/gecatsim/pyfiles/Phantom_Analytic.py @@ -648,10 +648,12 @@ def Phantom_Analytic_SetObjects(objs=None, debug=False): phantObject[i]['s'] = objs['clip'][i][:,3].T else: #TODO: consider change to [] - #phantObject[i]['eta'] = [] - #phantObject[i]['s'] = [] - phantObject[i]['eta'] = None - phantObject[i]['s'] = None + # phantObject[i]['eta'] = [] + # phantObject[i]['s'] = [] + # phantObject[i]['eta'] = None + # phantObject[i]['s'] = None + phantObject[i]['eta'] = np.empty((3, 0)) + phantObject[i]['s'] = np.empty(0) p1=pars[6]*np.pi / 180 p2=pars[7]*np.pi / 180 p3=pars[8]*np.pi / 180 diff --git a/gecatsim/pyfiles/Xray_Filter.py b/gecatsim/pyfiles/Xray_Filter.py index dadf7b5..a50dbd3 100644 --- a/gecatsim/pyfiles/Xray_Filter.py +++ b/gecatsim/pyfiles/Xray_Filter.py @@ -20,49 +20,59 @@ def Xray_Filter(cfg): def flat_filter(cfg): ''' - Apply the transmittance of flat filter at source side, dim: [Ebin, totalNumCells] + Apply the transmittance of flat filter at source side, optimized version ''' - cosineFactors = 1/np.cos(cfg.det.gammas)/np.cos(cfg.det.alphas) + cosineFactors = 1 / (np.cos(cfg.det.gammas) * np.cos(cfg.det.alphas)) Evec = cfg.sim.Evec - trans = np.ones([cfg.det.totalNumCells, cfg.spec.nEbin], dtype = np.single) + trans = np.ones([cfg.det.totalNumCells, cfg.spec.nEbin], dtype=np.single) if hasattr(cfg.protocol, "flatFilter"): - for ii in range(0, round(len(cfg.protocol.flatFilter)/2)): + # Pre-allocate array for mu values + mu_total = np.zeros_like(Evec, dtype=np.single) + + for ii in range(0, len(cfg.protocol.flatFilter) // 2): material = cfg.protocol.flatFilter[2*ii] depth = cfg.protocol.flatFilter[2*ii+1] mu = GetMu(material, Evec) - trans *= np.exp(-depth*0.1*cosineFactors @ mu.reshape(1, mu.size)) - cfg.src.filterTrans *= trans + mu_total += depth * 0.1 * mu + + # Single exponential calculation instead of multiple + trans *= np.exp(-np.outer(cosineFactors, mu_total)) + cfg.src.filterTrans *= trans return cfg def bowtie_filter(cfg): ''' - Apply the transmittance of bowtie filter, dim: [Ebin, totalNumCells] + Apply the transmittance of bowtie filter, optimized version ''' if not cfg.protocol.bowtie: return cfg - # find bowtie file bowtieFile = my_path.find("bowtie", cfg.protocol.bowtie, ".txt") - - # read bowtie file data = np.loadtxt(bowtieFile, dtype=np.single, comments=['#', '%']) gammas0 = data[:, 0] - t0 = data[:, 1:] # thickness in cm + t0 = data[:, 1:] bowtieMaterials = ['Al', 'graphite', 'Cu', 'Ti'] Evec = cfg.spec.Evec gammas1 = cfg.det.gammas + cos_alphas_inv = 1.0 / np.cos(cfg.det.alphas) + + # Pre-compute all mu values + mu_array = np.array([GetMu(material, Evec) for material in bowtieMaterials], dtype=np.single) + + # Vectorized interpolation for all materials at once + f_interps = [interpolate.interp1d(gammas0, t0[:, i], kind='linear', fill_value='extrapolate', assume_sorted=True) + for i in range(len(bowtieMaterials))] - muT = 0 - for i in range(len(bowtieMaterials)): - mu = GetMu(bowtieMaterials[i], Evec) - f = interpolate.interp1d(gammas0, t0[:, i], kind='linear', fill_value='extrapolate') - t1 = f(gammas1)/np.cos(cfg.det.alphas) - muT += t1 @ mu.reshape(1, mu.size) + # Compute total attenuation in one operation + muT = np.zeros((len(gammas1), len(Evec)), dtype=np.single) + for i, f in enumerate(f_interps): + t1 = f(gammas1) * cos_alphas_inv + muT += np.outer(t1, mu_array[i]) trans = np.exp(-muT) cfg.src.filterTrans *= trans From 1cdeefc2b4a3b85f374f55f2d8f5bd277ca042ad Mon Sep 17 00:00:00 2001 From: "Sung, Chi-Hsun" Date: Thu, 16 Apr 2026 11:39:21 +0200 Subject: [PATCH 2/4] changes in c++ --- gecatsim/clib_build/src/analytic_projector.c | 276 +++++++++---------- gecatsim/clib_build/src/analytic_projector.h | 10 +- gecatsim/clib_build/src/phantom.tab.cpp | 4 +- gecatsim/clib_build/src/phantom.tab.hpp | 2 +- 4 files changed, 146 insertions(+), 146 deletions(-) diff --git a/gecatsim/clib_build/src/analytic_projector.c b/gecatsim/clib_build/src/analytic_projector.c index bb9b5ad..4ccd5d4 100644 --- a/gecatsim/clib_build/src/analytic_projector.c +++ b/gecatsim/clib_build/src/analytic_projector.c @@ -32,7 +32,7 @@ int n_col_oversample_add_xtalk = 1; int n_row_oversample_add_xtalk = 1; -struct module_info +struct analytic_module_info { double *Height; double *Width; @@ -48,7 +48,7 @@ struct module_info int moduleOverlapType; }; -struct module_info modules = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0,0,0,0}; +struct analytic_module_info analytic_modules = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0,0,0,0}; struct pnt { @@ -67,7 +67,7 @@ struct bounding_info struct bounding_info bounding = {NULL,NULL}; -struct phantom_info +struct analytic_phantom_info { int numObjects; int *objectType; @@ -83,16 +83,16 @@ struct phantom_info double xbounds[2]; }; -struct phantom_info phantom = {0,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,{0,0}}; +struct analytic_phantom_info analytic_phantom = {0,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,{0,0}}; -struct material_info +struct analytic_material_info { int materialCount; int eBinCount; double *muTable; }; -struct material_info materials = {0,0,NULL}; +struct analytic_material_info analytic_materials = {0,0,NULL}; struct projector_args { @@ -463,12 +463,12 @@ int quartic_intersect_C(double *a0,double *Qlp,double *Qrp,double shp,double *al int i,out; double scale,displ,a0t[3],tmp[3],aa,b,c,d,e,f,C[5],*Ql,*Qr,sh; - Ql = &phantom.Qmatrix[obj*18]; - // printf("%1.12lf %1.12lf",phantom.Qmatrix[obj); + Ql = &analytic_phantom.Qmatrix[obj*18]; + // printf("%1.12lf %1.12lf",analytic_phantom.Qmatrix[obj); // for(i = 1;i<9;i++) printf("L %1.12lf %1.12lf diff: %1.12lf\r\n",Ql[i],Qlp[i],Ql[i]-Qlp[i]); - Qr = &phantom.Qmatrix[obj*18+9]; + Qr = &analytic_phantom.Qmatrix[obj*18+9]; // for(i = 1;i<9;i++) printf("R %1.12lf %1.12lf diff: %1.12lf\r\n",Qr[i],Qrp[i],Qr[i]-Qrp[i]); - sh = phantom.shape[obj]; + sh = analytic_phantom.shape[obj]; // printf("s diff: %1.12lf\r\n",sh-shp); scale = 1.0; @@ -505,9 +505,9 @@ int quartic_intersect(double *a0,double *alpha,double *tc,int obj) int i,out; double scale,displ,a0t[3],tmp[3],aa,b,c,d,e,f,C[5],*Ql,*Qr,sh; - Ql = &phantom.Qmatrix[obj*18]; - Qr = &phantom.Qmatrix[obj*18+9]; - sh = phantom.shape[obj]; + Ql = &analytic_phantom.Qmatrix[obj*18]; + Qr = &analytic_phantom.Qmatrix[obj*18+9]; + sh = analytic_phantom.shape[obj]; scale = 1.0; displ = 0; for(i = 0;i<3;i++) displ += a0[i]*a0[i]; @@ -541,12 +541,12 @@ DLLEXPORT int clip_all(double *a,double *alpha,double rayLength,double *tc2,int double b[3], s1, s2, tcrit, tmin, tmax, *eta, *s, den; int j, cp, firstPlane, materialIndex; - firstPlane = phantom.clipStartIndex[i]; - eta = &phantom.clipNormalVector[firstPlane*3]; - s = &phantom.clipDistance[firstPlane]; - cp = phantom.nClipPlanes[i]; - den = phantom.density[i]; - materialIndex = phantom.materialInd[i]; + firstPlane = analytic_phantom.clipStartIndex[i]; + eta = &analytic_phantom.clipNormalVector[firstPlane*3]; + s = &analytic_phantom.clipDistance[firstPlane]; + cp = analytic_phantom.nClipPlanes[i]; + den = analytic_phantom.density[i]; + materialIndex = analytic_phantom.materialInd[i]; tmin = -VERY_BIG;tmax = VERY_BIG; for(j = 0;j<3;j++) b[j] = a[j]+rayLength*alpha[j]; @@ -604,8 +604,8 @@ DLLEXPORT int quadratic_intersect(double *a0,double *alpha,int pars11,double *tc int out; double A,B,C,tmp,*Q,k; // The quadratic formula is of the form A t^2 + B t + C = 0 - Q = &phantom.Qmatrix[obj*18]; - k = phantom.shape[obj]; + Q = &analytic_phantom.Qmatrix[obj*18]; + k = analytic_phantom.shape[obj]; C = quadratic_form(a0,Q,a0)-k; B = 2.0*quadratic_form(alpha,Q,a0); A = quadratic_form(alpha,Q,alpha); @@ -665,29 +665,29 @@ double* my_memcpyd(double *src, double *dest, int bytes) DLLEXPORT void set_src_info(double *sourceWeights, int nSubSources) { - modules.sourceWeights = my_memcpyd(sourceWeights, modules.sourceWeights, nSubSources*sizeof(double)); - modules.nSubSources = nSubSources; + analytic_modules.sourceWeights = my_memcpyd(sourceWeights, analytic_modules.sourceWeights, nSubSources*sizeof(double)); + analytic_modules.nSubSources = nSubSources; } DLLEXPORT void set_module_info(double *Height, double *Width, int *Pix, double *Coords, int *Sub, double *Sampling, double *Weight, int nModuleTypes, int maxPix, int maxSubDets, int moduleOverlapType) { int i; - modules.Height = my_memcpyd(Height, modules.Height, nModuleTypes*sizeof(double)); - modules.Width = my_memcpyd(Width, modules.Width, nModuleTypes*sizeof(double)); + analytic_modules.Height = my_memcpyd(Height, analytic_modules.Height, nModuleTypes*sizeof(double)); + analytic_modules.Width = my_memcpyd(Width, analytic_modules.Width, nModuleTypes*sizeof(double)); for(i = 0; i < nModuleTypes; i++) { - modules.Height[i] = MAX(1e-7,modules.Height[i]); - modules.Width[i] = MAX(1e-7,modules.Width[i]); + analytic_modules.Height[i] = MAX(1e-7,analytic_modules.Height[i]); + analytic_modules.Width[i] = MAX(1e-7,analytic_modules.Width[i]); } - modules.Pix = my_memcpyi(Pix, modules.Pix, nModuleTypes*sizeof(int)); - modules.Coords = my_memcpyd(Coords, modules.Coords, maxPix*2*nModuleTypes*sizeof(double)); - modules.Sub = my_memcpyi(Sub, modules.Sub, nModuleTypes*sizeof(int)); - modules.Sampling = my_memcpyd(Sampling, modules.Sampling, 2*maxSubDets*nModuleTypes*sizeof(double)); - modules.Weight = my_memcpyd(Weight, modules.Weight, maxSubDets*nModuleTypes*sizeof(double)); - modules.maxPixPerModule = maxPix; - modules.maxSubDets = maxSubDets; - modules.moduleOverlapType = moduleOverlapType; + analytic_modules.Pix = my_memcpyi(Pix, analytic_modules.Pix, nModuleTypes*sizeof(int)); + analytic_modules.Coords = my_memcpyd(Coords, analytic_modules.Coords, maxPix*2*nModuleTypes*sizeof(double)); + analytic_modules.Sub = my_memcpyi(Sub, analytic_modules.Sub, nModuleTypes*sizeof(int)); + analytic_modules.Sampling = my_memcpyd(Sampling, analytic_modules.Sampling, 2*maxSubDets*nModuleTypes*sizeof(double)); + analytic_modules.Weight = my_memcpyd(Weight, analytic_modules.Weight, maxSubDets*nModuleTypes*sizeof(double)); + analytic_modules.maxPixPerModule = maxPix; + analytic_modules.maxSubDets = maxSubDets; + analytic_modules.moduleOverlapType = moduleOverlapType; } DLLEXPORT void set_bounding_info(int numObjs, int *vertexStartInd, double *vertLocs, int numVerts) @@ -697,38 +697,38 @@ DLLEXPORT void set_bounding_info(int numObjs, int *vertexStartInd, double *vertL bounding.vertexStartIndex = my_memcpyi(vertexStartInd, bounding.vertexStartIndex, (numObjs+1)*sizeof(int)); bounding.vertexLocations = my_memcpyd(vertLocs, bounding.vertexLocations, sizeof(double)*numVerts*3); - phantom.xbounds[0]=VERY_BIG; - phantom.xbounds[1]=-VERY_BIG; + analytic_phantom.xbounds[0]=VERY_BIG; + analytic_phantom.xbounds[1]=-VERY_BIG; for(i=0;iphantom.xbounds[1]) phantom.xbounds[1]=bounding.vertexLocations[i*3]; - if (bounding.vertexLocations[i*3]analytic_phantom.xbounds[1]) analytic_phantom.xbounds[1]=bounding.vertexLocations[i*3]; + if (bounding.vertexLocations[i*3]maxVertexPoints) maxVertexPoints = vertexPoints[i]; @@ -1695,7 +1695,7 @@ void compute_object_projections(double *srcHullPoints, int nSrcHullPoints, int m indexList = malloc(sizeof(int)*(maxVertexPoints*nSrcHullPoints+moduleEdges*2)); //dbug(3,"COP1\r\n"); - for(i = 0;i0) dbug(2,"listLength: %d\r\n",listLength); @@ -1788,7 +1788,7 @@ int any_objects_1(int moduleNumber, void *boundaries) int any_objs = 0, i; struct height_lims *heightLims = (struct height_lims *) boundaries; - for(i = 0;i < phantom.numObjects;i++) + for(i = 0;i < analytic_phantom.numObjects;i++) { //dbug(3,"Module number: %d minHeight[%d] = %lf maxHeight[%d] = %lf Different? %d",moduleNumber,i,heightLims.min[i],i,heightLims.max[i],(int)(heightLims.min[i] != heightLims.max[i])); if (heightLims[i].min != heightLims[i].max) @@ -1803,7 +1803,7 @@ int any_objects_2(int moduleNumber, void *boundaries) int any_objs = 0, i; struct box_lims *boxLims = (struct box_lims *) boundaries; - for(i = 0;i < phantom.numObjects;i++) + for(i = 0;i < analytic_phantom.numObjects;i++) { if (boxLims[i].minR != boxLims[i].maxR) {any_objs++;break;} @@ -1817,7 +1817,7 @@ void build_object_list1(double *pix_vlims, int *objectList, int *n_objlist, int int i; struct height_lims *heightLims = (struct height_lims *) boundaries; - for(i = 0;i= v+pix_vlims[0])) {objectList[n_objlist[0]] = i;n_objlist[0]++;} } @@ -1828,7 +1828,7 @@ void build_object_list2(double *pix_ulims, double *pix_vlims, int *objectList, i int i; struct box_lims *boxLims = (struct box_lims *) boundaries; - for(i = 0;i= uv[1]+pix_vlims[0]),(int) (boxLims[i].minR <= uv[0]+pix_ulims[1]),(int) (boxLims[i].maxR >= uv[0]+pix_ulims[0])); dbug(3,"parts for third one: %1.5lf %1.5lf %1.5lf %1.5lf \r\n ",boxLims[0].minR,uv[0],pix_ulims[1],uv[0]+pix_ulims[1]); @@ -1853,13 +1853,13 @@ void intersections_loop(double *detCenter, double *right, double *up, double *sa double *materialBuffer = NULL; double *pValueSpectrum = NULL; - materialBuffer = malloc(sizeof(double)*materials.materialCount); - pValueSpectrum = malloc(sizeof(double)*materials.eBinCount); + materialBuffer = malloc(sizeof(double)*analytic_materials.materialCount); + pValueSpectrum = malloc(sizeof(double)*analytic_materials.eBinCount); for(l = 0;l0.0) dbug(1,"materialBuffer[0]: %4.4lf\r\n",materialBuffer[0]); - for(m = 0;m0.0) dbug(1,"pValueSpectrum[0]: %4.4lf\r\n",pValueSpectrum[0]); //------------------- Accurate Detector Model, Mingye if(Accurate_Detector_Model_is_ON) { - for(m = 0;m < materials.eBinCount;m++) + for(m = 0;m < analytic_materials.eBinCount;m++) { - int the_index = (detIndex*n_col_oversample + (int)(l/n_row_oversample_add_xtalk))*materials.eBinCount + m; - thisView[the_index] += subviewWeight*detWeights[l]*modules.sourceWeights[k]*exp(-pValueSpectrum[m]); + int the_index = (detIndex*n_col_oversample + (int)(l/n_row_oversample_add_xtalk))*analytic_materials.eBinCount + m; + thisView[the_index] += subviewWeight*detWeights[l]*analytic_modules.sourceWeights[k]*exp(-pValueSpectrum[m]); } } //------------------- else { - for(m = 0;m < materials.eBinCount;m++) - thisView[detIndex*materials.eBinCount+m] += subviewWeight*detWeights[l]*modules.sourceWeights[k]*exp(-pValueSpectrum[m]); + for(m = 0;m < analytic_materials.eBinCount;m++) + thisView[detIndex*analytic_materials.eBinCount+m] += subviewWeight*detWeights[l]*analytic_modules.sourceWeights[k]*exp(-pValueSpectrum[m]); } } } @@ -1915,24 +1915,24 @@ void set_Accurate_Detector_Model(double *Paras) } if(Paras[5]==1 && Paras[1]<3 && Paras[2]==1) //col_xtalk>0, first subview in the first two views - if(modules.Sub[0]>n_col_oversample*n_row_oversample_add_xtalk) + if(analytic_modules.Sub[0]>n_col_oversample*n_row_oversample_add_xtalk) { int n1; - modules.maxSubDets-=2*n_row_oversample_add_xtalk; - modules.Sub[0]-=2*n_row_oversample_add_xtalk; + analytic_modules.maxSubDets-=2*n_row_oversample_add_xtalk; + analytic_modules.Sub[0]-=2*n_row_oversample_add_xtalk; double sum_weights=0; - for(n1=0;n1pix_vlims[1]) pix_vlims[1] = sampling[2*i+1]; if(sampling[2*i] phantom.numObjects) - //dbug(2,"nListObjects = %d, phantom.numObjects = %d\r\n", nListObjects, phantom.numObjects); + //if (nListObjects > analytic_phantom.numObjects) + //dbug(2,"nListObjects = %d, analytic_phantom.numObjects = %d\r\n", nListObjects, analytic_phantom.numObjects); //if (nListObjects > 0) //dbug(2,"nListObjects = %d, k = %d\r\n", nListObjects, k); break; @@ -2078,21 +2078,21 @@ DLLEXPORT void Projector(double *Paras, double subviewWeight, double *thisView, //------------------- Accurate Detector Model, Mingye if(Accurate_Detector_Model_is_ON) { - for(m = 0;m0) { min_val = VERY_BIG; - for(k = 0;k Date: Thu, 16 Apr 2026 12:01:57 +0200 Subject: [PATCH 3/4] path fix --- gecatsim/pyfiles/CommonTools.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gecatsim/pyfiles/CommonTools.py b/gecatsim/pyfiles/CommonTools.py index 2b206e6..fe1d933 100644 --- a/gecatsim/pyfiles/CommonTools.py +++ b/gecatsim/pyfiles/CommonTools.py @@ -30,7 +30,7 @@ def feval(funcName, *args): try: md = __import__(funcName) except: - md = __import__("xcist.src.app.pyfiles."+funcName, fromlist=[funcName]) # equal to: from gecatsim.foo import foo + md = __import__("gecatsim.src.app.pyfiles."+funcName, fromlist=[funcName]) # equal to: from gecatsim.foo import foo strip_leading_module = '.'.join(funcName.split('.')[1:]) func_name_only = funcName.split('.')[-1] From 85417e3bbb99917ad3bc231cfd0e665abef7108d Mon Sep 17 00:00:00 2001 From: Jiayong Zhang <100655819+zhangjy-ge@users.noreply.github.com> Date: Tue, 5 May 2026 16:50:42 -0400 Subject: [PATCH 4/4] support dynamic phantom with recalcPht with additional changes to be compatible with master branch --- gecatsim/pyfiles/C_Projector_SetData.py | 2 + gecatsim/pyfiles/C_Projector_Voxelized.py | 3 +- gecatsim/pyfiles/CommonTools.py | 43 ++- gecatsim/pyfiles/PhantomProjectorWrapper.py | 362 ++++---------------- 4 files changed, 94 insertions(+), 316 deletions(-) diff --git a/gecatsim/pyfiles/C_Projector_SetData.py b/gecatsim/pyfiles/C_Projector_SetData.py index bc8e77f..e5fa8c5 100644 --- a/gecatsim/pyfiles/C_Projector_SetData.py +++ b/gecatsim/pyfiles/C_Projector_SetData.py @@ -31,6 +31,8 @@ def C_Projector_SetData(cfg, viewId, subViewId): def get_projector_id(cfg): if not isinstance(cfg.phantom.callback, list): phantom_callbacks = [cfg.phantom.callback] + elif isinstance(cfg.phantom.callback, list) and isinstance(cfg.phantom.callback[0], list): + phantom_callbacks = cfg.phantom.callback[cfg.viewId] else: phantom_callbacks = cfg.phantom.callback projectorIDs = [] diff --git a/gecatsim/pyfiles/C_Projector_Voxelized.py b/gecatsim/pyfiles/C_Projector_Voxelized.py index 13ea98b..46f55ed 100644 --- a/gecatsim/pyfiles/C_Projector_Voxelized.py +++ b/gecatsim/pyfiles/C_Projector_Voxelized.py @@ -51,7 +51,8 @@ def C_Projector_Voxelized(cfg, viewId, subViewId): unused7 = 7 freeTheMemory = 0 - if viewId==cfg.sim.stopViewId and subViewId==cfg.sim.subViewCount-1 \ + # ZJY recalcPht is helpful in reducing memory usage in dynamic phantom + if (cfg.physics.recalcPht or (viewId==cfg.sim.stopViewId and subViewId==cfg.sim.subViewCount-1))\ and srcId==src.nSamples-1 and matId==cfg.phantom.numberOfMaterials-1: freeTheMemory = 1 diff --git a/gecatsim/pyfiles/CommonTools.py b/gecatsim/pyfiles/CommonTools.py index fe1d933..59b6a9d 100644 --- a/gecatsim/pyfiles/CommonTools.py +++ b/gecatsim/pyfiles/CommonTools.py @@ -27,18 +27,33 @@ def make_col(a): def feval(funcName, *args): + module_paths = [ + funcName, + f"gecatsim.pyfiles.{funcName}", + f"gecatsim.pyfiles.FlatPanel.{funcName}", + f"gecatsim.src.app.pyfiles.{funcName}" + ] + + module = None + for path in module_paths: + try: + module = importlib.import_module(path) + break + except ModuleNotFoundError: + continue + + if module is None: + raise ImportError(f"Could not import module: {funcName}. Tried paths: {module_paths}") + + # Extract the function name (the part after the last dot) + funcName_only = funcName.split('.')[-1] + try: - md = __import__(funcName) - except: - md = __import__("gecatsim.src.app.pyfiles."+funcName, fromlist=[funcName]) # equal to: from gecatsim.foo import foo - strip_leading_module = '.'.join(funcName.split('.')[1:]) - func_name_only = funcName.split('.')[-1] + func = getattr(module, funcName_only) + except AttributeError: + raise AttributeError(f"Function '{funcName_only}' not found in module '{module.__name__}'") - if len(strip_leading_module) > 0: - eval_name = f"md.{strip_leading_module}.{func_name_only}" - else: - eval_name = f"md.{func_name_only}" - return eval(eval_name)(*args) + return func(*args) def load_C_lib(): @@ -265,10 +280,14 @@ def source_cfg(*para): # initialize structs in cfg and structs attrList = ['sim', 'det', 'detNew', 'src', 'srcNew', 'spec', 'protocol', 'scanner', 'phantom', 'physics', 'recon', 'dose'] for attr in attrList: + # Ensure cfg has the attribute if not hasattr(cfg, attr): setattr(cfg, attr, emptyCFG()) - if not attr in dir(): - exec("%s = emptyCFG()" % attr) + + # Define the variable in the local scope if not already defined + if attr not in locals(): + globals()[attr] = emptyCFG() + # exec("%s = emptyCFG()" % attr) # globals()[attr] = ... is more explicit and safer for dynamic variable creation # execute scripts in cfg file exec(open(cfg_file).read()) diff --git a/gecatsim/pyfiles/PhantomProjectorWrapper.py b/gecatsim/pyfiles/PhantomProjectorWrapper.py index 088aeba..ed519fb 100644 --- a/gecatsim/pyfiles/PhantomProjectorWrapper.py +++ b/gecatsim/pyfiles/PhantomProjectorWrapper.py @@ -1,3 +1,5 @@ +gecatsim\pyfiles\PhantomProjectorWrapper.py +@@ -1,466 +1,94 @@ """ PhantomProjectorWrapper.py - Phantom and Projector Configuration Management @@ -109,7 +111,6 @@ def CopyCfgPhantom(cfgfrom, cfgto=None): return cfgto - def SplitCfgPhantom(cfg): """ Split a multi-phantom configuration into a list of single-phantom CFGs. @@ -143,324 +144,79 @@ def SplitCfgPhantom(cfg): for i in range(len(cfg.phantom.filename)): thiscfg = CFG() for thisattr in cfgphantom_attrs: - setattr(thiscfg.phantom, thisattr, deepcopy(getattr(cfg.phantom, thisattr)[i])) + # broadcast in this case + if not isinstance(getattr(cfg.phantom, thisattr), list) or len(getattr(cfg.phantom, thisattr))==1: + setattr(thiscfg.phantom, thisattr, deepcopy(getattr(cfg.phantom, thisattr))) + else: + setattr(thiscfg.phantom, thisattr, deepcopy(getattr(cfg.phantom, thisattr)[i])) cfg_list.append(thiscfg) return cfg_list - def PhantomWrapper(cfg): - """ - Prepare phantom configuration for the current view (initialization phase). - - ========================================================================== - WHEN THIS IS CALLED - ========================================================================== - Called by RunModels() during the FIRST view/subview only (initialization), - via the 'Phantom' callback in callback_map. The recalc_map typically has: - 'Phantom': ['viewId', 'subViewId'] - meaning it only recalculates when view or subview changes. - - ========================================================================== - PURPOSE - ========================================================================== - Handles DYNAMIC PHANTOM configuration where each view may have different - phantom data (useful for motion simulation, contrast agent wash-in, etc.). - - Detection: Checks if projectorCallback is a NESTED LIST: - - Nested list (dynamic): [["proj_view0"], ["proj_view1"], ...] - - Simple list (multiple phantoms): ["proj1", "proj2", ...] - - String (single phantom): "C_Projector_Voxelized" - - ========================================================================== - WHAT IT DOES FOR DYNAMIC PHANTOMS - ========================================================================== - 1. Stores original phantom config in cfg._original_dynamic_phantom - (needed to restore full config for subsequent views) - - 2. Extracts phantom attributes at index [viewId] from all list attributes - Example: cfg.phantom.filename[viewId] → temp_cfg.phantom.filename - - 3. Copies all other cfg attributes to temp_cfg (detector, source, etc.) - - 4. Returns temp_cfg which now has single-view phantom configuration - - ========================================================================== - FOR NON-DYNAMIC PHANTOMS - ========================================================================== - Simply returns cfg unchanged - no processing needed. - - Args: - cfg: Main configuration object with cfg.viewId set - - Returns: - cfg (possibly modified temp_cfg for dynamic phantoms) - - Note: The actual phantom data loading (Phantom_Voxelized) and projection - happen later in ProjectorWrapper, not here. This only prepares the config. - """ - # ======================================================================== - # DYNAMIC PHANTOM DETECTION - # ======================================================================== - # Check if projectorCallback is a nested list (list of lists) - # This indicates different phantom configurations per view - # Example: [["C_Projector"], ["C_Projector", "Analytic"], ...] - # ↑ view 0 ↑ view 1 (multiple phantoms) - if (isinstance(cfg.phantom.projectorCallback, list) and - len(cfg.phantom.projectorCallback) > 0 and - isinstance(cfg.phantom.projectorCallback[0], list)): - - # ==================================================================== - # PRESERVE ORIGINAL CONFIG - # ==================================================================== - # Store the full dynamic phantom config so ProjectorWrapper can - # restore it for subsequent views (each view needs its own config) - cfg._original_dynamic_phantom = deepcopy(cfg.phantom) - - # ==================================================================== - # EXTRACT VIEW-SPECIFIC PHANTOM ATTRIBUTES - # ==================================================================== - # Get all phantom attributes (filename, projectorCallback, materialList, etc.) - cfgphantom_attrs = [x for x in dir(cfg.phantom) if not x.startswith('__')] - - # Create a temporary cfg to hold the extracted single-view config - temp_cfg = CFG() - - # For each phantom attribute, extract the value at index [viewId] - # if it's a list with enough elements, otherwise keep as-is - for thisattr in cfgphantom_attrs: - phantom_attr = getattr(cfg.phantom, thisattr) - if isinstance(phantom_attr, list) and len(phantom_attr) > cfg.viewId: - # Extract view-specific value: e.g., filename[viewId] - setattr(temp_cfg.phantom, thisattr, phantom_attr[cfg.viewId]) + # nested loop, like [[pro1, proj2], [proj1], [proj2], ...] + is_dynamic_phantom = False + if isinstance(cfg.phantom.callback, list) and isinstance(cfg.phantom.callback[0], list): + is_dynamic_phantom = True + _cfg_list = SplitCfgPhantom(cfg) + _origCfgPhantom = CopyCfgPhantom(cfg) + cfg = CopyCfgPhantom(_cfg_list[cfg.viewId], cfg) + + if isinstance(cfg.phantom.callback, list): + # if duplicate phantom callback, raise an error + if (len(cfg.phantom.callback) != len(set(cfg.phantom.callback))) and len(cfg.phantom.callback) > 1: + print("Error! XCIST does not support repeated phantom types.\n") + sys.exit(1) + + cfg_list = SplitCfgPhantom(cfg) + origCfgPhantom = CopyCfgPhantom(cfg) + cfg_nom_list = [] + for thiscfg in cfg_list: + # use cfg here because we want to use cfg's other attr than phantom without the heavy copy overhead + cfg = CopyCfgPhantom(thiscfg, cfg) + cfg = feval(cfg.phantom.callback, cfg) + if hasattr(cfg.phantom, 'numberOfMaterials'): + cfg_nom_list.append(cfg.phantom.numberOfMaterials) + delattr(cfg.phantom, 'numberOfMaterials') else: - # Keep non-list or short-list attributes unchanged - setattr(temp_cfg.phantom, thisattr, phantom_attr) - - # ==================================================================== - # COPY NON-PHANTOM CFG ATTRIBUTES - # ==================================================================== - # Preserve all other configuration (detector, source, spectrum, etc.) - for attr in dir(cfg): - if not attr.startswith('__') and attr != 'phantom': - setattr(temp_cfg, attr, getattr(cfg, attr)) - - # Return the view-specific configuration - cfg = temp_cfg + # ncat does not have such property + cfg_nom_list.append(None) + + # restore original CFG + if is_dynamic_phantom: + cfg = CopyCfgPhantom(_origCfgPhantom, cfg) + cfg.phantom.numberOfMaterials = [None]*len(cfg.phantom.callback) + cfg.phantom.numberOfMaterials[cfg.viewId] = cfg_nom_list + else: + cfg = CopyCfgPhantom(origCfgPhantom, cfg) + cfg.phantom.numberOfMaterials = cfg_nom_list + else: + cfg = feval(cfg.phantom.callback, cfg) + - # For non-dynamic phantoms, return cfg unchanged return cfg def ProjectorWrapper(cfg, viewId, subViewId): - """ - Execute phantom projection for the current view and subview. - - ========================================================================== - WHEN THIS IS CALLED - ========================================================================== - Called by RunModels() for EVERY view/subview during phantom_scan, via the - 'Projector' callback in callback_map. The recalc_map typically has: - 'Projector': ['viewId', 'subViewId'] - - This is called AFTER all other RunModels callbacks (detector, source, - gantry, spectrum, filter, flux) have updated cfg for the current geometry. - - ========================================================================== - INPUT STATE - ========================================================================== - At entry, cfg contains: - - cfg.thisSubView: [totalNumCells, nEbin] float32 array - Already populated with detection flux (from Mammo_Detection_Flux) - This represents the X-ray intensity reaching each detector cell - before passing through the phantom. - - ========================================================================== - OUTPUT STATE - ========================================================================== - At exit, cfg.thisSubView has been MULTIPLIED by phantom transmission: - cfg.thisSubView *= transmission - - Where transmission = exp(-sum of path integrals through all materials) - - ========================================================================== - THREE MODES OF OPERATION - ========================================================================== - - MODE 1: DYNAMIC PHANTOM (nested list projectorCallback) - -------------------------------------------------------- - projectorCallback = [["proj_view0"], ["proj_view1"], ...] - - - Restores original phantom config from _original_dynamic_phantom - - Extracts view-specific phantom attributes at index [viewId] - - Then continues with Mode 2 or Mode 3 logic - - MODE 2: MULTIPLE PHANTOMS PER VIEW (simple list projectorCallback) - ------------------------------------------------------------------- - projectorCallback = ["C_Projector_Voxelized", "Analytic_Projector", ...] - - - Allocates partial_subview: [nSamples, totalNumCells, nEbin] - ⚠️ WARNING: This can be ~450 GB for full mammo resolution! - - - For CPU projector: pre-multiplies by source weights - - Iterates over each phantom: - 1. SplitCfgPhantom → get individual phantom configs - 2. CopyCfgPhantom → update cfg with this phantom - 3. feval(phantom.callback) → load phantom data (e.g., Phantom_Voxelized) - 4. feval(projectorCallback) → run projector (e.g., C_Projector_Voxelized) - The projector accumulates into partial_subview - - - Combines results: - CPU: thisSubView *= sum(partial_subview, axis=0) # sum over samples - GPU: thisSubView *= partial_subview[0,:,:] # GPU handles samples internally - - - Restores original phantom config for next view - - MODE 3: SINGLE PHANTOM (string projectorCallback) - ------------------------------------------------- - projectorCallback = "C_Projector_Voxelized" - - - Simplest case, no partial_subview allocation needed - - feval(phantom.callback) → load phantom data - - feval(projectorCallback) → run projector - - Projector directly modifies cfg.thisSubView - - ========================================================================== - PROJECTOR CALLBACK DETAILS (C_Projector_Voxelized) - ========================================================================== - The projector callback performs: - - 1. Double loop: for srcId in nSamples, for matId in nMaterials - 2. Calls C library: clib.voxelized_projector(phantom, rayOrigins, rayDirections) - 3. Computes path integrals through voxelized phantom - 4. Accumulates transmission: trans += weight * exp(-path_integral) - 5. Final: thisSubView *= trans (or partial_subview for multiple phantoms) - - Memory in projector (per call, float64): - - matPVS: [totalNumCells, nEbin] ~ 100 GB - - trans: [totalNumCells, nEbin] ~ 100 GB - - pValueSpectrum: [totalNumCells, nEbin] ~ 100 GB - These are freed when the projector function returns. - - ========================================================================== - MEMORY IMPACT - ========================================================================== - Mode 2 (multiple phantoms) is the most memory-intensive due to partial_subview. - Peak memory during multiple-phantom projection: - partial_subview (450 GB) + projector buffers (300 GB) = ~750 GB - - Mode 3 (single phantom) avoids partial_subview: - Only projector buffers (~300 GB) - - Args: - cfg: Main configuration with thisSubView already containing detection flux - viewId: Current view index (0 to protocol.viewCount-1) - subViewId: Current subview index (0 to subViewCount-1) - - Returns: - cfg with cfg.thisSubView multiplied by phantom transmission - """ - # ========================================================================== - # MODE 1: DYNAMIC PHANTOM HANDLING - # ========================================================================== - # Check if projectorCallback is a nested list (different config per view) - # This check runs EVERY view to extract the correct phantom for this viewId - if (isinstance(cfg.phantom.projectorCallback, list) and - len(cfg.phantom.projectorCallback) > 0 and - isinstance(cfg.phantom.projectorCallback[0], list)): - - # Restore the full dynamic phantom config (may have been modified by previous view) - if hasattr(cfg, '_original_dynamic_phantom'): - cfg.phantom = deepcopy(cfg._original_dynamic_phantom) - - # Extract phantom attributes at index [viewId] - cfgphantom_attrs = [x for x in dir(cfg.phantom) if not x.startswith('__')] - - temp_cfg = CFG() - - for thisattr in cfgphantom_attrs: - phantom_attr = getattr(cfg.phantom, thisattr) - if isinstance(phantom_attr, list) and len(phantom_attr) > viewId: - # Extract: filename[viewId], projectorCallback[viewId], etc. - setattr(temp_cfg.phantom, thisattr, phantom_attr[viewId]) - else: - setattr(temp_cfg.phantom, thisattr, phantom_attr) - - # Preserve non-phantom cfg attributes - for attr in dir(cfg): - if not attr.startswith('__') and attr != 'phantom': - setattr(temp_cfg, attr, getattr(cfg, attr)) - - # Continue with view-specific phantom config - cfg = temp_cfg - - # ========================================================================== - # MODE 2: MULTIPLE PHANTOMS PER VIEW - # ========================================================================== - # Check if projectorCallback is a simple list (not nested) = multiple phantoms + is_dynamic_phantom = False + if isinstance(cfg.phantom.callback, list) and isinstance(cfg.phantom.callback[0], list): + is_dynamic_phantom = True + _cfg_list = SplitCfgPhantom(cfg) + _origCfgPhantom = CopyCfgPhantom(cfg) + cfg = CopyCfgPhantom(_cfg_list[cfg.viewId], cfg) + if isinstance(cfg.phantom.projectorCallback, list): - # Split the multi-phantom config into individual phantom configs cfg_list = SplitCfgPhantom(cfg) - - # ====================================================================== - # ALLOCATE PARTIAL_SUBVIEW - LARGEST MEMORY ALLOCATION IN PIPELINE - # ====================================================================== - # Shape: [nSamples, totalNumCells, nEbin] - # This stores transmission for each source sample separately because - # multiple phantoms may have different spatial relationships to each sample - # ⚠️ WARNING: ~450 GB for full mammo (9 samples × 170M cells × 74 bins × 4B) - cfg.partial_subview = np.ones([cfg.src.nSamples, cfg.det.totalNumCells, cfg.spec.nEbin], dtype=np.single) - - if not cfg.protocol.gpu_optimization: - # ================================================================== - # CPU PROJECTOR PATH - # ================================================================== - # Pre-multiply by source sample weights for proper flux weighting - # Shape broadcast: [nSamples,1,1] * [nSamples, nCells, nEbin] - cfg.partial_subview *= cfg.src.weights[:, np.newaxis, np.newaxis] - - # ====================================================================== - # ITERATE OVER EACH PHANTOM - # ====================================================================== + origCfgPhantom = CopyCfgPhantom(cfg) for thiscfg in cfg_list: - # Copy this phantom's config into main cfg cfg = CopyCfgPhantom(thiscfg, cfg) - - # Load phantom data into C library (e.g., Phantom_Voxelized) - # This reads the .raw file and sets up material attenuation tables - cfg = feval(cfg.phantom.callback, cfg) - - # Execute projection (e.g., C_Projector_Voxelized) - # This multiplies partial_subview by this phantom's transmission cfg = feval(cfg.phantom.projectorCallback, cfg, viewId, subViewId) - - # ====================================================================== - # COMBINE PARTIAL RESULTS - # ====================================================================== - # After all phantoms processed, combine into thisSubView - if not cfg.protocol.gpu_optimization: - # CPU: Sum over source samples (weighted sum already applied above) - # Physics: Total transmission = weighted sum of sample transmissions - cfg.thisSubView *= np.sum(cfg.partial_subview, axis=0) - else: - # GPU: GPU projector handles source sample weighting internally - # Only need the first "sample" which contains combined result - cfg.thisSubView *= cfg.partial_subview[0,:,:] - - # Restore original phantom config for next view (important for dynamic phantom) - cfg.phantom = deepcopy(cfg._original_dynamic_phantom) - - # ========================================================================== - # MODE 3: SINGLE PHANTOM - # ========================================================================== + # restore original CFG + cfg = CopyCfgPhantom(origCfgPhantom, cfg) else: - # Simplest case: one phantom, one projector - # No partial_subview needed - projector writes directly to thisSubView - - # Load phantom data into C library - cfg = feval(cfg.phantom.callback, cfg) - - # Execute projection - directly modifies cfg.thisSubView cfg = feval(cfg.phantom.projectorCallback, cfg, viewId, subViewId) - return cfg \ No newline at end of file + if is_dynamic_phantom: + cfg = CopyCfgPhantom(_origCfgPhantom, cfg) + + return cfg