diff --git a/pytheranostics/data/output.json b/pytheranostics/data/output.json index db61af9..ee54f14 100644 --- a/pytheranostics/data/output.json +++ b/pytheranostics/data/output.json @@ -4,8 +4,8 @@ "ClinicalTrial": "MyClinicalTrial", "Radionuclide": "NA", "PatientID": "NA", - "Gender": "NA", - + "Gender": "NA", + "No_of_completed_cycles": "NA", "Cycle_01": [ { @@ -16,13 +16,14 @@ "InjectionTime": "NA", "InjectedActivity": "NA", "Weight_g": "NA", + "Height_cm": "NA", "Level": "NA", "Method": "NA", "OutputFormat": "NA", "ScaleDoseByDensity": "NA", "ReferenceTimePoint": "NA", "TimePoints_h": "NA", - "rois": {}, + "VOIs": {}, "Organ-level_AD": {} } ], @@ -35,13 +36,14 @@ "InjectionTime": "NA", "InjectedActivity": "NA", "Weight_g": "NA", + "Height_cm": "NA", "Level": "NA", "Method": "NA", "OutputFormat": "NA", "ScaleDoseByDensity": "NA", "ReferenceTimePoint": "NA", "TimePoints_h": "NA", - "rois": {}, + "VOIs": {}, "Organ-level_AD": {} } ], @@ -54,13 +56,14 @@ "InjectionTime": "NA", "InjectedActivity": "NA", "Weight_g": "NA", + "Height_cm": "NA", "Level": "NA", "Method": "NA", "OutputFormat": "NA", "ScaleDoseByDensity": "NA", "ReferenceTimePoint": "NA", "TimePoints_h": "NA", - "rois": {}, + "VOIs": {}, "Organ-level_AD": {} } ], @@ -73,13 +76,14 @@ "InjectionTime": "NA", "InjectedActivity": "NA", "Weight_g": "NA", + "Height_cm": "NA", "Level": "NA", "Method": "NA", "OutputFormat": "NA", "ScaleDoseByDensity": "NA", "ReferenceTimePoint": "NA", "TimePoints_h": "NA", - "rois": {}, + "VOIs": {}, "Organ-level_AD": {} } ], @@ -92,6 +96,7 @@ "InjectionTime": "NA", "InjectedActivity": "NA", "Weight_g": "NA", + "Height_cm": "NA", "ApplyBiokineticsFromPreviousCycle": "NA", "Level": "NA", "Method": "NA", @@ -99,7 +104,7 @@ "ScaleDoseByDensity": "NA", "ReferenceTimePoint": "NA", "TimePoints_h": "NA", - "rois": {}, + "VOIs": {}, "Organ-level_AD": {} } ], @@ -112,6 +117,7 @@ "InjectionTime": "NA", "InjectedActivity": "NA", "Weight_g": "NA", + "Height_cm": "NA", "ApplyBiokineticsFromPreviousCycle": "NA", "Level": "NA", "Method": "NA", @@ -119,7 +125,7 @@ "ScaleDoseByDensity": "NA", "ReferenceTimePoint": "NA", "TimePoints_h": "NA", - "rois": {}, + "VOIs": {}, "Organ-level_AD": {} } ] diff --git a/pytheranostics/data/phantom/human/ICRP_mass_male_source.csv b/pytheranostics/data/phantom/human/ICRP_mass_male_source.csv new file mode 100644 index 0000000..3c71886 --- /dev/null +++ b/pytheranostics/data/phantom/human/ICRP_mass_male_source.csv @@ -0,0 +1,28 @@ +Organ,Mass_g +Adrenals,14 +Brain,1450 +Esophagus,40 +Eyes,15 +Gallbladder Contents,58 +LLI Cont,75 +Small Intestine,350 +Stomach Contents,250 +ULI Cont,150 +Rectum,75 +Heart Contents,510 +Heart Wall,330 +Kidneys,310 +Liver,1800 +Lungs,1200 +Pancreas,140 +Prostate,17 +Salivary Glands,85 +Red Marrow,1170 +Cortical Bone,4400 +Trabecular Bone,1100 +Spleen,150 +Testes,35 +Thymus,25 +Thyroid,20 +Urinary Bladder Contents,211 +Total Body,73000 diff --git a/pytheranostics/data/phantom/human/ICRP_mass_male_target.csv b/pytheranostics/data/phantom/human/ICRP_mass_male_target.csv new file mode 100644 index 0000000..498bc57 --- /dev/null +++ b/pytheranostics/data/phantom/human/ICRP_mass_male_target.csv @@ -0,0 +1,27 @@ +Organ,Mass_g +Adrenals,14 +Brain,1450 +Esophagus,40 +Eyes,15 +Gallbladder Wall,10 +LLI Wall,150 +Small Intestine,650 +Stomach Wall,150 +ULI Wall,150 +Rectum,70 +Heart Wall,330 +Kidneys,310 +Liver,1800 +Lungs,1200 +Pancreas,140 +Prostate,17 +Salivary Glands,85 +Red Marrow,1170 +Cortical Bone,60 +Trabecular Bone,60 +Spleen,150 +Testes,35 +Thymus,25 +Thyroid,20 +Urinary Bladder Wall,50 +Total Body,73000 diff --git a/pytheranostics/dosimetry/base_dosimetry.py b/pytheranostics/dosimetry/base_dosimetry.py index f599f86..8af3cd5 100644 --- a/pytheranostics/dosimetry/base_dosimetry.py +++ b/pytheranostics/dosimetry/base_dosimetry.py @@ -123,9 +123,9 @@ def __init__( # DataFrame storing results self.results = self.initialize() - self.results_lesions = pandas.DataFrame() - self.results_salivaryglands = pandas.DataFrame() - self.df_ad = pandas.DataFrame() + self.results_dosimetry_lesions = pandas.DataFrame() + self.results_dosimetry_salivaryglands = pandas.DataFrame() + self.results_dosimetry_organs = pandas.DataFrame() # Sanity Checks: self.sanity_checks(metric="Volume_CT_mL") @@ -151,12 +151,13 @@ def default_config(self) -> None: "with_uptake": False, "fit_order": 1, "bounds": None, + "washout_ratio": None, } for key, value in defaults.items(): for region, _ in self.results.iterrows(): - if key not in self.config["rois"][region]: - self.config["rois"][region][key] = value + if key not in self.config["VOIs"][region]: + self.config["VOIs"][region][key] = value print( f"For {region}, the parameter '{key}' was not defined by the user, set to {value}." ) @@ -165,7 +166,7 @@ def extract_masks_and_correct_overlaps(self) -> None: """Extract masks and correct overlaps between regions.""" # Inform the user if some masks are unused and therefore excluded. for roi_name in self.nm_data.masks[0]: - if roi_name not in self.config["rois"] and roi_name != "BoneMarrow": + if roi_name not in self.config["VOIs"] and roi_name != "BoneMarrow": print( f"Although mask for {roi_name} is present, we are ignoring it because this region was not included in the" " configuration input file.\n" @@ -176,7 +177,7 @@ def extract_masks_and_correct_overlaps(self) -> None: time_id: extract_masks( time_id=time_id, mask_dataset=self.nm_data.masks, - requested_rois=list(self.config["rois"].keys()), + requested_rois=list(self.config["VOIs"].keys()), ) for time_id in self.nm_data.masks.keys() } @@ -185,13 +186,13 @@ def extract_masks_and_correct_overlaps(self) -> None: time_id: extract_masks( time_id=time_id, mask_dataset=self.ct_data.masks, - requested_rois=list(self.config["rois"].keys()), + requested_rois=list(self.config["VOIs"].keys()), ) for time_id in self.ct_data.masks.keys() } # Check availability of requested rois in existing masks - for roi_name in self.config["rois"]: + for roi_name in self.config["VOIs"]: if roi_name not in self.nm_data.masks[0] and roi_name != "BoneMarrow": raise AssertionError(f"The following mask was NOT found: {roi_name}\n") @@ -231,11 +232,20 @@ def check_mandatory_fields(self) -> None: self.config["ReferenceTimePoint"] = 0 if "Organ" in self.config["Level"]: - if "WholeBody" not in self.config["rois"]: - raise ValueError("Missing 'WholeBody' region parameters.") - - if "RemainderOfBody" not in self.config["rois"]: - raise ValueError("Missing 'RemainderOfBody' region parameters.") + if "WholeBody" not in self.config["VOIs"]: + if "No" in self.config["OrganLevel"]["AdditionalOptions"]["WholeBody"]: + pass + else: + raise ValueError("Missing 'WholeBody' region parameters.") + + if "RemainderOfBody" not in self.config["VOIs"]: + if ( + "No" + in self.config["OrganLevel"]["AdditionalOptions"]["RemainderOfBody"] + ): + pass + else: + raise ValueError("Missing 'RemainderOfBody' region parameters.") return None @@ -244,7 +254,7 @@ def initialize(self) -> pandas.DataFrame: tmp_results: Dict[str, List[float]] = { roi_name: [] for roi_name in self.nm_data.masks[0].keys() - if roi_name in self.config["rois"] + if roi_name in self.config["VOIs"] } cols: List[str] = ["Time_hr", "Volume_CT_mL", "Activity_MBq", "Density_HU"] @@ -294,7 +304,7 @@ def initialize_bone_marrow( ) -> Dict[str, List[float]]: """Initialize activity and times for Bone-Marrow blood-based measurements.""" if ( - "BoneMarrow" in self.config["rois"] + "BoneMarrow" in self.config["VOIs"] and self.clinical_data is not None and "BoneMarrow" not in self.nm_data.masks[0] ): @@ -387,8 +397,6 @@ def sanity_checks(self, metric: str) -> None: def compute_tia(self) -> None: """Compute Time-Integrated Activity over each source-organ.""" - # decay_constant = math.log(2) / (self.radionuclide["half_life"]) # 1/h # TODO: Check how to incorporate into bounds? (flake8) - if self.radionuclide["half_life_units"] != "hours": raise AssertionError( "Radionuclide Half-Life in Database should be in hours." @@ -403,7 +411,6 @@ def compute_tia(self) -> None: } for region, region_data in self.results.iterrows(): - fit_results = self.smart_fit_selection( region_data=region_data, region=region ) @@ -443,7 +450,7 @@ def compute_tia(self) -> None: tmp_tia_data["Lambda_eff"].append( [ fit_params[exp_params[i]] - for i in range(self.config["rois"][region]["fit_order"]) + for i in range(self.config["VOIs"][region]["fit_order"]) ] ) @@ -460,33 +467,19 @@ def compute_tia(self) -> None: def smart_fit_selection( self, region_data: pandas.Series, region: str ) -> lmfit.model.ModelResult: - """Select the best fit based on Akaike Information Criterion. - - If enabled in self.config, iterates through different valid fits orders and select best fit based on Akaike Information Criterion. - If a fit order is specified by the user, then the method will just perform fit following user's selected order and configuration. - - Parameters - ---------- - region_data : pandas.Series - Series containing Time and Activity. - region : str - Region of Interest - - Returns - ------- - lmfit.model.ModelResult - The best fit model based on Akaike Information Criterion. - """ + """Select the best fit based on Akaike Information Criterion.""" # If fit_order is defined by user: - if self.config["rois"][region]["fit_order"] is not None: + if self.config["VOIs"][region]["fit_order"] is not None: + print(region) fit_results, _ = exponential_fit_lmfit( x_data=numpy.array(region_data["Time_hr"]), y_data=numpy.array(region_data["Activity_MBq"]), - fixed_params=self.config["rois"][region]["fixed_parameters"], - num_exponentials=self.config["rois"][region]["fit_order"], - bounds=self.config["rois"][region]["bounds"], - params_init=self.config["rois"][region]["param_init"], - with_uptake=self.config["rois"][region]["with_uptake"], + fixed_params=self.config["VOIs"][region]["fixed_parameters"], + num_exponentials=self.config["VOIs"][region]["fit_order"], + bounds=self.config["VOIs"][region]["bounds"], + params_init=self.config["VOIs"][region]["param_init"], + with_uptake=self.config["VOIs"][region]["with_uptake"], + washout_ratio=self.config["VOIs"][region]["washout_ratio"], ) return fit_results @@ -515,7 +508,7 @@ def smart_fit_selection( y_data=numpy.array(region_data["Activity_MBq"]), fixed_params=None, num_exponentials=order, - bounds=self.config["rois"][region]["bounds"], + bounds=self.config["VOIs"][region]["bounds"], params_init={"A1": activity_init}, with_uptake=with_uptake, ) @@ -529,8 +522,8 @@ def smart_fit_selection( # If only one model fit, that is the winner. if len(aic_results) == 1: - self.config["rois"][region]["with_uptake"] = fit_config[0][0] - self.config["rois"][region]["fit_order"] = fit_config[0][1] + self.config["VOIs"][region]["with_uptake"] = fit_config[0][0] + self.config["VOIs"][region]["fit_order"] = fit_config[0][1] return all_fits[0] # If there are two more models, we check the top two models and compare their AIC. If the difference @@ -544,8 +537,8 @@ def smart_fit_selection( ): best_model_idx = aic_results[1][0] - self.config["rois"][region]["with_uptake"] = fit_config[best_model_idx][0] - self.config["rois"][region]["fit_order"] = fit_config[best_model_idx][1] + self.config["VOIs"][region]["with_uptake"] = fit_config[best_model_idx][0] + self.config["VOIs"][region]["fit_order"] = fit_config[best_model_idx][1] return all_fits[best_model_idx] @@ -636,8 +629,8 @@ def calculate_bed(self, kinetic: str) -> None: RADIOBIOLOGY_DATA_FILE = Path(this_dir, "data", "radiobiology.json") with open(RADIOBIOLOGY_DATA_FILE) as f: self.radiobiology_dic = json.load(f) - bed_df = self.df_ad[ - self.df_ad.index.isin(list(self.radiobiology_dic.keys())) + bed_df = self.results_dosimetry_organs[ + self.results_dosimetry_organs.index.isin(list(self.radiobiology_dic.keys())) ] # only organs that we know the radiobiology parameters organs = numpy.array(bed_df.index.unique()) bed = {} @@ -646,7 +639,7 @@ def calculate_bed(self, kinetic: str) -> None: t_repair = self.radiobiology_dic[organ]["t_repair"] alpha_beta = self.radiobiology_dic[organ]["alpha_beta"] AD = ( - float(self.df_ad.loc[bed_df.index == organ]["AD[Gy/GBq]"].values[0]) + float(bed_df.loc[bed_df.index == organ]["AD_total[Gy/GBq]"].values[0]) * float(self.config["InjectedActivity"]) / 1000 ) # Gy @@ -702,7 +695,9 @@ def calculate_bed(self, kinetic: str) -> None: ) print(f"{organ}", bed[organ]) - self.df_ad["BED[Gy]"] = self.df_ad.index.map(bed) + self.results_dosimetry_organs["BED[Gy]"] = ( + self.results_dosimetry_organs.index.map(bed) + ) def save_images_and_masks_at(self, time_id: int) -> None: """Save CT, NM and masks for a specific time point. @@ -718,7 +713,7 @@ def save_images_and_masks_at(self, time_id: int) -> None: time_id=time_id, out_path=self.db_dir, name="SPECT" ) self.nm_data.save_masks_to_nii_at( - time_id=time_id, out_path=self.db_dir, regions=self.config["rois"] + time_id=time_id, out_path=self.db_dir, regions=self.config["VOIs"] ) return None @@ -774,18 +769,21 @@ def write_json_data( cycle["InjectionTime"] = self.config["InjectionTime"] cycle["InjectedActivity"] = self.config["InjectedActivity"] cycle["Weight_g"] = self.config["PatientWeight_g"] + cycle["Height_cm"] = self.config["PatientHeight_cm"] cycle["Level"] = self.config["Level"] - cycle["Method"] = self.config["Method"] - cycle["OutputFormat"] = self.config["OutputFormat"] - cycle["ScaleDoseByDensity"] = self.config.get( - "ScaleDoseByDensity", cycle.get("ScaleDoseByDensity", "NA") - ) + if cycle["Level"] == "Organ": + cycle["Method"] = self.config["OrganLevel"] + elif cycle["Level"] == "Voxel": + cycle["Method"] = self.config["VoxelLevel"] + cycle["ScaleDoseByDensity"] = self.config.get( + "ScaleDoseByDensity", cycle.get("ScaleDoseByDensity", "NA") + ) cycle["ReferenceTimePoint"] = self.config["ReferenceTimePoint"] cycle["TimePoints_h"] = self.results["Time_hr"][0] - for organ in self.config["rois"].keys(): - if organ not in cycle["rois"]: - cycle["rois"][organ] = { + for organ in self.config["VOIs"].keys(): + if organ not in cycle["VOIs"]: + cycle["VOIs"][organ] = { "volumes_mL": {}, "activity_MBq": {}, "timepoints_h": {}, @@ -817,32 +815,32 @@ def write_json_data( "BED_Gy_uncertainty": {}, } - cycle["rois"][organ]["volumes_mL"]["different_tps"] = self.results.loc[ + cycle["VOIs"][organ]["volumes_mL"]["different_tps"] = self.results.loc[ organ, "Volume_CT_mL" ] - cycle["rois"][organ]["volumes_mL"]["uncertainty"] = "NA" - cycle["rois"][organ]["volumes_mL"]["mean"] = numpy.mean( + cycle["VOIs"][organ]["volumes_mL"]["uncertainty"] = "NA" + cycle["VOIs"][organ]["volumes_mL"]["mean"] = numpy.mean( self.results.loc[organ, "Volume_CT_mL"] ) - cycle["rois"][organ]["volumes_mL"]["mean_uncertainty"] = "NA" - cycle["rois"][organ]["activity_MBq"]["values"] = self.results.loc[ - organ, "Activity_MBq" + cycle["VOIs"][organ]["volumes_mL"]["mean_uncertainty"] = "NA" + cycle["VOIs"][organ]["activity_MBq"]["values"] = [ + float(x) for x in self.results.loc[organ, "Activity_MBq"] ] - cycle["rois"][organ]["activity_MBq"]["uncertainty"] = "NA" - cycle["rois"][organ]["timepoints_h"]["values"] = self.results.loc[ + cycle["VOIs"][organ]["activity_MBq"]["uncertainty"] = "NA" + cycle["VOIs"][organ]["timepoints_h"]["values"] = self.results.loc[ organ, "Time_hr" ] - cycle["rois"][organ]["doserate_MBq_per_h"]["values"] = "NA" - cycle["rois"][organ]["doserate_MBq_per_h"]["uncertainty"] = "NA" + cycle["VOIs"][organ]["doserate_MBq_per_h"]["values"] = "NA" + cycle["VOIs"][organ]["doserate_MBq_per_h"]["uncertainty"] = "NA" try: - cycle["rois"][organ]["density_HU"]["different_tps"] = self.results.loc[ + cycle["VOIs"][organ]["density_HU"]["different_tps"] = self.results.loc[ organ, "Density_HU" ] except (KeyError, AttributeError): # TODO: Handle errors explicitly pass - cycle["rois"][organ]["density_HU"]["uncertainty"] = "NA" + cycle["VOIs"][organ]["density_HU"]["uncertainty"] = "NA" try: - cycle["rois"][organ]["density_HU"]["mean"] = numpy.mean( + cycle["VOIs"][organ]["density_HU"]["mean"] = numpy.mean( self.results.loc[organ, "Density_HU"] ) except ( @@ -851,96 +849,107 @@ def write_json_data( TypeError, ): # TODO: Handle errors explicitly pass - cycle["rois"][organ]["density_HU"]["mean_uncertainty"] = "NA" - cycle["rois"][organ]["density_gml"]["different_tps"] = "NA" - cycle["rois"][organ]["density_gml"]["uncertainty"] = "NA" - cycle["rois"][organ]["density_gml"]["mean"] = "NA" - cycle["rois"][organ]["density_gml"]["mean_uncertainty"] = "NA" - cycle["rois"][organ]["mass_g"]["different_tps"] = "NA" - cycle["rois"][organ]["mass_g"]["uncertainty"] = "NA" - cycle["rois"][organ]["mass_g"]["mean"] = "NA" - cycle["rois"][organ]["mass_g"]["mean_uncertainty"] = "NA" - cycle["rois"][organ]["fitting_eq"] = self.config["rois"][organ]["fit_order"] - cycle["rois"][organ]["no_of_fit_params"] = "NA" - cycle["rois"][organ]["fit_params"] = list( + cycle["VOIs"][organ]["density_HU"]["mean_uncertainty"] = "NA" + cycle["VOIs"][organ]["density_gml"]["different_tps"] = "NA" + cycle["VOIs"][organ]["density_gml"]["uncertainty"] = "NA" + cycle["VOIs"][organ]["density_gml"]["mean"] = "NA" + cycle["VOIs"][organ]["density_gml"]["mean_uncertainty"] = "NA" + cycle["VOIs"][organ]["mass_g"]["different_tps"] = "NA" + cycle["VOIs"][organ]["mass_g"]["uncertainty"] = "NA" + cycle["VOIs"][organ]["mass_g"]["mean"] = "NA" + cycle["VOIs"][organ]["mass_g"]["mean_uncertainty"] = "NA" + cycle["VOIs"][organ]["fitting_eq"] = self.config["VOIs"][organ]["fit_order"] + cycle["VOIs"][organ]["no_of_fit_params"] = "NA" + cycle["VOIs"][organ]["fit_params"] = list( self.results.loc[organ, "Fit_params"] ) - cycle["rois"][organ]["fit_params_uncertainty"] = "NA" - cycle["rois"][organ]["R_2"] = self.results.loc[organ, "R_squared_AIC"][0] - cycle["rois"][organ]["AIC"] = self.results.loc[organ, "R_squared_AIC"][1] - cycle["rois"][organ]["TIA_MBqh"] = self.results.loc[organ, "TIA_MBq_h"] - cycle["rois"][organ]["TIA_MBqh_uncertainty"] = "NA" - cycle["rois"][organ]["TIA_h"] = self.results.loc[organ, "TIA_h"] - cycle["rois"][organ]["TIA_h_uncertainty"] = "NA" - cycle["rois"][organ]["mean_AD_Gy"] = "NA" - cycle["rois"][organ]["mean_AD_Gy_uncertainty"] = "NA" - cycle["rois"][organ]["min_AD_Gy"] = "NA" - cycle["rois"][organ]["max_AD_Gy"] = "NA" - cycle["rois"][organ]["peak_AD_Gy"] = "NA" - cycle["rois"][organ]["repair_halflife"] = "NA" - cycle["rois"][organ]["alpha_beta"] = "NA" - cycle["rois"][organ]["composition"] = "NA" - cycle["rois"][organ]["total_s_value"] = "NA" - cycle["rois"][organ]["total_s_value_uncertainty"] = "NA" + cycle["VOIs"][organ]["washout_ratio"] = self.config["VOIs"][organ][ + "washout_ratio" + ] + cycle["VOIs"][organ]["fit_params_uncertainty"] = "NA" + cycle["VOIs"][organ]["R_2"] = ( + "NA" + if pandas.isna(self.results.loc[organ, "R_squared_AIC"][0]) + else self.results.loc[organ, "R_squared_AIC"][0] + ) + cycle["VOIs"][organ]["AIC"] = ( + "NA" + if pandas.isna(self.results.loc[organ, "R_squared_AIC"][1]) + else self.results.loc[organ, "R_squared_AIC"][1] + ) + cycle["VOIs"][organ]["TIA_MBqh"] = self.results.loc[organ, "TIA_MBq_h"] + cycle["VOIs"][organ]["TIA_MBqh_uncertainty"] = "NA" + cycle["VOIs"][organ]["TIA_h"] = self.results.loc[organ, "TIA_h"] + cycle["VOIs"][organ]["TIA_h_uncertainty"] = "NA" + cycle["VOIs"][organ]["mean_AD_Gy"] = "NA" + cycle["VOIs"][organ]["mean_AD_Gy_uncertainty"] = "NA" + cycle["VOIs"][organ]["min_AD_Gy"] = "NA" + cycle["VOIs"][organ]["max_AD_Gy"] = "NA" + cycle["VOIs"][organ]["peak_AD_Gy"] = "NA" + cycle["VOIs"][organ]["repair_halflife"] = "NA" + cycle["VOIs"][organ]["alpha_beta"] = "NA" + cycle["VOIs"][organ]["composition"] = "NA" + cycle["VOIs"][organ]["total_s_value"] = "NA" + cycle["VOIs"][organ]["total_s_value_uncertainty"] = "NA" if "Lesion" in organ or "TTB" in organ: - cycle["rois"][organ]["density_gml"]["different_tps"] = "NA" - cycle["rois"][organ]["density_gml"]["uncertainty"] = "NA" - cycle["rois"][organ]["density_gml"]["mean"] = self.results_lesions.loc[ - organ, "Density_g_per_mL" - ] - cycle["rois"][organ]["density_gml"]["mean_uncertainty"] = "NA" - cycle["rois"][organ]["mass_g"]["different_tps"] = "NA" - cycle["rois"][organ]["mass_g"]["uncertainty"] = "NA" - cycle["rois"][organ]["mass_g"]["mean"] = self.results_lesions.loc[ - organ, "Mass_g" - ] - cycle["rois"][organ]["mass_g"]["mean_uncertainty"] = "NA" - cycle["rois"][organ]["composition"] = self.results_lesions.loc[ - organ, "Composition" - ] - cycle["rois"][organ]["total_s_value"] = self.results_lesions.loc[ - organ, "Total_S_Value" - ] - cycle["rois"][organ]["total_s_value_uncertainty"] = "NA" - cycle["rois"][organ]["mean_AD_Gy"] = self.results_lesions.loc[ + cycle["VOIs"][organ]["density_gml"]["different_tps"] = "NA" + cycle["VOIs"][organ]["density_gml"]["uncertainty"] = "NA" + cycle["VOIs"][organ]["density_gml"]["mean"] = ( + self.results_dosimetry_lesions.loc[organ, "Density_g_per_mL"] + ) + cycle["VOIs"][organ]["density_gml"]["mean_uncertainty"] = "NA" + cycle["VOIs"][organ]["mass_g"]["different_tps"] = "NA" + cycle["VOIs"][organ]["mass_g"]["uncertainty"] = "NA" + cycle["VOIs"][organ]["mass_g"]["mean"] = ( + self.results_dosimetry_lesions.loc[organ, "Mass_g"] + ) + cycle["VOIs"][organ]["mass_g"]["mean_uncertainty"] = "NA" + cycle["VOIs"][organ]["composition"] = ( + self.results_dosimetry_lesions.loc[organ, "Composition"] + ) + cycle["VOIs"][organ]["total_s_value"] = ( + self.results_dosimetry_lesions.loc[organ, "Total_S_Value"] + ) + cycle["VOIs"][organ]["total_s_value_uncertainty"] = "NA" + cycle["VOIs"][organ]["mean_AD_Gy"] = self.results_dosimetry_lesions.loc[ organ, "AD_Gy" ] - cycle["rois"][organ]["mean_AD_Gy_uncertainty"] = "NA" + cycle["VOIs"][organ]["mean_AD_Gy_uncertainty"] = "NA" if "BoneMarrow" in organ: - cycle["rois"][organ]["volumes_mL"]["different_tps"] = 1170 - cycle["rois"][organ]["volumes_mL"]["uncertainty"] = "NA" - cycle["rois"][organ]["volumes_mL"]["mean"] = 1170 + cycle["VOIs"][organ]["volumes_mL"]["different_tps"] = 1170 + cycle["VOIs"][organ]["volumes_mL"]["uncertainty"] = "NA" + cycle["VOIs"][organ]["volumes_mL"]["mean"] = 1170 if "Gland" in organ: - cycle["rois"][organ]["density_gml"]["different_tps"] = "NA" - cycle["rois"][organ]["density_gml"]["uncertainty"] = "NA" - cycle["rois"][organ]["density_gml"]["mean"] = ( - self.results_salivaryglands.loc[organ, "Density_g_per_mL"] + cycle["VOIs"][organ]["density_gml"]["different_tps"] = "NA" + cycle["VOIs"][organ]["density_gml"]["uncertainty"] = "NA" + cycle["VOIs"][organ]["density_gml"]["mean"] = ( + self.results_dosimetry_salivaryglands.loc[organ, "Density_g_per_mL"] ) - cycle["rois"][organ]["density_gml"]["mean_uncertainty"] = "NA" - cycle["rois"][organ]["mass_g"]["different_tps"] = "NA" - cycle["rois"][organ]["mass_g"]["uncertainty"] = "NA" - cycle["rois"][organ]["mass_g"]["mean"] = ( - self.results_salivaryglands.loc[organ, "Mass_g"] + cycle["VOIs"][organ]["density_gml"]["mean_uncertainty"] = "NA" + cycle["VOIs"][organ]["mass_g"]["different_tps"] = "NA" + cycle["VOIs"][organ]["mass_g"]["uncertainty"] = "NA" + cycle["VOIs"][organ]["mass_g"]["mean"] = ( + self.results_dosimetry_salivaryglands.loc[organ, "Mass_g"] ) - cycle["rois"][organ]["mass_g"]["mean_uncertainty"] = "NA" - cycle["rois"][organ]["composition"] = self.results_salivaryglands.loc[ - organ, "Composition" - ] - cycle["rois"][organ]["total_s_value"] = self.results_salivaryglands.loc[ - organ, "Total_S_Value" - ] - cycle["rois"][organ]["total_s_value_uncertainty"] = "NA" - cycle["rois"][organ]["mean_AD_Gy"] = self.results_salivaryglands.loc[ - organ, "AD_Gy" - ] - cycle["rois"][organ]["mean_AD_Gy_uncertainty"] = "NA" + cycle["VOIs"][organ]["mass_g"]["mean_uncertainty"] = "NA" + cycle["VOIs"][organ]["composition"] = ( + self.results_dosimetry_salivaryglands.loc[organ, "Composition"] + ) + cycle["VOIs"][organ]["total_s_value"] = ( + self.results_dosimetry_salivaryglands.loc[organ, "Total_S_Value"] + ) + cycle["VOIs"][organ]["total_s_value_uncertainty"] = "NA" + cycle["VOIs"][organ]["mean_AD_Gy"] = ( + self.results_dosimetry_salivaryglands.loc[organ, "AD_Gy"] + ) + cycle["VOIs"][organ]["mean_AD_Gy_uncertainty"] = "NA" if self.config["Level"] == "Organ": - for organ in self.df_ad.index: - if organ in self.df_ad.index: + for organ in self.results_dosimetry_organs.index: + if organ in self.results_dosimetry_organs.index: cycle["Organ-level_AD"][organ] = { "AD[Gy/GBq]": {}, "AD[Gy/GBq]_uncertianty": {}, @@ -949,19 +958,21 @@ def write_json_data( "BED[Gy]": {}, "BED[Gy]_uncertianty": {}, } - cycle["Organ-level_AD"][organ]["AD[Gy/GBq]"] = self.df_ad.loc[ - organ, "AD[Gy/GBq]" - ] + cycle["Organ-level_AD"][organ]["AD[Gy/GBq]"] = ( + self.results_dosimetry_organs.loc[organ, "AD_total[Gy/GBq]"] + ) cycle["Organ-level_AD"][organ]["AD[Gy/GBq]_uncertainty"] = "NA" - cycle["Organ-level_AD"][organ]["AD[Gy]"] = self.df_ad.loc[ - organ, "AD[Gy]" - ] + cycle["Organ-level_AD"][organ]["AD[Gy]"] = ( + self.results_dosimetry_organs.loc[organ, "AD_total[Gy]"] + ) cycle["Organ-level_AD"][organ]["AD[Gy]_uncertainty"] = "NA" - if "BED[Gy]" in self.df_ad.columns: + if "BED[Gy]" in self.results_dosimetry_organs.columns: cycle["Organ-level_AD"][organ]["BED[Gy]"] = ( - self.df_ad.loc[organ, "BED[Gy]"] - if pandas.notna(self.df_ad.loc[organ, "BED[Gy]"]) + self.results_dosimetry_organs.loc[organ, "BED[Gy]"] + if pandas.notna( + self.results_dosimetry_organs.loc[organ, "BED[Gy]"] + ) else "NA" ) else: @@ -969,7 +980,9 @@ def write_json_data( cycle["Organ-level_AD"][organ]["BED[Gy]_uncertianty"] = "NA" - if "Yes" in self.config["LesionDosimetry"]: + if "Yes" in self.config["OrganLevel"]["AdditionalOptions"].get( + "LesionDosimetry" + ): cycle["Organ-level_AD"]["TTB"] = { "mass_g": {}, "volumes_mL": {}, @@ -979,22 +992,23 @@ def write_json_data( "AD[Gy/GBq]": {}, "AD[Gy/GBq]_uncertianty": {}, } - cycle["Organ-level_AD"]["TTB"]["mass_g"] = self.results_lesions.loc[ - "TTB", "Mass_g" - ] - cycle["Organ-level_AD"]["TTB"]["volumes_mL"] = self.results_lesions.loc[ - "TTB", "Volume_CT_mL" - ] - cycle["Organ-level_AD"]["TTB"]["TIA_h"] = self.results_lesions.loc[ - "TTB", "TIA_h" - ] - cycle["Organ-level_AD"]["TTB"]["AD[Gy]"] = self.results_lesions.loc[ - "TTB", "AD_Gy" - ] + cycle["Organ-level_AD"]["TTB"]["mass_g"] = ( + self.results_dosimetry_lesions.loc["TTB", "Mass_g"] + ) + cycle["Organ-level_AD"]["TTB"]["volumes_mL"] = ( + self.results_dosimetry_lesions.loc["TTB", "Volume_CT_mL"] + ) + cycle["Organ-level_AD"]["TTB"]["TIA_h"] = ( + self.results_dosimetry_lesions.loc["TTB", "TIA_h"] + ) + cycle["Organ-level_AD"]["TTB"]["AD[Gy]"] = ( + self.results_dosimetry_lesions.loc["TTB", "AD_Gy"] + ) cycle["Organ-level_AD"]["TTB"]["AD[Gy]_uncertainty"] = "NA" - cycle["Organ-level_AD"]["TTB"]["AD[Gy/GBq]"] = self.results_lesions.loc[ - "TTB", "AD_Gy" - ] / (float(self.config["InjectedActivity"]) / 1000) + cycle["Organ-level_AD"]["TTB"]["AD[Gy/GBq]"] = ( + self.results_dosimetry_lesions.loc["TTB", "AD_Gy"] + / (float(self.config["InjectedActivity"]) / 1000) + ) cycle["Organ-level_AD"]["TTB"]["AD[Gy/GBq]_uncertianty"] = "NA" with open(file_path, "w") as file: diff --git a/pytheranostics/dosimetry/organ_s_dosimetry.py b/pytheranostics/dosimetry/organ_s_dosimetry.py index 062cbec..0f460d5 100644 --- a/pytheranostics/dosimetry/organ_s_dosimetry.py +++ b/pytheranostics/dosimetry/organ_s_dosimetry.py @@ -55,13 +55,21 @@ def check_mandatory_fields_organ(self) -> None: return None @staticmethod - def _load_human_mass_table() -> pandas.DataFrame: + def _load_human_mass_target_organs_table(self) -> pandas.DataFrame: """Load the reference human phantom masses.""" with resource_path( - "pytheranostics.data", "phantom/human/human_phantom_masses.csv" + "pytheranostics.data", + f"phantom/human/ICRP_mass_{self.config['Gender'].lower()}_target.csv", + ) as masses_path: + self.target_organ_masses = pandas.read_csv(masses_path, index_col=0) + + def _load_human_mass_source_organs_table(self) -> pandas.DataFrame: + """Load the reference human phantom masses.""" + with resource_path( + "pytheranostics.data", + f"phantom/human/ICRP_mass_{self.config['Gender'].lower()}_source.csv", ) as masses_path: - masses = pandas.read_csv(masses_path, index_col=0) - return masses + self.source_organ_masses = pandas.read_csv(masses_path, index_col=0) def composition_and_density_from_HU(self, density: float) -> Tuple[str, float]: """Determine composition and density for a given CT HU value.""" @@ -123,22 +131,27 @@ def apply_sphere_method(self, df: pandas.DataFrame) -> pandas.DataFrame: def calculate_ttb(self): """Compute Total Tumor Burden (TTB) metrics and append to results_lesions.""" metrics = { - "Mass_g": self.results_lesions["Mass_g"].sum(), - "Volume_CT_mL": self.results_lesions["Volume_CT_mL"].sum(), - "TIA_h": self.results_lesions["TIA_h"].sum(), + "Mass_g": self.results_dosimetry_lesions["Mass_g"].sum(), + "Volume_CT_mL": self.results_dosimetry_lesions["Volume_CT_mL"].sum(), + "TIA_h": self.results_dosimetry_lesions["TIA_h"].sum(), "AD_Gy": ( - (self.results_lesions["Mass_g"] * self.results_lesions["AD_Gy"]).sum() + ( + self.results_dosimetry_lesions["Mass_g"] + * self.results_dosimetry_lesions["AD_Gy"] + ).sum() ) / ( - self.results_lesions["Mass_g"].sum() - if self.results_lesions["Mass_g"].sum() > 0 + self.results_dosimetry_lesions["Mass_g"].sum() + if self.results_dosimetry_lesions["Mass_g"].sum() > 0 else 0 ), } TTB = pandas.DataFrame(metrics, index=["TTB"]) - self.results_lesions = pandas.concat([self.results_lesions, TTB], axis=0) - return self.results_lesions + self.results_dosimetry_lesions = pandas.concat( + [self.results_dosimetry_lesions, TTB], axis=0 + ) + return self.results_dosimetry_lesions def prepare_data(self) -> None: """ @@ -219,6 +232,7 @@ def prepare_data(self) -> None: self.results_fitting.loc["Red Marrow"][ "Volume_CT_mL" ] = 1170 # TODO volume hardcoded, think about alternatives + self.results_fitting.loc["RemainderOfBody"]["Volume_CT_mL"] = ( self.config["PatientWeight_g"] - self.results_fitting.loc[ @@ -237,23 +251,25 @@ def prepare_data(self) -> None: self.results_fitting_organs = self.results_fitting[ ~lesion_mask ].copy() # all non-lesion entries - self.results_fitting_lesions = self.results_fitting[ - lesion_mask - ].copy() # only lesion entries - + self.results_fitting_lesions = self.results[ + self.results.index.str.contains("Lesion") + ] + elif "No" in self.config["OrganLevel"]["AdditionalOptions"].get( + "LesionDosimetry" + ): + self.results_fitting_organs = self.results_fitting.copy() if "TotalTumorBurden" in self.results_fitting.index: self.results_fitting.drop("TotalTumorBurden", axis=0, inplace=True) - if output_type == "Export": fmt = organ_conf["Output"]["ExportFormat"] if fmt.lower() == "olinda": - self.results_fitting = self.results_fitting.rename( + self.results_fitting_organs = self.results_fitting_organs.rename( index={"RemainderOfBody": "Total Body"} ) - self.results_fitting.loc["Total Body"]["Volume_CT_mL"] = ( + self.results_fitting_organs.loc["Total Body"]["Volume_CT_mL"] = ( self.config["PatientWeight_g"] - - self.results_fitting.loc[ - ~self.results_fitting.index.isin( + - self.results_fitting_organs.loc[ + ~self.results_fitting_organs.index.isin( ["Total Body", "RemainderOfBody"] ), "Volume_CT_mL", @@ -263,6 +279,10 @@ def prepare_data(self) -> None: print("Not Implemented yet.") else: print(f"Export format {fmt} not recognized.") + elif "Lesion" in self.config["Level"]: + self.results_fitting_lesions = self.results[ + self.results.index.str.contains("Lesion") + ] return None @@ -299,7 +319,6 @@ def process_dosimetry(self) -> None: output_type = organ_conf["Output"]["Type"] # 'Export' or 'Calculate' self.prepare_data() - print("Processing dosimetry at the organ level of OARs") if output_type == "Export": fmt = organ_conf["Output"]["ExportFormat"] @@ -314,7 +333,8 @@ def process_dosimetry(self) -> None: print(f"Export format {fmt} not recognized.") elif output_type == "Calculate": - self.calculate_absorbed_dose() + if "Organ" in self.config["Level"]: + self.calculate_absorbed_dose() else: raise ValueError(f"Unknown output type: {output_type}") @@ -323,7 +343,7 @@ def process_dosimetry(self) -> None: "LesionDosimetry" ): self.results_dosimetry_lesions = self.apply_sphere_method( - self.results_fitting_lesions.index.str.contains("Lesion") + self.results_fitting_lesions ) self.results_dosimetry_lesions = self.calculate_ttb() if "Yes" in self.config["OrganLevel"]["AdditionalOptions"].get( @@ -338,12 +358,12 @@ def calculate_absorbed_dose(self) -> pandas.DataFrame: """Calculate absorbed dose per target organ based on model and disintegration data.""" model_files = { "Female": { - "beta": f'177Lu_S_values_female_{self.config["OrganLevel"]["Calculation"]["SValueSource"].lower()}_BETA.csv', - "gamma": f'177Lu_S_values_female_{self.config["OrganLevel"]["Calculation"]["SValueSource"].lower()}_GAMMA.csv', + "gamma": f'{self.config["Radionuclide"]}_S_values_female_{self.config["OrganLevel"]["Calculation"]["SValueSource"].lower()}_GAMMA.csv', + "beta": f'{self.config["Radionuclide"]}_S_values_female_{self.config["OrganLevel"]["Calculation"]["SValueSource"].lower()}_BETA.csv', }, "Male": { - "beta": f'177Lu_S_values_male_{self.config["OrganLevel"]["Calculation"]["SValueSource"].lower()}_BETA.csv', - "gamma": f'177Lu_S_values_male_{self.config["OrganLevel"]["Calculation"]["SValueSource"].lower()}_GAMMA.csv', + "gamma": f'{self.config["Radionuclide"]}_S_values_male_{self.config["OrganLevel"]["Calculation"]["SValueSource"].lower()}_GAMMA.csv', + "beta": f'{self.config["Radionuclide"]}_S_values_male_{self.config["OrganLevel"]["Calculation"]["SValueSource"].lower()}_BETA.csv', }, } @@ -359,7 +379,10 @@ def calculate_absorbed_dose(self) -> pandas.DataFrame: svalues_gamma = self.load_svalues(gamma_path) print("Source organs available in the model:", svalues_beta.columns.tolist()) - print("Source organs present :", self.results_fitting.index.tolist()) + print("Source organs present :", self.results_fitting_organs.index.tolist()) + + self._load_human_mass_target_organs_table() + self._load_human_mass_source_organs_table() self.source_organs_missing = set(svalues_beta.columns) - set( self.results_fitting.index @@ -391,7 +414,7 @@ def calculate_absorbed_dose(self) -> pandas.DataFrame: ) # Apply mass scaling - dose_df = self.perform_mass_scaling(dose_df, self.config["Gender"]) + dose_df = self.perform_mass_scaling(dose_df) # Calculate absorbed dose in Gy for injected activity dose_df["AD_total[Gy/GBq]"] = ( @@ -403,7 +426,10 @@ def calculate_absorbed_dose(self) -> pandas.DataFrame: dose_df["AD_total[Gy]"] = dose_df["AD_total[Gy/GBq]"] / 1000 * injected_activity dose_df = dose_df.reset_index(drop=True) - self.df_ad = dose_df.copy() + self.results_dosimetry_organs = dose_df.copy() + self.results_dosimetry_organs = self.results_dosimetry_organs.set_index( + "Target organ" + ) return dose_df def hollow_organ_correction(self, df: pandas.DataFrame) -> pandas.DataFrame: @@ -439,14 +465,14 @@ def redistribute_ROB_into_source_organs_missing( if "RemainderOfBody" not in tia_series.index: return tia_series - missing_organs_df = self.organ_masses.loc[ - self.organ_masses.index.isin(self.source_organs_missing) + missing_organs_df = self.source_organ_masses.loc[ + self.source_organ_masses.index.isin(self.source_organs_missing) ] print(f"Missing organs DataFrame:\n{missing_organs_df.index.tolist()}") missing_organs_df = missing_organs_df.drop(index="Total Body") - mass_source_organs = self.organ_masses.loc[ + mass_source_organs = self.source_organ_masses.loc[ [org for org in tia_series.index if org != "RemainderOfBody"], "Mass_g" ].sum() @@ -475,13 +501,7 @@ def redistribute_ROB_into_source_organs_missing( def apply_s_value(self, tia_df, s_values, radiation_type) -> pandas.DataFrame: """Multiply S-values by TIA to compute dose matrix for radiation type.""" - # Path to organ masses - masses = self._load_human_mass_table() - gender = self.config.get("Gender", "Male") - if gender not in masses.columns: - raise ValueError(f"Unknown gender '{gender}' for mass table.") - self.organ_masses = masses[[gender]].rename(columns={gender: "Mass_g"}) - + print(f"Applying S-values for {radiation_type} radiation...") # Handle remainder of the body # Redistribute ROB TIA into missing source organs if needed - approach consistent with MIRDcalc software if "RemainderOfBody" in tia_df.index: @@ -548,39 +568,49 @@ def apply_s_value(self, tia_df, s_values, radiation_type) -> pandas.DataFrame: dose_df = s_values_subset.multiply(tia_series, axis=1) # ROB + print("Calculating Remainder of Body contributions...") dose_df["RemainderOfBody"] = 0 - total_body_mass = self.config["PatientWeight_g"] - ROB_mass = tia_df.loc["RemainderOfBody", "Volume_CT_mL"] + total_body_mass = self.source_organ_masses.loc["Total Body", "Mass_g"] + source_organs = tia_series.index if radiation_type == "beta": - organs = s_values_subset.index.difference( + target_organs = s_values_subset.index.difference( tia_series.index.difference(["Red Marrow", "Osteogenic Cells"]) ) + ROB_mass = tia_df.loc["RemainderOfBody", "Volume_CT_mL"] if radiation_type == "gamma": - organs = s_values_subset.index - # adjust source organs so that they are different for gamma and beta radiation (in beta it is total - source organs + bone + skeleton) + target_organs = s_values_subset.index + total_mass_source_organs = self.source_organ_masses.loc[ + self.source_organ_masses.index.intersection(tia_series.index), + "Mass_g", + ].sum() - for target_organ in organs: + ROB_mass = total_body_mass - total_mass_source_organs + + for target_organ in target_organs: contribution_from_sources = 0 - for source_organ in s_values.columns.difference(tia_series.index): - if target_organ.split()[0] == source_organ.split()[0]: - continue - if target_organ == "Osteogenic Cells" and source_organ in [ - "Cortical Bone", - "Trabecular Bone", - "Red Marrow", - ]: - continue - if target_organ == "Red Marrow" and source_organ in [ - "Trabecular Bone" - ]: + for source_organ in source_organs: + if radiation_type == "beta": + if target_organ.split()[0] == source_organ.split()[0]: + continue + if target_organ == "Red Marrow" and source_organ in [ + "Trabecular Bone" + ]: + continue + if target_organ == "Osteogenic Cells" and source_organ in [ + "Cortical Bone", + "Trabecular Bone", + "Red Marrow", + ]: + continue + + if source_organ == "Total Body": continue + if target_organ == "Total Body": continue - if source_organ == "Total Body": - continue - source_organ_mass = self.organ_masses.loc[ + source_organ_mass = self.source_organ_masses.loc[ source_organ, "Mass_g" ] @@ -591,14 +621,20 @@ def apply_s_value(self, tia_df, s_values, radiation_type) -> pandas.DataFrame: source_organ_mass / ROB_mass ) - s_value_ROB_to_target = ( - s_values.loc[target_organ, "Total Body"] - * (total_body_mass / ROB_mass) - ) - contribution_from_sources - dose_df.at[target_organ, "RemainderOfBody"] = ( + if target_organ == "Total Body": + s_value_ROB_to_target = ( + s_values.loc[target_organ, "Total Body"] + # * (total_body_mass / ROB_mass) since it is a target organ, the scaling will happen at mass scaling method, so not to have it twice its not performed here + ) - contribution_from_sources + else: + s_value_ROB_to_target = ( + s_values.loc[target_organ, "Total Body"] + * (total_body_mass / ROB_mass) + ) - contribution_from_sources + + dose_df.at[target_organ, "RemainderOfBody"] = float( s_value_ROB_to_target * tia_df.loc["RemainderOfBody", "TIA_s"] ) - else: # TODO print("No ROB option selected. Proceeding without ROB handling.") @@ -606,21 +642,28 @@ def apply_s_value(self, tia_df, s_values, radiation_type) -> pandas.DataFrame: return dose_df def perform_mass_scaling( - self, df: pandas.DataFrame, gender: str + self, + df: pandas.DataFrame, ) -> pandas.DataFrame: """Apply mass scaling to absorbed dose calculations based on patient-specific organ masses.""" - masses = self._load_human_mass_table() - if gender not in masses.columns: - raise ValueError(f"Unknown gender '{gender}' for mass table.") - model_masses_df = masses[[gender]].rename(columns={gender: "Mass_g"}) - print("Performing mass scaling...") for organ in df["Target organ"]: + mass_key = ( + "Total Body" if organ in ["Total Body", "RemainderOfBody"] else organ + ) + fit_key = ( + "RemainderOfBody" + if organ in ["Total Body", "RemainderOfBody"] + else organ + ) - if organ in model_masses_df.index and organ in self.results_fitting.index: - model_mass = model_masses_df.loc[organ, "Mass_g"] - patient_mass = self.results_fitting.loc[organ, "Volume_CT_mL"] + if ( + mass_key in self.target_organ_masses.index + and fit_key in self.results_fitting.index + ): + model_mass = self.target_organ_masses.loc[mass_key, "Mass_g"] + patient_mass = self.results_fitting.loc[fit_key, "Volume_CT_mL"] if ( pandas.notna(patient_mass) @@ -630,7 +673,6 @@ def perform_mass_scaling( if "AD_beta[Gy/GBq]" in df.columns: scaling_factor_beta = model_mass / patient_mass - print(f"Scaling factor for {organ} [β]: {scaling_factor_beta}") df.loc[ df["Target organ"] == organ, "AD_beta[Gy/GBq]" ] *= scaling_factor_beta @@ -674,13 +716,13 @@ def create_Olinda_file(self, dirname: str, savefile: bool = False) -> None: ind = template[template["Data"] == "[BEGIN NUCLIDES]"].index template.loc[ind[0] + 1, "Data"] = formatted_radionuclide + "|" - for organ in self.results_fitting.index: + for organ in self.results_fitting_organs.index: indices = template[template["Data"].str.contains(organ)].index source_organ = template.iloc[indices[0]].str.split("|")[0][0] mass_phantom = template.iloc[indices[0]].str.split("|")[0][1] - kinetic_data = self.results_fitting.loc[organ]["TIA_h"] - mass_data = round(self.results_fitting.loc[organ]["Volume_CT_mL"], 1) + kinetic_data = self.results_fitting_organs.loc[organ]["TIA_h"] + mass_data = round(self.results_fitting_organs.loc[organ]["Volume_CT_mL"], 1) # Update the template DataFrame template.iloc[indices[0]] = ( @@ -763,7 +805,7 @@ def read_Olinda_results(self, olinda_results_path: str) -> None: df_ad["Total"].fillna(0, inplace=True) df_ad["AD[Gy]"] = df_ad["Total"] * float(self.config["InjectedActivity"]) / 1000 df_ad = df_ad.rename(columns={"Total": "AD[Gy/GBq]"}) - self.df_ad = df_ad + self.results_dosimetry_organs = df_ad def compute_dose(self): """Compute Time Integrated Activity.""" diff --git a/pytheranostics/fits/fits.py b/pytheranostics/fits/fits.py index a5cc94c..e8874dc 100644 --- a/pytheranostics/fits/fits.py +++ b/pytheranostics/fits/fits.py @@ -23,6 +23,7 @@ def exponential_fit_lmfit( bounds: Optional[Dict[str, Tuple[float, float]]] = None, params_init: Optional[Dict[str, float]] = None, with_uptake: bool = False, + washout_ratio: Optional[int] = None, ) -> Tuple[lmfit.model.ModelResult, Callable]: """Fit data to a sum of exponentials with flexible parameter fixing using lmfit. @@ -105,6 +106,8 @@ def model_func(x, **params): will_have_expr = True if num_exponentials == 3 and with_uptake and amp_name == "C1": will_have_expr = True + if num_exponentials == 2 and washout_ratio and amp_name == "B2": + will_have_expr = True if amp_name in fixed_params: params.add(amp_name, value=fixed_params[amp_name], vary=False) @@ -158,6 +161,18 @@ def model_func(x, **params): "Parameters 'A1', 'B1', and 'C1' must be present to apply the constraint 'C1 = -(A1 + B1)'." ) + if num_exponentials == 2 and washout_ratio is not None: + if "B2" in params and "A2" in params: + print( + f"Applying constant ratio between decay constants of washout phases: B2 = {washout_ratio} * A2" + ) + params["B2"].set(expr=f"{washout_ratio} * A2") + else: + raise ValueError( + "Parameters 'A2' and 'B2' must be present to apply the constraint 'B2 = washout_ratio * A2'." + ) + + print(y_data, params) # Perform the fit result = model.fit( y_data, params, x=x_data, weights=None if sigma is None else 1.0 / sigma diff --git a/pytheranostics/imaging_ds/cycle_loader.py b/pytheranostics/imaging_ds/cycle_loader.py index b3ebee8..90d8a36 100644 --- a/pytheranostics/imaging_ds/cycle_loader.py +++ b/pytheranostics/imaging_ds/cycle_loader.py @@ -171,12 +171,19 @@ def extract_injection_from_first_tp_spect( weight_g = int(round(float(pw) * 1000.0)) except Exception: pass - + height_cm: Optional[int] = None + pw = getattr(ds, "PatientSize", None) # cm + if pw is not None: + try: + height_cm = int(round(float(pw) * 100.0)) + except Exception: + pass return { "InjectionDate": inj_date or "", "InjectionTime": inj_time or "", "InjectedActivity": injected_activity, # Bq (int) or None "PatientWeight_g": weight_g, + "PatientHeight_cm": height_cm, } diff --git a/pytheranostics/misc_tools/report_generator.py b/pytheranostics/misc_tools/report_generator.py index 759540f..4c8bcf6 100644 --- a/pytheranostics/misc_tools/report_generator.py +++ b/pytheranostics/misc_tools/report_generator.py @@ -11,6 +11,7 @@ from reportlab.lib.units import inch from reportlab.platypus import ( Image, + PageBreak, Paragraph, SimpleDocTemplate, Spacer, @@ -25,21 +26,33 @@ def signature_block(person, styles, width=2.5 * inch, height=0.6 * inch): Placeholder for signature, line, and text. Returns a small table that can be inserted side-by-side with others. """ - # Empty row for signature space - sig_space = Spacer(1, height) + elements = [] + + # If a signature path is provided, insert the image + if person.get("signature"): + sig_img = Image(person["signature"]) + # Scale the image to fit width and maintain aspect ratio + sig_img.drawHeight = height + sig_img.drawWidth = width + elements.append(sig_img) + else: + # Empty space if no signature image + elements.append(Spacer(1, height)) # Line row line = Table([[""]], colWidths=[width]) line.setStyle(TableStyle([("LINEABOVE", (0, 0), (-1, -1), 1, colors.black)])) + elements.append(line) # Text row text = Paragraph( f"{person['name']}
{person['title']}
{person['affiliation']}", styles["Normal"], ) + elements.append(text) - # Wrap everything in a column table (1 col, 3 rows) - block = Table([[sig_space], [line], [text]], colWidths=[width]) + # Wrap everything in a column table (1 col, stacked) + block = Table([[e] for e in elements], colWidths=[width]) block.setStyle(TableStyle([("ALIGN", (0, 0), (-1, -1), "CENTER")])) return block @@ -69,23 +82,26 @@ def create_dosimetry_pdf(json_file, output_file, calculated_by=None, approved_by # Title title = Paragraph( - "DOSIMETRY REPORT", styles["Title"] + "DOSIMETRY ASSESSMENT", styles["Title"] ) elements.append(title) elements.append(Spacer(1, 0.5 * inch)) # Subject Information Section elements.append(Paragraph("Subject Information", styles["Heading2"])) - # Subject Information Table subject_data = [ - ["Clinical Trial", "PR.21"], + ["Clinical Trial", data.get("ClinicalTrial")], + ["Radiopharmaceutical", "177Lu-PSMA-617"], + ["Mode of administration", "I.V."], ["ID", data.get("PatientID")], ["Sex", data.get("Gender")], + ["Weight kg", data.get("Cycle_01", {})[0].get("Weight_g", "") / 1000], + ["Height cm", data.get("Cycle_01", {})[0].get("Height_cm", "")], ["Number of cycles ", data.get("No_of_completed_cycles")], ] - subject_table = Table(subject_data, colWidths=[1.5 * inch, 4 * inch]) + subject_table = Table(subject_data, colWidths=[2 * inch, 3.5 * inch]) subject_table.setStyle( TableStyle( [ @@ -146,11 +162,89 @@ def create_dosimetry_pdf(json_file, output_file, calculated_by=None, approved_by elements.append(caption) elements.append(Spacer(1, 0.2 * inch)) + elements.append(PageBreak()) + for i in range(1, data.get("No_of_completed_cycles") + 1): cycle_info(i, elements, styles, data) + elements.append(PageBreak()) + elements.append(Paragraph("Patient Summary", styles["Heading2"])) + # =============================== + # CUMULATIVE ORGAN TABLE + # =============================== + elements.append( + Paragraph("Cumulative Absorbed Dose Summary", styles["Heading3"]) + ) + + total_tia_kidneys = 0 + total_tia_salivary = 0 + total_tia_marrow = 0 + + total_ad_kidneys = 0 + total_ad_salivary = 0 + total_ad_marrow = 0 + + total_bed_kidneys = 0 + + for i in range(1, data.get("No_of_completed_cycles") + 1): + therapy_info = data.get(f"Cycle_0{i}", {})[0] + + # Kidneys + total_tia_kidneys += ( + therapy_info["VOIs"]["Kidney_Left"]["TIA_h"] + + therapy_info["VOIs"]["Kidney_Right"]["TIA_h"] + ) + total_ad_kidneys += therapy_info["Organ-level_AD"]["Kidneys"]["AD[Gy]"] + total_bed_kidneys += therapy_info["Organ-level_AD"]["Kidneys"]["BED[Gy]"] + + # Red marrow + total_tia_marrow += therapy_info["VOIs"]["BoneMarrow"]["TIA_h"] + total_ad_marrow += therapy_info["Organ-level_AD"]["Red Marrow"]["AD[Gy]"] + + # Salivary glands + total_tia_salivary += ( + therapy_info["VOIs"]["ParotidGland_Left"]["TIA_h"] + + therapy_info["VOIs"]["ParotidGland_Right"]["TIA_h"] + + therapy_info["VOIs"]["SubmandibularGland_Left"]["TIA_h"] + + therapy_info["VOIs"]["SubmandibularGland_Right"]["TIA_h"] + ) + total_ad_salivary += therapy_info["Organ-level_AD"]["Salivary Glands"]["AD[Gy]"] + + # Build the cumulative table + cumulative_data = [ + ["Organ", "Cumulative TIA (h)", "Cumulative AD (Gy)", "Cumulative BED (Gy)"], + [ + "Kidneys", + round(total_tia_kidneys, 2), + round(total_ad_kidneys, 2), + round(total_bed_kidneys, 2), + ], + ["Red Marrow", round(total_tia_marrow, 2), round(total_ad_marrow, 2), "-"], + [ + "Salivary glands", + round(total_tia_salivary, 2), + round(total_ad_salivary, 2), + "-", + ], + ] + + cumulative_table = Table(cumulative_data, colWidths=[1.5 * inch, 1.7 * inch]) + cumulative_table.setStyle( + TableStyle( + [ + ("ALIGN", (0, 0), (-1, -1), "LEFT"), + ("VALIGN", (0, 0), (-1, -1), "MIDDLE"), + ("FONTNAME", (0, 0), (-1, -1), "Helvetica"), + ("FONTSIZE", (0, 0), (-1, -1), 12), + ("GRID", (0, 0), (-1, -1), 1, colors.black), + ("BACKGROUND", (0, 0), (0, -1), colors.lightblue), + ] + ) + ) + + elements.append(cumulative_table) # Paths to your three images image_paths = [ calling_folder / "TestDoseDB/Gy_cummulative.png", @@ -162,7 +256,7 @@ def create_dosimetry_pdf(json_file, output_file, calculated_by=None, approved_by imgs = [] for path in image_paths: img = Image(str(path)) - scale = min(max_width / img.imageWidth, max_height / img.imageHeight) / 2.8 + scale = min(max_width / img.imageWidth, max_height / img.imageHeight) / 3 img.drawWidth = img.imageWidth * scale img.drawHeight = img.imageHeight * scale imgs.append(img) @@ -186,7 +280,7 @@ def create_dosimetry_pdf(json_file, output_file, calculated_by=None, approved_by # Add caption elements.append(Spacer(1, 0.2 * inch)) caption = Paragraph( - "Figure 2: Cumulative absorbed dose, absorbed dose per cycle, and absorbed dose per GBq per cycle for target organs.", + "Figure 2: Cumulative AD, AD per cycle, and AD per GBq per cycle for target organs.", styles["Normal"], ) elements.append(caption) @@ -205,7 +299,7 @@ def create_dosimetry_pdf(json_file, output_file, calculated_by=None, approved_by img = Image(str(path)) # Scale to fit 2×2 layout scale = min( - (max_width / 2.2) / img.imageWidth, (max_height / 2.2) / img.imageHeight + (max_width / 2.5) / img.imageWidth, (max_height / 2.5) / img.imageHeight ) img.drawWidth = img.imageWidth * scale img.drawHeight = img.imageHeight * scale @@ -214,7 +308,7 @@ def create_dosimetry_pdf(json_file, output_file, calculated_by=None, approved_by # Arrange in 2×2 structure trend_table_data = [[trend_imgs[0], trend_imgs[1]], [trend_imgs[2], trend_imgs[3]]] - trend_table = Table(trend_table_data, colWidths=[max_width / 2.2, max_width / 2.2]) + trend_table = Table(trend_table_data, colWidths=[max_width / 2.5, max_width / 2.5]) trend_table.setStyle( TableStyle( @@ -226,18 +320,14 @@ def create_dosimetry_pdf(json_file, output_file, calculated_by=None, approved_by ] ) ) - - elements.append(Spacer(1, 0.3 * inch)) elements.append(trend_table) # Add caption - elements.append(Spacer(1, 0.2 * inch)) caption = Paragraph( "Figure 3: Trends of hematological and renal function, and PSA.", styles["Normal"], ) elements.append(caption) - elements.append(Spacer(1, 0.2 * inch)) # =============================== # Signatures Section @@ -272,11 +362,27 @@ def create_dosimetry_pdf(json_file, output_file, calculated_by=None, approved_by ) elements.append(app_table) + # =============================== + # Appendix + # =============================== + + elements.append(PageBreak()) + title = Paragraph("Appendix", styles["Heading2"]) + elements.append(title) + fig_title = Paragraph( + "Biodistribution and kinetic analysis", + styles["Heading3"], + ) + elements.append(fig_title) + + for i in range(1, data.get("No_of_completed_cycles") + 1): + biodistribution_per_cycle(i, elements, styles, data) + # Build PDF doc.build(elements) -def cycle_info(cycle_n, elements, styles, data): +def biodistribution_per_cycle(cycle_n, elements, styles, data): """Add cycle information to the PDF report elements. Parameters @@ -290,54 +396,18 @@ def cycle_info(cycle_n, elements, styles, data): data : dict Patient data dictionary. """ - # Therapy Information Section + if cycle_n == 1: + pass + else: + elements.append(PageBreak()) + therapy_title = Paragraph(f"Cycle {cycle_n}", styles["Heading2"]) elements.append(therapy_title) - # Therapy Information Table - therapy_info = data.get(f"Cycle_0{cycle_n}", {}) - therapy_data = [ - ["Radiopharmaceutical", "177Lu-PSMA-617"], - ["Mode of administration", "I.V."], - ["Administered Activity (MBq)", therapy_info[0].get("InjectedActivity", "")], - [ - "Date of injection", - ( - datetime.strptime( - therapy_info[0].get("InjectionDate", ""), "%Y%m%d" - ).strftime("%Y-%m-%d") - if therapy_info[0].get("InjectionDate", "") - else "" - ), - ], - ] - - therapy_table = Table(therapy_data, colWidths=[2.5 * inch, 3 * inch]) - therapy_table.setStyle( - TableStyle( - [ - ("ALIGN", (0, 0), (-1, -1), "LEFT"), - ("VALIGN", (0, 0), (-1, -1), "MIDDLE"), - ("FONTNAME", (0, 0), (-1, -1), "Helvetica"), - ("FONTSIZE", (0, 0), (-1, -1), 12), - ("GRID", (0, 0), (-1, -1), 1, colors.black), - ("BACKGROUND", (0, 0), (0, -1), colors.lightgrey), - ] - ) - ) - - elements.append(therapy_table) - page_width, page_height = letter max_width = page_width - 2 * 72 # 1-inch margins max_height = page_height - 2 * 72 - fig_title = Paragraph( - "Time-activity curves, fit functions, and fit parameters", - styles["Heading3"], - ) - elements.append(fig_title) - calling_folder = Path().absolute() # notebook folder pattern = calling_folder / f"TestDoseDB/*_fit_Cycle_0{cycle_n}.png" @@ -357,19 +427,62 @@ def cycle_info(cycle_n, elements, styles, data): elements.append(Spacer(1, 0.3 * inch)) elements.append(img) + +def cycle_info(cycle_n, elements, styles, data): + """Add cycle information to the PDF report elements. + + Parameters + ---------- + cycle_n : int + The cycle number. + elements : list + List of reportlab elements to append to. + styles : dict + ReportLab styles dictionary. + data : dict + Patient data dictionary. + """ + # Therapy Information Section + therapy_title = Paragraph(f"Cycle {cycle_n}", styles["Heading2"]) + elements.append(therapy_title) + + # Therapy Information Table + therapy_info = data.get(f"Cycle_0{cycle_n}", {}) + + # Format injection date nicely + inj_date_raw = therapy_info[0].get("InjectionDate", "") + inj_date = ( + datetime.strptime(inj_date_raw, "%Y%m%d").strftime("%Y-%m-%d") + if inj_date_raw + else "" + ) + + # Build a clean paragraph instead of a table + therapy_info_para = Paragraph( + f"" + f"Administered activity: {therapy_info[0].get('InjectedActivity', '')} MBq
" + f"Date of injection: {inj_date}" + f"", + styles["Normal"], + ) + + # Add to document + elements.append(therapy_info_para) + elements.append(Spacer(1, 0.089 * inch)) + fig_title = Paragraph( "Absorbed dose results for the organs at risk", styles["Heading3"] ) elements.append(fig_title) organ_data_Gy_GBq = [ - ["Organ", "TIA (h)", "AD (Gy/GBq)", "AD(Gy)", "BED (Gy)"], + ["Organ", "TIA (h)", "AD (Gy/GBq)", "AD (Gy)", "BED (Gy)"], [ "Kidneys", round( ( - therapy_info[0]["rois"]["Kidney_Left"]["TIA_h"] - + therapy_info[0]["rois"]["Kidney_Right"]["TIA_h"] + therapy_info[0]["VOIs"]["Kidney_Left"]["TIA_h"] + + therapy_info[0]["VOIs"]["Kidney_Right"]["TIA_h"] ), 2, ), @@ -379,7 +492,7 @@ def cycle_info(cycle_n, elements, styles, data): ], [ "Red Marrow", - round((therapy_info[0]["rois"]["BoneMarrow"]["TIA_h"]), 2), + round((therapy_info[0]["VOIs"]["BoneMarrow"]["TIA_h"]), 2), round(therapy_info[0]["Organ-level_AD"]["Red Marrow"]["AD[Gy/GBq]"], 2), round(therapy_info[0]["Organ-level_AD"]["Red Marrow"]["AD[Gy]"], 2), "-", @@ -388,10 +501,10 @@ def cycle_info(cycle_n, elements, styles, data): "Salivary glands", round( ( - therapy_info[0]["rois"]["ParotidGland_Left"]["TIA_h"] - + therapy_info[0]["rois"]["ParotidGland_Right"]["TIA_h"] - + therapy_info[0]["rois"]["SubmandibularGland_Left"]["TIA_h"] - + therapy_info[0]["rois"]["SubmandibularGland_Right"]["TIA_h"] + therapy_info[0]["VOIs"]["ParotidGland_Left"]["TIA_h"] + + therapy_info[0]["VOIs"]["ParotidGland_Right"]["TIA_h"] + + therapy_info[0]["VOIs"]["SubmandibularGland_Left"]["TIA_h"] + + therapy_info[0]["VOIs"]["SubmandibularGland_Right"]["TIA_h"] ), 2, ), diff --git a/pytheranostics/misc_tools/tools.py b/pytheranostics/misc_tools/tools.py index 44b447a..a9c6c40 100644 --- a/pytheranostics/misc_tools/tools.py +++ b/pytheranostics/misc_tools/tools.py @@ -3,7 +3,6 @@ from datetime import datetime from typing import Dict, Tuple -import matplotlib.pyplot as plt import numpy from scipy.ndimage import median_filter @@ -199,9 +198,9 @@ def extract_exponential_params_from_json( and all the parameters of the fit. """ with_uptake = False - # Determine order: - parameters = json_data[cycle][0]["rois"][region]["fit_params"] + parameters = json_data[cycle][0]["VOIs"][region]["fit_params"] + washout_ratio = json_data[cycle][0]["VOIs"][region]["washout_ratio"] # Handle Legacy: if len(parameters) in [3, 5]: @@ -214,7 +213,6 @@ def extract_exponential_params_from_json( fixed_parameters: Dict[str, float] = {} all_parameters: Dict[str, float] = {} - for order in range(fit_order): fixed_parameters[f"{param_name_base[order]}2"] = parameters[ exponential_idxs[order] @@ -226,14 +224,16 @@ def extract_exponential_params_from_json( all_parameters[f"{param_name_base[order]}2"] = parameters[ exponential_idxs[order] ] - + if washout_ratio is not None: + # Fix A1 ONLY + if "B1" in all_parameters: + fixed_parameters["B1"] = all_parameters["B1"] if fit_order == 2 and parameters[0] == -parameters[2]: with_uptake = True if fit_order == 3 and parameters[4] == -(parameters[0] + parameters[2]): with_uptake = True - - return fixed_parameters, with_uptake, all_parameters + return fixed_parameters, with_uptake, all_parameters, washout_ratio def extract_exponential_params_from_json_legacy( @@ -266,7 +266,7 @@ def extract_exponential_params_from_json_legacy( When the parameter configuration is incompatible with new format. """ # Read Parameters: - parameters = json_data[cycle][0]["rois"][region]["fit_params"] + parameters = json_data[cycle][0]["VOIs"][region]["fit_params"] if len(parameters) not in [3, 5]: raise AssertionError( @@ -321,7 +321,7 @@ def initialize_biokinetics_from_prior_cycle( dict Updated configuration dictionary. """ - for roi, roi_info in config["rois"].items(): + for roi, roi_info in config["VOIs"].items(): if ( "biokinectics_from_previous_cycle" in roi_info @@ -329,15 +329,18 @@ def initialize_biokinetics_from_prior_cycle( ): # Get previous cycle parameters: - fixed_param, with_uptake, all_params = extract_exponential_params_from_json( - json_data=prior_treatment_data, cycle=cycle, region=roi + fixed_param, with_uptake, all_params, washout_ratio = ( + extract_exponential_params_from_json( + json_data=prior_treatment_data, cycle=cycle, region=roi + ) ) - config["rois"][roi] = { + config["VOIs"][roi] = { "fixed_parameters": fixed_param, "fit_order": len(all_params) // 2, "param_init": all_params, "with_uptake": with_uptake, + "washout_ratio": washout_ratio, } print( @@ -347,45 +350,3 @@ def initialize_biokinetics_from_prior_cycle( print("") return config - - -def plot_MIP(SPECT=None, vmax=300000, figsize=(10, 5), ax=None): - """Plot Maximum Intensity Projection (MIP) of SPECT data. - - Parameters - ---------- - ax : matplotlib.axes.Axes, optional - Axes object to plot on. If None, a new figure and axes will be created. - SPECT : numpy.ndarray - SPECT data array. - vmax : int, optional - Maximum value for display, by default 300000. - figsize : tuple, optional - Figure size (width, height) in inches, by default (10, 5). - Only used if ax is None. - - Returns - ------- - fig : matplotlib.figure.Figure - Figure object. - ax : matplotlib.axes.Axes - Axes object with the plot. - """ - if ax is None: - fig, ax = plt.subplots(figsize=figsize) - else: - fig = ax.get_figure() - - plt.sca(ax) - plt.imshow( - SPECT.max(axis=0).T, cmap="Greys", interpolation="Gaussian", vmax=vmax, vmin=0 - ) - plt.xlim(30, 100) - plt.ylim(0, 234) - - plt.axis("off") - - plt.xticks([]) - plt.yticks([]) - - return fig, ax diff --git a/pytheranostics/plots/plots.py b/pytheranostics/plots/plots.py index 1d9f1ec..4d9307b 100644 --- a/pytheranostics/plots/plots.py +++ b/pytheranostics/plots/plots.py @@ -86,7 +86,7 @@ def plot_tac_residuals( else: weights = None # Generate x-values for plotting the fitted model starting from x=0 - x_fit = numpy.linspace(0, x_data[-1] * 3, 500) + x_fit = numpy.linspace(0, x_data.max() * 3, 500) y_fit = result.eval(x=x_fit) # First subplot: Linear scale plot @@ -96,7 +96,7 @@ def plot_tac_residuals( # Plot fitted model ax1.plot(x_fit, y_fit, color="red") ax1.set_xlim(left=0) # Start x-axis from zero - ax1.set_xlim(right=x_data[-1] * 2) # Start y-axis from zero + ax1.set_xlim(right=x_data.max() * 2) # Start y-axis from zero ax1.set_ylim(bottom=0) # Start y-axis from zero ax1.set_title(region) ax1.set_xlabel(x_label) @@ -119,7 +119,7 @@ def plot_tac_residuals( # Plot fitted model ax2.plot(x_fit, y_fit, color="red") ax2.set_xlim(left=0) # Start x-axis from zero - ax2.set_xlim(right=x_data[-1] * 2) # Start y-axis from zero + ax2.set_xlim(right=x_data.max() * 2) # Start y-axis from zero ax2.set_yscale("log") ax2.set_title(title_text) ax2.set_xlabel(x_label) @@ -147,3 +147,78 @@ def plot_tac_residuals( plt.show() return None + + +def plot_MIP_with_mask_outlines(ax, SPECT, masks=None, vmax=300000, label=None): + """Plot Maximum Intensity Projection (MIP) of SPECT data with masks outlines. + + Parameters + ---------- + ax : _type_ + _description_ + SPECT : _type_ + _description_ + masks : _type_, optional + _description_, by default None + vmax : int, optional + _description_, by default 300000 + """ + plt.sca(ax) + spect_mip = SPECT.max(axis=0) + plt.imshow(spect_mip.T, cmap="Greys", interpolation="Gaussian", vmax=vmax, vmin=0) + + if masks is not None: + for organ, mask in masks.items(): + organ_lower = organ.lower() + print(organ_lower) + if "peak" in organ_lower: + continue + else: + if "kidney" in organ_lower: + color = "lime" + elif "parotid" in organ_lower: + color = "red" + elif "submandibular" in organ_lower: + color = "red" + elif "lesion" in organ_lower: + color = "m" + else: + continue + + mip_mask = mask.max(axis=0) + if mip_mask.shape != spect_mip.shape: + mip_mask = mip_mask.T + + plt.contour( + numpy.transpose(mip_mask, (1, 0)), + levels=[0.5], + colors=[color], + linewidths=1.5, + alpha=0.5, + ) + + # --- Add label at mask center --- + if label is True: + ys, xs = numpy.where(mip_mask > 0) + + if len(xs) > 0: + # Corrected for transpose in contour + x_center = ys.mean() - 0.2 * ys.mean() # was xs + y_center = xs.mean() # was ys + + plt.text( + x_center, + y_center, + organ, + color=color, + fontsize=8, + ha="center", + va="center", + alpha=0.7, + ) + + plt.xlim(30, 100) + plt.ylim(0, 234) + plt.axis("off") + plt.xticks([]) + plt.yticks([])