diff --git a/.gitignore b/.gitignore index 58c2889..61c596f 100644 --- a/.gitignore +++ b/.gitignore @@ -100,3 +100,8 @@ ENV/ # mypy .mypy_cache/ +ogo/util/.DS_Store +ogo/util/graphcuts/.DS_Store +ogo/.vscode/launch.json +.vscode/launch.json +.DS_Store diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..54c917e --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,5 @@ +{ + "idf.pythonInstallPath": "/Users/skboyd/miniconda3/bin/python", + "python-envs.defaultEnvManager": "ms-python.python:conda", + "python-envs.defaultPackageManager": "ms-python.python:conda" +} \ No newline at end of file diff --git a/ogo/cli/AnalyzeBMD.py b/ogo/cli/AnalyzeBMD.py index 0c7cd4d..abe1b3f 100644 --- a/ogo/cli/AnalyzeBMD.py +++ b/ogo/cli/AnalyzeBMD.py @@ -22,12 +22,15 @@ # Start Script -def AnalyzeBMD(image_filename, mask_filename, labels, output_filename, noheader, overwrite=False): +def AnalyzeBMD(image_filename, mask_filename, labels, output_filename, custom_labels=None, noheader=False, overwrite=False): if output_filename: ogo.message("Start of Script...") - # Master list of labels - labelsDict = lb.labels_dict + if custom_labels: + [labelsDict,labelsHdr] = ogo.custom_labelsDict(custom_labels) + else: + # Master list of labels + labelsDict = lb.labels_dict # Check if output exists and should overwrite ogo.pass_check_if_output_exists(output_filename,overwrite) @@ -62,7 +65,7 @@ def AnalyzeBMD(image_filename, mask_filename, labels, output_filename, noheader, if labels[0] == 0: if output_filename: ogo.message('No labels specified. Building list...') - labels = np.unique(array) + labels = np.unique(array).tolist() # Remove label 0 from list (background) if labels[0] == 0: @@ -78,7 +81,7 @@ def AnalyzeBMD(image_filename, mask_filename, labels, output_filename, noheader, os.sys.exit() # if (str(labels[idx]) not in labelsDict): if labels[idx] not in labelsDict: - print("ERROR: Label {:d} does not exist in master list of labels.".format(labels[idx])) + print("ERROR: Label {:d} does not exist in list of labels.".format(labels[idx])) os.sys.exit() if output_filename: ogo.message('!> label {:3d} {:12s} --> {:d} voxels'.format(labels[idx], labelsDict[labels[idx]]['LABEL'], @@ -108,10 +111,13 @@ def AnalyzeBMD(image_filename, mask_filename, labels, output_filename, noheader, parameters_dict['Label'] = lab parameters_dict['LabelDesc'] = labelsDict[lab]['LABEL'] - parameters_dict['Integral BMD [mg/cc]'] = '{:.3f}'.format(bmd_outcomes['Integral BMD [mg/cc]']) - parameters_dict['Integral BMC [mg]'] = '{:.3f}'.format(bmd_outcomes['Integral BMC [mg]']) - parameters_dict['Bone Volume [mm^3]'] = '{:.3f}'.format(bmd_outcomes['Bone Volume [mm^3]']) - parameters_dict['Bone Volume [cm^3]'] = '{:.3f}'.format(bmd_outcomes['Bone Volume [cm^3]']) + parameters_dict['Integral BMD [mg/cc]'] = '{:.4f}'.format(bmd_outcomes['Integral BMD [mg/cc]']) + parameters_dict['Integral BMC [mg]'] = '{:.4f}'.format(bmd_outcomes['Integral BMC [mg]']) + parameters_dict['Bone Volume [mm^3]'] = '{:.4f}'.format(bmd_outcomes['Bone Volume [mm^3]']) + parameters_dict['Bone Volume [cm^3]'] = '{:.4f}'.format(bmd_outcomes['Bone Volume [cm^3]']) + parameters_dict['mean_signal'] = '{:.4f}'.format(bmd_outcomes['mean_signal']) + parameters_dict['std_noise'] = '{:.4f}'.format(bmd_outcomes['std_noise']) + parameters_dict['SNR'] = '{:.4f}'.format(bmd_outcomes['SNR']) if output_filename: if idx == 0 and not noheader: @@ -144,6 +150,8 @@ def main(): If no output text file is defined then results are output to the screen. +Outputs the mean and standard deviation of the ROI to allow SNR calculation. + WARNING: If your input image is not calibrated then the results here will be incorrect. There is no calibration done as part of this application. @@ -171,9 +179,11 @@ def main(): help='Specify list of labels; space separated (e.g. 1 2 3) (Default: report all available)') parser.add_argument('--output_filename', default=None, metavar='TEXTFILE', help='Output file name (*.txt) (Default: dump to screen)') + parser.add_argument('--custom_labels', default=None, metavar='TEXTFILE', + help='Custom labels of type ITKSnap (*.txt) (Default: OgoMasterLabels)') parser.add_argument('--noheader', action='store_true', help='Suppress printing header (default: %(default)s)') parser.add_argument('--overwrite', action='store_true', help='Overwrite output without asking') - print() + #print() # Parse and display args = parser.parse_args() diff --git a/ogo/cli/AssignPhantomLabels.py b/ogo/cli/AssignPhantomLabels.py new file mode 100644 index 0000000..b4ce887 --- /dev/null +++ b/ogo/cli/AssignPhantomLabels.py @@ -0,0 +1,278 @@ +# /------------------------------------------------------------------------------+ +# | 24-APR-2026 | +# | Copyright (c) Bone Imaging Laboratory | +# | All rights reserved | +# | bonelab@ucalgary.ca | +# +------------------------------------------------------------------------------+ + +# Post-processing script for ML-based calibration phantom rod segmentation. +# Takes the raw CT scan and the ML model's binary rod segmentation (all rods +# labeled 1) and automatically identifies the correct calibration phantom by +# fitting the measured rod HU values against known phantom densities from ogo. +# Rods are then relabeled with their correct phantom-specific labels. + +# Imports +import argparse +import inspect +import re + +import numpy as np +import nibabel as nib +import scipy as sp +from ogo.calib.standard_calibration import StandardCalibration +from ogo.calib.mindways_calibration import MindwaysCalibration +import ogo.util.Helper as ogo +from ogo.util.echo_arguments import echo_arguments +from ogo.util.cli_tools import ( + check_overwrite, + validate_input_file +) + + +#------------------------------------------------------------------------------+ +# Helper Functions +def get_known_phantoms(): + """Return a list of all known phantom names.""" + source = inspect.getsource(ogo.get_phantom) + known_phantoms = re.findall(r"phantom_type in '([^']+)'", source) + return known_phantoms + + +def get_all_phantoms(): + """Return a list of all phantom dicts.""" + known_phantoms = get_known_phantoms() + return [ogo.get_phantom(name) for name in known_phantoms] + + +def extract_rod_statistics(raw_array, rod_array, quiet): + """Extract connected components and compute mean HU for each rod.""" + labeled_array, num_features = sp.ndimage.label(rod_array) + n_rods = num_features + + if not quiet: + ogo.message(f'Found {n_rods} rods') + + rod_stats = [] + + for rod_id in range(1, n_rods + 1): + rod_mask = (labeled_array == rod_id) + mean_hu = raw_array[rod_mask].mean() + + rod_stats.append({ + "rod_id": rod_id, + "mean_hu": mean_hu, + "mask": rod_mask + }) + + if not quiet: + ogo.message(f' Rod {rod_id}: mean HU = {mean_hu:.2f}') + + # Sort rods by mean HU ascending to match to phantom densities later + rod_stats.sort(key=lambda x: x["mean_hu"]) + + return rod_stats + + +def fit_phantom(rod_stats, phantom_dict, quiet): + """Fit phantom calibration to rod statistics.""" + hu_values = [rod["mean_hu"] for rod in rod_stats] + densities = phantom_dict["densities"] + + if len(hu_values) != len(densities): + return None + + # Sort densities ascending and reorder HU to match + paired = sorted(zip(densities, hu_values), key=lambda x: x[0]) + densities_sorted = [d for d, h in paired] + hu_sorted = [h for d, h in paired] + + try: + if phantom_dict["type"] == "k2hpo4": + water = phantom_dict["h2o_densities"] + + # Sort all three lists together by density ascending + triples = sorted(zip(densities, hu_values, water), key=lambda x: x[0]) + densities_sorted = [d for d, h, w in triples] + hu_sorted = [h for d, h, w in triples] + water_sorted = [w for d, h, w in triples] + + calibrator = MindwaysCalibration() + calibrator.fit(hu_sorted, densities_sorted, water_sorted) + else: + calibrator = StandardCalibration() + calibrator.fit(hu_sorted, densities_sorted) + + return calibrator + + except Exception as e: + if not quiet: + ogo.message(f' Fit failed for {phantom_dict["name"]}: {e}') + return None + + +def rank_phantoms(rod_stats, all_phantoms, quiet): + """Try fitting every known phantom and rank by R^2.""" + n_rods = len(rod_stats) + results = [] + + for phantom_dict in all_phantoms: + if phantom_dict["number_rods"] != n_rods: + continue + + calibrator = fit_phantom(rod_stats, phantom_dict, quiet) + + if calibrator is not None: + r2 = calibrator.r_value ** 2 + if not quiet: + ogo.message(f' {phantom_dict["name"]}: R^2 = {r2:.4f}') + results.append((phantom_dict, calibrator, r2)) + + results.sort(key=lambda x: x[2], reverse=True) + return results + + +def assign_labels(rod_stats, best_phantom, quiet): + """Assign labels by matching HU rank to density rank.""" + density_label_pairs = sorted( + zip(best_phantom["densities"], best_phantom["rod_labels"]), + key=lambda x: x[0] + ) + + if not quiet: + ogo.message('\nLabel assignments:') + + for rod, (density, label) in zip(rod_stats, density_label_pairs): + rod["assigned_label"] = label + rod["assigned_density"] = density + if not quiet: + ogo.message(f' Rod {rod["rod_id"]}: mean HU = {rod["mean_hu"]:.2f} -> label {label} ({density} mg/cc)') + + return rod_stats + + +def save_labeled_nifti(rod_stats, reference_img, output_file, quiet): + """Build output label image with correct rod labels and save.""" + output_array = np.zeros(reference_img.shape, dtype=np.int16) + + for rod in rod_stats: + output_array[rod["mask"]] = rod["assigned_label"] + + output_img = nib.Nifti1Image(output_array, affine=reference_img.affine, header=reference_img.header) + nib.save(output_img, output_file) + + if not quiet: + ogo.message(f'\nSaved labeled output to: {output_file}') + + +#------------------------------------------------------------------------------+ +# Main function +def AssignPhantomLabels(raw_nifti_file, rod_nifti_file, output_file, overwrite, quiet): + """Assign correct labels to segmented phantom rods based on HU values. + + This program automatically identifies the calibration phantom from segmented + phantom rods by fitting measured rod HU values against known phantom densities. + Rods are then relabeled with their correct phantom-specific labels. + """ + + # Check if output exists and should overwrite + if not check_overwrite(output_file, overwrite): + return 1 + + # Load images - raw NIfTI image and segmented image of rods + if not quiet: + ogo.message(f'Reading raw image: {raw_nifti_file}') + + raw_img = nib.load(raw_nifti_file) + raw_array = raw_img.get_fdata() + + if not quiet: + ogo.message(f'Reading rod segmentation: {rod_nifti_file}') + + rod_img = nib.load(rod_nifti_file) + rod_array = rod_img.get_fdata() + + # Extract rod statistics + rod_stats = extract_rod_statistics(raw_array, rod_array, quiet) + + if len(rod_stats) == 0: + ogo.message('ERROR: No rods found in segmentation.') + return 1 + + # Get all known phantoms and rank by fit + all_phantoms = get_all_phantoms() + ranked = rank_phantoms(rod_stats, all_phantoms, quiet) + + if not ranked: + ogo.message(f'ERROR: No known phantom matched {len(rod_stats)} rods. Check segmentation.') + return 1 + + # Select best matching phantom + best_phantom, best_calibrator, best_r2 = ranked[0] + if not quiet: + ogo.message(f'\nBest match: {best_phantom["name"]} (R^2 = {best_r2:.4f})') + + if best_r2 < 0.99: + ogo.message(f' WARNING: R^2 = {best_r2:.4f} is low — check segmentation') + + # Assign labels and save output + rod_stats = assign_labels(rod_stats, best_phantom, quiet) + save_labeled_nifti(rod_stats, rod_img, output_file, quiet) + + if not quiet: + ogo.message('[DONE]') + + return 0 + + +def main(): + """Setup argument parsing and run program.""" + + description = '''Assign correct labels to segmented calibration phantom rods. + +Post-processing script for ML-based calibration phantom rod segmentation. +Takes the raw CT scan and the ML model's binary rod segmentation (all rods +labeled 1) and automatically identifies the correct calibration phantom by +fitting the measured rod HU values against known phantom densities. Rods are +then relabeled with their correct phantom-specific labels. + +Output: _labeled.nii.gz — rod mask with correct phantom labels +''' + + epilog = '''Example calls: +ogoAssignPhantomLabels raw.nii.gz rods.nii.gz output.nii.gz +ogoAssignPhantomLabels raw.nii.gz rods.nii.gz output.nii.gz --quiet +''' + + # Setup argument parsing + parser = argparse.ArgumentParser( + formatter_class=argparse.RawTextHelpFormatter, + prog="ogoAssignPhantomLabels", + description=description, + epilog=epilog + ) + + parser.add_argument('raw_nifti_file', type=validate_input_file(['.nii', '.nii.gz']), + help='Raw CT image file (*.nii, *.nii.gz)') + parser.add_argument('rod_nifti_file', type=validate_input_file(['.nii', '.nii.gz']), + help='Rod segmentation file (*.nii, *.nii.gz)') + parser.add_argument('output_file', type=str, + help='Output labeled rod segmentation file (*.nii, *.nii.gz)') + parser.add_argument('--overwrite', '-ow', action='store_true', + help='Overwrite output without asking') + parser.add_argument('--quiet', '-q', action='store_true', + help='Suppress informational output to stdout') + + print() + + # Parse and display + args = parser.parse_args() + + if not args.quiet: + ogo.message(echo_arguments('AssignPhantomLabels', vars(args))) + + # Run program + return AssignPhantomLabels(**vars(args)) + + +if __name__ == '__main__': + exit(main()) diff --git a/ogo/cli/CalculateLinearAttenuation.py b/ogo/cli/CalculateLinearAttenuation.py new file mode 100644 index 0000000..efc5763 --- /dev/null +++ b/ogo/cli/CalculateLinearAttenuation.py @@ -0,0 +1,297 @@ +# /------------------------------------------------------------------------------+ +# | 7-NOV-2025 | +# | Copyright (c) Bone Imaging Laboratory | +# | All rights reserved | +# | bonelab@ucalgary.ca | +# +------------------------------------------------------------------------------+ + +# Imports +import sys +import os +import argparse +from ogo.util.echo_arguments import echo_arguments +import ogo.util.spectral_util as md + +def print_dict(d): + keys = d.keys() + for k in keys: + print(' {:10s}: {:30s}'.format('material',d[k]['material'])) + print(' {:>30s}: {:12.4f} {}'.format('energy',d[k]['energy'],'keV')) + print(' {:>30s}: {:12.4f} {}'.format('volume_fraction',d[k]['volume_fraction'],'')) + print(' {:>30s}: {:12.4f} {}'.format('mass_concentration',d[k]['mass_concentration'],'g/cm3')) + print(' {:>30s}: {:12.4f} {}'.format('mass_attenuation',d[k]['mass_attenuation'],'cm2/g')) + print(' {:>30s}: {:12.4f} {}'.format('icru_mass_density',d[k]['icru_mass_density'],'g/cm3')) + print(' {:>30s}: {:12.4f} {}'.format('mu',d[k]['mu'],'/cm')) + print(' {:>30s}: {:12.4f} {}'.format('ctn',d[k]['ctn'],'HU')) + +def CalculateLinearAttenuation(energy, material, volume_fraction, icru_mass_density,ctn_measured, header=False, delimiter=',', description='', quiet=False): + + # Initialize variables + entry = [] + inferred_volume_fraction = False + choices = list(md.mass_density().keys()) + linear_attenuation_dict = {} + if ctn_measured is not None: + calculate_error = True + else: + calculate_error = False + + # Check inputs + if energy<30 or energy>200: + os.sys.exit('[ERROR] Energy is expected to be between 30 and 200 keV.') + if not material: + os.sys.exit('[ERROR] At least one ICRU material must be defined.') + if volume_fraction: + if sum(volume_fraction)>1.0: + os.sys.exit('[ERROR] Total volume fraction defined cannot be greater than 1.0.') + if icru_mass_density: + if len(icru_mass_density) != len(material): + os.sys.exit('[ERROR] N_material = {} and N_icru_mass_density = {}.\n'.format(len(material), len(icru_mass_density)) + + ' Every material needs an ICRU mass density defined.') + for rho in icru_mass_density: + if rho <= 0.0: + os.sys.exit('[ERROR] Invalid ICRU mass density. Must be positive: {}'.format(rho)) + if len(material) != len(volume_fraction): + if (len(material) == 1) and (not volume_fraction): + volume_fraction = [1.0] # only one material defined, so we know volume_fraction is 1.0 + elif len(material) > 1 and (len(material)-len(volume_fraction))==1: + volume_fraction.append(1.0 - sum(volume_fraction)) # one less volume_fraction defined, so we assume its value + inferred_volume_fraction = True + else: + os.sys.exit('[ERROR] N_material = {} and N_volume_fraction = {}.\n'.format(len(material),len(volume_fraction))+ + ' Every ICRU material needs a volume fraction defined.') + for idx,mat in enumerate(material): + if mat not in choices: + os.sys.exit('[ERROR] Material is not valid: {}'.format(mat)) + linear_attenuation_dict[idx] = { + 'material': mat, + 'energy': energy, + 'volume_fraction': 0.0, + 'mass_concentration': 0.0, + 'mass_attenuation': 0.0, + 'icru_mass_density': 0.0, + 'mu': 0.0, + 'ctn': 0.0 + } + + # Determine mass_concentration from known volume fractions of ICRU materials + for idx,input_volume_fraction in enumerate(volume_fraction): + mat = material[idx] + if input_volume_fraction<0: + os.sys.exit('[ERROR] Invalid volume fraction. Cannot be negative: {}'.format(input_volume_fraction)) + if icru_mass_density: + icru_density = icru_mass_density[idx] + else: + icru_density = md.mass_density()[material[idx]] + mass_concentration = (input_volume_fraction * icru_density) + linear_attenuation_dict[idx]['volume_fraction'] = input_volume_fraction + linear_attenuation_dict[idx]['icru_mass_density'] = icru_density + linear_attenuation_dict[idx]['mass_concentration'] = mass_concentration + if mass_concentration > icru_density: + os.sys.exit('[ERROR] Material = {}, ICRU density = {:.3f}, volume_fraction = {:.3f}\n'.format(mat,icru_density,input_volume_fraction) + + ' Mass concentration {} cannot be greater than its ICRU density.'.format(mass_concentration)) + + # Look up mass attenuations from NIST + for idx,mat in enumerate(material): + mass_attenuation = md.interpolate_mass_attenuation(mat,energy,method='log-log',smoothing_factor=1e-6) + linear_attenuation_dict[idx]['mass_attenuation'] = mass_attenuation + + # Calculate linear attenuation + mu_final = 0.0 + for idx,mat in enumerate(material): + mu = linear_attenuation_dict[idx]['mass_attenuation'] * linear_attenuation_dict[idx]['icru_mass_density'] * linear_attenuation_dict[idx]['volume_fraction'] + linear_attenuation_dict[idx]['mu'] = mu + mu_final += mu + + mu_water = md.interpolate_mass_attenuation('water',energy,method='log-log',smoothing_factor=1e-6) # we need this to convert to CT numbers + + # Calculate the CT number (HU) + for idx,mat in enumerate(material): + mu = linear_attenuation_dict[idx]['mu'] + ctn = round(md.mu_to_ct_number(mu,mu_water)) + linear_attenuation_dict[idx]['ctn'] = ctn + + # CT number based on total linear attenuation + ctn_final = round(md.mu_to_ct_number(mu_final,mu_water)) + + # Calculate error if a CT value is provided? + if calculate_error: + mu_measured = md.ct_number_to_mu(ctn_measured,mu_water) + mu_delta = mu_final - mu_measured + mu_percent_error = (mu_delta)/mu_measured * 100.0 + ctn_delta = ctn_final - ctn_measured + if ctn_measured != 0: + ctn_percent_error = (ctn_delta)/ctn_measured * 100.0 + + # Print results to screen + if not quiet: + print('--------------------------------------------------------------------------------') + print('Materials:') + print('--------------------------------------------------------------------------------') + print_dict(linear_attenuation_dict) + if inferred_volume_fraction: + print(' {:>44s}'.format('Note that volume fraction of this materials is inferred.')) + print('--------------------------------------------------------------------------------') + print('Summed material results:') + print('--------------------------------------------------------------------------------') + print(' {:>30s}: {:12.4f} [1]'.format('total_volume_fraction',sum(volume_fraction))) + print(' {:>30s}: {:12.4f} [/cm]'.format('mu',mu_final)) + print(' {:>30s}: {:12.4f} [HU]'.format('ctn',ctn_final)) + print(' {:>30s}: {:s}'.format('description',description)) + if calculate_error: + print('--------------------------------------------------------------------------------') + print('Error calculations:') + print('--------------------------------------------------------------------------------') + print(' {:>30s}: {:12.1f} [HU]'.format('CT number',ctn_measured)) + print(' {:>30s}: {:12.4f} [HU]'.format('delta',ctn_delta)) + print(' {:>30s}: {:12.4f} [%]'.format('error',ctn_percent_error)) + print('') + print(' {:>30s}: {:12.4f} [/cm]'.format('mu from CT number',mu_measured)) + print(' {:>30s}: {:12.4f} [/cm]'.format('delta',mu_delta)) + print(' {:>30s}: {:12.4f} [%]'.format('error',mu_percent_error)) + print('--------------------------------------------------------------------------------') + + # Primary outcomes only + entry.append( ['description', '{:s}'.format(description),'TXT'] ) + entry.append( ['energy', '{:.4f}'.format(energy), '[keV]'] ) + for idx,mat in enumerate(material): + entry.append( ['material'+'_{}'.format(idx), '{:s}'.format(linear_attenuation_dict[idx]['material']), '[]'] ) + entry.append( ['volume_fraction'+'_{}'.format(idx), '{:.4f}'.format(100.0 * linear_attenuation_dict[idx]['volume_fraction']), '[%]'] ) + entry.append( ['mass_concentration'+'_{}'.format(idx), '{:.4f}'.format(linear_attenuation_dict[idx]['mass_concentration']), '[g/cm3]'] ) + entry.append( ['mass_attenuation'+'_{}'.format(idx), '{:.4f}'.format(linear_attenuation_dict[idx]['mass_attenuation']), '[cm2/g]'] ) + entry.append( ['icru_mass_density'+'_{}'.format(idx), '{:.4f}'.format(linear_attenuation_dict[idx]['icru_mass_density']), '[g/cm3]'] ) + entry.append( ['mu'+'_{}'.format(idx), '{:.4f}'.format(linear_attenuation_dict[idx]['mu']), '[/cm]'] ) + entry.append( ['ctn'+'_{}'.format(idx), '{:.4f}'.format(linear_attenuation_dict[idx]['ctn']), '[HU]'] ) + entry.append( ['mu_final', '{:.4f}'.format(mu_final), '[/cm]'] ) + entry.append( ['ctn_final', '{:d}'.format(ctn_final), '[/cm]'] ) + if calculate_error: + entry.append( ['ctn_measured', '{:.1f}'.format(ctn_measured), '[HU]'] ) + entry.append( ['ctn_delta', '{:.4f}'.format(ctn_delta), '[H]'] ) + if ctn_measured != 0: + entry.append( ['ctn_percent_error', '{:.4f}'.format(ctn_percent_error), '[%]'] ) + else: + entry.append( ['ctn_percent_error', 'n/a', '[%]'] ) + entry.append( ['mu_measured', '{:.4f}'.format(mu_measured), '[/cm]'] ) + entry.append( ['mu_delta', '{:.4f}'.format(mu_delta), '[/cm]'] ) + entry.append( ['mu_percent_error', '{:.4f}'.format(mu_percent_error), '[%]'] ) + + # Output results for stdout + if quiet: + # Print the output + entry = list(zip(*entry)) + out = sys.stdout + if header: + out.write (delimiter.join(entry[0]) + "\n") + out.write (delimiter.join(entry[2]) + "\n") + out.write (delimiter.join(entry[1]) + "\n") + + +def main(): + # Build description including dynamic ICRU mass densities + base_description = ''' + Given the volume fractions of ICRU-defined materials and the X-ray + energy of the virtual monoenergetic image (VMI) calculate the expected + linear attenuation. + + If multiple materials are defined (e.g., HA and water) then the + calculated linear attenuation sum of the parts. + + The total volume fractions must sum to 1. For an HA insert, the volume + fraction is based on its partial density of pure HA and the ICRU total + density. Here are some examples: + + 101.7 mg HA/cm3 insert: HA volume fraction is 101.7 / 3160.0 = 0.0322 + 202.0 mg HA/cm3 insert: HA volume fraction is 202.0 / 3160.0 = 0.0639 + + Knowing the volume fraction of the inserts from the certifications, + we can calculate the volume fraction of the water materials as: + + 101.7 mg HA/cm3 insert: water volume fraction is 1 - 0.0322 = 0.9678 + 202.0 mg HA/cm3 insert: water volume fraction is 1 - 0.0639 = 0.9361 + + If the number of volume fractions is 1 less than the number of materials + then we assume that the remaining volume fraction is for the material that + was undefined. Usually this would be the water component of an HA phantom. + + If the CT number measured in Hounsfield units is provided from the + scan then the error relative to the theoretical value will be output. + + The ICRU-defined mass density can be overridden by specifying with the + ICRU argument. There must be a mass density defined for each material. +''' + materials_list = list(md.mass_density().keys()) + # Query ICRU mass densities + densities = md.mass_density() + mat_lines = ['\n Valid ICRU materials are (mass_density in g/cm3):'] + for m in materials_list: + d = densities.get(m, float('nan')) + mat_lines.append(f' {m:14s}: {d:.6f}') + mat_lines.append('\n Use quiet mode for output to STDOUT.\n') + description = base_description + '\n'.join(mat_lines) + epilog = ''' +Example calls: + ogoCalculateLinearAttenuation 140 --material cha --quiet + ogoCalculateLinearAttenuation 140 --material cha water \\ + --volume_fraction 0.202 0.972 + ogoCalculateLinearAttenuation 140 --material cha water \\ + --volume_fraction 0.202 \\ + --ctn_measured 135 \\ + -d "PCD_ESP_120kV_140keV_Br40_04Th +''' + + # Setup argument parsing + parser = argparse.ArgumentParser( + formatter_class=argparse.RawTextHelpFormatter, + prog="ogoCalculateLinearAttenuation", + description=description, + epilog=epilog + ) + parser.add_argument('energy', + type=float, + default=50.0, + metavar='Energy', + help='Energy of the scan (keV)') + parser.add_argument('--material', + nargs='*', + default=[], + metavar='MAT1 MAT2', + help='Specify ICRU materials (bone, muscle, cha, etc)') + parser.add_argument('--volume_fraction', + type=float, nargs='*', + default=[], + metavar='f1 f2', + help='Specify volume fraction for each material (sum to 1)') + parser.add_argument('--icru_mass_density', + type=float, nargs='*', + default=[], + metavar='rho1 rho2', + help='Override ICRU mass density for each material (g/cm3)') + parser.add_argument('--ctn_measured', + type=float, + default=None, + metavar='CTN', + help='Measured CT number in HU') + parser.add_argument('--header', '-H', + action="store_true", + help="""Print a header line first (default: %(default)s)""") + parser.add_argument ("--delimiter", + metavar='VAL', + default = "\t", + help="""Delimiter character (default: tab, '\\t')""") + parser.add_argument ("--description", "-d", + metavar='TXT', + default = '', + help="""Description to pass to output (e.g., filename) (default: %(default)s)""") + parser.add_argument("--quiet", + action='store_true', + help='Overwrite output without asking') + # Parse and display + args = parser.parse_args() + if not args.quiet: + print(echo_arguments('CalculateLinearAttenuation', vars(args))) + + # Run program + CalculateLinearAttenuation(**vars(args)) + +if __name__ == '__main__': + main() diff --git a/ogo/cli/ImageChangeExtent.py b/ogo/cli/ImageChangeExtent.py new file mode 100644 index 0000000..e3e3ecf --- /dev/null +++ b/ogo/cli/ImageChangeExtent.py @@ -0,0 +1,311 @@ +# /------------------------------------------------------------------------------+ +# | 11-APR-2026 | +# | Copyright (c) Bone Imaging Laboratory | +# | All rights reserved | +# | bonelab@ucalgary.ca | +# +------------------------------------------------------------------------------+ + +# Imports +import argparse +import SimpleITK as sitk +from ogo.util.echo_arguments import echo_arguments +import ogo.util.Helper as ogo +from ogo.util.cli_tools import ( + check_overwrite, + validate_integer_min, + validate_input_file, + validate_output_file, + load_config_from_json, + save_config_to_json, + print_image_info +) + +#------------------------------------------------------------------------------+ +# Main function +def ImageChangeExtent(input_file, output_file, offset, extent, filler, overwrite, quiet): + + # Check if output exists and should overwrite + if not check_overwrite(output_file, overwrite): + return 1 + + # Check that offset and extent are mutually exclusive + if offset is not None and extent is not None: + ogo.message('[ERROR] Cannot specify both --offset and --extent.') + ogo.message(' Please specify one or the other.') + return 1 + + # Check that at least one is specified + if offset is None and extent is None: + ogo.message('[ERROR] Must specify either --offset or --extent.') + return 1 + + # Read input image + if not quiet: + ogo.message(f'Reading input file: {input_file}') + + input_image = sitk.ReadImage(input_file) + + if not quiet: + print_image_info(input_file, input_image) + + # Get image properties + input_size = input_image.GetSize() + input_origin = input_image.GetOrigin() + input_spacing = input_image.GetSpacing() + pixel_type = input_image.GetPixelID() + + # Process based on offset or extent + if offset is not None: + # Validate offset dimensions + if len(offset) != 3: + ogo.message(f'[ERROR] Offset must have exactly 3 values, got {len(offset)}') + return 1 + + ox, oy, oz = offset + sx, sy, sz = input_size + + # Validate negative offsets don't equal or exceed half dimension + if ox < 0 and abs(ox) >= sx // 2: + ogo.message(f'[ERROR] Negative X offset ({ox}) equals or exceeds half the X dimension ({sx // 2})') + return 1 + if oy < 0 and abs(oy) >= sy // 2: + ogo.message(f'[ERROR] Negative Y offset ({oy}) equals or exceeds half the Y dimension ({sy // 2})') + return 1 + if oz < 0 and abs(oz) >= sz // 2: + ogo.message(f'[ERROR] Negative Z offset ({oz}) equals or exceeds half the Z dimension ({sz // 2})') + return 1 + + # Calculate new extent based on offset + # Offset is applied symmetrically on both sides of each dimension + xmin, ymin, zmin = -ox, -oy, -oz + xmax = sx + ox + ymax = sy + oy + zmax = sz + oz + + if not quiet: + if ox >= 0 and oy >= 0 and oz >= 0: + ogo.message(f'Expanding image by offset: [{ox}, {oy}, {oz}]') + elif ox <= 0 and oy <= 0 and oz <= 0: + ogo.message(f'Cropping image by offset: [{ox}, {oy}, {oz}]') + else: + ogo.message(f'Adjusting image by offset: [{ox}, {oy}, {oz}]') + ogo.message(f'New extent (voxels): [{xmin}, {xmax}, {ymin}, {ymax}, {zmin}, {zmax}]') + + else: # extent is not None + # Validate extent + if len(extent) != 6: + ogo.message(f'[ERROR] Extent must have exactly 6 values, got {len(extent)}') + return 1 + + xmin, xmax, ymin, ymax, zmin, zmax = extent + + # Validate min < max for each dimension + if xmin >= xmax: + ogo.message(f'[ERROR] xmin ({xmin}) must be less than xmax ({xmax})') + return 1 + if ymin >= ymax: + ogo.message(f'[ERROR] ymin ({ymin}) must be less than ymax ({ymax})') + return 1 + if zmin >= zmax: + ogo.message(f'[ERROR] zmin ({zmin}) must be less than zmax ({zmax})') + return 1 + + if not quiet: + ogo.message(f'Using explicit extent: [{xmin}, {xmax}, {ymin}, {ymax}, {zmin}, {zmax}]') + + # Calculate new size + new_size = [xmax - xmin, ymax - ymin, zmax - zmin] + + # Create output image filled with filler value + if not quiet: + ogo.message(f'Creating output image with size: {new_size}') + ogo.message(f'Filler value: {filler}') + + output_image = sitk.Image(new_size, pixel_type) + output_image.SetSpacing(input_spacing) + output_image.SetDirection(input_image.GetDirection()) + + # Calculate new origin + # Origin shifts when we change the extent + index_offset = [float(xmin), float(ymin), float(zmin)] + new_origin = input_image.TransformContinuousIndexToPhysicalPoint(index_offset) + output_image.SetOrigin(new_origin) + + # Fill with filler value + if filler != 0: + output_image = output_image + filler + + # Determine the region to paste from input image + # Calculate overlap between input image extent [0, sx) x [0, sy) x [0, sz) + # and new extent [xmin, xmax) x [ymin, ymax) x [zmin, zmax) + + paste_dest_start = [max(0, -xmin), max(0, -ymin), max(0, -zmin)] + crop_source_start = [max(0, xmin), max(0, ymin), max(0, zmin)] + + overlap_size = [ + min(input_size[0], xmax) - max(0, xmin), + min(input_size[1], ymax) - max(0, ymin), + min(input_size[2], zmax) - max(0, zmin) + ] + + # Check if there's any overlap + if overlap_size[0] <= 0 or overlap_size[1] <= 0 or overlap_size[2] <= 0: + ogo.message('[WARNING] No overlap between input image and new extent.') + ogo.message(' Output will contain only filler values.') + else: + # Extract the region from input image + if not quiet: + ogo.message(f'Extracting region from input: start={crop_source_start}, size={overlap_size}') + ogo.message(f'Pasting into output at: {paste_dest_start}') + + region = sitk.RegionOfInterest(input_image, overlap_size, crop_source_start) + + # Paste into output image + output_image = sitk.Paste(output_image, region, region.GetSize(), [0, 0, 0], paste_dest_start) + + # Display extent information + if not quiet: + # Input image extent + input_extent_voxels = [0, input_size[0], 0, input_size[1], 0, input_size[2]] + input_phys_min = input_image.TransformContinuousIndexToPhysicalPoint([0.0, 0.0, 0.0]) + input_phys_max = input_image.TransformContinuousIndexToPhysicalPoint([float(input_size[0]), float(input_size[1]), float(input_size[2])]) + + ogo.message(f'Input image extent:') + ogo.message(f' Voxels: [{input_extent_voxels[0]}, {input_extent_voxels[1]}, {input_extent_voxels[2]}, {input_extent_voxels[3]}, {input_extent_voxels[4]}, {input_extent_voxels[5]}]') + ogo.message(f' Physical: [{input_phys_min[0]:.2f}, {input_phys_max[0]:.2f}, {input_phys_min[1]:.2f}, {input_phys_max[1]:.2f}, {input_phys_min[2]:.2f}, {input_phys_max[2]:.2f}] mm') + + # Output image extent - use input image transformation with extent coordinates + output_phys_min = input_image.TransformContinuousIndexToPhysicalPoint([float(xmin), float(ymin), float(zmin)]) + output_phys_max = input_image.TransformContinuousIndexToPhysicalPoint([float(xmax), float(ymax), float(zmax)]) + + ogo.message(f'Output image extent:') + ogo.message(f' Voxels: [{xmin}, {xmax}, {ymin}, {ymax}, {zmin}, {zmax}]') + ogo.message(f' Physical: [{output_phys_min[0]:.2f}, {output_phys_max[0]:.2f}, {output_phys_min[1]:.2f}, {output_phys_max[1]:.2f}, {output_phys_min[2]:.2f}, {output_phys_max[2]:.2f}] mm') + + # Calculate voxel statistics + total_input_voxels = input_size[0] * input_size[1] * input_size[2] + total_output_voxels = new_size[0] * new_size[1] * new_size[2] + + if overlap_size[0] > 0 and overlap_size[1] > 0 and overlap_size[2] > 0: + pasted_voxels = overlap_size[0] * overlap_size[1] * overlap_size[2] + percent_pasted = (pasted_voxels / total_input_voxels) * 100.0 + else: + pasted_voxels = 0 + percent_pasted = 0.0 + + ogo.message(f'Voxel statistics:') + ogo.message(f' Input voxels: {total_input_voxels:,}') + ogo.message(f' Output voxels: {total_output_voxels:,}') + ogo.message(f' Pasted voxels: {pasted_voxels:,} ({percent_pasted:.1f}% of input)') + + # Write output + if not quiet: + ogo.message(f'Writing output file: {output_file}') + + sitk.WriteImage(output_image, output_file) + + if not quiet: + print_image_info(output_file, output_image) + ogo.message('[DONE]') + + return 0 + +def main(): + # Setup description + description=''' +Adjust the extent of an image by either cropping or expanding. + +This tool modifies the spatial extent of an image using either relative offsets or +absolute extent values. Offsets are applied symmetrically to both sides of each dimension. + +OFFSET MODE (--offset X Y Z): + - Positive values: Expand the image (pad with filler value) + - Negative values: Crop the image + - Values are applied symmetrically on both sides of each dimension + - Negative offsets cannot equal or exceed half the dimension size + + Example: --offset 10 10 0 adds 10 voxels to each side in X and Y dimensions + Example: --offset -5 -5 -5 removes 5 voxels from each side in all dimensions + +EXTENT MODE (--extent xmin xmax ymin ymax zmin zmax): + - Define absolute voxel coordinates for the new extent + - Can specify both expansion (negative min or max beyond image) and cropping + - Coordinates are in voxel index space + + Example: --extent -10 100 -10 100 0 50 for a 90x90x50 image starting at [-10,-10,0] + +The two modes are mutually exclusive - use one or the other, not both. +Areas outside the original image are filled with the filler value (default: 0). +''' + epilog=''' +Example calls: +ogoImageChangeExtent input.nii.gz output.nii.gz --offset 10 10 10 +ogoImageChangeExtent input.nii.gz output.nii.gz --offset -5 -5 0 +ogoImageChangeExtent input.nii.gz output.nii.gz --extent -10 200 -10 200 0 100 +ogoImageChangeExtent input.nii.gz output.nii.gz -off 20 20 20 --filler -1000 +ogoImageChangeExtent input.nii.gz output.nii.gz -ext 0 150 0 150 10 90 -ow -q +''' + + # Setup argument parsing + parser = argparse.ArgumentParser( + formatter_class=argparse.RawTextHelpFormatter, + prog="ogoImageChangeExtent", + description=description, + epilog=epilog + ) + + # JSON configuration argument + parser.add_argument('--json', nargs='?', const='', metavar='FILE', + help='Output arguments as JSON to stdout (no file), or load arguments from JSON file') + + parser.add_argument('input_file', nargs='?', type=validate_input_file(['.nii', '.nii.gz']), + help='Input image file (*.nii, *.nii.gz)') + parser.add_argument('output_file', nargs='?', type=validate_output_file(['.nii', '.nii.gz']), + help='Output image file (*.nii, *.nii.gz)') + parser.add_argument('--offset', '-off', type=int, nargs=3, default=None, metavar=('X', 'Y', 'Z'), + help='Relative offset for each dimension (3 integers); positive=expand, negative=crop') + parser.add_argument('--extent', '-ext', type=int, nargs=6, default=None, metavar=('xmin', 'xmax', 'ymin', 'ymax', 'zmin', 'zmax'), + help='Absolute extent in voxel coordinates (6 integers)') + parser.add_argument('--filler', '-f', type=int, default=0, metavar='VALUE', + help='Value to use for padding/filling (default: %(default)s)') + parser.add_argument('--overwrite', '-ow', action='store_true', help='Overwrite output without asking') + parser.add_argument('--quiet', '-q', action='store_true', help='Suppress informational output to stdout') + + print() + + # Parse and display + args = parser.parse_args() + + # Handle --json flag + if args.json is not None: + if args.json == '': + # --json with no argument: output to stdout + save_config_to_json(args, None) + return 0 + else: + # --json file.json: load from file + if not args.json.endswith('.json'): + parser.error('JSON file must have .json extension') + json_config = load_config_from_json(args.json) + # Update args with JSON values (command line args override JSON) + for key, value in json_config.items(): + if hasattr(args, key): + # Only update if argument was not explicitly provided on command line + if getattr(args, key) == parser.get_default(key): + setattr(args, key, value) + + # Check that required positional arguments are provided + if args.input_file is None or args.output_file is None: + parser.error('input_file and output_file are required') + + if not args.quiet: + ogo.message(echo_arguments('ImageChangeExtent', vars(args))) + + # Run program (exclude json argument) + run_args = vars(args).copy() + run_args.pop('json', None) + result = ImageChangeExtent(**run_args) + return result + +if __name__ == '__main__': + main() diff --git a/ogo/cli/ImageCombine.py b/ogo/cli/ImageCombine.py new file mode 100644 index 0000000..0fdc0d7 --- /dev/null +++ b/ogo/cli/ImageCombine.py @@ -0,0 +1,424 @@ +# /------------------------------------------------------------------------------+ +# | 11-APR-2026 | +# | Copyright (c) Bone Imaging Laboratory | +# | All rights reserved | +# | bonelab@ucalgary.ca | +# +------------------------------------------------------------------------------+ + +# Imports +import argparse +import SimpleITK as sitk +from ogo.util.echo_arguments import echo_arguments +import ogo.util.Helper as ogo +from ogo.util.cli_tools import ( + check_overwrite, + validate_integer_min, + validate_input_file, + validate_integer_range, + validate_output_file, + validate_extent, + load_config_from_json, + save_config_to_json, + print_image_info +) + +#------------------------------------------------------------------------------+ +# Main function +def ImageCombine(input_file1, input_file2, output_file, mode, position, side_position, spacer, filler, set_extent, overwrite, quiet): + + # Check if output exists and should overwrite + if not check_overwrite(output_file, overwrite): + return 1 + + # Read input image 1 first to get dimensions if side_position is used + if not quiet: + ogo.message(f'Reading input file 1: {input_file1}') + + input_image1 = sitk.ReadImage(input_file1) + + if not quiet: + print_image_info(input_file1, input_image1) + + # Check that position and side_position are not both set + if position != [0, 0, 0] and side_position is not None: + ogo.message('[ERROR] Cannot set both --position and --side_position. Use one or the other.') + return 1 + + # If side_position is set, calculate position from image1 dimensions + if side_position is not None: + dim = input_image1.GetSize() + if side_position == 'X': + position = [dim[0] + spacer, 0, 0] + elif side_position == 'Y': + position = [0, dim[1] + spacer, 0] + elif side_position == 'Z': + position = [0, 0, dim[2] + spacer] + if not quiet: + ogo.message(f'Side position {side_position} with spacer {spacer} set position to: {position}') + + # Display parameters + if not quiet: + ogo.message(f'Mode: {mode}') + ogo.message(f'Position: {position}') + ogo.message(f'Filler value: {filler}') + + if not quiet: + ogo.message(f'Reading input file 2: {input_file2}') + + input_image2 = sitk.ReadImage(input_file2) + + if not quiet: + print_image_info(input_file2, input_image2) + + # Check that both images have the same spacing + spacing1 = input_image1.GetSpacing() + spacing2 = input_image2.GetSpacing() + + if spacing1 != spacing2: + ogo.message('[ERROR] Image spacings do not match.') + ogo.message(f' input_file1 spacing: {spacing1}') + ogo.message(f' input_file2 spacing: {spacing2}') + return 1 + + # Determine pixel type for combined image + pixel_type1 = input_image1.GetPixelID() + pixel_type2 = input_image2.GetPixelID() + + # Check for unusual pixel types + if pixel_type1 != sitk.sitkUInt8 and pixel_type1 != sitk.sitkInt16: + ogo.message(f'[WARNING] input_file1 has unusual pixel type: {sitk.GetPixelIDValueAsString(pixel_type1)}') + + if pixel_type2 != sitk.sitkUInt8 and pixel_type2 != sitk.sitkInt16: + ogo.message(f'[WARNING] input_file2 has unusual pixel type: {sitk.GetPixelIDValueAsString(pixel_type2)}') + + # Use UInt8 if both images are UInt8, otherwise use Int16 + if pixel_type1 == sitk.sitkUInt8 and pixel_type2 == sitk.sitkUInt8: + combined_pixel_type = sitk.sitkUInt8 + else: + combined_pixel_type = sitk.sitkInt16 + + if not quiet: + ogo.message(f'Combined image pixel type: {sitk.GetPixelIDValueAsString(combined_pixel_type)}') + + # Validate filler value for UInt8 and Int8 output + if combined_pixel_type == sitk.sitkUInt8: + if filler < 0: + ogo.message(f'[WARNING] Filler value {filler} is out of range for UInt8 (0-255). Setting filler to 0.') + filler = 0 + elif filler > 255: + ogo.message(f'[WARNING] Filler value {filler} is out of range for UInt8 (0-255). Setting filler to 255.') + filler = 255 + elif combined_pixel_type == sitk.sitkInt8: + if filler < -128: + ogo.message(f'[WARNING] Filler value {filler} is out of range for Int8 (-128-127). Setting filler to -128.') + filler = -128 + elif filler > 127: + ogo.message(f'[WARNING] Filler value {filler} is out of range for Int8 (-128-127). Setting filler to 127.') + filler = 127 + + # Get dimensions of both images + dim1 = input_image1.GetSize() + dim2 = input_image2.GetSize() + + # Calculate extent of combined image + if set_extent is not None: + # Use user-provided extent + minX, maxX, minY, maxY, minZ, maxZ = set_extent + if not quiet: + ogo.message(f'Using user-provided extent: X[{minX}:{maxX}], Y[{minY}:{maxY}], Z[{minZ}:{maxZ}]') + + # Validate that both images fit within the user-provided extent + # Image1 is at [0, 0, 0] + if 0 < minX or dim1[0] > maxX: + ogo.message(f'[ERROR] Image1 (X extent [0:{dim1[0]}]) does not fit within specified extent X[{minX}:{maxX}]') + return 1 + if 0 < minY or dim1[1] > maxY: + ogo.message(f'[ERROR] Image1 (Y extent [0:{dim1[1]}]) does not fit within specified extent Y[{minY}:{maxY}]') + return 1 + if 0 < minZ or dim1[2] > maxZ: + ogo.message(f'[ERROR] Image1 (Z extent [0:{dim1[2]}]) does not fit within specified extent Z[{minZ}:{maxZ}]') + return 1 + + # Image2 is at position + if position[0] < minX or (position[0] + dim2[0]) > maxX: + ogo.message(f'[ERROR] Image2 (X extent [{position[0]}:{position[0] + dim2[0]}]) does not fit within specified extent X[{minX}:{maxX}]') + return 1 + if position[1] < minY or (position[1] + dim2[1]) > maxY: + ogo.message(f'[ERROR] Image2 (Y extent [{position[1]}:{position[1] + dim2[1]}]) does not fit within specified extent Y[{minY}:{maxY}]') + return 1 + if position[2] < minZ or (position[2] + dim2[2]) > maxZ: + ogo.message(f'[ERROR] Image2 (Z extent [{position[2]}:{position[2] + dim2[2]}]) does not fit within specified extent Z[{minZ}:{maxZ}]') + return 1 + else: + # Calculate extent automatically based on image positions + # Image1 is at [0, 0, 0], Image2 is at position + minX = min(0, position[0]) + maxX = max(dim1[0], position[0] + dim2[0]) + minY = min(0, position[1]) + maxY = max(dim1[1], position[1] + dim2[1]) + minZ = min(0, position[2]) + maxZ = max(dim1[2], position[2] + dim2[2]) + + # Dimensions of combined image + combined_dim = [maxX - minX, maxY - minY, maxZ - minZ] + + if not quiet: + ogo.message(f'Combined image dimensions: {combined_dim}') + ogo.message(f'Combined image extent (voxels): X[{minX}:{maxX}], Y[{minY}:{maxY}], Z[{minZ}:{maxZ}]') + + # Calculate physical extent using image1's coordinate system + phys_min = input_image1.TransformContinuousIndexToPhysicalPoint([float(minX), float(minY), float(minZ)]) + phys_max = input_image1.TransformContinuousIndexToPhysicalPoint([float(maxX), float(maxY), float(maxZ)]) + ogo.message(f'Combined image extent (physical): X[{phys_min[0]:.2f}:{phys_max[0]:.2f}], Y[{phys_min[1]:.2f}:{phys_max[1]:.2f}], Z[{phys_min[2]:.2f}:{phys_max[2]:.2f}] mm') + + # Create combined image initialized with filler value + combined_image = sitk.Image(combined_dim, combined_pixel_type) + combined_image.SetSpacing(input_image1.GetSpacing()) + + # Adjust origin to account for extent that may extend beyond image1 + index_offset = [float(minX), float(minY), float(minZ)] + new_origin = input_image1.TransformContinuousIndexToPhysicalPoint(index_offset) + combined_image.SetOrigin(new_origin) + + combined_image.SetDirection(input_image1.GetDirection()) + + # Initialize with filler value + combined_image = combined_image + filler + + if not quiet: + ogo.message(f'Created combined image initialized with filler value: {filler}') + + # Paste image1 into combined_image + # Image1 is at [0, 0, 0] in original coordinate system + # In combined image, it's at position [0-minX, 0-minY, 0-minZ] + paste_position1 = [0 - minX, 0 - minY, 0 - minZ] + + if not quiet: + ogo.message(f'Pasting image1 at position: {paste_position1}') + + combined_image = sitk.Paste( + combined_image, + input_image1, + input_image1.GetSize(), + [0, 0, 0], + paste_position1 + ) + + # Paste image2 into combined_image based on mode + # Image2 is at position in original coordinate system + # In combined image, it's at position [position[0]-minX, position[1]-minY, position[2]-minZ] + paste_position2 = [position[0] - minX, position[1] - minY, position[2] - minZ] + + if not quiet: + ogo.message(f'Pasting image2 at position: {paste_position2}') + + if mode == 'overlay': + # Overlay mode: Paste image2, overwriting existing values except where image2 is 0 + # Extract the region where image2 will go + region = sitk.RegionOfInterest( + combined_image, + input_image2.GetSize(), + paste_position2 + ) + # Cast image2 to match combined_image pixel type + image2_casted = sitk.Cast(input_image2, combined_pixel_type) + # Copy physical space information from region to match + image2_casted.CopyInformation(region) + + # Create mask where image2 is non-zero + mask = sitk.Cast(image2_casted != 0, combined_pixel_type) + + # Combine: use image2 where mask is non-zero, otherwise keep region + result = mask * image2_casted + (1 - mask) * region + + # Paste the result back + combined_image = sitk.Paste( + combined_image, + result, + result.GetSize(), + [0, 0, 0], + paste_position2 + ) + elif mode == 'add': + # Add mode: Add image2 values to existing combined_image values + # Extract the region where image2 will go + region = sitk.RegionOfInterest( + combined_image, + input_image2.GetSize(), + paste_position2 + ) + # Cast image2 to match combined_image pixel type for arithmetic operation + image2_casted = sitk.Cast(input_image2, combined_pixel_type) + # Copy physical space information from region to match for arithmetic operation + image2_casted.CopyInformation(region) + # Add image2 to this region + region = sitk.Add(region, image2_casted) + # Paste the result back + combined_image = sitk.Paste( + combined_image, + region, + region.GetSize(), + [0, 0, 0], + paste_position2 + ) + elif mode == 'subtract': + # Subtract mode: Subtract image2 values from existing combined_image values + # Extract the region where image2 will go + region = sitk.RegionOfInterest( + combined_image, + input_image2.GetSize(), + paste_position2 + ) + # Cast image2 to match combined_image pixel type for arithmetic operation + image2_casted = sitk.Cast(input_image2, combined_pixel_type) + # Copy physical space information from region to match for arithmetic operation + image2_casted.CopyInformation(region) + # Subtract image2 from this region + region = sitk.Subtract(region, image2_casted) + # Paste the result back + combined_image = sitk.Paste( + combined_image, + region, + region.GetSize(), + [0, 0, 0], + paste_position2 + ) + + if not quiet: + ogo.message(f'Writing output file: {output_file}') + + sitk.WriteImage(combined_image, output_file) + + if not quiet: + print_image_info(output_file, combined_image) + + # Display final extent of output image in both voxels and physical coordinates + ogo.message(f'Output image extent (voxels): X[{minX}:{maxX}], Y[{minY}:{maxY}], Z[{minZ}:{maxZ}]') + + output_origin = combined_image.GetOrigin() + output_size = combined_image.GetSize() + output_spacing = combined_image.GetSpacing() + output_phys_max = combined_image.TransformContinuousIndexToPhysicalPoint( + [float(output_size[0]), float(output_size[1]), float(output_size[2])] + ) + ogo.message(f'Output image extent (physical): X[{output_origin[0]:.2f}:{output_phys_max[0]:.2f}], Y[{output_origin[1]:.2f}:{output_phys_max[1]:.2f}], Z[{output_origin[2]:.2f}:{output_phys_max[2]:.2f}] mm') + + if not quiet: + ogo.message('[DONE]') + + return 0 + +def main(): + # Setup description + description=''' +Utility to combine two NIFTI image files with flexible positioning and blending modes. +This performs the combination in voxel units, not physical space (i.e., Origin). + +This tool takes two input NIFTI images and combines them according to the specified mode. +It is not required that the two input images have compatible dimensions. The output image +will have dimensions that cover the extent of the two input images. The filler variable +will be used in regions where neither image provides data within the extent. + +The mode determines how overlapping voxels are handled: + - add: Add voxel values where images overlap (image 1 + image 2) + - subtract: Subtract voxel values where images overlap (image 1 - image 2) + - overlay: Prioritize image 2 voxels (overwrite image 1 voxels) + +The position parameter defines the relative position of input_file2 with respect to +input_file1 in voxel coordinates. Alternatively, use --side_position for convenient +positioning along dimension boundaries of input_file1 (X, Y, or Z axis). The --spacer +parameter can be used with --side_position to add additional spacing (in voxels) between +the images. + +Optionally, use --set_extent to manually specify the extent of the output combined image. +''' + epilog=''' +Example calls: +ogoImageCombine input1.nii.gz input2.nii.gz output.nii.gz +ogoImageCombine input1.nii.gz input2.nii.gz output.nii.gz --mode overlay --position 10 20 30 +ogoImageCombine input1.nii.gz input2.nii.gz output.nii.gz -m add --filler 255 +ogoImageCombine input1.nii.gz input2.nii.gz output.nii.gz -m subtract -p -5 0 10 -ow -q +ogoImageCombine input1.nii.gz input2.nii.gz output.nii.gz --side_position X +ogoImageCombine input1.nii.gz input2.nii.gz output.nii.gz -s Z -m overlay -ow +ogoImageCombine input1.nii.gz input2.nii.gz output.nii.gz -s X --spacer 10 +ogoImageCombine input1.nii.gz input2.nii.gz output.nii.gz -e -10 100 0 200 0 150 -m add +''' + + # Setup argument parsing + parser = argparse.ArgumentParser( + formatter_class=argparse.RawTextHelpFormatter, + prog="ogoImageCombine", + description=description, + epilog=epilog + ) + + # JSON configuration argument + parser.add_argument('--json', nargs='?', const='', metavar='FILE', + help='Output arguments as JSON to stdout (no file), or load arguments from JSON file') + + parser.add_argument('input_file1', nargs='?', type=validate_input_file(['.nii', '.nii.gz']), + help='First input image file (*.nii, *.nii.gz)') + parser.add_argument('input_file2', nargs='?', type=validate_input_file(['.nii', '.nii.gz']), + help='Second input image file (*.nii, *.nii.gz)') + parser.add_argument('output_file', nargs='?', type=validate_output_file(['.nii', '.nii.gz']), + help='Output combined image file (*.nii, *.nii.gz)') + parser.add_argument('--mode', '-m', default='overlay', choices=['add', 'subtract', 'overlay'], + help='Mode for combining overlapping voxels (default: %(default)s)') + parser.add_argument('--position', '-p', type=int, nargs=3, default=[0, 0, 0], metavar='INT', + help='Relative position of input_file2 [x, y, z] in voxels (default: 0 0 0)') + parser.add_argument('--side_position', '-s', default=None, choices=['X', 'Y', 'Z'], + help='Position input_file2 at dimension boundary (X, Y, or Z) of input_file1. Alternative to --position') + parser.add_argument('--spacer', '-sp', type=validate_integer_min(0), default=0, metavar='INT', + help='Additional spacing (in voxels) when using --side_position (must be >= 0, default: %(default)s)') + parser.add_argument('--filler', '-f', type=validate_integer_range(-32768, 32767), default=0, metavar='VALUE', + help='Filler value for regions without data (default: %(default)s)') + parser.add_argument('--set_extent', '-e', type=int, nargs=6, default=None, metavar=('xmin', 'xmax', 'ymin', 'ymax', 'zmin', 'zmax'), + help='Manually set extent of combined image [xmin xmax ymin ymax zmin zmax] in voxels') + parser.add_argument('--overwrite','-ow', action='store_true', help='Overwrite output without asking') + parser.add_argument('--quiet','-q', action='store_true', help='Suppress informational output to stdout') + + print() + + # Parse and display + args = parser.parse_args() + + # Handle --json flag + if args.json is not None: + if args.json == '': + # --json with no argument: output to stdout + save_config_to_json(args, None) + return 0 + else: + # --json file.json: load from file + if not args.json.endswith('.json'): + parser.error('JSON file must have .json extension') + json_config = load_config_from_json(args.json) + # Update args with JSON values (command line args override JSON) + for key, value in json_config.items(): + if hasattr(args, key): + # Only update if argument was not explicitly provided on command line + if getattr(args, key) == parser.get_default(key): + setattr(args, key, value) + + # Check that required positional arguments are provided + if args.input_file1 is None or args.input_file2 is None or args.output_file is None: + parser.error('input_file1, input_file2, and output_file are required') + + # Validate set_extent if provided + if args.set_extent is not None: + try: + args.set_extent = validate_extent(args.set_extent) + except argparse.ArgumentTypeError as e: + parser.error(str(e)) + + if not args.quiet: + ogo.message(echo_arguments('ImageCombine', vars(args))) + + # Run program (exclude json argument) + run_args = vars(args).copy() + run_args.pop('json', None) + result = ImageCombine(**run_args) + return result + +if __name__ == '__main__': + main() diff --git a/ogo/cli/ImageCrop.py b/ogo/cli/ImageCrop.py index 3cd627d..82f31db 100644 --- a/ogo/cli/ImageCrop.py +++ b/ogo/cli/ImageCrop.py @@ -14,82 +14,66 @@ from ogo.util.echo_arguments import echo_arguments import ogo.util.Helper as ogo import ogo.dat.OgoMasterLabels as lb +from ogo.util.cli_tools import ( + check_overwrite, + check_image_dimensions_match, + validate_integer_range, + validate_input_file, + validate_output_file, + print_image_info +) -def printImageInfo(im): - dim = im.GetSize() - spacing = im.GetSpacing() - origin = im.GetOrigin() - phys_dim = [x * y for x, y in zip(dim, spacing)] - position = [math.floor(x / y) for x, y in zip(origin, spacing)] +def ImageCrop(input_file, output_file, input_label_file, output_label_file, voi, labels, slab, offset, overwrite, quiet): - guard = '!-------------------------------------------------------------------------------' - print(guard) - print('!> dim {:>8} {:>8} {:>8}'.format(*dim)) - print('!> off {:>8} {:>8} {:>8}'.format('-', '-', '-')) - print('!> pos {:>8} {:>8} {:>8}'.format(*position)) - print('!> element size in mm {:>8.4f} {:>8.4f} {:>8.4f}'.format(*spacing)) - print('!> phys dim in mm {:>8.4f} {:>8.4f} {:>8.4f}'.format(*phys_dim)) - print(guard) - -def ImageCrop(input_image, output_image, output_label, label_image, voi, labels, slab, overwrite): - - # Define output_image if not set explicitly - if output_image is None: - basename = os.path.basename(input_image) - name, ext = os.path.splitext(input_image) - if 'gz' in ext: - name = os.path.splitext(name)[0] # Manages files with double extension - ext = '.nii' + ext - output_image = name + '_CROP.nii.gz' - ogo.message('Output image will be:') - ogo.message(' \"{}\"'.format(output_image)) - - # Define output_label if not set explicitly - if output_label is None and label_image is not None: - basename = os.path.basename(label_image) - name, ext = os.path.splitext(label_image) - if 'gz' in ext: - name = os.path.splitext(name)[0] # Manages files with double extension - ext = '.nii' + ext - output_label = name + '_CROP.nii.gz' - if label_image is not None: + if not quiet: + ogo.message('Output image will be:') + ogo.message(' "{}"'.format(output_file)) + + if output_label_file is not None and not quiet: ogo.message('Output label will be:') - ogo.message(' \"{}\"'.format(output_label)) + ogo.message(' "{}"'.format(output_label_file)) - # Check if output exists and should overwrite - if os.path.isfile(output_image) and not overwrite: - result = input('File \"{}\" already exists. Overwrite? [y/n]: '.format(output_image)) - if result.lower() not in ['y', 'yes']: - ogo.message('Not overwriting. Exiting...') - os.sys.exit() + # Check if outputs exist and should overwrite + # Collect all output files that need to be checked + output_files = [output_file] + if output_label_file is not None: + output_files.append(output_label_file) + + # Check all output files together - only prompt once if needed + if not overwrite: + existing_files = [f for f in output_files if os.path.isfile(f)] + if existing_files: + files_list = '\n '.join(['"{}"'.format(f) for f in existing_files]) + result = input(f'The following file(s) already exist:\n {files_list}\n Overwrite? [y/n]: ') + if result.lower() not in ['y', 'yes']: + print('Not overwriting. Exiting...') + return 1 # Read input - if not os.path.isfile(input_image): - os.sys.exit('[ERROR] Cannot find file \"{}\"'.format(input_image)) - - ogo.message('Reading input CT image to be cropped:') - ogo.message(' \"{}\"'.format(input_image)) - ct = sitk.ReadImage(input_image) - printImageInfo(ct) + if not quiet: + ogo.message('Reading input CT image to be cropped:') + ogo.message(' "{}"'.format(input_file)) + ct = sitk.ReadImage(input_file) + if not quiet: + print_image_info(input_file, ct, brief=True) # If using labels to define bounds - if label_image is not None: - ogo.message('Reading label image') - ogo.message(' \"{}\"'.format(label_image)) - if os.path.isfile(label_image): - ct_labels = sitk.ReadImage(label_image, sitk.sitkUInt8) - printImageInfo(ct_labels) + if input_label_file is not None: + if not quiet: + ogo.message('Reading label image') + ogo.message(' "{}"'.format(input_label_file)) + if os.path.isfile(input_label_file): + ct_labels = sitk.ReadImage(input_label_file, sitk.sitkUInt8) + if not quiet: + print_image_info(input_label_file, ct_labels, brief=True) # check dimensions of label and ct images are the same - if not ((ct.GetSize()[0] == ct_labels.GetSize()[0]) and \ - (ct.GetSize()[1] == ct_labels.GetSize()[1]) and \ - (ct.GetSize()[2] == ct_labels.GetSize()[2])): - ogo.message('[ERROR] Label image must be same dimensions as input image.') - os.sys.exit() + if not check_image_dimensions_match(ct, ct_labels, "input_file", "input_label_file"): + return 1 if not labels: - ogo.message('Labels must be defined to crop using labels.') - os.sys.exit() + ogo.message('[ERROR] Labels must be defined to crop using labels.') + return 1 # create mask of labels ct_mask = ct_labels==-1 # we want these to be all zeros @@ -103,67 +87,153 @@ def ImageCrop(input_image, output_image, output_label, label_image, voi, labels, desc = 'unknown label' bounding_box = filt.GetBoundingBox(label) - ogo.message(' label {:2d} ({:s}):'.format(label,desc)) + if not quiet: + ogo.message(' label {:2d} ({:s}):'.format(label,desc)) if label in labels: - ogo.message(' [' + ', '.join('{:3d}'.format(i) for i in bounding_box) + '] - used to define VOI') + if not quiet: + ogo.message(' [' + ', '.join('{:3d}'.format(i) for i in bounding_box) + '] - used to define VOI') bin_part = ct_labels==label mask = 1 - bin_part ct_mask = sitk.Mask(ct_mask, mask) ct_mask = ct_mask + bin_part else: - ogo.message(' [' + ', '.join('{:3d}'.format(i) for i in bounding_box) + ']') + if not quiet: + ogo.message(' [' + ', '.join('{:3d}'.format(i) for i in bounding_box) + ']') filt.Execute(ct_mask) - voi = filt.GetBoundingBox(1) - ogo.message('') - ogo.message(' VOI based on labels:') - ogo.message(' [' + ', '.join('{:3d}'.format(i) for i in voi) + ']') + voi = list(filt.GetBoundingBox(1)) # Convert tuple to list + if not quiet: + ogo.message('') + ogo.message(' VOI based on labels:') + ogo.message(' [' + ', '.join('{:3d}'.format(i) for i in voi) + ']') + + # Apply offset if specified + if offset is not None: + dim = ct.GetSize() + # Calculate new VOI bounds with offset, clamped to image dimensions + requested_start_x = voi[0] - offset[0] + requested_end_x = voi[0] + voi[3] + offset[0] + requested_start_y = voi[1] - offset[1] + requested_end_y = voi[1] + voi[4] + offset[1] + requested_start_z = voi[2] - offset[2] + requested_end_z = voi[2] + voi[5] + offset[2] + + new_start_x = max(0, requested_start_x) + new_end_x = min(dim[0], requested_end_x) + new_start_y = max(0, requested_start_y) + new_end_y = min(dim[1], requested_end_y) + new_start_z = max(0, requested_start_z) + new_end_z = min(dim[2], requested_end_z) + + # Check if clamping occurred and warn the user + clamped_dims = [] + if requested_start_x < 0 or requested_end_x > dim[0]: + clamped_dims.append('X') + if requested_start_y < 0 or requested_end_y > dim[1]: + clamped_dims.append('Y') + if requested_start_z < 0 or requested_end_z > dim[2]: + clamped_dims.append('Z') + + if clamped_dims and not quiet: + ogo.message('[WARNING] Offset expansion exceeded image boundaries in dimension(s): {}'.format(', '.join(clamped_dims))) + ogo.message('[WARNING] VOI has been clamped to fit within image bounds [0-{}, 0-{}, 0-{}]'.format(dim[0]-1, dim[1]-1, dim[2]-1)) + + # Update VOI with new values + voi[0] = new_start_x + voi[1] = new_start_y + voi[2] = new_start_z + voi[3] = new_end_x - new_start_x + voi[4] = new_end_y - new_start_y + voi[5] = new_end_z - new_start_z + + if not quiet: + ogo.message('') + ogo.message(' VOI after applying offset [{}, {}, {}]:'.format(*offset)) + ogo.message(' [' + ', '.join('{:3d}'.format(i) for i in voi) + ']') # Reset VOI if slab dim = ct.GetSize() if slab: - ogo.message('Using Z bounds only.') + if not quiet: + ogo.message('Using Z bounds only.') voi = [0, 0, voi[2], dim[0], dim[1], voi[5]] - ogo.message('') - ogo.message(' Final VOI:') - ogo.message(' [' + ', '.join('{:3d}'.format(i) for i in voi) + ']') + if not quiet: + ogo.message('') + ogo.message(' Final VOI:') + ogo.message(' [' + ', '.join('{:3d}'.format(i) for i in voi) + ']') else: - ogo.message('') - ogo.message(' Final VOI:') - ogo.message(' [' + ', '.join('{:3d}'.format(i) for i in voi) + ']') + if not quiet: + ogo.message('') + ogo.message(' Final VOI:') + ogo.message(' [' + ', '.join('{:3d}'.format(i) for i in voi) + ']') # Take the subvolume if voi[0]<0 or voi[1]<0 or voi[2]<0 or voi[3]>dim[0] or voi[4]>dim[1] or voi[5]>dim[2]: - os.sys.exit('[ERROR] Defined VOI is out of range: {} {} {} {} {} {}'.format(*voi)) + ogo.message('[ERROR] Defined VOI is out of range: {} {} {} {} {} {}'.format(*voi)) + return 1 ct_out = ct[voi[0]:(voi[0]+voi[3]), voi[1]:(voi[1]+voi[4]), voi[2]:(voi[2]+voi[5])] - ogo.message('Writing output image to file {}'.format(output_image)) - sitk.WriteImage(ct_out, output_image) - printImageInfo(ct_out) + if not quiet: + ogo.message('Writing output image to file {}'.format(output_file)) + sitk.WriteImage(ct_out, output_file) + if not quiet: + print_image_info(output_file, ct_out) - if label_image is not None: + if input_label_file is not None and output_label_file is not None: ct_label_out = ct_labels[voi[0]:(voi[0]+voi[3]), voi[1]:(voi[1]+voi[4]), voi[2]:(voi[2]+voi[5])] - ogo.message('Writing output labels to file {}'.format(output_label)) - sitk.WriteImage(ct_label_out, output_label) - printImageInfo(ct_label_out) - - ogo.message('Done ogoImageCrop!') + if not quiet: + ogo.message('Writing output labels to file {}'.format(output_label_file)) + sitk.WriteImage(ct_label_out, output_label_file) + if not quiet: + print_image_info(output_label_file, ct_label_out, brief=True) + + if not quiet: + ogo.message('Done ogoImageCrop!') + + return 0 def main(): # Setup description description = ''' -Utility to crop a NIFTI file to a subvolume. +Utility to crop a NIfTI file to a subvolume. + +There are two main methods of cropping. One is to define the volume of interest (VOI) +and the other is to use the labels in a mask of the image to be cropped. + +Define the VOI using six integers:origX, origY, origZ, dimX, dimY, dimZ +For example: --voi 50 100 150 200 250 300 +This crops from voxel (50,100,150) with dimensions (200,250,300). + +Define an input file with labels (usually from ML or similar) and specify which labels +to use as a basis for cropping. It will automatically find the VOI that includes all +the specified labels. + +For example: --input_label_file labels.nii.gz --labels 5 10 15 +This finds the smallest VOI containing labels 5, 10, and 15, then crops. + +Optional flags: + +--slab: When using label-based cropping, this option keeps the full X and Y + dimensions and only crops in the Z direction based on the label bounds. + Useful for extracting vertebral slabs or similar structures. + +--offset: When using label-based cropping, this option adds a symmetric offset to the + VOI defined by the labels. For example, --offset 10 10 5 adds 10 voxels on + each side in X (20 total), 10 on each side in Y (20 total), and 5 on each + side in Z (10 total). -Cropping can be performed based on the labels within an image. A list -of labels can be defined. +Outputs are the cropped image and optionally the cropped label file if an input label +file was used. The output label file will contain the same labels as the input label +file, but only within the cropped region. ''' epilog = ''' Example calls: ogoImageCrop input.nii.gz output.nii.gz --voi 0 512 0 512 137 448 -ogoImageCrop input.nii.gz output.nii.gz --label_image labels.nii.gz --labels 8 9 10 --slab +ogoImageCrop input.nii.gz output.nii.gz --input_label_file labels.nii.gz -l 8 9 10 --slab +ogoImageCrop input.nii.gz output_cropped.nii.gz --input_label_file labels.nii.gz -l 5 10 15 --offset 10 10 5 -ow -q ''' # Setup argument parsing @@ -173,21 +243,25 @@ def main(): description=description, epilog=epilog ) - parser.add_argument('input_image', help='Input image file (*.nii, *.nii.gz)') - parser.add_argument('--output_image', metavar='FILE', default=None, help='Output CT image file (*.nii, *.nii.gz)') - parser.add_argument('--output_label', metavar='FILE', default=None, help='Output label file (*.nii, *.nii.gz)') - parser.add_argument('--label_image', metavar='FILE', default=None, help='Label image used to define crop (*.nii, *.nii.gz)') - parser.add_argument('--voi', type=int, nargs=6, default=[0,1,0,1,0,1], metavar='0', help='[origX, origY, origZ, dimX, dimY, dimZ') - parser.add_argument('--labels', type=int, nargs='*', default=[8], metavar='LABEL', help='List of labels to define crop region (default: %(default)s)') + parser.add_argument('input_file', type=validate_input_file(['.nii', '.nii.gz']), help='Input image file (*.nii, *.nii.gz)') + parser.add_argument('output_file', type=validate_output_file(['.nii', '.nii.gz']), help='Output CT image file (*.nii, *.nii.gz)') + parser.add_argument('--input_label_file','-ilf', type=validate_input_file(['.nii', '.nii.gz']), metavar='FILE', default=None, help='Label image used to define crop (*.nii, *.nii.gz)') + parser.add_argument('--output_label_file', '-olf', type=validate_output_file(['.nii', '.nii.gz']), metavar='FILE', default=None, help='Output label file (*.nii, *.nii.gz)') + parser.add_argument('--voi', type=int, nargs=6, default=[0,1,0,1,0,1], metavar='0', help='[origX, origY, origZ, dimX, dimY, dimZ]') + parser.add_argument('--labels', '-l', type=validate_integer_range(0, 255), nargs='*', default=[], metavar='LABEL', help='List of labels to define crop region (default: %(default)s)') + parser.add_argument('--offset','-os', type=int, nargs=3, default=None, metavar=('X', 'Y', 'Z'), help='Symmetric offset to expand VOI in label-based cropping [offsetX, offsetY, offsetZ]') parser.add_argument('--slab', action='store_true', help='Crop a Z slab (overrides X and Y VOI)') - parser.add_argument('--overwrite', action='store_true', help='Overwrite output without asking') + parser.add_argument('--overwrite', '-ow', action='store_true', help='Overwrite output without asking') + parser.add_argument('--quiet', '-q', action='store_true', help='Suppress informational output to stdout') # Parse and display args = parser.parse_args() - print(echo_arguments('ImageCrop', vars(args))) + if not args.quiet: + ogo.message(echo_arguments('ImageCrop', vars(args))) # Run program - ImageCrop(**vars(args)) + result = ImageCrop(**vars(args)) + return result if __name__ == '__main__': diff --git a/ogo/cli/ImageExam.py b/ogo/cli/ImageExam.py index 08378bd..a18bf65 100644 --- a/ogo/cli/ImageExam.py +++ b/ogo/cli/ImageExam.py @@ -8,58 +8,117 @@ # Imports import argparse import os -import vtk import numpy as np -import math -from vtk.util.numpy_support import vtk_to_numpy from ogo.util.echo_arguments import echo_arguments import ogo.util.Helper as ogo import SimpleITK as sitk +from ogo.util.cli_tools import ( + validate_integer_min, + validate_input_file, + print_image_info +) -def ImageExam(input_filename, header, histogram, bins): +def ImageExam(input_file, header, histogram, bins, quiet): # Check valid number of histogram bins if bins < 1: - os.sys.exit('[ERROR] Number of histogram bins must be greater than zero: bins = {}'.format(bins)) + ogo.message('[ERROR] Number of histogram bins must be greater than zero: bins = {}'.format(bins)) + return 1 - # ogo.message('Start ogoImageExam!') - - # Read input - if not os.path.isfile(input_filename): - os.sys.exit('[ERROR] Cannot find file \"{}\"'.format(input_filename)) - - if input_filename.lower().endswith('.nii'): - reader = vtk.vtkNIFTIImageReader() - elif input_filename.lower().endswith('.nii.gz'): - reader = vtk.vtkNIFTIImageReader() - else: - os.sys.exit('[ERROR] Cannot find reader for file \"{}\"'.format(input_filename)) - - print() - ogo.message('Reading input image ' + input_filename) - reader.SetFileName(input_filename) - reader.Update() - - # ogo.message('Image dimensions and type') - ogo.aix(input_filename, reader.GetOutput()) + # Read input using SimpleITK + if not quiet: + ogo.message('Reading input image ' + input_file) + + image = sitk.ReadImage(input_file) + + if not quiet: + print_image_info(input_file, image) if header: - ogo.message('Header from NIFTI file (VTK reader)') - ogo.infoNIFTI(reader) + if not quiet: + ogo.message('Header from NIFTI file (SimpleITK reader)') - ogo.message('Header from NIFTI file –-> Using SimpleITK reader') - ct = sitk.ReadImage(input_filename) - for k in ct.GetMetaDataKeys(): - v = ct.GetMetaData(k) - print(' {:20s} = {}'.format(k,v)) + for k in image.GetMetaDataKeys(): + v = image.GetMetaData(k) + if not quiet: + print(' {:20s} = {}'.format(k, v)) if histogram: - ogo.message('Histogram of NIFTI image data') - ogo.histogram(reader.GetOutput(), bins) - - + if not quiet: + ogo.message('Histogram of NIFTI image data') + + # Get numpy array from SimpleITK image + array = sitk.GetArrayFromImage(image).ravel() + + # Check if data type is integer + is_integer_type = np.issubdtype(array.dtype, np.integer) + + # Compute histogram with integer-aligned bins for integer data + if is_integer_type: + # For integer data, align bin edges to integer values + # Use strict [a, b) intervals by extending beyond max + data_min = int(np.min(array)) + data_max = int(np.max(array)) + + # Create integer bin edges from min to max+1 + # The last edge at max+1 ensures that max value falls in [max, max+1) + # This enforces strict [a, b) for all bins (no special case for last bin) + bin_edges = np.linspace(data_min, data_max + 1, bins + 1, dtype=float) + bin_edges = np.round(bin_edges).astype(int) + + # Ensure we don't have duplicate edges + bin_edges = np.unique(bin_edges) + if len(bin_edges) < bins + 1: + # If unique edges are fewer than requested bins, use what we have + if not quiet: + ogo.message(f'[INFO] Adjusted to {len(bin_edges)-1} bins for integer alignment') + + counts, bin_edges = np.histogram(array, bins=bin_edges) + else: + # For floating-point data, use default binning + counts, bin_edges = np.histogram(array, bins=bins) + + # Display histogram + guard = '!-------------------------------------------------------------------------------' + + # Print histogram bins differently for integer vs float data + actual_bins = len(counts) + + if is_integer_type: + # For integer data, show bin range instead of center + # All bins use [min, max) notation (right edge excluded) + line_format_header = '!> {:>18s} {:>8s} {:>12s}' + line_format_data = '!> [{:>7.0f}, {:>7.0f}) {:>8d} {:>12d}' + line_format_summary = '!> {:>18.0f} {:>8s} {:>12s}' + + print(guard) + print(line_format_header.format('Range [min, max)', 'Bin', 'Count')) + print(guard) + + for i in range(actual_bins): + print(line_format_data.format(bin_edges[i], bin_edges[i + 1], i, counts[i])) + else: + # For float data, show bin center + line_format_header = '!> {:>8s} {:>8s} {:>12s}' + line_format_data = '!> {:>8.3f} {:>8d} {:>12d}' + line_format_summary = '!> {:>8.3f} {:>8s} {:>12s}' + + print(guard) + print(line_format_header.format('Value', 'Bin', 'Count')) + print(guard) + + for i in range(actual_bins): + bin_center = (bin_edges[i] + bin_edges[i + 1]) / 2.0 + print(line_format_data.format(bin_center, i, counts[i])) + + print(guard) + print(line_format_summary.format(np.min(array), 'Min', '')) + print(line_format_summary.format(np.max(array), 'Max', '')) + print(line_format_summary.format(np.mean(array), 'Mean', '')) + print(line_format_summary.format(np.std(array), 'Std', '')) + print(guard) - # ogo.message('Done ogoImageExam!') + return 0 def main(): @@ -71,7 +130,8 @@ def main(): epilog = ''' Example calls: ogoImageExam input.nii.gz -ogoImageExam --histogram input.nii.gz +ogoImageExam input.nii.gz --histogram +ogoImageExam input.nii.gz -H --histogram -b 256 -q ''' # Setup argument parsing @@ -81,18 +141,21 @@ def main(): description=description, epilog=epilog ) - parser.add_argument('input_filename', help='Input image file (*.nii, *.nii.gz)') - parser.add_argument('--header', action='store_true', help='Show NIFTI header (default: %(default)s)') + parser.add_argument('input_file', type=validate_input_file(['.nii', '.nii.gz']), help='Input image file (*.nii, *.nii.gz)') + parser.add_argument('--header', '-H', action='store_true', help='Show NIFTI header (default: %(default)s)') parser.add_argument('--histogram', action='store_true', help='Show histogram (default: %(default)s)') - parser.add_argument('--bins', type=int, default=128, metavar='#', + parser.add_argument('--bins', '-b', type=validate_integer_min(1), default=128, metavar='#', help='Number of bins in histogram (default: %(default)s)') + parser.add_argument('--quiet', '-q', action='store_true', help='Suppress informational output to stdout') # Parse and display args = parser.parse_args() - # print(echo_arguments('ImageExam', vars(args))) + if not args.quiet: + print(echo_arguments('ImageExam', vars(args))) # Run program - ImageExam(**vars(args)) + result = ImageExam(**vars(args)) + return result if __name__ == '__main__': diff --git a/ogo/cli/ImageExtent.py b/ogo/cli/ImageExtent.py new file mode 100644 index 0000000..ec01897 --- /dev/null +++ b/ogo/cli/ImageExtent.py @@ -0,0 +1,269 @@ +# /------------------------------------------------------------------------------+ +# | 11-APR-2026 | +# | Copyright (c) Bone Imaging Laboratory | +# | All rights reserved | +# | bonelab@ucalgary.ca | +# +------------------------------------------------------------------------------+ + +# Imports +import argparse +import SimpleITK as sitk +from ogo.util.echo_arguments import echo_arguments +import ogo.util.Helper as ogo +from ogo.util.cli_tools import ( + check_overwrite, + validate_integer_min, + validate_input_file, + validate_output_file, + validate_extent, + load_config_from_json, + save_config_to_json, + print_image_info +) + +#------------------------------------------------------------------------------+ +# Main function +def ImageChangeExtent(input_file, output_file, offset, extent, filler, overwrite, quiet): + + # Check if output exists and should overwrite + if not check_overwrite(output_file, overwrite): + return 1 + + # Check that offset and extent are mutually exclusive + if offset is not None and extent is not None: + ogo.message('[ERROR] Cannot specify both --offset and --extent.') + ogo.message(' Please specify one or the other.') + return 1 + + # Check that at least one is specified + if offset is None and extent is None: + ogo.message('[ERROR] Must specify either --offset or --extent.') + return 1 + + # Read input image + if not quiet: + ogo.message(f'Reading input file: {input_file}') + + input_image = sitk.ReadImage(input_file) + + if not quiet: + print_image_info(input_file, input_image) + + # Get image properties + input_size = input_image.GetSize() + input_origin = input_image.GetOrigin() + input_spacing = input_image.GetSpacing() + pixel_type = input_image.GetPixelID() + + # Process based on offset or extent + if offset is not None: + # Validate offset dimensions + if len(offset) != 3: + ogo.message(f'[ERROR] Offset must have exactly 3 values, got {len(offset)}') + return 1 + + ox, oy, oz = offset + sx, sy, sz = input_size + + # Validate negative offsets don't exceed half dimension + if ox < 0 and abs(ox) > sx // 2: + ogo.message(f'[ERROR] Negative X offset ({ox}) exceeds half the X dimension ({sx // 2})') + return 1 + if oy < 0 and abs(oy) > sy // 2: + ogo.message(f'[ERROR] Negative Y offset ({oy}) exceeds half the Y dimension ({sy // 2})') + return 1 + if oz < 0 and abs(oz) > sz // 2: + ogo.message(f'[ERROR] Negative Z offset ({oz}) exceeds half the Z dimension ({sz // 2})') + return 1 + + # Calculate new extent based on offset + # Offset is applied symmetrically on both sides of each dimension + xmin, ymin, zmin = -ox, -oy, -oz + xmax = sx + ox + ymax = sy + oy + zmax = sz + oz + + if not quiet: + if ox >= 0 and oy >= 0 and oz >= 0: + ogo.message(f'Expanding image by offset: [{ox}, {oy}, {oz}]') + elif ox <= 0 and oy <= 0 and oz <= 0: + ogo.message(f'Cropping image by offset: [{ox}, {oy}, {oz}]') + else: + ogo.message(f'Adjusting image by offset: [{ox}, {oy}, {oz}]') + ogo.message(f'New extent (voxels): [{xmin}, {xmax}, {ymin}, {ymax}, {zmin}, {zmax}]') + + else: # extent is not None + # Extent already validated by validate_extent in main() + xmin, xmax, ymin, ymax, zmin, zmax = extent + + if not quiet: + ogo.message(f'Using explicit extent: [{xmin}, {xmax}, {ymin}, {ymax}, {zmin}, {zmax}]') + + # Calculate new size + new_size = [xmax - xmin, ymax - ymin, zmax - zmin] + + # Create output image filled with filler value + if not quiet: + ogo.message(f'Creating output image with size: {new_size}') + ogo.message(f'Filler value: {filler}') + + output_image = sitk.Image(new_size, pixel_type) + output_image.SetSpacing(input_spacing) + output_image.SetDirection(input_image.GetDirection()) + + # Calculate new origin + # Origin shifts when we change the extent + index_offset = [float(xmin), float(ymin), float(zmin)] + new_origin = input_image.TransformContinuousIndexToPhysicalPoint(index_offset) + output_image.SetOrigin(new_origin) + + # Fill with filler value + if filler != 0: + output_image = output_image + filler + + # Determine the region to paste from input image + # Calculate overlap between input image extent [0, sx) x [0, sy) x [0, sz) + # and new extent [xmin, xmax) x [ymin, ymax) x [zmin, zmax) + + paste_dest_start = [max(0, -xmin), max(0, -ymin), max(0, -zmin)] + crop_source_start = [max(0, xmin), max(0, ymin), max(0, zmin)] + + overlap_size = [ + min(input_size[0], xmax) - max(0, xmin), + min(input_size[1], ymax) - max(0, ymin), + min(input_size[2], zmax) - max(0, zmin) + ] + + # Check if there's any overlap + if overlap_size[0] <= 0 or overlap_size[1] <= 0 or overlap_size[2] <= 0: + ogo.message('[WARNING] No overlap between input image and new extent.') + ogo.message(' Output will contain only filler values.') + else: + # Extract the region from input image + if not quiet: + ogo.message(f'Extracting region from input: start={crop_source_start}, size={overlap_size}') + ogo.message(f'Pasting into output at: {paste_dest_start}') + + region = sitk.RegionOfInterest(input_image, overlap_size, crop_source_start) + + # Paste into output image + output_image = sitk.Paste(output_image, region, region.GetSize(), [0, 0, 0], paste_dest_start) + + # Write output + if not quiet: + ogo.message(f'Writing output file: {output_file}') + + sitk.WriteImage(output_image, output_file) + + if not quiet: + print_image_info(output_file, output_image) + ogo.message('[DONE]') + + return 0 + +def main(): + # Setup description + description=''' +Adjust the extent of an image by either cropping or expanding. + +This tool modifies the spatial extent of an image using either relative offsets or +absolute extent values. Offsets are applied symmetrically to both sides of each dimension. + +OFFSET MODE (--offset X Y Z): + - Positive values: Expand the image (pad with filler value) + - Negative values: Crop the image + - Values are applied symmetrically on both sides of each dimension + - Negative offsets cannot exceed half the dimension size + + Example: --offset 10 10 0 adds 10 voxels to each side in X and Y dimensions + Example: --offset -5 -5 -5 removes 5 voxels from each side in all dimensions + +EXTENT MODE (--extent xmin xmax ymin ymax zmin zmax): + - Define absolute voxel coordinates for the new extent + - Can specify both expansion (negative min or max beyond image) and cropping + - Coordinates are in voxel index space + + Example: --extent -10 100 -10 100 0 50 for a 90x90x50 image starting at [-10,-10,0] + +The two modes are mutually exclusive - use one or the other, not both. +Areas outside the original image are filled with the filler value (default: 0). +''' + epilog=''' +Example calls: +ogoImageChangeExtent input.nii.gz output.nii.gz --offset 10 10 10 +ogoImageChangeExtent input.nii.gz output.nii.gz --offset -5 -5 0 +ogoImageChangeExtent input.nii.gz output.nii.gz --extent -10 200 -10 200 0 100 +ogoImageChangeExtent input.nii.gz output.nii.gz -off 20 20 20 --filler -1000 +ogoImageChangeExtent input.nii.gz output.nii.gz -ext 0 150 0 150 10 90 -ow -q +''' + + # Setup argument parsing + parser = argparse.ArgumentParser( + formatter_class=argparse.RawTextHelpFormatter, + prog="ogoImageChangeExtent", + description=description, + epilog=epilog + ) + + # JSON configuration argument + parser.add_argument('--json', nargs='?', const='', metavar='FILE', + help='Output arguments as JSON to stdout (no file), or load arguments from JSON file') + + parser.add_argument('input_file', nargs='?', type=validate_input_file(['.nii', '.nii.gz']), + help='Input image file (*.nii, *.nii.gz)') + parser.add_argument('output_file', nargs='?', type=validate_output_file(['.nii', '.nii.gz']), + help='Output image file (*.nii, *.nii.gz)') + parser.add_argument('--offset', '-off', type=int, nargs=3, default=None, metavar=('X', 'Y', 'Z'), + help='Relative offset for each dimension (3 integers); positive=expand, negative=crop') + parser.add_argument('--extent', '-ext', type=int, nargs=6, default=None, metavar=('xmin', 'xmax', 'ymin', 'ymax', 'zmin', 'zmax'), + help='Absolute extent in voxel coordinates (6 integers)') + parser.add_argument('--filler', '-f', type=int, default=0, metavar='VALUE', + help='Value to use for padding/filling (default: %(default)s)') + parser.add_argument('--overwrite', '-ow', action='store_true', help='Overwrite output without asking') + parser.add_argument('--quiet', '-q', action='store_true', help='Suppress informational output to stdout') + + print() + + # Parse and display + args = parser.parse_args() + + # Handle --json flag + if args.json is not None: + if args.json == '': + # --json with no argument: output to stdout + save_config_to_json(args, None) + return 0 + else: + # --json file.json: load from file + if not args.json.endswith('.json'): + parser.error('JSON file must have .json extension') + json_config = load_config_from_json(args.json) + # Update args with JSON values (command line args override JSON) + for key, value in json_config.items(): + if hasattr(args, key): + # Only update if argument was not explicitly provided on command line + if getattr(args, key) == parser.get_default(key): + setattr(args, key, value) + + # Check that required positional arguments are provided + if args.input_file is None or args.output_file is None: + parser.error('input_file and output_file are required') + + # Validate extent if provided + if args.extent is not None: + try: + validate_extent(args.extent) + except argparse.ArgumentTypeError as e: + parser.error(str(e)) + + if not args.quiet: + ogo.message(echo_arguments('ImageChangeExtent', vars(args))) + + # Run program (exclude json argument) + run_args = vars(args).copy() + run_args.pop('json', None) + result = ImageChangeExtent(**run_args) + return result + +if __name__ == '__main__': + main() diff --git a/ogo/cli/ImageFilterByMask.py b/ogo/cli/ImageFilterByMask.py new file mode 100644 index 0000000..9009442 --- /dev/null +++ b/ogo/cli/ImageFilterByMask.py @@ -0,0 +1,253 @@ +# /------------------------------------------------------------------------------+ +# | 11-APR-2026 | +# | Copyright (c) Bone Imaging Laboratory | +# | All rights reserved | +# | bonelab@ucalgary.ca | +# +------------------------------------------------------------------------------+ + +# Imports +import argparse +import SimpleITK as sitk +from ogo.util.echo_arguments import echo_arguments +import ogo.util.Helper as ogo +from ogo.util.cli_tools import ( + check_overwrite, + validate_integer_min, + validate_input_file, + validate_output_file, + load_config_from_json, + save_config_to_json, + print_image_info +) + +#------------------------------------------------------------------------------+ +# Main function +def ImageFilterByMask(input_file, mask_file, output_file, mask_dilate, mask_erode, morphological, overwrite, quiet): + + # Check if output exists and should overwrite + if not check_overwrite(output_file, overwrite): + return 1 + + # Read input image + if not quiet: + ogo.message(f'Reading input file: {input_file}') + + input_image = sitk.ReadImage(input_file) + + if not quiet: + print_image_info(input_file, input_image) + + # Check pixel type + if input_image.GetPixelID() != sitk.sitkInt16: + ogo.message(f'[ERROR] Input file must be Int16, got {sitk.GetPixelIDValueAsString(input_image.GetPixelID())}') + return 1 + + # Read mask + if not quiet: + ogo.message(f'Reading mask file: {mask_file}') + + mask_image = sitk.ReadImage(mask_file) + + if not quiet: + print_image_info(mask_file, mask_image) + + # Check mask pixel type + if mask_image.GetPixelID() not in [sitk.sitkUInt8, sitk.sitkUInt16]: + ogo.message(f'[ERROR] Mask file must be UInt8 or UInt16, got {sitk.GetPixelIDValueAsString(mask_image.GetPixelID())}') + return 1 + + # Check that images have matching dimensions + if input_image.GetSize() != mask_image.GetSize(): + ogo.message('[ERROR] Input and mask dimensions do not match.') + ogo.message(f' input_file: {input_image.GetSize()}') + ogo.message(f' mask_file: {mask_image.GetSize()}') + return 1 + + # Convert mask to binary (0 or 1) for morphological operations + processed_mask = sitk.Cast(mask_image != 0, sitk.sitkUInt8) + + # Warn if both dilate and erode are specified without explicit morphological parameter + if mask_dilate > 0 and mask_erode > 0 and morphological is None: + ogo.message('[WARNING] Both dilate and erode specified. Dilation will be applied first, then erosion.') + ogo.message(' This is equivalent to a morphological closing operation.') + ogo.message(' To explicitly control operation order, use --morphological opening or --morphological closing.') + + # Warn if morphological operation selected but no dilate/erode specified + if morphological in ['opening', 'closing'] and mask_dilate == 0 and mask_erode == 0: + ogo.message(f'[WARNING] Morphological {morphological} selected but neither --mask_dilate nor --mask_erode specified.') + ogo.message(' No morphological operations will be applied to the mask.') + + # Determine operation order based on morphological parameter + if morphological == 'opening': + # Opening: erode then dilate (removes small objects, smooths boundaries from outside) + if not quiet: + ogo.message(f'Using morphological opening (erode: {mask_erode} voxels, dilate: {mask_dilate} voxels)') + + if mask_erode > 0: + if not quiet: + ogo.message(f'Eroding mask by {mask_erode} voxels (opening: step 1/2)') + erode_filter = sitk.BinaryErodeImageFilter() + erode_filter.SetKernelRadius(mask_erode) + processed_mask = erode_filter.Execute(processed_mask) + + if mask_dilate > 0: + if not quiet: + ogo.message(f'Dilating mask by {mask_dilate} voxels (opening: step 2/2)') + dilate_filter = sitk.BinaryDilateImageFilter() + dilate_filter.SetKernelRadius(mask_dilate) + processed_mask = dilate_filter.Execute(processed_mask) + + elif morphological == 'closing': + # Closing: dilate then erode (fills small holes, smooths boundaries from inside) + if not quiet: + ogo.message(f'Using morphological closing (dilate: {mask_dilate} voxels, erode: {mask_erode} voxels)') + + if mask_dilate > 0: + if not quiet: + ogo.message(f'Dilating mask by {mask_dilate} voxels (closing: step 1/2)') + dilate_filter = sitk.BinaryDilateImageFilter() + dilate_filter.SetKernelRadius(mask_dilate) + processed_mask = dilate_filter.Execute(processed_mask) + + if mask_erode > 0: + if not quiet: + ogo.message(f'Eroding mask by {mask_erode} voxels (closing: step 2/2)') + erode_filter = sitk.BinaryErodeImageFilter() + erode_filter.SetKernelRadius(mask_erode) + processed_mask = erode_filter.Execute(processed_mask) + + else: + # Default behavior: dilate then erode + if mask_dilate > 0: + if not quiet: + ogo.message(f'Dilating mask by {mask_dilate} voxels') + dilate_filter = sitk.BinaryDilateImageFilter() + dilate_filter.SetKernelRadius(mask_dilate) + processed_mask = dilate_filter.Execute(processed_mask) + + if mask_erode > 0: + if not quiet: + ogo.message(f'Eroding mask by {mask_erode} voxels') + erode_filter = sitk.BinaryErodeImageFilter() + erode_filter.SetKernelRadius(mask_erode) + processed_mask = erode_filter.Execute(processed_mask) + + # Create binary mask as Int16: 1 where mask is non-zero, 0 where mask is zero + binary_mask = sitk.Cast(processed_mask, sitk.sitkInt16) + + # Apply mask: where mask=0, output=0; where mask!=0, output=input value + if not quiet: + ogo.message('Applying mask to input image') + + output_image = input_image * binary_mask + + # Ensure output is Int16 + output_image = sitk.Cast(output_image, sitk.sitkInt16) + + # Copy metadata from input image + output_image.CopyInformation(input_image) + + # Write output + if not quiet: + ogo.message(f'Writing output file: {output_file}') + + sitk.WriteImage(output_image, output_file) + + if not quiet: + print_image_info(output_file, output_image) + ogo.message('[DONE]') + + return 0 + +def main(): + # Setup description + description=''' +Utility to apply a binary mask to an input image. + +This tool applies a mask to set voxels to zero where the mask is non-zero, and passes +through the input image values where the mask is zero. The input image must be +Int16 and the mask must be UInt8 or UInt16. The output will always be Int16. + +Optional morphological operations (dilation and erosion) can be applied to the +mask before applying it to the input image. + +Morphological operation types: + - opening: Erode then dilate (removes small objects, smooths from outside) + - closing: Dilate then erode (fills small holes, smooths from inside) + - Default (no --morphological): Dilate then erode if both specified +''' + epilog=''' +Example calls: +ogoImageFilterByMask input.nii.gz mask.nii.gz output.nii.gz +ogoImageFilterByMask input.nii.gz mask.nii.gz output.nii.gz --mask_dilate 2 +ogoImageFilterByMask input.nii.gz mask.nii.gz output.nii.gz --mask_erode 1 -ow +ogoImageFilterByMask input.nii.gz mask.nii.gz output.nii.gz -d 3 -e 1 -q +ogoImageFilterByMask input.nii.gz mask.nii.gz output.nii.gz -d 2 -e 2 --morphological opening +ogoImageFilterByMask input.nii.gz mask.nii.gz output.nii.gz -d 3 -e 3 -mo closing +''' + + # Setup argument parsing + parser = argparse.ArgumentParser( + formatter_class=argparse.RawTextHelpFormatter, + prog="ogoImageFilterByMask", + description=description, + epilog=epilog + ) + + # JSON configuration argument + parser.add_argument('--json', nargs='?', const='', metavar='FILE', + help='Output arguments as JSON to stdout (no file), or load arguments from JSON file') + + parser.add_argument('input_file', nargs='?', type=validate_input_file(['.nii', '.nii.gz']), + help='Input image file (*.nii, *.nii.gz) - must be Int16') + parser.add_argument('mask_file', nargs='?', type=validate_input_file(['.nii', '.nii.gz']), + help='Mask image file (*.nii, *.nii.gz) - must be UInt8 or UInt16') + parser.add_argument('output_file', nargs='?', type=validate_output_file(['.nii', '.nii.gz']), + help='Output masked image file (*.nii, *.nii.gz)') + parser.add_argument('--mask_dilate', '-d', type=validate_integer_min(0), default=0, metavar='VOXELS', + help='Dilate mask by specified number of voxels (default: %(default)s)') + parser.add_argument('--mask_erode', '-e', type=validate_integer_min(0), default=0, metavar='VOXELS', + help='Erode mask by specified number of voxels (default: %(default)s)') + parser.add_argument('--morphological', '-mo', choices=['opening', 'closing'], default=None, metavar='TYPE', + help='Morphological operation type: opening (erode then dilate) or closing (dilate then erode)') + parser.add_argument('--overwrite','-ow', action='store_true', help='Overwrite output without asking') + parser.add_argument('--quiet','-q', action='store_true', help='Suppress informational output to stdout') + + print() + + # Parse and display + args = parser.parse_args() + + # Handle --json flag + if args.json is not None: + if args.json == '': + # --json with no argument: output to stdout + save_config_to_json(args, None) + return 0 + else: + # --json file.json: load from file + if not args.json.endswith('.json'): + parser.error('JSON file must have .json extension') + json_config = load_config_from_json(args.json) + # Update args with JSON values (command line args override JSON) + for key, value in json_config.items(): + if hasattr(args, key): + # Only update if argument was not explicitly provided on command line + if getattr(args, key) == parser.get_default(key): + setattr(args, key, value) + + # Check that required positional arguments are provided + if args.input_file is None or args.mask_file is None or args.output_file is None: + parser.error('input_file, mask_file, and output_file are required') + + if not args.quiet: + ogo.message(echo_arguments('ImageFilterByMask', vars(args))) + + # Run program (exclude json argument) + run_args = vars(args).copy() + run_args.pop('json', None) + result = ImageFilterByMask(**run_args) + return result + +if __name__ == '__main__': + main() diff --git a/ogo/cli/ImageMirror.py b/ogo/cli/ImageMirror.py new file mode 100644 index 0000000..4e45971 --- /dev/null +++ b/ogo/cli/ImageMirror.py @@ -0,0 +1,173 @@ +# /------------------------------------------------------------------------------+ +# | 11-APR-2026 | +# | Copyright (c) Bone Imaging Laboratory | +# | All rights reserved | +# | bonelab@ucalgary.ca | +# +------------------------------------------------------------------------------+ + +# Imports +import argparse +import SimpleITK as sitk +from ogo.util.echo_arguments import echo_arguments +import ogo.util.Helper as ogo +from ogo.util.cli_tools import ( + check_overwrite, + validate_input_file, + validate_output_file, + load_config_from_json, + save_config_to_json, + print_image_info +) + +#------------------------------------------------------------------------------+ +# Main function +def ImageMirror(input_file, output_file, mirror_plane, flip_about_origin, overwrite, quiet): + + # Check if output exists and should overwrite + if not check_overwrite(output_file, overwrite): + return 1 + + # Read input image + if not quiet: + ogo.message(f'Reading input file: {input_file}') + + input_image = sitk.ReadImage(input_file) + + if not quiet: + print_image_info(input_file, input_image) + + # Determine which axis to flip based on mirror plane + flip_axes = [False, False, False] + + if mirror_plane == 'X': + flip_axes[0] = True # Flip X axis + axis_name = 'X' + elif mirror_plane == 'Y': + flip_axes[1] = True # Flip Y axis + axis_name = 'Y' + elif mirror_plane == 'Z': + flip_axes[2] = True # Flip Z axis + axis_name = 'Z' + + if not quiet: + ogo.message(f'Flipping image along {axis_name} axis') + if flip_about_origin: + ogo.message(f'Flip about origin: True (flipping across coordinate system origin)') + else: + ogo.message(f'Flip about origin: False (flipping within coordinate grid)') + + # Apply flip filter + # Note: SimpleITK automatically updates origin and direction matrix + # to preserve physical space. Voxel indices are reversed. + output_image = sitk.Flip(input_image, flip_axes, flip_about_origin) + + if not quiet: + ogo.message(f'Voxel order reversed along {axis_name} axis') + if not flip_about_origin: + ogo.message(f'Physical location of voxels preserved') + + # Write output + if not quiet: + ogo.message(f'Writing output file: {output_file}') + + sitk.WriteImage(output_image, output_file) + + if not quiet: + print_image_info(output_file, output_image) + ogo.message('[DONE]') + + return 0 + +def main(): + # Setup description + description=''' +Mirror an image by flipping along a specified axis. + +This tool flips an image along one of the three principal axes: + - X: Flip along the X axis + - Y: Flip along the Y axis + - Z: Flip along the Z axis + +The mirroring operation reverses the order of voxels along the specified axis. + +IMPORTANT: SimpleITK automatically updates the origin and direction matrix of the +image so that the physical location of the pixels remains consistent. If you view +the flipped image in a physical-space-aware tool like ITK-SNAP or 3D Slicer, it +may look the same as the original, but the voxel indices (i, j, k) will be reversed. + +Flip modes: + - flip_about_origin=False (default): Flips within the coordinate grid, + preserving physical location of voxels + - flip_about_origin=True: Flips across the coordinate system's origin, + which may change the physical location +''' + epilog=''' +Example calls: +ogoImageMirror input.nii.gz output.nii.gz X +ogoImageMirror input.nii.gz output.nii.gz Y -ow +ogoImageMirror input.nii.gz output.nii.gz Z --quiet +ogoImageMirror input.nii.gz mirrored.nii.gz X --flip_about_origin +ogoImageMirror input.nii.gz mirrored.nii.gz Y -fao -q +''' + + # Setup argument parsing + parser = argparse.ArgumentParser( + formatter_class=argparse.RawTextHelpFormatter, + prog="ogoImageMirror", + description=description, + epilog=epilog + ) + + # JSON configuration argument + parser.add_argument('--json', nargs='?', const='', metavar='FILE', + help='Output arguments as JSON to stdout (no file), or load arguments from JSON file') + + parser.add_argument('input_file', nargs='?', type=validate_input_file(['.nii', '.nii.gz']), + help='Input image file (*.nii, *.nii.gz)') + parser.add_argument('output_file', nargs='?', type=validate_output_file(['.nii', '.nii.gz']), + help='Output image file (*.nii, *.nii.gz)') + parser.add_argument('mirror_plane', nargs='?', choices=['X', 'Y', 'Z'], + help='Axis to flip: X, Y, or Z') + parser.add_argument('--flip_about_origin', '-fao', action='store_true', default=False, + help='Flip across coordinate system origin (default: flip within coordinate grid)') + parser.add_argument('--overwrite', '-ow', action='store_true', help='Overwrite output without asking') + parser.add_argument('--quiet', '-q', action='store_true', help='Suppress informational output to stdout') + + print() + + # Parse and display + args = parser.parse_args() + + # Handle --json flag + if args.json is not None: + if args.json == '': + # --json with no argument: output to stdout + save_config_to_json(args, None) + return 0 + else: + # --json file.json: load from file + if not args.json.endswith('.json'): + parser.error('JSON file must have .json extension') + json_config = load_config_from_json(args.json) + # Update args with JSON values (command line args override JSON) + for key, value in json_config.items(): + if hasattr(args, key): + # Only update if argument was not explicitly provided on command line + if getattr(args, key) == parser.get_default(key): + setattr(args, key, value) + + # Check that required positional arguments are provided + if args.input_file is None or args.output_file is None or args.mirror_plane is None: + parser.error('input_file, output_file, and mirror_plane are required') + + if not args.quiet: + ogo.message(echo_arguments('ImageMirror', vars(args))) + + # Run program (exclude json argument) + run_args = vars(args).copy() + run_args.pop('json', None) + result = ImageMirror(**run_args) + return result + +if __name__ == '__main__': + main() diff --git a/ogo/cli/ImageThreshold.py b/ogo/cli/ImageThreshold.py index 94d08d3..7d783ca 100644 --- a/ogo/cli/ImageThreshold.py +++ b/ogo/cli/ImageThreshold.py @@ -14,168 +14,182 @@ from ogo.util.echo_arguments import echo_arguments import ogo.util.Helper as ogo import ogo.dat.OgoMasterLabels as lb +from ogo.util.cli_tools import ( + check_overwrite, + validate_input_file, + validate_output_file, + validate_integer_range, + load_config_from_json, + save_config_to_json, + print_image_info +) -def ImageThreshold(input_image, output_image, lower, upper, inside, outside, no_output, overwrite): - +# Pixel type value ranges (min, max) +PIXEL_TYPE_RANGES = { + sitk.sitkInt8: (-128, 127), + sitk.sitkUInt8: (0, 255), + sitk.sitkInt16: (-32768, 32767), + sitk.sitkUInt16: (0, 65535), + sitk.sitkInt32: (-2147483648, 2147483647), + sitk.sitkUInt32: (0, 4294967295), +} + +def ImageThreshold(input_file, output_file, lower, upper, inside, outside, overwrite, quiet): + # Check if output exists and should overwrite - if os.path.isfile(output_image) and not overwrite and not no_output: - result = input('File \"{}\" already exists. Overwrite? [y/n]: '.format(output_image)) - if result.lower() not in ['y', 'yes']: - print('Not overwriting. Exiting...') - os.sys.exit() - + if not check_overwrite(output_file, overwrite): + return 1 + + # Validate threshold parameters + if lower > upper: + ogo.message('[ERROR] Lower threshold cannot be greater than upper threshold.') + ogo.message(f' lower={lower}, upper={upper}') + return 1 + + if inside == outside: + if not quiet: + ogo.message('[WARNING] Setting inside equal to outside makes no sense!') + # Read input - if not os.path.isfile(input_image): - os.sys.exit('[ERROR] Cannot find file \"{}\"'.format(input_image)) - - ogo.message('Reading input CT image to be cropped:') - ogo.message(' \"{}\"'.format(input_image)) - ct = sitk.ReadImage(input_image) - - # Checking image type - ogo.message('Input image type is {}'.format(ct.GetPixelIDTypeAsString())) - image_type = ct.GetPixelIDValue() - if image_type is sitk.sitkInt8: - type_min = -128 - type_max = 127 - elif image_type is sitk.sitkUInt8: - type_min = 0 - type_max = 255 - elif image_type is sitk.sitkInt16: - type_min = -32768 - type_max = 32767 - elif image_type is sitk.sitkUInt16: - type_min = 0 - type_max = 65535 - elif image_type is sitk.sitkInt32: - type_min = -2147483648 - type_max = 2147483647 - elif image_type is sitk.sitkUInt32: - type_min = 0 - type_max = 4294967295 - else: + if not quiet: + ogo.message(f'Reading input image:') + ogo.message(f' {input_file}') + + input_image = sitk.ReadImage(input_file) + + if not quiet: + print_image_info(input_file, input_image) + + # Checking image type and validate thresholds + image_type = input_image.GetPixelIDValue() + if not quiet: + ogo.message(f'Input image type is {input_image.GetPixelIDTypeAsString()}') + + # Get type min/max from constant dictionary + if image_type not in PIXEL_TYPE_RANGES: ogo.message('[ERROR] Unexpected image input type.') + return 1 + + type_min, type_max = PIXEL_TYPE_RANGES[image_type] + + # Adjust thresholds if they exceed type limits + lower_adjusted = lower + upper_adjusted = upper + if lower < type_min: - lower = type_min - ogo.message('[WARNING] Invalid lower threshold. Changed to {}'.format(lower)) + lower_adjusted = type_min + if not quiet: + ogo.message(f'[WARNING] Lower threshold {lower} < type minimum {type_min}. Adjusted to {type_min}') + if upper > type_max: - upper = type_max - ogo.message('[WARNING] Invalid upper threshold. Changed to {}'.format(upper)) + upper_adjusted = type_max + if not quiet: + ogo.message(f'[WARNING] Upper threshold {upper} > type maximum {type_max}. Adjusted to {type_max}') - # Image statistics + # Get image statistics stats = sitk.StatisticsImageFilter() - stats.Execute(ct) - - if lower < stats.GetMinimum(): - ogo.message('Lower ') - # Report information about input image - dim = ct.GetSize() - spacing = ct.GetSpacing() - origin = ct.GetOrigin() - phys_dim = [x * y for x, y in zip(dim, spacing)] - position = [math.floor(x / y) for x, y in zip(origin, spacing)] - - guard = '!-------------------------------------------------------------------------------' - print(guard) - print('!> dim {:>8} {:>8} {:>8}'.format(*dim)) - print('!> off {:>8} {:>8} {:>8}'.format('-', '-', '-')) - print('!> pos {:>8} {:>8} {:>8}'.format(*position)) - print('!> element size in mm {:>8.4f} {:>8.4f} {:>8.4f}'.format(*spacing)) - print('!> phys dim in mm {:>8.4f} {:>8.4f} {:>8.4f}'.format(*phys_dim)) - print('!> ') - print('!> min {:>8.0f}'.format(stats.GetMinimum())) - print('!> max {:>8.0f}'.format(stats.GetMaximum())) - print(guard) - - # Threshold - print(guard) - print('!> lower {:>8}'.format(lower)) - print('!> upper {:>8}'.format(upper)) - print('!> inside {:>8}'.format(inside)) - print('!> outside {:>8}'.format(outside)) - print(guard) + stats.Execute(input_image) - if lower > upper: - os.sys.exit('[ERROR] Lower threshold cannot be greater than upper threshold.') - if inside == outside: - os.message('[WARNING] Setting inside equal to outside makes no sense!') - - ct_thres = sitk.BinaryThresholdImageFilter() - ct_thres.SetInsideValue(inside) - ct_thres.SetOutsideValue(outside) - ct_thres.SetLowerThreshold(lower) - ct_thres.SetUpperThreshold(upper) - ct_out = ct_thres.Execute(ct) - - # Report stats on resulting image - stats = sitk.LabelIntensityStatisticsImageFilter() - stats.Execute(ct_out,ct) - - n_labels = stats.GetNumberOfLabels() - basename = os.path.basename(input_image) - - line_hdr = 'line_{},{},'.format('hdr','name') - line_units = 'line_{},{},'.format('units','[text]') - line_data = 'line_{},{},'.format('data',basename) - line_hdr += '{},{},{},{},{},{},{}'.format('desc','label','total_vol','n_voxels','cx','cy','cz') - line_units += '{},{},{},{},{},{},{}'.format('[text]','[#]','[mm3]','[N]','[i]','[j]','[k]') - - report = '' - report += ' {:>20s}\n'.format('_______________________________________________________________________Output') - report += ' {:>20s} {:>10s} {:>10s} {:>19s}\n'.format('Label','Volume','N Voxels','Centroid') - - if n_labels > 0: - for label in stats.GetLabels(): - try: - desc = lb.labels_dict[label]['LABEL'] - except KeyError: - desc = 'unknown label' - - centroid = stats.GetCentroid(label) - bounding_box = stats.GetBoundingBox(label) - bb=[0]*3 - bb[0] = bounding_box[0] + int(math.ceil(bounding_box[3]/2)) - bb[1] = bounding_box[1] + int(math.ceil(bounding_box[4]/2)) - bb[2] = bounding_box[2] + int(math.ceil(bounding_box[5]/2)) - line_data += '{},{},{:.1f},{},{},{},{}'.format(desc,label,stats.GetPhysicalSize(label),stats.GetNumberOfPixels(label),bb[0],bb[1],bb[2]) - - report += ' {:>15s}{:>5s} {:10.1f} {:10d} ({:5d},{:5d},{:5d})\n'\ - .format('('+desc+')',str(label),stats.GetPhysicalSize(label),stats.GetNumberOfPixels(label),bb[0],bb[1],bb[2]) - else: - line_data += '{},{},{:.1f},{},{},{},{}'.format('empty','empty',0,0,0,0,0) - report += ' {:>20s}\n'.format('-- No data --') - - print(report) - print(line_hdr) - print(line_units) - print(line_data) - print() - - if not no_output: - ogo.message('Writing output to file {}'.format(output_image)) - sitk.WriteImage(ct_out, output_image) - else: - ogo.message('Writing output to file {} is SUPPRESSED'.format(output_image)) - - ## Report information about input image - #dim = ct_out.GetSize() - #spacing = ct_out.GetSpacing() - #origin = ct_out.GetOrigin() - #phys_dim = [x * y for x, y in zip(dim, spacing)] - #position = [math.floor(x / y) for x, y in zip(origin, spacing)] - - ogo.message('Done ogoImageThreshold!') + # Display threshold information + if not quiet: + guard = '!-------------------------------------------------------------------------------' + print(guard) + print(f'!> Threshold Parameters:') + print(f'!> lower {lower_adjusted:>12}') + print(f'!> upper {upper_adjusted:>12}') + print(f'!> inside {inside:>12}') + print(f'!> outside {outside:>12}') + print(guard) + + # Apply threshold + threshold_filter = sitk.BinaryThresholdImageFilter() + threshold_filter.SetInsideValue(inside) + threshold_filter.SetOutsideValue(outside) + threshold_filter.SetLowerThreshold(lower_adjusted) + threshold_filter.SetUpperThreshold(upper_adjusted) + output_image = threshold_filter.Execute(input_image) + + # Calculate label statistics on the result + label_stats = sitk.LabelIntensityStatisticsImageFilter() + label_stats.Execute(output_image, input_image) + + n_labels = label_stats.GetNumberOfLabels() + + # Display label statistics + if not quiet: + if n_labels > 0: + ogo.message('') + ogo.message('Label Statistics:') + ogo.message(f' {"Label":<20s} {"Volume":<12s} {"N Voxels":<12s} {"Centroid":>19s}') + + for label in label_stats.GetLabels(): + # Get label description + desc = lb.labels_dict.get(label, {}).get('LABEL', 'unknown label') + + # Calculate bounding box center + bbox = label_stats.GetBoundingBox(label) + center = [ + bbox[0] + bbox[3] // 2, + bbox[1] + bbox[4] // 2, + bbox[2] + bbox[5] // 2 + ] + + phys_size = label_stats.GetPhysicalSize(label) + n_pixels = label_stats.GetNumberOfPixels(label) + + ogo.message(f' {desc:<15s}({label:<3d}) ' + f'{phys_size:>10.1f} mm³ ' + f'{n_pixels:>10d} ' + f'({center[0]:>5d},{center[1]:>5d},{center[2]:>5d})') + else: + ogo.message('') + ogo.message('[INFO] No labels found in output image.') + + # Write output + if not quiet: + ogo.message(f'Writing output to: {output_file}') + + sitk.WriteImage(output_image, output_file) + + if not quiet: + print_image_info(output_file, output_image) + ogo.message('[DONE]') + + return 0 def main(): # Setup description description = ''' Utility to threshold a NIFTI file. + +This tool applies binary thresholding to an image, setting pixels within the +threshold range [lower, upper] to the 'inside' value, and all other pixels to +the 'outside' value. + +The threshold operation uses SimpleITK's BinaryThresholdImageFilter, which +includes both lower and upper bounds (i.e., lower ≤ pixel ≤ upper). + +THRESHOLD BEHAVIOR: + - Pixels where (lower ≤ value ≤ upper) are set to 'inside' value + - All other pixels are set to 'outside' value + +TYPICAL USAGE: + - Create binary mask: --inside 1 --outside 0 + - Extract specific range: --lower 250 --upper 3000 + - Bone segmentation: --lower 250 --inside 127 --outside 0 + +JSON CONFIGURATION: + - To output current arguments as JSON: --json + - To load arguments from JSON file: --json my_config.json ''' epilog = ''' Example calls: -ogoImageThreshold input.nii.gz --output_image output.nii.gz --upper 1500 -ogoImageThreshold input.nii.gz --output_image output.nii.gz --no_output --lower 2500 --upper 5000 --inside 81 +ogoImageThreshold input.nii.gz output.nii.gz --upper 1500 +ogoImageThreshold input.nii.gz output.nii.gz --lower 250 --upper 3000 --inside 1 --outside 0 +ogoImageThreshold input.nii.gz output.nii.gz --lower 2500 --upper 5000 --inside 81 --quiet +ogoImageThreshold input.nii.gz output.nii.gz --lower 500 --json > my_config.json +ogoImageThreshold --json my_config.json ''' # Setup argument parsing @@ -185,21 +199,61 @@ def main(): description=description, epilog=epilog ) - parser.add_argument('input_image', help='Input image file (*.nii, *.nii.gz)') - parser.add_argument('output_image', help='Output image file (*.nii, *.nii.gz)') - parser.add_argument('--lower', type=int, default=250, metavar='VAL', help='Lower threshold (min=-32768, default: %(default)s)') - parser.add_argument('--upper', type=int, default=3000, metavar='VAL', help='Upper threshold (max=32767, default: %(default)s)') - parser.add_argument('--inside', type=int, default=127, metavar='VAL', help='Inside value (default: %(default)s)') - parser.add_argument('--outside', type=int, default=0, metavar='VAL', help='Outside value (default: %(default)s)') - parser.add_argument('--no_output', action='store_true', help='Suppress writing output (useful if just looking for statistics)') - parser.add_argument('--overwrite', action='store_true', help='Overwrite output without asking') + + # JSON configuration argument + parser.add_argument('--json', nargs='?', const='', metavar='FILE', + help='Output arguments as JSON to stdout (no file), or load arguments from JSON file') + + parser.add_argument('input_file', nargs='?', type=validate_input_file(['.nii', '.nii.gz']), + help='Input image file (*.nii, *.nii.gz)') + parser.add_argument('output_file', nargs='?', type=validate_output_file(['.nii', '.nii.gz']), + help='Output image file (*.nii, *.nii.gz)') + parser.add_argument('--lower','-l', type=int, default=250, metavar='VAL', + help='Lower threshold (default: %(default)s)') + parser.add_argument('--upper', '-u', type=int, default=3000, metavar='VAL', + help='Upper threshold (default: %(default)s)') + parser.add_argument('--inside', '-i', type=validate_integer_range(0, 255), default=127, metavar='VAL', + help='Inside value (0-255, default: %(default)s)') + parser.add_argument('--outside', '-o', type=validate_integer_range(0, 255), default=0, metavar='VAL', + help='Outside value (0-255, default: %(default)s)') + parser.add_argument('--overwrite', '-ow', action='store_true', + help='Overwrite output without asking') + parser.add_argument('--quiet', '-q', action='store_true', + help='Suppress informational output to stdout') # Parse and display args = parser.parse_args() - print(echo_arguments('ImageThreshold', vars(args))) + + # Handle --json flag + if args.json is not None: + if args.json == '': + # --json with no argument: output to stdout + save_config_to_json(args, None) + return 0 + else: + # --json file.json: load from file + if not args.json.endswith('.json'): + parser.error('JSON file must have .json extension') + json_config = load_config_from_json(args.json) + # Update args with JSON values (command line args override JSON) + for key, value in json_config.items(): + if hasattr(args, key): + # Only update if argument was not explicitly provided on command line + if getattr(args, key) == parser.get_default(key): + setattr(args, key, value) + + # Check that required positional arguments are provided + if args.input_file is None or args.output_file is None: + parser.error('input_file and output_file are required') + + if not args.quiet: + ogo.message(echo_arguments('ImageThreshold', vars(args))) - # Run program - ImageThreshold(**vars(args)) + # Run program (exclude json argument) + run_args = vars(args).copy() + run_args.pop('json', None) + result = ImageThreshold(**run_args) + return result if __name__ == '__main__': diff --git a/ogo/cli/MaterialDecomposition.py b/ogo/cli/MaterialDecomposition.py new file mode 100644 index 0000000..0252123 --- /dev/null +++ b/ogo/cli/MaterialDecomposition.py @@ -0,0 +1,431 @@ +# /------------------------------------------------------------------------------+ +# | 5-DEC-2025 | +# | Copyright (c) Bone Imaging Laboratory | +# | All rights reserved | +# | bonelab@ucalgary.ca | +# +------------------------------------------------------------------------------+ + +# Imports +import sys +import os +import argparse +import math +import copy +import SimpleITK as sitk +import numpy as np +from ogo.dat.MassAttenuationTables import mass_attenuation_tables +from ogo.util.echo_arguments import echo_arguments +import ogo.util.Helper as ogo +import ogo.util.spectral_util as md + +# Just used to do some debugging +def print_pattern(pattern): + basename,dirname,name,ext = md.parse_filename(pattern) + print('{:>20s}: {:s}'.format('INPUT', pattern)) + print('{:>20s}: {:s}'.format('basename', basename)) + print('{:>20s}: {:s}'.format('dirname', dirname)) + print('{:>20s}: {:s}'.format('name', name)) + print('{:>20s}: {:s}'.format('ext', ext)) + +def MaterialDecomposition(images, energies, materials, pattern, suppress, quiet, overwrite): + + # Initialize variables + choices = list(md.mass_density().keys()) + + no_print_str = '--suppressed--' + + # Check that input images exist and are type NIfTI + for image in images: + if not os.path.isfile(image): + os.sys.exit('[ERROR] Cannot find file \"{}\"'.format(image)) + if not (image.lower().endswith('.nii') or image.lower().endswith('.nii.gz')): + os.sys.exit('[ERROR] Input must be type NIFTI file: \"{}\"'.format(image)) + + # Check that the energies are in a valid range and not equal + for energy in energies: + if energy<40 or energy>190: + os.sys.exit('[ERROR] Energies must be between 40 and 190 keV.') + if energies[0] == energies[1]: + os.sys.exit('[ERROR] Energies cannot be equal.') + + # Check that the materials are valid + if len(materials)<2 or len(materials)>3: + os.sys.exit('[ERROR] Only 2- or 3-material decomposition is valid.') + for material in materials: + if material not in choices: + os.sys.exit('[ERROR] Material is not valid: {}'.format(material)) + + # Look up mass attenuations of materials + mass_attenuations=[] + for energy in energies: + for material in materials: + mass_attenuations.append(md.interpolate_mass_attenuation(material,energy)) + + # Look up ICRU mass densities of materials + mass_densities=[] + for material in materials: + mass_densities.append(md.mass_density()[material]) + + if not quiet: + ogo.message('Low energy image at {:.1f} keV:'.format(energies[0])) + basename,dirname,name,ext = md.parse_filename(images[0]) + ogo.message(' {:s}'.format(name+ext)) + ogo.message('') + ogo.message('High energy image at {:.1f} keV:'.format(energies[1])) + basename,dirname,name,ext = md.parse_filename(images[1]) + ogo.message(' {:s}'.format(name+ext)) + ogo.message('') + ogo.message('Materials: [' + ', '.join('{:s}'.format(i) for i in materials) + ']') + ogo.message('') + ogo.message('Mass attenuations: [' + ', '.join('{:.4f}'.format(i) for i in mass_attenuations) + ']') + ogo.message('') + ogo.message('Mass densities: [' + ', '.join('{:.4f}'.format(i) for i in mass_densities) + ']') + ogo.message('') + + # Define output filenames and check for overwrite + ofiles = [] + if pattern is not None: + basename,dirname,name,ext = md.parse_filename(pattern) + if not ext: + ext='.nii.gz' + if dirname and not os.path.isdir(os.path.dirname(pattern)): + os.sys.exit('[ERROR] Output directory does not exist: {}'.format(os.path.dirname(pattern))) + else: + basename,dirname,name,ext = md.parse_filename(images[0]) + + for material in materials: + if material not in suppress: + ofiles.append(os.path.join(dirname, name + '_' + '{}'.format(len(materials)) + material + ext)) + else: + ofiles.append(no_print_str) + + for ofile in ofiles: + if os.path.isfile(ofile) and not overwrite: + result = input('File \"{}\" already exists. Overwrite? [y/n]: '.format(ofile)) + if result.lower() not in ['y', 'yes']: + ogo.message('Not overwriting. Exiting...') + os.sys.exit() + + # Read in CT data --------------------------------------------------------------------------------- + if not quiet: + ogo.message('Input files') + ct1 = sitk.ReadImage(images[0]) + if not quiet: + ogo.message(' low image: {:s}'.format(os.path.basename(images[0]))) + ogo.message(' energy: {:.1f}'.format(energies[0])) + ogo.message(' directory: {:s}'.format(os.path.dirname(images[0]))) + ogo.message(' dimension: {:d}'.format(ct1.GetDimension())) + ogo.message(' type: {:s}'.format(ct1.GetPixelIDTypeAsString())) + ogo.message(' size: [' + ', '.join('{:d}'.format(i) for i in ct1.GetSize()) + ']') + ogo.message(' origin: [' + ', '.join('{:.2f}'.format(i) for i in ct1.GetOrigin()) + ']') + ogo.message(' spacing: [' + ', '.join('{:.2f}'.format(i) for i in ct1.GetSpacing()) + ']') + ogo.message(' direction: [' + ', '.join('{:.1f}'.format(i) for i in ct1.GetDirection()) + ']') + ogo.message('') + + ct2 = sitk.ReadImage(images[1]) + if not quiet: + ogo.message(' high image: {:s}'.format(os.path.basename(images[1]))) + ogo.message(' directory: {:s}'.format(os.path.dirname(images[1]))) + ogo.message(' dimension: {:d}'.format(ct2.GetDimension())) + ogo.message(' type: {:s}'.format(ct2.GetPixelIDTypeAsString())) + ogo.message(' size: [' + ', '.join('{:d}'.format(i) for i in ct2.GetSize()) + ']') + ogo.message(' origin: [' + ', '.join('{:.2f}'.format(i) for i in ct2.GetOrigin()) + ']') + ogo.message(' spacing: [' + ', '.join('{:.2f}'.format(i) for i in ct2.GetSpacing()) + ']') + ogo.message(' direction: [' + ', '.join('{:.1f}'.format(i) for i in ct2.GetDirection()) + ']') + ogo.message('') + + # Check dimensions of images are the same and pixel type + if ct1.GetDimension() != ct2.GetDimension(): + os.sys.exit('[ERROR] CT input files are not the same dimensions.') + if ct1.GetSize() != ct2.GetSize(): + os.sys.exit('[ERROR] CT input files are not the same sizes.') + if ct1.GetPixelID() != sitk.sitkInt16: + ogo.message('[WARNING] Unexpected CT data type: {}'.format(ct1.GetPixelIDTypeAsString())) + if ct2.GetPixelID() != sitk.sitkInt16: + ogo.message('[WARNING] Unexpected CT data type: {}'.format(ct2.GetPixelIDTypeAsString())) + if (ct1.GetOrigin() != ct2.GetOrigin() or + ct1.GetSpacing() != ct2.GetSpacing() or + ct1.GetDirection() != ct2.GetDirection()): + ogo.message('[WARNING] Origin, Spacing or Direction are different for input images.') + ogo.message(' Output settings are taken from image 1 input.') + + origin = ct1.GetOrigin() + spacing = ct1.GetSpacing() + direction = ct1.GetDirection() + + if not quiet: + odir = os.path.dirname(ofiles[0]) + if not odir: odir='.' + ogo.message('Output files') + for idx, ofile in enumerate(ofiles): + ogo.message('{:>15s}: {:s}'.format(materials[idx],os.path.basename(ofile))) + ogo.message('{:>15s}: {:s}'.format('DIR',odir)) + ogo.message('') + + # Convert SimpleITK images to numpy arrays and flatten them to 1D + np_low = sitk.GetArrayFromImage(ct1) + np_low_flattened=np_low.reshape(-1) + + np_high = sitk.GetArrayFromImage(ct2) + np_high_flattened=np_high.reshape(-1) + + original_shape = np_low.shape + + # Convert images from CT number to linear attenuation + mu_water_low = md.interpolate_mass_attenuation('water',energies[0]) + mu_water_high = md.interpolate_mass_attenuation('water',energies[1]) + + np_low_flattened_mu = md.ct_number_to_mu(np_low_flattened,mu_water_low) + np_high_flattened_mu = md.ct_number_to_mu(np_high_flattened,mu_water_high) + + # 2-material decomposition ------------------------------------------------------------------------ + if len(materials)==2: + if not quiet: + ogo.message('Performing 2-material decomposition without calibration phantom.') + ogo.message(' Materials are: [' + ', '.join('{}'.format(i) for i in materials) + ']') + ogo.message(' No calibration phantom used.') + ogo.message('') + + # Ax = b, where + # + # ⎡ a b ⎤⎡ rho_mat1 ⎤ = ⎡ e ⎤ + # ⎣ c d ⎦⎣ rho_mat2 ⎦ ⎣ f ⎦ + + A_matrix = np.array(mass_attenuations).reshape(2, 2) + + if not quiet: + condition_A = np.linalg.cond(A_matrix) + ogo.message(' Condition of A: {:12.4e}'.format(condition_A)) + if condition_A>10: + ogo.message(' [WARNING] Condition of A is high!') + ogo.message('') + + B_matrix = np.vstack((np_low_flattened_mu,np_high_flattened_mu)) + + solution = md.solve_system_equations(A_matrix,B_matrix,True) + + # Reshape + np_ct_material1 = solution[0,:].reshape(original_shape) + np_ct_material2 = solution[1,:].reshape(original_shape) + + # Create SimpleITK images + ct_material1 = sitk.GetImageFromArray(np_ct_material1) + ct_material1.SetOrigin(origin) + ct_material1.SetSpacing(spacing) + ct_material1.SetDirection(direction) + + ct_material2 = sitk.GetImageFromArray(np_ct_material2) + ct_material2.SetOrigin(origin) + ct_material2.SetSpacing(spacing) + ct_material2.SetDirection(direction) + + # Write output + if not quiet: + ogo.message('Writing output') + if ofiles[0] is not no_print_str: + if not quiet: + ogo.message(' {}'.format(ofiles[0])) + sitk.WriteImage(ct_material1, ofiles[0]) + if ofiles[1] is not no_print_str: + if not quiet: + ogo.message(' {}'.format(ofiles[1])) + sitk.WriteImage(ct_material2, ofiles[1]) + + # 3-material decomposition ------------------------------------------------------------------------ + if len(materials)==3: + if not quiet: + ogo.message('Performing 3-material decomposition without calibration phantom.') + ogo.message(' Materials are: [' + ', '.join('{}'.format(i) for i in materials) + ']') + ogo.message(' No calibration phantom used.') + ogo.message('') + + # Ax = b, where + # + # ⎡ a b c ⎤⎡ rho_mat1 ⎤ ⎡ j ⎤ + # ⎜ d e f ⎟⎜ rho_mat2 ⎟ = ⎜ k ⎟ + # ⎣ g h i ⎦⎣ rho_mat3 ⎦ ⎣ 1 ⎦ + + A_matrix = np.array(mass_attenuations).reshape(2, 3) # Mass attenuations + eq3 = np.array([[1.0/mass_densities[0], 1.0/mass_densities[1], 1.0/mass_densities[2]]]) # Inverse mass density + A_matrix = np.append(A_matrix, eq3, axis=0) + + if not quiet: + condition_A = np.linalg.cond(A_matrix) + ogo.message(' Condition of A: {:12.4e}'.format(condition_A)) + if condition_A>10: + ogo.message(' [WARNING] Condition of A is high!') + ogo.message('') + + B_matrix = np.vstack((np_low_flattened_mu,np_high_flattened_mu,np.ones(np_low_flattened_mu.size))) # Add 1's for third equation + + solution = md.solve_system_equations(A_matrix,B_matrix,True) + + # Reshape + np_ct_material1 = solution[0,:].reshape(original_shape) + np_ct_material2 = solution[1,:].reshape(original_shape) + np_ct_material3 = solution[2,:].reshape(original_shape) + + # Create SimpleITK images + ct_material1 = sitk.GetImageFromArray(np_ct_material1) + ct_material1.SetOrigin(origin) + ct_material1.SetSpacing(spacing) + ct_material1.SetDirection(direction) + + ct_material2 = sitk.GetImageFromArray(np_ct_material2) + ct_material2.SetOrigin(origin) + ct_material2.SetSpacing(spacing) + ct_material2.SetDirection(direction) + + ct_material3 = sitk.GetImageFromArray(np_ct_material3) + ct_material3.SetOrigin(origin) + ct_material3.SetSpacing(spacing) + ct_material3.SetDirection(direction) + + # Write output + if not quiet: + ogo.message('Writing output') + if ofiles[0] is not no_print_str: + if not quiet: + ogo.message(' {}'.format(ofiles[0])) + sitk.WriteImage(ct_material1, ofiles[0]) + if ofiles[1] is not no_print_str: + if not quiet: + ogo.message(' {}'.format(ofiles[1])) + sitk.WriteImage(ct_material2, ofiles[1]) + if ofiles[2] is not no_print_str: + if not quiet: + ogo.message(' {}'.format(ofiles[2])) + sitk.WriteImage(ct_material3, ofiles[2]) + + if not quiet: + ogo.message('') + ogo.message('Done.') + + +def main(): + description = ''' + Performs image-space material decomposition based on a pair of + images representing LOW and HIGH acquisitions. These can be + from VMIs or any other LOW/HIGH pair. + + Typical use is 2-material decomposition into water and HA or a + 3-material decomposition into HA, water and adipose. However + any combinations can be employed, although do this at your own + risk. + + Input images are expected to be CT numbers in Hounsfield units. + + Valid ICRU materials are: + adipose + air + blood + bone + calcium + cha + k2hpo4 + muscle + water + softtissue + redmarrow + yellowmarrow + spongiosa + iodine + + An output image is generated for two or three materials + depending on type type of composition. The default output + filename has the suffix decomposition type (2 or 3) and + material added to it. For example, _2cha, _2water or possibly + _3cha, _3water, _3adipose. + + Use the PATTERN option if a new directory and base filename + is preferred, otherwise output is written to same directory + as input images and same base filename with suffix. + + Suppress writing an output image by specifying which + material(s). + + For 3-material decomposition that includes adipose we can + set a limit because of the effects of noise. Using Liu et + al. Medical Physics, 2009 we can choose 600 mg/cm3. + +''' + epilog = ''' +Example calls: + ogoMaterialDecomposition im1.nii im2.nii 40 80 cha water + + ogoMaterialDecomposition \\ + PCD_HIPQC_120kV_060KEV_LOW_Qr40_04Th_10Slices.nii.gz \\ + PCD_HIPQC_120kV_084KEV_HGH_Qr40_04Th_10Slices.nii.gz \\ + 60 84 cha water adipose --pattern PCD_HIPQC_120kV.nii.gz + + ogoMaterialDecomposition \\ + PCD_HIPQC_120kV_060KEV_LOW_Qr40_04Th_10Slices.nii.gz \\ + PCD_HIPQC_120kV_084KEV_HGH_Qr40_04Th_10Slices.nii.gz \\ + 60 84 cha water --pattern PCD_HIPQC_120kV.nii.gz + + ogoMaterialDecomposition \\ + PCD_ESP_140kV_040keV_Br40_04Th_10Slices.nii.gz \\ + PCD_ESP_140kV_090keV_Br40_04Th_10Slices.nii.gz \\ + 40 90 cha water --pattern ./PCD_ESP_140kV.nii.gz + + ogoMaterialDecomposition \\ + PCD_ESP_140kV_040keV_Br40_04Th_10Slices.nii.gz \\ + PCD_ESP_140kV_090keV_Br40_04Th_10Slices.nii.gz \\ + 40 90 cha water adipose --pattern ./PCD_ESP_140kV.nii.gz + +''' + + # Setup argument parsing + parser = argparse.ArgumentParser( + formatter_class=argparse.RawTextHelpFormatter, + prog="ogoMaterialDecomposition", + description=description, + epilog=epilog + ) + parser.add_argument('images', + metavar='Filename Filename', + nargs=2, + help='CT low and high images (*.nii.gz)') + parser.add_argument('energies', + type=float, + default=[], + nargs=2, + metavar='keV keV', + help='Energies of low and high images (40 to 190 keV)') + parser.add_argument('materials', + nargs='*', + default=[], + metavar='MAT1 MAT2 (MAT3)', + help='Specify 2 or 3 ICRU materials (eg. cha, water, adipose)') + parser.add_argument("--pattern", + default=None, + metavar='FILE', + help='Directory and filename pattern.') + parser.add_argument('--suppress', + nargs='*', + default=[], + metavar='MAT1 MAT2 (MAT3)', + help='File outputs to suppress (e.g., HA, water, adipose)') + #parser.add_argument('--max_adipose', + # type=float, + # default=[], + # nargs=1, + # metavar='mg/cm3', + # help='Maximum allowable adipose density (df. 600 mg/cm3)') + parser.add_argument("--quiet", + action='store_true', + help='Overwrite output without asking') + parser.add_argument("--overwrite", + action='store_true', + default=False, + help='Overwrite output files without asking (default: %(default)s)') + # Parse and display + args = parser.parse_args() + if not args.quiet: + print(echo_arguments('MaterialDecomposition', vars(args))) + + # Run program + MaterialDecomposition(**vars(args)) + +if __name__ == '__main__': + main() diff --git a/ogo/cli/MorphologicalOperation.py b/ogo/cli/MorphologicalOperation.py index 65c309c..2c97207 100644 --- a/ogo/cli/MorphologicalOperation.py +++ b/ogo/cli/MorphologicalOperation.py @@ -5,162 +5,363 @@ # | bonelab@ucalgary.ca | # +------------------------------------------------------------------------------+ -script_version = 1.0 - # Imports import argparse -import os -import sys -import vtk -import math -import numpy as np import SimpleITK as sitk -from scipy.spatial import procrustes -from datetime import date -from vtk.util.numpy_support import vtk_to_numpy, numpy_to_vtk from ogo.util.echo_arguments import echo_arguments import ogo.util.Helper as ogo -import ogo.dat.OgoMasterLabels as lb +from ogo.util.cli_tools import ( + check_overwrite, + validate_input_file, + validate_output_file, + validate_integer_range, + validate_integer_min, + load_config_from_json, + save_config_to_json, + print_image_info +) -def get_labels(ct): - filt = sitk.LabelShapeStatisticsImageFilter() - filt.Execute(ct) - labels = filt.GetLabels() - return labels -# +------------------------------------------------------------------------------+ -def MorphologicalOperation(input_image, output_image, operation_labels, operation, kernel, overwrite, kernels=None): +def MorphologicalOperation(input_file, output_file, operation_labels, dilate, erode, morphological, overwrite, quiet): + # Check if output exists and should overwrite - if os.path.isfile(output_image) and not overwrite: - result = input('File \"{}\" already exists. Overwrite? [y/n]: '.format(output_image)) - if result.lower() not in ['y', 'yes']: - ogo.message('Not overwriting. Exiting...') - os.sys.exit() - + if not check_overwrite(output_file, overwrite): + return 1 + + # Validate dilate and erode parameters + if not dilate: + dilate = [] + if not erode: + erode = [] + + # Check that at least one operation is specified + if not dilate and not erode: + ogo.message('[ERROR] At least one of --dilate or --erode must be specified.') + return 1 + + # Validate all dilate values are positive + for d in dilate: + if d <= 0: + ogo.message(f'[ERROR] All dilate kernel values must be positive. Found: {d}') + return 1 + + # Validate all erode values are positive + for e in erode: + if e <= 0: + ogo.message(f'[ERROR] All erode kernel values must be positive. Found: {e}') + return 1 + # Read input - if not os.path.isfile(input_image): - os.sys.exit('[ERROR] Cannot find file \"{}\"'.format(input_image)) - - if not (input_image.lower().endswith('.nii') or input_image.lower().endswith('.nii.gz')): - os.sys.exit('[ERROR] Input must be type NIFTI file: \"{}\"'.format(input_image)) + if not quiet: + ogo.message(f'Reading input image:') + ogo.message(f' {input_file}') - ogo.message('Reading image: ') - ogo.message('\"{}\"'.format(input_image)) - ct = sitk.ReadImage(input_image, sitk.sitkUInt8) + input_image = sitk.ReadImage(input_file, sitk.sitkUInt8) - # Get all the available labels in the input image - labels = get_labels(ct) - n_labels = len(labels) - ogo.message('Input image contains the following labels:') - ogo.message(' [' + ', '.join('{:d}'.format(i) for i in labels) + ']') + if not quiet: + print_image_info(input_file, input_image, brief=True) + + # Get all available labels in the input image + label_filter = sitk.LabelShapeStatisticsImageFilter() + label_filter.Execute(input_image) + available_labels = list(label_filter.GetLabels()) + + if not quiet: + ogo.message(f'Available labels in input image: {available_labels}') # Set the labels on which we will perform morphological operations if not operation_labels: - operation_labels = labels + operation_labels = available_labels + if not quiet: + ogo.message('No operation_labels specified - will process all labels') else: - for label in operation_labels: - if label not in labels: - os.sys.exit('[ERROR] Label {} is not in input file.'.format(label)) - - # Validate kernels if mulitple provided - if kernels: - if len(kernels) != len(operation_labels): - os.sys.exit('[ERROR] Number of kernels must match number of operation labels.') + # Filter operation_labels to only include labels that exist in the image + original_operation_labels = operation_labels.copy() + valid_operation_labels = [] + invalid_labels = [] + + # Track which indices have valid labels (for per-label kernels) + valid_indices = [] + + for idx, label in enumerate(original_operation_labels): + if label in available_labels: + valid_operation_labels.append(label) + valid_indices.append(idx) + else: + invalid_labels.append(label) + + # Warn about invalid labels + if invalid_labels: + if not quiet: + ogo.message(f'[WARNING] Label(s) {invalid_labels} not found in input image and will be skipped.') + ogo.message(f' Available labels: {available_labels}') + + # Update operation_labels to only valid labels + operation_labels = valid_operation_labels + + # If no valid labels remain, nothing to do + if not operation_labels: + ogo.message('[ERROR] None of the specified operation_labels exist in the input image.') + ogo.message(f' Available labels: {available_labels}') + return 1 + + # Filter dilate and erode to match valid labels if they were specified per-label + if len(dilate) == len(original_operation_labels): + dilate = [dilate[i] for i in valid_indices] + + if len(erode) == len(original_operation_labels): + erode = [erode[i] for i in valid_indices] + + # Expand dilate and erode to match number of labels + num_labels = len(operation_labels) + + if len(dilate) == 1: + dilate_kernels = dilate * num_labels + elif len(dilate) == num_labels: + dilate_kernels = dilate + elif len(dilate) == 0: + dilate_kernels = [0] * num_labels else: - kernels = [kernel] * len(operation_labels) - - # Set the morphological operation - ogo.message('Setting morphological operation to: [{}]'.format(operation)) - if operation == 'dilate': - morphFilt = sitk.BinaryDilateImageFilter() - elif operation == 'erode': - morphFilt = sitk.BinaryErodeImageFilter() - elif operation == 'opening': - morphFilt = sitk.BinaryMorphologicalOpeningImageFilter() - elif operation == 'closing': - morphFilt = sitk.BinaryMorphologicalClosingImageFilter() + ogo.message('[ERROR] Number of --dilate values must be either 1 or match number of operation labels.') + ogo.message(f' operation_labels: {num_labels}, dilate values: {len(dilate)}') + return 1 + + if len(erode) == 1: + erode_kernels = erode * num_labels + elif len(erode) == num_labels: + erode_kernels = erode + elif len(erode) == 0: + erode_kernels = [0] * num_labels else: - os.sys.exit('[ERROR] Unknown morphological operation: {}'.format(operation)) + ogo.message('[ERROR] Number of --erode values must be either 1 or match number of operation labels.') + ogo.message(f' operation_labels: {num_labels}, erode values: {len(erode)}') + return 1 + + # Warn if both dilate and erode are specified without explicit morphological parameter + if dilate and erode and morphological is None: + if not quiet: + ogo.message('[WARNING] Both --dilate and --erode specified. Dilation will be applied first, then erosion.') + ogo.message(' This is equivalent to a morphological closing operation.') + ogo.message(' To explicitly control operation order, use --morphological opening or --morphological closing.') - # Cycle through the labels - seg = ct<0 # Zero all labels + # Warn if morphological operation selected but parameters not aligned + if morphological == 'opening' and not erode: + if not quiet: + ogo.message('[WARNING] Morphological opening selected but --erode not specified.') + ogo.message(' Only dilation will be applied.') + + if morphological == 'closing' and not dilate: + if not quiet: + ogo.message('[WARNING] Morphological closing selected but --dilate not specified.') + ogo.message(' Only erosion will be applied.') + + # Display operation info + if not quiet: + if morphological == 'opening': + ogo.message('Morphological operation: opening (erode then dilate)') + elif morphological == 'closing': + ogo.message('Morphological operation: closing (dilate then erode)') + else: + ogo.message('Morphological operation: dilate then erode (default)') + + # Initialize output as empty segmentation + output_seg = input_image < 0 # All zeros + + # Process each label + for label in available_labels: + # Create binary mask for this label + label_mask = input_image == label - for idx, label in enumerate(labels): - ct_thres = ct==label if label in operation_labels: - kernel_idx = operation_labels.index(label) - current_kernel = kernels[kernel_idx] + # Find the kernels for this label + label_idx = operation_labels.index(label) + current_dilate = dilate_kernels[label_idx] + current_erode = erode_kernels[label_idx] - if current_kernel <= 0: - os.sys.exit('[ERROR] Kernel must be positive: {}'.format(current_kernel)) + # Apply morphological operations based on morphological parameter + processed_mask = label_mask - morphFilt.SetKernelRadius(current_kernel) - ogo.message(' label {:3d}: {:s}, kernel={:d}'.format(label, operation, current_kernel)) - ct_part = morphFilt.Execute(ct_thres) - else: - ogo.message(' label {:3d}: {:s}'.format(label, 'skip')) - ct_part = ct_thres + if morphological == 'opening': + # Opening: erode then dilate + if current_erode > 0: + erode_filter = sitk.BinaryErodeImageFilter() + erode_filter.SetKernelRadius(current_erode) + processed_mask = erode_filter.Execute(processed_mask) + + if current_dilate > 0: + dilate_filter = sitk.BinaryDilateImageFilter() + dilate_filter.SetKernelRadius(current_dilate) + processed_mask = dilate_filter.Execute(processed_mask) + + if not quiet: + ogo.message(f' Label {label:3d}: opening (erode={current_erode}, dilate={current_dilate})') + + elif morphological == 'closing': + # Closing: dilate then erode + if current_dilate > 0: + dilate_filter = sitk.BinaryDilateImageFilter() + dilate_filter.SetKernelRadius(current_dilate) + processed_mask = dilate_filter.Execute(processed_mask) + + if current_erode > 0: + erode_filter = sitk.BinaryErodeImageFilter() + erode_filter.SetKernelRadius(current_erode) + processed_mask = erode_filter.Execute(processed_mask) + + if not quiet: + ogo.message(f' Label {label:3d}: closing (dilate={current_dilate}, erode={current_erode})') - bin_seg = seg>0 - bin_part = ct_part>0 - overlap = (bin_seg + bin_part)==2 - mask = 1 - overlap - this_label = sitk.Mask(ct_part, mask) - seg = seg + label*(this_label>0) + else: + # Default: dilate then erode + if current_dilate > 0: + dilate_filter = sitk.BinaryDilateImageFilter() + dilate_filter.SetKernelRadius(current_dilate) + processed_mask = dilate_filter.Execute(processed_mask) + + if current_erode > 0: + erode_filter = sitk.BinaryErodeImageFilter() + erode_filter.SetKernelRadius(current_erode) + processed_mask = erode_filter.Execute(processed_mask) + + if not quiet: + ogo.message(f' Label {label:3d}: dilate={current_dilate}, erode={current_erode}') + + morphed_mask = processed_mask + else: + # No operation for this label - keep as is + morphed_mask = label_mask + if not quiet: + ogo.message(f' Label {label:3d}: unchanged') + + # Handle overlaps - existing labels have priority + existing_seg = output_seg > 0 + new_label_mask = morphed_mask > 0 - ct_final = seg + # Cast to integers and find overlap + existing_seg_int = sitk.Cast(existing_seg, sitk.sitkUInt8) + new_label_mask_int = sitk.Cast(new_label_mask, sitk.sitkUInt8) + overlap = (existing_seg_int + new_label_mask_int) == 2 + + # Mask out overlapping regions from new label + no_overlap_mask = sitk.Cast(overlap == 0, sitk.sitkUInt8) + this_label_clean = sitk.Mask(morphed_mask, no_overlap_mask) + + # Add this label to the output + output_seg = output_seg + sitk.Cast(label * (this_label_clean > 0), sitk.sitkUInt8) + + # Write output + if not quiet: + ogo.message(f'Writing output to: {output_file}') - ogo.message('Writing output file:') - ogo.message(' {}'.format(output_image)) - sitk.WriteImage(ct_final, output_image) + sitk.WriteImage(output_seg, output_file) - ogo.message('Done.') + if not quiet: + print_image_info(output_file, output_seg, brief=True) + ogo.message('[DONE]') + return 0 + def main(): description = ''' -Function performs morphological operations. Although it is possible to -apply this to images, it is meant for application with labels. The operations -include: - -- dilate -- erode -- closing (erosion of the dilation) -- opening (dilation of the erosion) - -A list of specific labels can be defined to restrict the operation to those -labels only. The labels are processed in the order they are defined in the -list. +Perform morphological operations on labeled images. + +This tool applies binary morphological operations (dilation and erosion) to specific +labels in a segmented image. Each label is processed independently, and overlaps are +resolved by giving priority to labels processed first (in the order they appear in the image). + +MORPHOLOGICAL OPERATIONS: + - opening: Erode then dilate (removes small protrusions, smooths from outside) + - closing: Dilate then erode (fills small holes, smooths from inside) + - Default: Dilate then erode (if both --dilate and --erode specified) + +DILATE AND ERODE: + Specify kernel radii for dilation and erosion operations: + - Single value: applies to all processed labels + - Multiple values: must match number of operation_labels (one kernel per label) + - Kernel radius determines structuring element size (1=3x3x3, 2=5x5x5, etc.) + +LABEL PROCESSING: + - If no operation_labels specified, all labels in the image are processed + - Labels are processed in the order they appear in the image + - Overlaps are resolved by keeping the first label (earlier labels have priority) + +JSON CONFIGURATION: + - To output current arguments as JSON: --json + - To load arguments from JSON file: --json my_config.json ''' epilog = ''' -Example call: - -ogoMorphologicalOperation image_in.nii.gz image_out.nii.gz --operation erode --operation_labels 1 2 3 --kernels 2 3 4 -ogoMorphologicalOperation image_in.nii.gz image_out.nii.gz --operation erode --operation_labels 1 2 3 --kernel 3 - +Example calls: +ogoMorphologicalOperation input.nii.gz output.nii.gz --dilate 2 +ogoMorphologicalOperation input.nii.gz output.nii.gz --erode 3 --operation_labels 1 2 3 +ogoMorphologicalOperation input.nii.gz output.nii.gz --dilate 2 --erode 1 --morphological closing +ogoMorphologicalOperation input.nii.gz output.nii.gz --dilate 2 3 4 --erode 1 1 1 --operation_labels 1 2 3 +ogoMorphologicalOperation input.nii.gz output.nii.gz --dilate 2 --morphological opening --json > config.json +ogoMorphologicalOperation --json config.json ''' + parser = argparse.ArgumentParser( formatter_class=argparse.RawTextHelpFormatter, prog="ogoMorphologicalOperation", description=description, epilog=epilog ) + + # JSON configuration argument + parser.add_argument('--json', nargs='?', const='', metavar='FILE', + help='Output arguments as JSON to stdout (no file), or load arguments from JSON file') + + parser.add_argument('input_file', nargs='?', type=validate_input_file(['.nii', '.nii.gz']), + help='Input label image file (*.nii, *.nii.gz)') + parser.add_argument('output_file', nargs='?', type=validate_output_file(['.nii', '.nii.gz']), + help='Output label image file (*.nii, *.nii.gz)') + parser.add_argument('--operation_labels','-l', type=validate_integer_range(0, 255), nargs='*', default=[], + metavar='LABEL', help='List of labels to process (default: all labels in image)') + parser.add_argument('--dilate', '-d', type=validate_integer_min(1), nargs='*', default=[], metavar='KERNEL', + help='Kernel radius/radii for dilation (single value or one per label)') + parser.add_argument('--erode', '-e', type=validate_integer_min(1), nargs='*', default=[], metavar='KERNEL', + help='Kernel radius/radii for erosion (single value or one per label)') + parser.add_argument('--morphological', '-mo', choices=['opening', 'closing'], default=None, metavar='TYPE', + help='Morphological operation type: opening (erode then dilate) or closing (dilate then erode)') + parser.add_argument('--overwrite', '-ow', action='store_true', + help='Overwrite output without asking') + parser.add_argument('--quiet', '-q', action='store_true', + help='Suppress informational output to stdout') - parser.add_argument('input_image', help='Input image file (*.nii, *.nii.gz)') - parser.add_argument('output_image', help='Output image file (*.nii, *.nii.gz)') - parser.add_argument('--operation_labels', type=int, nargs='*', default=[], metavar='LABEL', help='List of labels to apply operation (default: all)') - parser.add_argument('--operation', default='dilate', choices=['dilate', 'erode', 'opening', 'closing'], - help='Select morphological operation (default: %(default)s)') - parser.add_argument('--kernel', type=int, default=1, metavar='KERNEL', help='Default kernel for morphological operation (default: %(default)s)') - parser.add_argument('--kernels', type=int, nargs='*', metavar='KERNELS', help='List of kernels for each label (overrides --kernel)') - parser.add_argument('--overwrite', action='store_true', help='Overwrite output file without asking') - + # Parse and display args = parser.parse_args() - print(echo_arguments('MorphologicalOperation', vars(args))) + + # Handle --json flag + if args.json is not None: + if args.json == '': + # --json with no argument: output to stdout + save_config_to_json(args, None) + return 0 + else: + # --json file.json: load from file + if not args.json.endswith('.json'): + parser.error('JSON file must have .json extension') + json_config = load_config_from_json(args.json) + # Update args with JSON values (command line args override JSON) + for key, value in json_config.items(): + if hasattr(args, key): + # Only update if argument was not explicitly provided on command line + if getattr(args, key) == parser.get_default(key): + setattr(args, key, value) + + # Check that required positional arguments are provided + if args.input_file is None or args.output_file is None: + parser.error('input_file and output_file are required') + + if not args.quiet: + ogo.message(echo_arguments('MorphologicalOperation', vars(args))) - # Run program - MorphologicalOperation(**vars(args)) + # Run program (exclude json argument) + run_args = vars(args).copy() + run_args.pop('json', None) + result = MorphologicalOperation(**run_args) + return result if __name__ == '__main__': diff --git a/ogo/cli/StyleGuide.py b/ogo/cli/StyleGuide.py new file mode 100644 index 0000000..1b54fc4 --- /dev/null +++ b/ogo/cli/StyleGuide.py @@ -0,0 +1,220 @@ +# /------------------------------------------------------------------------------+ +# | 29-MAR-2026 | +# | Copyright (c) Bone Imaging Laboratory | +# | All rights reserved | +# | bonelab@ucalgary.ca | +# +------------------------------------------------------------------------------+ + +# Imports +import argparse +import SimpleITK as sitk +from ogo.util.echo_arguments import echo_arguments +import ogo.util.Helper as ogo +from ogo.util.cli_tools import ( + check_overwrite, + check_label_lengths, + check_image_dimensions_match, + validate_integer_range, + validate_integer_min, + validate_float_range, + validate_float_min, + validate_input_file, + validate_output_file, + load_config_from_json, + save_config_to_json, + print_image_info +) + +#------------------------------------------------------------------------------+ +# Main function +def StyleGuide(input_file, mask_file, output_file, from_labels, to_labels, labels_to_keep, collection, threshold, slice_percent, bins, header, overwrite, quiet): + + # Check if output exists and should overwrite + if not check_overwrite(output_file, overwrite): + return 1 + + # Check that from_labels and to_labels have the same length + if not check_label_lengths(from_labels, to_labels): + return 1 + + # Check that at least some labels are defined + if not from_labels and not to_labels and not labels_to_keep: + ogo.message('ERROR: No labels defined.') + ogo.message(' Define --labels_to_keep or --from_labels & --to_labels') + ogo.message(' There is nothing to do.') + return 1 + + # Check if threshold was explicitly set + if threshold is not None: + ogo.message(f' user set threshold: {threshold}') + else: + threshold = 1.0 # Set default value programmatically + + # Check if slice_percent was explicitly set + if slice_percent is not None: + ogo.message(f' user set slice_percent: {slice_percent}') + else: + slice_percent = 1.0 # Set default value programmatically + + # Check if bins was explicitly set + if bins is not None: + ogo.message(f' user set bins: {bins}') + else: + bins = 10 # Set default value programmatically + + # Display label information + if not quiet: + ogo.message(f' from_labels: {from_labels}') + ogo.message(f' to_labels: {to_labels}') + ogo.message(f' labels_to_keep: {labels_to_keep}') + ogo.message(f' threshold: {threshold}') + ogo.message(f' slice_percent: {slice_percent}') + ogo.message(f' bins: {bins}') + + # Read input file based on extension + if not quiet: + ogo.message(f'Reading input file: {input_file}') + + input_image = sitk.ReadImage(input_file) + + if not quiet: + print_image_info(input_file, input_image) + + # Read mask file based on extension + if not quiet: + ogo.message(f'Reading mask file: {mask_file}') + + mask_image = sitk.ReadImage(mask_file, sitk.sitkUInt8) + + if not quiet: + print_image_info(mask_file, mask_image) + + # Check that images have matching dimensions + if not check_image_dimensions_match(input_image, mask_image, "input_file", "mask_file"): + return 1 + + # Function logic goes here + + if not quiet: + ogo.message('[DONE]') + + return 0 + +def main(): + # Setup description + description=''' +This program demonstrates the coding styles and best practices for Ogo CLI tools. + +Basically it does nothing... but it shows how to set up a CLI with argparse, +how to validate inputs, and how to structure the code for readability and +maintainability. + +STYLE GUIDE RULES: + 1. Use the TYPE option in argparse to check valid inputs whenever possible. + 2. Minimize use of positional arguments for mandatory inputs only. + 3. Use flags with long names, but include short names in addition if needed (see below). + 4. Return 0 for successful completion, return 1 for early exit. + 5. Formatting: + – Use 'snake_case' for variables and function names (e.g., input_file). + – Use 'PascalCase' for class names (e.g., ImageProcessor). + - Use 'UPPER_SNAKE_CASE' for constants (e.g., MAX_ITERATIONS). + - Do NOT use 'camelCase' (e.g., inputLabels, outputLabels). + +EXAMPLE VARIABLE NAMES: + input_file / output_file / mask_file / calib_file + input_path / output_path + from_labels / to_labels / labels_to_keep + overwrite + +COMMON ARGUMENT OPTIONS: + -h, --help Shows user information about CLI + -ow, --overwrite Overwrite output if it exists + -o, --output Output file (e.g., sort, gcc) + -q, --quiet Display less output (useful for scripts) + -d, --debug Show debugging output (verbose) + -f, --force Force destructive actions without confirmation + --json Display JSON output + -n, --dry_run Describe changes without executing + -a, --all All (e.g., ps, fetchmail) + -u, --user User (e.g., ps, ssh) + -v, --version Version + +We can use JSON files instead of command line arguments for more complex configurations. + - To output current arguments as JSON: --json + - To load arguments from JSON file: --json my_config.json +''' + epilog=''' +Example calls: +ogoStyleGuide input.nii.gz output.nii.gz -fl 1 2 3 -tl 10 20 30 +ogoStyleGuide input.nii.gz output.nii.gz -lk 5 10 15 --quiet +ogoStyleGuide input.nii.gz output.nii.gz -lk 12 3 --json > my_config.json +ogoStyleGuide --json my_config.json +''' + + # Setup argument parsing + parser = argparse.ArgumentParser( + formatter_class=argparse.RawTextHelpFormatter, + prog="ogoStyleGuide", + description=description, + epilog=epilog + ) + + # JSON configuration argument + parser.add_argument('--json', nargs='?', const='', metavar='FILE', + help='Output arguments as JSON to stdout (no file), or load arguments from JSON file') + + parser.add_argument('input_file', nargs='?', type=validate_input_file(['.nii', '.nii.gz']), + help='Input image file (*.nii, *.nii.gz)') + parser.add_argument('mask_file', nargs='?', type=validate_input_file(['.nii', '.nii.gz']), + help='Mask image file (*.nii, *.nii.gz)') + parser.add_argument('output_file', nargs='?', type=validate_output_file(['.csv', '.txt']), + help='Output file (*.csv, *.txt)') + parser.add_argument('--from_labels','-fl', type=validate_integer_range(0, 255), nargs='*', default=[], metavar='ID', help='Target label input IDs; space separated (e.g. 1 2 3)') + parser.add_argument('--to_labels','-tl', type=validate_integer_range(0, 255), nargs='*', default=[], metavar='ID', help='Target label output IDs; space separated (e.g. 4 5 6)') + parser.add_argument('--labels_to_keep','-lk', type=validate_integer_range(0, 255), nargs='*', default=[], metavar='ID', help='Labels to keep (others removed); (e.g. 7 8 9)') + parser.add_argument('--collection', default='all', choices=['all', 'ossai', 'skeleton', 'cardio', 'organs', 'muscle', 'spine', 'gastro'], + help='Select collection (default: %(default)s)') + parser.add_argument('--threshold', '-t', type=validate_float_min(1.0), default=None, help='Threshold value (must be >= 10.0, default: 1.0)') + parser.add_argument('--slice_percent', '-s', type=validate_float_range(1.0, 5.0), default=None, help='Slice percentage (must be between 1.0 and 5.0, default: 1.0)') + parser.add_argument('--bins', '-b', type=validate_integer_min(1), default=None, help='Number of bins (must be >= 1, default: 10)') + parser.add_argument('--header','-H', action='store_true', help='Include header in output') + parser.add_argument('--overwrite','-ow', action='store_true', help='Overwrite output without asking') + parser.add_argument('--quiet','-q', action='store_true', help='Suppress informational output to stdout') + + print() + + # Parse and display + args = parser.parse_args() + + # Handle --json flag + if args.json is not None: + if args.json == '': + # --json with no argument: output to stdout + save_config_to_json(args, None) + return 0 + else: + # --json file.json: load from file + if not args.json.endswith('.json'): + parser.error('JSON file must have .json extension') + json_config = load_config_from_json(args.json) + # Update args with JSON values (command line args override JSON) + for key, value in json_config.items(): + if hasattr(args, key): + # Only update if argument was not explicitly provided on command line + if getattr(args, key) == parser.get_default(key): + setattr(args, key, value) + + # Check that required positional arguments are provided + if args.input_file is None or args.mask_file is None or args.output_file is None: + parser.error('input_file, mask_file, and output_file are required') + + if not args.quiet: + ogo.message(echo_arguments('StyleGuide', vars(args))) + + # Run program (exclude json argument) + run_args = vars(args).copy() + run_args.pop('json', None) + StyleGuide(**run_args) + +if __name__ == '__main__': + main() diff --git a/ogo/cli/Visualize.py b/ogo/cli/Visualize.py index d8859e1..0e7c445 100644 --- a/ogo/cli/Visualize.py +++ b/ogo/cli/Visualize.py @@ -305,7 +305,7 @@ def vis3d(input_filename, select, collection, gaussian, radius, isosurface, elev if offscreen: renderWindow.Finalize() -def vis2d(input_filename, outfile, mask_image, window, level, nThreads, image_orientation, slice_percent, offscreen, overwrite, func): +def vis2d(input_filename, outfile, mask_image, window, level, nThreads, image_orientation, slice_percent, offscreen, overwrite, show_labels, func): # Check if output exists and should overwrite if outfile: @@ -335,6 +335,12 @@ def vis2d(input_filename, outfile, mask_image, window, level, nThreads, image_or if not (input_filename.lower().endswith('.nii') or input_filename.lower().endswith('.nii.gz')): os.sys.exit('[ERROR] Input must be type NIFTI file: \"{}\"'.format(input_filename)) + # Check if show_labels is set without mask_image + if show_labels and mask_image is None: + ogo.message('[WARNING]: show_labels is only possible if mask_image is used.') + ogo.message(' Ignoring show_labels option.') + show_labels = False + ogo.message('Reading input image ' + input_filename) reader = vtk.vtkNIFTIImageReader() reader.SetFileName(input_filename) @@ -378,6 +384,9 @@ def vis2d(input_filename, outfile, mask_image, window, level, nThreads, image_or renderer = vtk.vtkRenderer() renderer.AddActor(inputSlice) + # Initialize labels list + labels = [] + # Read mask if not mask_image is None: if not os.path.isfile(mask_image): @@ -391,6 +400,34 @@ def vis2d(input_filename, outfile, mask_image, window, level, nThreads, image_or mask_reader.SetFileName(mask_image) mask_reader.Update() + ogo.message('Scanning for available labels:') + labels = getLabels(mask_reader.GetOutput()) + + # Get image data for centroid calculation for placing label text onto the rendered image + mask_output = mask_reader.GetOutput() + dims = mask_output.GetDimensions() + spacing = mask_output.GetSpacing() + origin = mask_output.GetOrigin() + mask_array = vtk_to_numpy(mask_output.GetPointData().GetScalars()).reshape(dims, order='F') + + for idx,label in enumerate(labels): + rgb = lb.labels_dict[label]['RGB'] + desc = lb.labels_dict[label]['LABEL'] + + # Calculate centroid for this label + label_indices = np.argwhere(mask_array == label) + if len(label_indices) > 0: + # Compute centroid in index space + centroid_ijk = np.mean(label_indices, axis=0) + # Convert to physical coordinates + centroid_xyz = origin + centroid_ijk * spacing + # Store in labels_dict + lb.labels_dict[label]['CENTROID'] = centroid_xyz + ogo.message('{:5d}: ({:3d},{:3d},{:3d}) – {} – Centroid: ({:.1f}, {:.1f}, {:.1f})'.format( + label, rgb[0], rgb[1], rgb[2], desc, centroid_xyz[0], centroid_xyz[1], centroid_xyz[2])) + else: + ogo.message('{:5d}: ({:3d},{:3d},{:3d}) – {}'.format(label, rgb[0], rgb[1], rgb[2], desc)) + lut = get_lookup_table() mask_inputMapper = vtk.vtkOpenGLImageSliceMapper() @@ -404,7 +441,7 @@ def vis2d(input_filename, outfile, mask_image, window, level, nThreads, image_or mask_imageProperty = vtk.vtkImageProperty() mask_imageProperty.SetInterpolationTypeToNearest() mask_imageProperty.SetLookupTable(lut) - mask_imageProperty.SetOpacity(0.5) + mask_imageProperty.SetOpacity(0.25) mask_inputSlice = vtk.vtkImageSlice() mask_inputSlice.SetMapper(mask_inputMapper) @@ -438,6 +475,37 @@ def vis2d(input_filename, outfile, mask_image, window, level, nThreads, image_or os.sys.exit() renderer.ResetCamera() + # Add text labels at centroids if show_labels is True + if show_labels and mask_image is not None: + for idx, label in enumerate(labels): + if 'CENTROID' in lb.labels_dict[label]: + centroid = lb.labels_dict[label]['CENTROID'] + desc = lb.labels_dict[label]['LABEL'] + rgb = lb.labels_dict[label]['RGB'] + + # Transform centroid by the NIfTI transformation matrix + point = [centroid[0], centroid[1], centroid[2], 1.0] + transformed_point = [0.0, 0.0, 0.0, 0.0] + mat.MultiplyPoint(point, transformed_point) + + # Create billboard text actor for 3D position that always faces camera + textActor = vtk.vtkBillboardTextActor3D() + textActor.SetInput(desc) + textActor.SetPosition(transformed_point[0], transformed_point[1], transformed_point[2]) + + # Set text properties for readability + textProperty = textActor.GetTextProperty() + textProperty.SetFontSize(18) + textProperty.SetColor(rgb[0]/255.0, rgb[1]/255.0, rgb[2]/255.0) + textProperty.SetBold(True) + textProperty.SetFontFamilyToArial() + textProperty.SetJustificationToCentered() + textProperty.SetVerticalJustificationToCentered() + textProperty.SetBackgroundColor(0, 0, 0) + textProperty.SetBackgroundOpacity(0.5) + + renderer.AddActor(textActor) + # Set slice ogo.message('Percent slice set to {:.1f}%'.format(slice_percent)) camera = renderer.GetActiveCamera() @@ -584,6 +652,7 @@ def main(): parser_vis2d.add_argument('--image_orientation', default='coronal', choices=['sagittal', 'coronal', 'axial'], help='Initiate a particular orientation (default: %(default)s)') parser_vis2d.add_argument('--slice_percent', default=float(50), type=float, help='Set percent slice through image (default: %(default)s)') + parser_vis2d.add_argument('--show_labels', action='store_true', help='Show labels (default: %(default)s)') parser_vis2d.add_argument('--outfile', default=None, metavar='FN', help='Output image file (*.tif) (default: %(default)s)') parser_vis2d.add_argument('--offscreen', action='store_true', help='Set to offscreen rendering (default: %(default)s)') parser_vis2d.add_argument('--overwrite', action='store_true', help='Overwrite output without asking') diff --git a/ogo/dat/MassAttenuationTables.py b/ogo/dat/MassAttenuationTables.py index 0b559ab..3ed2a04 100644 --- a/ogo/dat/MassAttenuationTables.py +++ b/ogo/dat/MassAttenuationTables.py @@ -5,7 +5,7 @@ # | bonelab@ucalgary.ca | # +------------------------------------------------------------------------------+ -'''Mass attenuation tables for internal calibration +r'''Mass attenuation tables for internal calibration The variable ``mass_attenuation_tables`` is defined containing: @@ -20,6 +20,17 @@ cha_table triglyceride_table water_table + ctwater_table + softtissue_table + redmarrow_table + yellowmarrow_table + spongiosa_table + iodine_table + calcium_table + liver_table + brain_table + iron_table + Each of these are pandas dataframe with two columns: @@ -31,19 +42,29 @@ Energy ranges from 1 KeV to 20 MeV but is not equality spaced for all materials. -The following material definitions came from [1]: +The following material definitions came from [1] unless otherwise specified: .. table:: Mass Attenuation Table =============== ===================== ========================== Ogo Mass Attenuation NIST =============== ===================== ========================== - Adipose Adipose Adipose Tissue (ICRU-44) + Adipose Adipose [2] XCOM: Adipose Tissue (ICRU-44) Air Air Air, Dry (near sea level) Blood Blood Blood, Whole (ICRU-44) - Cortical Bone Bone Bone, Cortical (ICRU-44) - Skeletal Muscle Muscle Muscle, Skeletal (ICRU-44) - Water Water Water, Liquid + Cortical Bone Bone [2] XCOM: Bone, Cortical (ICRU-44) + Skeletal Muscle Muscle [2] XCOM: Muscle, Skeletal (ICRU-44) + Water Water Water, Liquid (ICRU-44) + CTWater CTWater CTWater, as defined by QRM + Soft Tissue SoftTissue Tissue, Soft (ICRU-44) + Red Marrow RedMarrow Skeleton Red Marrow Adult (ICRU-46) + Yellow Marrow YellowMarrow Skeleton Yellow Marrow Adult (ICRU-46) + Spongiosa Spongiosa [2] XCOM: Skeleton Spongiosa Adult (ICRU-46) + Iodine Iodine Iodine + Calcium Calcium Calcium + Liver Liver [2] XCOM: Liver Adult Health (ICRU-46) + Brain Brain [2] XCOM: Brain Adult (ICRU-46) + Iron Iron [2] XCOM: Iron =============== ===================== ========================== Triglyceride, K\ :sub:`2`\ HPO\ :sub:`4`, and CHA were calculated from @@ -68,96 +89,393 @@ import numpy as np -adipose_data = np.array([[ - 1.00000E-03, - 1.03542E-03, - 1.07210E-03, - 1.07210E-03, - 1.50000E-03, - 2.00000E-03, - 2.47200E-03, - 2.47200E-03, - 2.64140E-03, - 2.82240E-03, - 2.82240E-03, - 3.00000E-03, - 4.00000E-03, - 5.00000E-03, - 6.00000E-03, - 8.00000E-03, - 1.00000E-02, - 1.50000E-02, - 2.00000E-02, - 3.00000E-02, - 4.00000E-02, - 5.00000E-02, - 6.00000E-02, - 8.00000E-02, - 1.00000E-01, - 1.50000E-01, - 2.00000E-01, - 3.00000E-01, - 4.00000E-01, - 5.00000E-01, - 6.00000E-01, - 8.00000E-01, - 1.00000E+00, - 1.25000E+00, - 1.50000E+00, - 2.00000E+00, - 3.00000E+00, - 4.00000E+00, - 5.00000E+00, - 6.00000E+00, - 8.00000E+00, - 1.00000E+01, - 1.50000E+01, - 2.00000E+01 +# Liver, Adult Healthy +# +# https://physics.nist.gov/cgi-bin/Xcom/xcom2?Method=Mix&Output2=Hand +# +# H 10.2 +# C 13.9 +# N 3.0 +# O 71.6 +# Na 0.2 +# S 0.3 +# Cl 0.2 +# P 0.3 +# K 0.3 + +liver_data = np.array([[ + 5.000E-03, + 1.000E-02, + 1.500E-02, + 2.000E-02, + 2.500E-02, + 3.000E-02, + 3.500E-02, + 4.000E-02, + 4.500E-02, + 5.000E-02, + 5.500E-02, + 6.000E-02, + 6.500E-02, + 7.000E-02, + 7.500E-02, + 8.000E-02, + 8.500E-02, + 9.000E-02, + 9.500E-02, + 1.000E-01, + 1.050E-01, + 1.100E-01, + 1.150E-01, + 1.200E-01, + 1.250E-01, + 1.300E-01, + 1.350E-01, + 1.400E-01, + 1.450E-01, + 1.500E-01, + 1.550E-01, + 1.600E-01, + 1.650E-01, + 1.700E-01, + 1.750E-01, + 1.800E-01, + 1.850E-01, + 1.900E-01, + 1.950E-01, + 2.000E-01 + +], [ + 4.242E+01, + 5.402E+00, + 1.706E+00, + 8.258E-01, + 5.163E-01, + 3.798E-01, + 3.096E-01, + 2.691E-01, + 2.437E-01, + 2.265E-01, + 2.142E-01, + 2.049E-01, + 1.976E-01, + 1.917E-01, + 1.867E-01, + 1.824E-01, + 1.786E-01, + 1.752E-01, + 1.721E-01, + 1.693E-01, + 1.667E-01, + 1.643E-01, + 1.621E-01, + 1.600E-01, + 1.579E-01, + 1.560E-01, + 1.542E-01, + 1.525E-01, + 1.508E-01, + 1.492E-01, + 1.476E-01, + 1.461E-01, + 1.447E-01, + 1.433E-01, + 1.419E-01, + 1.406E-01, + 1.394E-01, + 1.381E-01, + 1.369E-01, + 1.358E-01 +]]) +liver_data[0, :] = liver_data[0, :] * 1000 + +# Brain +# +# https://physics.nist.gov/cgi-bin/Xcom/xcom2?Method=Mix&Output2=Hand +# +# H 10.7 +# C 14.5 +# N 2.2 +# O 71.2 +# Na 0.2 +# S 0.2 +# Cl 0.3 +# P 0.4 +# K 0.3 + +brain_data = np.array([[ + 5.000E-03, + 1.000E-02, + 1.500E-02, + 2.000E-02, + 2.500E-02, + 3.000E-02, + 3.500E-02, + 4.000E-02, + 4.500E-02, + 5.000E-02, + 5.500E-02, + 6.000E-02, + 6.500E-02, + 7.000E-02, + 7.500E-02, + 8.000E-02, + 8.500E-02, + 9.000E-02, + 9.500E-02, + 1.000E-01, + 1.050E-01, + 1.100E-01, + 1.150E-01, + 1.200E-01, + 1.250E-01, + 1.300E-01, + 1.350E-01, + 1.400E-01, + 1.450E-01, + 1.500E-01, + 1.550E-01, + 1.600E-01, + 1.650E-01, + 1.700E-01, + 1.750E-01, + 1.800E-01, + 1.850E-01, + 1.900E-01, + 1.950E-01, + 2.000E-01 + ], [ - 2.628E+03, - 2.392E+03, - 2.176E+03, - 2.182E+03, - 8.622E+02, - 3.800E+02, - 2.053E+02, - 2.072E+02, + 4.242E+01, + 5.411E+00, + 1.710E+00, + 8.283E-01, + 5.181E-01, + 3.811E-01, + 3.107E-01, + 2.702E-01, + 2.447E-01, + 2.275E-01, + 2.151E-01, + 2.058E-01, + 1.985E-01, + 1.925E-01, + 1.874E-01, + 1.831E-01, + 1.793E-01, + 1.759E-01, + 1.729E-01, + 1.701E-01, + 1.675E-01, + 1.650E-01, + 1.628E-01, + 1.606E-01, + 1.586E-01, + 1.567E-01, + 1.549E-01, + 1.531E-01, + 1.514E-01, + 1.498E-01, + 1.483E-01, + 1.468E-01, + 1.453E-01, + 1.439E-01, + 1.426E-01, + 1.412E-01, + 1.400E-01, + 1.387E-01, + 1.375E-01, + 1.364E-01 +]]) +brain_data[0, :] = brain_data[0, :] * 1000 + +# Iron +# +# https://physics.nist.gov/cgi-bin/Xcom/xcom2?Method=Mix&Output2=Hand +# +# Fe 100.0 + +iron_data = np.array([[ + 5.000E-03, + 1.000E-02, + 1.500E-02, + 2.000E-02, + 2.500E-02, + 3.000E-02, + 3.500E-02, + 4.000E-02, + 4.500E-02, + 5.000E-02, + 5.500E-02, + 6.000E-02, + 6.500E-02, + 7.000E-02, + 7.500E-02, + 8.000E-02, + 8.500E-02, + 9.000E-02, + 9.500E-02, + 1.000E-01, + 1.050E-01, + 1.100E-01, + 1.150E-01, + 1.200E-01, + 1.250E-01, + 1.300E-01, + 1.350E-01, + 1.400E-01, + 1.450E-01, + 1.500E-01, + 1.550E-01, + 1.600E-01, + 1.650E-01, + 1.700E-01, + 1.750E-01, + 1.800E-01, + 1.850E-01, + 1.900E-01, + 1.950E-01, + 2.000E-01 + +], [ + 1.399E+02, 1.707E+02, - 1.405E+02, - 1.420E+02, - 1.188E+02, - 5.054E+01, + 5.709E+01, + 2.568E+01, + 1.370E+01, + 8.176E+00, + 5.285E+00, + 3.629E+00, + 2.614E+00, + 1.957E+00, + 1.515E+00, + 1.205E+00, + 9.816E-01, + 8.164E-01, + 6.915E-01, + 5.952E-01, + 5.198E-01, + 4.597E-01, + 4.113E-01, + 3.717E-01, + 3.391E-01, + 3.118E-01, + 2.888E-01, + 2.693E-01, + 2.526E-01, + 2.381E-01, + 2.256E-01, + 2.146E-01, + 2.050E-01, + 1.964E-01, + 1.888E-01, + 1.820E-01, + 1.759E-01, + 1.704E-01, + 1.654E-01, + 1.608E-01, + 1.567E-01, + 1.528E-01, + 1.493E-01, + 1.460E-01 +]]) +iron_data[0, :] = iron_data[0, :] * 1000 + +# Adipose #2 +# +# https://physics.nist.gov/cgi-bin/Xcom/xcom2?Method=Mix&Output2=Hand +# +# H 11.4 +# C 59.8 +# N 0.7 +# O 27.8 +# Na 0.1 +# S 0.1 +# Cl 0.1 + +adipose_data = np.array([[ + 5.000E-03, + 1.000E-02, + 1.500E-02, + 2.000E-02, + 2.500E-02, + 3.000E-02, + 3.500E-02, + 4.000E-02, + 4.500E-02, + 5.000E-02, + 5.500E-02, + 6.000E-02, + 6.500E-02, + 7.000E-02, + 7.500E-02, + 8.000E-02, + 8.500E-02, + 9.000E-02, + 9.500E-02, + 1.000E-01, + 1.050E-01, + 1.100E-01, + 1.150E-01, + 1.200E-01, + 1.250E-01, + 1.300E-01, + 1.350E-01, + 1.400E-01, + 1.450E-01, + 1.500E-01, + 1.550E-01, + 1.600E-01, + 1.650E-01, + 1.700E-01, + 1.750E-01, + 1.800E-01, + 1.850E-01, + 1.900E-01, + 1.950E-01, + 2.000E-01 + +], [ 2.587E+01, - 1.494E+01, - 6.300E+00, 3.268E+00, 1.083E+00, - 5.677E-01, + 5.678E-01, + 3.868E-01, 3.063E-01, + 2.643E-01, 2.396E-01, + 2.235E-01, 2.123E-01, + 2.039E-01, 1.974E-01, + 1.920E-01, + 1.875E-01, + 1.835E-01, 1.800E-01, + 1.769E-01, + 1.740E-01, + 1.713E-01, 1.688E-01, + 1.665E-01, + 1.643E-01, + 1.622E-01, + 1.602E-01, + 1.584E-01, + 1.565E-01, + 1.548E-01, + 1.531E-01, + 1.515E-01, 1.500E-01, - 1.368E-01, - 1.187E-01, - 1.062E-01, - 9.696E-02, - 8.965E-02, - 7.873E-02, - 7.078E-02, - 6.330E-02, - 5.760E-02, - 4.940E-02, - 3.955E-02, - 3.377E-02, - 2.995E-02, - 2.725E-02, - 2.368E-02, - 2.145E-02, - 1.843E-02, - 1.698E-02 + 1.485E-01, + 1.470E-01, + 1.456E-01, + 1.442E-01, + 1.429E-01, + 1.416E-01, + 1.404E-01, + 1.392E-01, + 1.380E-01, + 1.368E-01 ]]) adipose_data[0, :] = adipose_data[0, :] * 1000 @@ -349,209 +667,201 @@ ]]) blood_data[0, :] = blood_data[0, :] * 1000 +# Cortical bone +# +# https://physics.nist.gov/cgi-bin/Xcom/xcom2?Method=Mix&Output2=Hand +# +# H 3.4 +# C 15.5 +# N 4.2 +# O 43.5 +# Ca 22.5 +# Na 0.1 +# S 0.3 +# Mg 0.2 +# P 10.3 + bone_data = np.array([[ - 1.00000E-03, - 1.03542E-03, - 1.07210E-03, - 1.07210E-03, - 1.18283E-03, - 1.30500E-03, - 1.30500E-03, - 1.50000E-03, - 2.00000E-03, - 2.14550E-03, - 2.14550E-03, - 2.30297E-03, - 2.47200E-03, - 2.47200E-03, - 3.00000E-03, - 4.00000E-03, - 4.03810E-03, - 4.03810E-03, - 5.00000E-03, - 6.00000E-03, - 8.00000E-03, - 1.00000E-02, - 1.50000E-02, - 2.00000E-02, - 3.00000E-02, - 4.00000E-02, - 5.00000E-02, - 6.00000E-02, - 8.00000E-02, - 1.00000E-01, - 1.50000E-01, - 2.00000E-01, - 3.00000E-01, - 4.00000E-01, - 5.00000E-01, - 6.00000E-01, - 8.00000E-01, - 1.00000E+00, - 1.25000E+00, - 1.50000E+00, - 2.00000E+00, - 3.00000E+00, - 4.00000E+00, - 5.00000E+00, - 6.00000E+00, - 8.00000E+00, - 1.00000E+01, - 1.50000E+01, - 2.00000E+01 + 5.000E-03, + 1.000E-02, + 1.500E-02, + 2.000E-02, + 2.500E-02, + 3.000E-02, + 3.500E-02, + 4.000E-02, + 4.500E-02, + 5.000E-02, + 5.500E-02, + 6.000E-02, + 6.500E-02, + 7.000E-02, + 7.500E-02, + 8.000E-02, + 8.500E-02, + 9.000E-02, + 9.500E-02, + 1.000E-01, + 1.050E-01, + 1.100E-01, + 1.150E-01, + 1.200E-01, + 1.250E-01, + 1.300E-01, + 1.350E-01, + 1.400E-01, + 1.450E-01, + 1.500E-01, + 1.550E-01, + 1.600E-01, + 1.650E-01, + 1.700E-01, + 1.750E-01, + 1.800E-01, + 1.850E-01, + 1.900E-01, + 1.950E-01, + 2.000E-01 ], [ - 3.781E+03, - 3.452E+03, - 3.150E+03, - 3.156E+03, - 2.434E+03, - 1.873E+03, - 1.883E+03, - 1.295E+03, - 5.869E+02, - 4.824E+02, - 7.114E+02, - 5.916E+02, - 4.907E+02, - 4.962E+02, - 2.958E+02, - 1.331E+02, - 1.296E+02, - 3.332E+02, - 1.917E+02, - 1.171E+02, - 5.323E+01, - 2.851E+01, - 9.032E+00, - 4.001E+00, - 1.331E+00, - 6.655E-01, - 4.242E-01, - 3.148E-01, - 2.229E-01, - 1.855E-01, - 1.480E-01, - 1.309E-01, - 1.113E-01, - 9.908E-02, - 9.022E-02, - 8.332E-02, - 7.308E-02, - 6.566E-02, - 5.871E-02, - 5.346E-02, - 4.607E-02, - 3.745E-02, - 3.257E-02, - 2.946E-02, - 2.734E-02, - 2.467E-02, - 2.314E-02, - 2.132E-02, - 2.068E-02 + 1.917E+02, + 2.852E+01, + 9.032E+00, + 4.001E+00, + 2.158E+00, + 1.331E+00, + 9.061E-01, + 6.656E-01, + 5.189E-01, + 4.242E-01, + 3.601E-01, + 3.148E-01, + 2.818E-01, + 2.570E-01, + 2.379E-01, + 2.229E-01, + 2.108E-01, + 2.008E-01, + 1.926E-01, + 1.855E-01, + 1.795E-01, + 1.743E-01, + 1.697E-01, + 1.656E-01, + 1.620E-01, + 1.587E-01, + 1.557E-01, + 1.529E-01, + 1.504E-01, + 1.480E-01, + 1.459E-01, + 1.438E-01, + 1.419E-01, + 1.401E-01, + 1.383E-01, + 1.367E-01, + 1.351E-01, + 1.337E-01, + 1.322E-01, + 1.309E-01 ]]) bone_data[0, :] = bone_data[0, :] * 1000 +# Muscle +# +# https://physics.nist.gov/cgi-bin/Xcom/xcom2?Method=Mix&Output2=Hand +# +# H 10.2 +# C 14.3 +# N 3.4 +# O 71 +# Na 0.1 +# S 0.3 +# Cl 0.1 +# P 0.2 +# K 0.4 + muscle_data = np.array([[ - 1.00000E-03, - 1.03542E-03, - 1.07210E-03, - 1.07210E-03, - 1.50000E-03, - 2.00000E-03, - 2.14550E-03, - 2.14550E-03, - 2.30297E-03, - 2.47200E-03, - 2.47200E-03, - 2.64140E-03, - 2.82240E-03, - 2.82240E-03, - 3.00000E-03, - 3.60740E-03, - 3.60740E-03, - 4.00000E-03, - 5.00000E-03, - 6.00000E-03, - 8.00000E-03, - 1.00000E-02, - 1.50000E-02, - 2.00000E-02, - 3.00000E-02, - 4.00000E-02, - 5.00000E-02, - 6.00000E-02, - 8.00000E-02, - 1.00000E-01, - 1.50000E-01, - 2.00000E-01, - 3.00000E-01, - 4.00000E-01, - 5.00000E-01, - 6.00000E-01, - 8.00000E-01, - 1.00000E+00, - 1.25000E+00, - 1.50000E+00, - 2.00000E+00, - 3.00000E+00, - 4.00000E+00, - 5.00000E+00, - 6.00000E+00, - 8.00000E+00, - 1.00000E+01, - 1.50000E+01, - 2.00000E+01 + 5.000E-03, + 1.000E-02, + 1.500E-02, + 2.000E-02, + 2.500E-02, + 3.000E-02, + 3.500E-02, + 4.000E-02, + 4.500E-02, + 5.000E-02, + 5.500E-02, + 6.000E-02, + 6.500E-02, + 7.000E-02, + 7.500E-02, + 8.000E-02, + 8.500E-02, + 9.000E-02, + 9.500E-02, + 1.000E-01, + 1.050E-01, + 1.100E-01, + 1.150E-01, + 1.200E-01, + 1.250E-01, + 1.300E-01, + 1.350E-01, + 1.400E-01, + 1.450E-01, + 1.500E-01, + 1.550E-01, + 1.600E-01, + 1.650E-01, + 1.700E-01, + 1.750E-01, + 1.800E-01, + 1.850E-01, + 1.900E-01, + 1.950E-01, + 2.000E-01 ], [ - 3.719E+03, - 3.393E+03, - 3.094E+03, - 3.100E+03, - 1.251E+03, - 5.594E+02, - 4.581E+02, - 4.626E+02, - 3.776E+02, - 3.085E+02, - 3.140E+02, - 2.597E+02, - 2.145E+02, - 2.160E+02, - 1.812E+02, - 1.057E+02, - 1.100E+02, - 8.127E+01, 4.206E+01, - 2.446E+01, - 1.037E+01, - 5.356E+00, + 5.357E+00, 1.693E+00, - 8.205E-01, + 8.206E-01, + 5.138E-01, 3.783E-01, + 3.087E-01, 2.685E-01, - 2.262E-01, + 2.433E-01, + 2.263E-01, + 2.140E-01, 2.048E-01, + 1.975E-01, + 1.916E-01, + 1.866E-01, 1.823E-01, + 1.785E-01, + 1.751E-01, + 1.721E-01, 1.693E-01, + 1.667E-01, + 1.643E-01, + 1.621E-01, + 1.599E-01, + 1.579E-01, + 1.560E-01, + 1.542E-01, + 1.525E-01, + 1.508E-01, 1.492E-01, - 1.358E-01, - 1.176E-01, - 1.052E-01, - 9.598E-02, - 8.874E-02, - 7.793E-02, - 7.007E-02, - 6.265E-02, - 5.701E-02, - 4.896E-02, - 3.931E-02, - 3.369E-02, - 3.000E-02, - 2.741E-02, - 2.401E-02, - 2.192E-02, - 1.915E-02, - 1.786E-02 + 1.476E-01, + 1.461E-01, + 1.447E-01, + 1.433E-01, + 1.419E-01, + 1.406E-01, + 1.394E-01, + 1.381E-01, + 1.369E-01, + 1.358E-01 ]]) muscle_data[0, :] = muscle_data[0, :] * 1000 @@ -862,6 +1172,13 @@ ]]) triglyceride_data[0, :] = triglyceride_data[0, :] * 1000 +# Water +# +# https://physics.nist.gov/cgi-bin/Xcom/xcom2?Method=Mix&Output2=Hand +# +# H 11.2 +# O 88.8 + water_data = np.array([[ 1.00000E-03, 1.50000E-03, @@ -939,6 +1256,805 @@ ]]) water_data[0, :] = water_data[0, :] * 1000 +# CTWater, QRM +# +# https://physics.nist.gov/cgi-bin/Xcom/xcom2?Method=Mix&Output2=Hand +# +# H 8.05 +# C 69.92 +# N 2.0 +# O 18.01 +# Ca 1.88 +# Cl 0.14 + +ctwater_data = np.array([[ + 5.000E-03, + 1.000E-02, + 1.500E-02, + 2.000E-02, + 2.500E-02, + 3.000E-02, + 3.500E-02, + 4.000E-02, + 4.500E-02, + 5.000E-02, + 5.500E-02, + 6.000E-02, + 6.500E-02, + 7.000E-02, + 7.500E-02, + 8.000E-02, + 8.500E-02, + 9.000E-02, + 9.500E-02, + 1.000E-01, + 1.050E-01, + 1.100E-01, + 1.150E-01, + 1.200E-01, + 1.250E-01, + 1.300E-01, + 1.350E-01, + 1.400E-01, + 1.450E-01, + 1.500E-01, + 1.550E-01, + 1.600E-01, + 1.650E-01, + 1.700E-01, + 1.750E-01, + 1.800E-01, + 1.850E-01, + 1.900E-01, + 1.950E-01, + 2.000E-01 +], [ + 3.453E+01, + 4.676E+00, + 1.535E+00, + 7.633E-01, + 4.862E-01, + 3.621E-01, + 2.976E-01, + 2.601E-01, + 2.364E-01, + 2.203E-01, + 2.086E-01, + 1.998E-01, + 1.928E-01, + 1.871E-01, + 1.823E-01, + 1.782E-01, + 1.746E-01, + 1.713E-01, + 1.684E-01, + 1.656E-01, + 1.631E-01, + 1.608E-01, + 1.586E-01, + 1.566E-01, + 1.546E-01, + 1.527E-01, + 1.510E-01, + 1.493E-01, + 1.476E-01, + 1.461E-01, + 1.446E-01, + 1.431E-01, + 1.417E-01, + 1.403E-01, + 1.390E-01, + 1.378E-01, + 1.365E-01, + 1.353E-01, + 1.342E-01, + 1.330E-01 +]]) +ctwater_data[0, :] = ctwater_data[0, :] * 1000 + +softtissue_data = np.array([[ + 1.00000E-03, + 1.03542E-03, + 1.07210E-03, + 1.07210E-03, + 1.50000E-03, + 2.00000E-03, + 2.14550E-03, + 2.14550E-03, + 2.30297E-03, + 2.47200E-03, + 2.47200E-03, + 2.64140E-03, + 2.82240E-03, + 2.82240E-03, + 3.00000E-03, + 3.60740E-03, + 3.60740E-03, + 4.00000E-03, + 5.00000E-03, + 6.00000E-03, + 8.00000E-03, + 1.00000E-02, + 1.50000E-02, + 2.00000E-02, + 3.00000E-02, + 4.00000E-02, + 5.00000E-02, + 6.00000E-02, + 8.00000E-02, + 1.00000E-01, + 1.50000E-01, + 2.00000E-01, + 3.00000E-01, + 4.00000E-01, + 5.00000E-01, + 6.00000E-01, + 8.00000E-01, + 1.00000E+00, + 1.25000E+00, + 1.50000E+00, + 2.00000E+00, + 3.00000E+00, + 4.00000E+00, + 5.00000E+00, + 6.00000E+00, + 8.00000E+00, + 1.00000E+01, + 1.50000E+01, + 2.00000E+01 +], [ + 3.71200E+03, + 3.38600E+03, + 3.08700E+03, + 3.09900E+03, + 1.25100E+03, + 5.59600E+02, + 4.58300E+02, + 4.65000E+02, + 3.80000E+02, + 3.10200E+02, + 3.15800E+02, + 2.61200E+02, + 2.15800E+02, + 2.18800E+02, + 1.83600E+02, + 1.07300E+02, + 1.10500E+02, + 8.16100E+01, + 4.22400E+01, + 2.45600E+01, + 1.04200E+01, + 5.37900E+00, + 1.69900E+00, + 8.23000E-01, + 3.79000E-01, + 2.68800E-01, + 2.26400E-01, + 2.04800E-01, + 1.82300E-01, + 1.69300E-01, + 1.49200E-01, + 1.35800E-01, + 1.17500E-01, + 1.05200E-01, + 9.59800E-02, + 8.87300E-02, + 7.79300E-02, + 7.00600E-02, + 6.26500E-02, + 5.70100E-02, + 4.89500E-02, + 3.93100E-02, + 3.36900E-02, + 3.00000E-02, + 2.74100E-02, + 2.40100E-02, + 2.19300E-02, + 1.91500E-02, + 1.78700E-02 +]]) +softtissue_data[0, :] = softtissue_data[0, :] * 1000 + +redmarrow_data = np.array([[ + 1.00000E-02, + 1.50000E-02, + 2.00000E-02, + 3.00000E-02, + 4.00000E-02, + 5.00000E-02, + 6.00000E-02, + 8.00000E-02, + 1.00000E-01, + 1.50000E-01, + 2.00000E-01, + 3.00000E-01, + 4.00000E-01, + 5.00000E-01, + 6.00000E-01, + 8.00000E-01, + 1.00000E+00, + 1.50000E+00, + 2.00000E+00, + 3.00000E+00, + 4.00000E+00, + 5.00000E+00, + 6.00000E+00, + 8.00000E+00, + 1.00000E+01, + 1.50000E+01, + 2.00000E+01, + 3.00000E+01, + 4.00000E+01, + 5.00000E+01, + 6.00000E+01, + 8.00000E+01, + 1.00000E+02 +], [ + 4.35000E+00, + 1.41000E+00, + 7.04000E-01, + 3.46000E-01, + 2.55000E-01, + 2.20000E-01, + 2.01000E-01, + 1.81000E-01, + 1.69000E-01, + 1.49000E-01, + 1.36000E-01, + 1.18000E-01, + 1.05000E-01, + 9.62000E-02, + 8.89000E-02, + 7.81000E-02, + 7.02000E-02, + 5.72000E-02, + 4.90000E-02, + 3.93000E-02, + 3.36000E-02, + 2.99000E-02, + 2.72000E-02, + 2.38000E-02, + 2.16000E-02, + 1.87000E-02, + 1.73000E-02, + 1.62000E-02, + 1.58000E-02, + 1.57000E-02, + 1.56000E-02, + 1.58000E-02, + 1.60000E-02 +]]) +redmarrow_data[0, :] = redmarrow_data[0, :] * 1000 + +yellowmarrow_data = np.array([[ + 1.00000E-02, + 1.50000E-02, + 2.00000E-02, + 3.00000E-02, + 4.00000E-02, + 5.00000E-02, + 6.00000E-02, + 8.00000E-02, + 1.00000E-01, + 1.50000E-01, + 2.00000E-01, + 3.00000E-01, + 4.00000E-01, + 5.00000E-01, + 6.00000E-01, + 8.00000E-01, + 1.00000E+00, + 1.50000E+00, + 2.00000E+00, + 3.00000E+00, + 4.00000E+00, + 5.00000E+00, + 6.00000E+00, + 8.00000E+00, + 1.00000E+01, + 1.50000E+01, + 2.00000E+01, + 3.00000E+01, + 4.00000E+01, + 5.00000E+01, + 6.00000E+01, + 8.00000E+01, + 1.00000E+02 +], [ + 3.10000E+00, + 1.03000E+00, + 5.48000E-01, + 3.01000E-01, + 2.37000E-01, + 2.11000E-01, + 1.97000E-01, + 1.80000E-01, + 1.69000E-01, + 1.50000E-01, + 1.37000E-01, + 1.19000E-01, + 1.06000E-01, + 9.70000E-02, + 8.97000E-02, + 7.88000E-02, + 7.08000E-02, + 5.77000E-02, + 4.94000E-02, + 3.96000E-02, + 3.38000E-02, + 2.99000E-02, + 2.72000E-02, + 2.36000E-02, + 2.14000E-02, + 1.84000E-02, + 1.69000E-02, + 1.56000E-02, + 1.51000E-02, + 1.49000E-02, + 1.49000E-02, + 1.50000E-02, + 1.51000E-02 +]]) +yellowmarrow_data[0, :] = yellowmarrow_data[0, :] * 1000 + +# Spongiosa: +# Source is https://physics.nist.gov/cgi-bin/Xcom/xcom2?Method=Mix&Output2=Hand +# H 0.085 +# C 0.404 +# N 0.028 +# O 0.367 +# Ca 0.074 +# Na 0.001 +# S 0.002 +# Cl 0.002 +# Mg 0.001 +# P 0.034 +# K 0.001 +# Fe 0.001 + +spongiosa_data = np.array([[ + 5.000E-03, + 1.000E-02, + 1.500E-02, + 2.000E-02, + 2.500E-02, + 3.000E-02, + 3.500E-02, + 4.000E-02, + 4.500E-02, + 5.000E-02, + 5.500E-02, + 6.000E-02, + 6.500E-02, + 7.000E-02, + 7.500E-02, + 8.000E-02, + 8.500E-02, + 9.000E-02, + 9.500E-02, + 1.000E-01, + 1.050E-01, + 1.100E-01, + 1.150E-01, + 1.200E-01, + 1.250E-01, + 1.300E-01, + 1.350E-01, + 1.400E-01, + 1.450E-01, + 1.500E-01, + 1.550E-01, + 1.600E-01, + 1.650E-01, + 1.700E-01, + 1.750E-01, + 1.800E-01, + 1.850E-01, + 1.900E-01, + 1.950E-01, + 2.000E-01 +], [ + 8.256E+01, + 1.201E+01, + 3.833E+00, + 1.755E+00, + 9.994E-01, + 6.606E-01, + 4.862E-01, + 3.868E-01, + 3.256E-01, + 2.854E-01, + 2.578E-01, + 2.378E-01, + 2.229E-01, + 2.114E-01, + 2.022E-01, + 1.947E-01, + 1.884E-01, + 1.831E-01, + 1.785E-01, + 1.745E-01, + 1.708E-01, + 1.676E-01, + 1.647E-01, + 1.619E-01, + 1.594E-01, + 1.571E-01, + 1.549E-01, + 1.529E-01, + 1.510E-01, + 1.491E-01, + 1.474E-01, + 1.457E-01, + 1.442E-01, + 1.426E-01, + 1.412E-01, + 1.398E-01, + 1.384E-01, + 1.371E-01, + 1.358E-01, + 1.346E-01 +]]) +spongiosa_data[0, :] = spongiosa_data[0, :] * 1000 + +iodine_data = np.array([[ + 1.00E-03, + 1.04E-03, + 1.07E-03, + 1.07E-03, + 1.50E-03, + 2.00E-03, + 3.00E-03, + 4.00E-03, + 4.56E-03, + 4.56E-03, + 4.70E-03, + 4.85E-03, + 4.85E-03, + 5.00E-03, + 5.19E-03, + 5.19E-03, + 6.00E-03, + 8.00E-03, + 1.00E-02, + 1.50E-02, + 2.00E-02, + 3.00E-02, + 3.32E-02, + 3.32E-02, + 4.00E-02, + 5.00E-02, + 6.00E-02, + 8.00E-02, + 1.00E-01, + 1.50E-01, + 2.00E-01, + 3.00E-01, + 4.00E-01, + 5.00E-01, + 6.00E-01, + 8.00E-01, + 1.00E+00, + 1.02E+00, + 1.25E+00, + 1.50E+00, + 2.00E+00, + 2.04E+00, + 3.00E+00, + 4.00E+00, + 5.00E+00, + 6.00E+00, + 7.00E+00, + 8.00E+00, + 9.00E+00, + 1.00E+01, + 1.10E+01, + 1.20E+01, + 1.30E+01, + 1.40E+01, + 1.50E+01, + 1.60E+01, + 1.80E+01, + 2.00E+01, + 2.20E+01, + 2.40E+01, + 2.60E+01, + 2.80E+01, + 3.00E+01, + 4.00E+01, + 5.00E+01, + 6.00E+01, + 8.00E+01, + 1.00E+02, + 1.50E+02, + 2.00E+02, + 3.00E+02, + 4.00E+02, + 5.00E+02, + 6.00E+02, + 8.00E+02, + 1.00E+03, + 1.50E+03, + 2.00E+03, + 3.00E+03, + 4.00E+03, + 5.00E+03, + 6.00E+03, + 8.00E+03, + 1.00E+04, + 1.50E+04, + 2.00E+04, + 3.00E+04, + 4.00E+04, + 5.00E+04, + 6.00E+04, + 8.00E+04, + 1.00E+05 +], [ + 9.10E+03, + 8.46E+03, + 7.87E+03, + 8.20E+03, + 3.92E+03, + 2.00E+03, + 7.42E+02, + 3.61E+02, + 2.59E+02, + 7.55E+02, + 7.08E+02, + 6.64E+02, + 8.93E+02, + 8.43E+02, + 7.67E+02, + 8.84E+02, + 6.18E+02, + 2.92E+02, + 1.63E+02, + 5.51E+01, + 2.54E+01, + 8.56E+00, + 6.55E+00, + 3.58E+01, + 2.21E+01, + 1.23E+01, + 7.58E+00, + 3.51E+00, + 1.94E+00, + 6.98E-01, + 3.66E-01, + 1.77E-01, + 1.22E-01, + 9.70E-02, + 8.31E-02, + 6.75E-02, + 5.84E-02, + 5.76E-02, + 5.11E-02, + 4.65E-02, + 4.12E-02, + 4.09E-02, + 3.72E-02, + 3.61E-02, + 3.61E-02, + 3.66E-02, + 3.73E-02, + 3.82E-02, + 3.91E-02, + 4.00E-02, + 4.10E-02, + 4.19E-02, + 4.28E-02, + 4.37E-02, + 4.46E-02, + 4.53E-02, + 4.68E-02, + 4.82E-02, + 4.95E-02, + 5.08E-02, + 5.19E-02, + 5.30E-02, + 5.40E-02, + 5.81E-02, + 6.13E-02, + 6.38E-02, + 6.76E-02, + 7.04E-02, + 7.49E-02, + 7.77E-02, + 8.09E-02, + 8.28E-02, + 8.40E-02, + 8.49E-02, + 8.61E-02, + 8.69E-02, + 8.80E-02, + 8.86E-02, + 8.93E-02, + 8.97E-02, + 8.99E-02, + 9.00E-02, + 9.03E-02, + 9.04E-02, + 9.06E-02, + 9.07E-02, + 9.08E-02, + 9.08E-02, + 9.09E-02, + 9.09E-02, + 9.09E-02, + 9.09E-02 +]]) +iodine_data[0, :] = iodine_data[0, :] * 1000 + +calcium_data = np.array([[ + 1.00E-03, + 1.50E-03, + 2.00E-03, + 3.00E-03, + 4.00E-03, + 4.04E-03, + 4.04E-03, + 5.00E-03, + 6.00E-03, + 8.00E-03, + 1.00E-02, + 1.50E-02, + 2.00E-02, + 3.00E-02, + 4.00E-02, + 5.00E-02, + 6.00E-02, + 8.00E-02, + 1.00E-01, + 1.50E-01, + 2.00E-01, + 3.00E-01, + 4.00E-01, + 5.00E-01, + 6.00E-01, + 8.00E-01, + 1.00E+00, + 1.02E+00, + 1.25E+00, + 1.50E+00, + 2.00E+00, + 2.04E+00, + 3.00E+00, + 4.00E+00, + 5.00E+00, + 6.00E+00, + 7.00E+00, + 8.00E+00, + 9.00E+00, + 1.00E+01, + 1.10E+01, + 1.20E+01, + 1.30E+01, + 1.40E+01, + 1.50E+01, + 1.60E+01, + 1.80E+01, + 2.00E+01, + 2.20E+01, + 2.40E+01, + 2.60E+01, + 2.80E+01, + 3.00E+01, + 4.00E+01, + 5.00E+01, + 6.00E+01, + 8.00E+01, + 1.00E+02, + 1.50E+02, + 2.00E+02, + 3.00E+02, + 4.00E+02, + 5.00E+02, + 6.00E+02, + 8.00E+02, + 1.00E+03, + 1.50E+03, + 2.00E+03, + 3.00E+03, + 4.00E+03, + 5.00E+03, + 6.00E+03, + 8.00E+03, + 1.00E+04, + 1.50E+04, + 2.00E+04, + 3.00E+04, + 4.00E+04, + 5.00E+04, + 6.00E+04, + 8.00E+04, + 1.00E+05 +], [ + 4.866E+03, + 1.713E+03, + 7.998E+02, + 2.676E+02, + 1.218E+02, + 1.187E+02, + 1.023E+03, + 6.025E+02, + 3.731E+02, + 1.727E+02, + 9.341E+01, + 2.979E+01, + 1.306E+01, + 4.079E+00, + 1.830E+00, + 1.019E+00, + 6.579E-01, + 3.655E-01, + 2.571E-01, + 1.674E-01, + 1.376E-01, + 1.116E-01, + 9.783E-02, + 8.851E-02, + 8.148E-02, + 7.122E-02, + 6.388E-02, + 6.320E-02, + 5.709E-02, + 5.207E-02, + 4.524E-02, + 4.478E-02, + 3.780E-02, + 3.395E-02, + 3.170E-02, + 3.035E-02, + 2.948E-02, + 2.892E-02, + 2.859E-02, + 2.839E-02, + 2.829E-02, + 2.825E-02, + 2.825E-02, + 2.831E-02, + 2.839E-02, + 2.849E-02, + 2.875E-02, + 2.903E-02, + 2.933E-02, + 2.965E-02, + 2.996E-02, + 3.029E-02, + 3.059E-02, + 3.201E-02, + 3.322E-02, + 3.424E-02, + 3.588E-02, + 3.712E-02, + 3.926E-02, + 4.062E-02, + 4.229E-02, + 4.330E-02, + 4.398E-02, + 4.448E-02, + 4.515E-02, + 4.559E-02, + 4.625E-02, + 4.661E-02, + 4.701E-02, + 4.723E-02, + 4.736E-02, + 4.745E-02, + 4.759E-02, + 4.766E-02, + 4.778E-02, + 4.784E-02, + 4.790E-02, + 4.794E-02, + 4.795E-02, + 4.797E-02, + 4.800E-02, + 4.800E-02 +]]) +calcium_data[0, :] = calcium_data[0, :] * 1000 + # Place in one table and clean up variables names = [ 'adipose', @@ -950,6 +2066,16 @@ 'cha', 'triglyceride', 'water', + 'ctwater', + 'softtissue', + 'redmarrow', + 'yellowmarrow', + 'spongiosa', + 'iodine', + 'calcium', + 'liver', + 'brain', + 'iron' ] mass_attenuation_tables = {} diff --git a/ogo/util/Helper.py b/ogo/util/Helper.py index 67ad5de..a312f2c 100644 --- a/ogo/util/Helper.py +++ b/ogo/util/Helper.py @@ -101,7 +101,7 @@ def get_phantom(phantom_type): calibration_dict['type'] = 'k2hpo4' calibration_dict['serial'] = '#4408' calibration_dict['number_rods'] = 4 - calibration_dict['rod_labels'] = [131, 132, 133, 134] + calibration_dict['rod_labels'] = [121, 122, 123, 124] calibration_dict['rod_names'] = ['A', 'B', 'C', 'D'] calibration_dict['densities'] = [58.88, -53.40, 157.05, 157.13] calibration_dict['h2o_densities'] = [1103.57, 1056.95, 1119.52, 1120.10] @@ -112,9 +112,9 @@ def get_phantom(phantom_type): calibration_dict['type'] = 'CHA' calibration_dict['serial'] = '#BDC-71' calibration_dict['number_rods'] = 3 - calibration_dict['rod_labels'] = [141, 142, 143] + calibration_dict['rod_labels'] = [125, 126, 127] calibration_dict['rod_names'] = ['A', 'B', 'C'] - calibration_dict['densities'] = [100, 400, 800] + calibration_dict['densities'] = [100.4, 399.9, 795.7] calibration_dict['h2o_densities'] = [None] elif phantom_type in 'QRM-BDC 6-rod': @@ -123,9 +123,9 @@ def get_phantom(phantom_type): calibration_dict['type'] = 'CHA' calibration_dict['serial'] = '#BDC-106' calibration_dict['number_rods'] = 6 - calibration_dict['rod_labels'] = [151, 152, 153, 154, 155, 156] + calibration_dict['rod_labels'] = [128, 129, 130, 131, 132, 133] calibration_dict['rod_names'] = ['A', 'B', 'C', 'D', 'E', 'F'] - calibration_dict['densities'] = [0, 100, 200, 400, 600, 800] + calibration_dict['densities'] = [0, 98.0, 200.0, 395.7, 596.4, 798.0] calibration_dict['h2o_densities'] = [None] elif phantom_type in 'Image Analysis QCT-3D Plus': @@ -145,7 +145,7 @@ def get_phantom(phantom_type): calibration_dict['type'] = 'CHA' calibration_dict['serial'] = 'unknown' calibration_dict['number_rods'] = 5 - calibration_dict['rod_labels'] = [121, 122, 123, 124, 125] + calibration_dict['rod_labels'] = [116, 117, 118, 119, 120] calibration_dict['rod_names'] = ['A', 'B', 'C', 'D', 'E'] calibration_dict['densities'] = [0, 50, 100, 150, 200] calibration_dict['h2o_densities'] = [None] @@ -399,6 +399,40 @@ def infoNIFTI(reader): print('!> {:30s} = {}'.format('NIFTIHeader', reader.GetNIFTIHeader())) print(guard) +def custom_labelsDict(file_path): + """Takes a labels file that is used from ITKSnap and returns a DICT + """ + labels_dict = {} + labels_hdr = '' + + # Get the header + try: + with open(file_path, 'r') as file: + for line in file: + # Process each line + # The 'line' variable will contain the current line, including the newline character '\n' + # You might want to strip whitespace, including the newline, using .strip() + #processed_line = line.strip() + if line.startswith("#") or line.strip() == "": # header information or blank line + labels_hdr += line + else: # labels + info = line.split() + ind = int(info[0]) + rgb = [int(info[1]), int(info[2]), int(info[3])] + alpha = float(info[4]) + vis = int(info[5]) + idx = int(info[6]) + lab = (line.split('\"')[1]).strip() # tricky as we need to include spaces in the label name + labels_dict[ind]={"RGB":rgb,"A":alpha,"VIS":vis,"IDX":idx,"LABEL":lab} + #print(ind,rgb,alpha,vis,idx,lab) + + except FileNotFoundError: + print(f"Error: The file '{file_path}' was not found.") + except Exception as e: + print(f"An error occurred: {e}") + + return labels_dict,labels_hdr + def applyMask(imageData, maskData): """Applies the mask to the image. @@ -649,9 +683,17 @@ def bmd_metrics(vtk_image): """ spacing = vtk_image.GetSpacing() numpy_image = vtk2numpy(vtk_image) + nonzero_data = numpy_image[numpy_image != 0] voxel_count = np.count_nonzero(numpy_image) voxel_volume = spacing[0] * spacing[1] * spacing[2] # [mm^3] voxel_volume2 = voxel_volume / 1000 # [cm^3] + + mean_signal = np.mean(nonzero_data) + std_noise = np.std(nonzero_data) + if mean_signal == 0: + signal_to_noise = 0 + else: + signal_to_noise = mean_signal / std_noise # BMD measures BMD_total = numpy_image.sum() # [mg/cc K2HPO4] @@ -667,7 +709,10 @@ def bmd_metrics(vtk_image): 'Integral BMD [mg/cc]': BMD_AVG, 'Integral BMC [mg]': BMC, 'Bone Volume [mm^3]': VOLUME_mm, - 'Bone Volume [cm^3]': VOLUME_cm + 'Bone Volume [cm^3]': VOLUME_cm, + 'mean_signal': mean_signal, + 'std_noise': std_noise, + 'SNR': signal_to_noise } @@ -1450,7 +1495,7 @@ def maskThreshold(imageData, threshold_values): combined_mask = None for i, threshold_value in enumerate(threshold_values, start=1): - message("Threshold Value: %d" % threshold_value) + #message("Threshold Value: %d" % threshold_value) # Threshold the image for the current value thres = vtk.vtkImageThreshold() diff --git a/ogo/util/cli_tools.py b/ogo/util/cli_tools.py new file mode 100644 index 0000000..0137e0b --- /dev/null +++ b/ogo/util/cli_tools.py @@ -0,0 +1,538 @@ +# /------------------------------------------------------------------------------+ +# | 29-MAR-2026 | +# | Copyright (c) Bone Imaging Laboratory | +# | All rights reserved | +# | bonelab@ucalgary.ca | +# +------------------------------------------------------------------------------+ + +""" +CLI tools and validation utility functions for argument parsing and file operations. +""" + +# Imports +import os +import math +import json +import argparse +import numpy as np +import SimpleITK as sitk + +def check_overwrite(output_file, overwrite): + """ + Check if output file can be overwritten. + + Args: + output_file: Path to the output file + overwrite: Whether to overwrite without asking + + Returns: + True if it's safe to proceed (file doesn't exist, overwrite is True, or user confirmed) + False if user declined to overwrite + """ + if not overwrite and os.path.isfile(output_file): + result = input(f'File "{output_file}" already exists. Overwrite? [y/n]: ') + if result.lower() not in ['y', 'yes']: + print('Not overwriting. Exiting...') + return False + return True + +def check_label_lengths(from_labels, to_labels): + """ + Check that from_labels and to_labels have the same length. + + Args: + from_labels: List of source labels + to_labels: List of target labels + + Returns: + True if lengths match or both are empty + False if lengths don't match (prints error and exits) + """ + if len(from_labels) != len(to_labels): + print(f'ERROR: Number of from_labels must equal number of to_labels.') + print(f' [from_labels #{len(from_labels)} != to_labels #{len(to_labels)}]') + return False + return True + +def check_image_dimensions_match(image1, image2, name1="image1", name2="image2"): + """ + Check if two images have the same dimensions. + + Args: + image1: First SimpleITK image + image2: Second SimpleITK image + name1: Name/description of first image for error message + name2: Name/description of second image for error message + + Returns: + True if dimensions match, False otherwise + """ + size1 = image1.GetSize() + size2 = image2.GetSize() + + if size1 != size2: + print(f'ERROR: Image dimensions do not match.') + print(f' {name1}: {size1}') + print(f' {name2}: {size2}') + return False + return True + +def validate_integer_range(min_val, max_val): + """ + Factory function that returns a validator for integers within a range. + + Args: + min_val: Minimum allowed value (inclusive) + max_val: Maximum allowed value (inclusive) + + Returns: + A validator function for use with argparse type parameter + """ + def validator(value): + """Validates that the integer is within the specified range.""" + try: + ivalue = int(value) + except ValueError: + raise argparse.ArgumentTypeError(f'Invalid integer value: "{value}"') + + if ivalue < min_val or ivalue > max_val: + raise argparse.ArgumentTypeError( + f'Value must be between {min_val} and {max_val}, got {ivalue}' + ) + return ivalue + + return validator + +def validate_integer_min(min_val): + """ + Factory function that returns a validator for integers greater than or equal to a minimum value. + + Args: + min_val: Minimum allowed value (inclusive) + + Returns: + A validator function for use with argparse type parameter + """ + def validator(value): + """Validates that the integer is greater than or equal to the minimum value.""" + try: + ivalue = int(value) + except ValueError: + raise argparse.ArgumentTypeError(f'Invalid integer value: "{value}"') + + if ivalue < min_val: + raise argparse.ArgumentTypeError( + f'Value must be greater than or equal to {min_val}, got {ivalue}' + ) + return ivalue + + return validator + +def validate_float_range(min_val, max_val): + """ + Factory function that returns a validator for floats within a range. + + Args: + min_val: Minimum allowed value (inclusive) + max_val: Maximum allowed value (inclusive) + + Returns: + A validator function for use with argparse type parameter + """ + def validator(value): + """Validates that the float is within the specified range.""" + try: + fvalue = float(value) + except ValueError: + raise argparse.ArgumentTypeError(f'Invalid float value: "{value}"') + + if fvalue < min_val or fvalue > max_val: + raise argparse.ArgumentTypeError( + f'Value must be between {min_val} and {max_val}, got {fvalue}' + ) + return fvalue + + return validator + +def validate_path_exists(value): + """ + Validator function to check if a path or its parent directory exists. + + This function validates two cases: + 1. The path itself exists (file or directory) + 2. The parent directory exists (useful for output files that don't exist yet) + + Args: + value: Path string to validate + + Returns: + The path if it or its parent directory exists + + Raises: + argparse.ArgumentTypeError: If neither the path nor its parent directory exists + """ + # Check if the path exists as-is + if os.path.exists(value): + return value + + # If not, check if the parent directory exists (for output files) + parent_dir = os.path.dirname(value) + + # Handle case where dirname returns empty string (current directory) + if parent_dir == '': + parent_dir = '.' + + if os.path.exists(parent_dir): + return value + + # Neither path nor parent directory exists + raise argparse.ArgumentTypeError( + f'Path does not exist: "{value}". Parent directory "{parent_dir}" also does not exist.' + ) + +def validate_input_file(allowed_extensions): + """ + Factory function that returns a validator for input files. + Checks both file existence and valid extensions. + + Args: + allowed_extensions: List or tuple of allowed file extensions (e.g., ['.nii', '.nii.gz', '.aim']) + + Returns: + A validator function for use with argparse type parameter + """ + def validator(filepath): + """Validates input file existence and extension.""" + # Check if file exists + if not os.path.isfile(filepath): + raise argparse.ArgumentTypeError( + f'File does not exist: "{filepath}"' + ) + + # Check file extension + filepath_lower = filepath.lower() + if not any(filepath_lower.endswith(ext) for ext in allowed_extensions): + # Get the actual extension + if '.' in filepath: + actual_ext = '.' + filepath.rsplit('.', 1)[1] + if filepath_lower.endswith('.nii.gz'): + actual_ext = '.nii.gz' + else: + actual_ext = 'no extension' + + raise argparse.ArgumentTypeError( + f'Invalid file extension: "{actual_ext}". ' + f'Allowed extensions: {", ".join(allowed_extensions)}' + ) + + return filepath + + return validator + +def validate_output_file(allowed_extensions): + """ + Factory function that returns a validator for output files. + Checks valid extensions only (does not check file existence). + + Args: + allowed_extensions: List or tuple of allowed file extensions (e.g., ['.nii', '.nii.gz', '.csv', '.txt']) + + Returns: + A validator function for use with argparse type parameter + """ + def validator(filepath): + """Validates output file extension.""" + # Check file extension + filepath_lower = filepath.lower() + if not any(filepath_lower.endswith(ext) for ext in allowed_extensions): + # Get the actual extension + if '.' in filepath: + actual_ext = '.' + filepath.rsplit('.', 1)[1] + if filepath_lower.endswith('.nii.gz'): + actual_ext = '.nii.gz' + else: + actual_ext = 'no extension' + + raise argparse.ArgumentTypeError( + f'Invalid output file extension: "{actual_ext}". ' + f'Allowed extensions: {", ".join(allowed_extensions)}' + ) + + return filepath + + return validator + +def validate_float_min(min_val): + """ + Factory function that returns a validator for floats greater than or equal to a minimum value. + + Args: + min_val: Minimum allowed value (inclusive) + + Returns: + A validator function for use with argparse type parameter + """ + def validator(value): + """Validates that the float is greater than or equal to the minimum value.""" + try: + fvalue = float(value) + except ValueError: + raise argparse.ArgumentTypeError(f'Invalid float value: "{value}"') + + if fvalue < min_val: + raise argparse.ArgumentTypeError( + f'Value must be greater than or equal to {min_val}, got {fvalue}' + ) + return fvalue + + return validator + +def validate_extent(extent_list): + """ + Validator for image extent specified as [xmin, xmax, ymin, ymax, zmin, zmax]. + Ensures that min < max for each dimension. + + Args: + extent_list: List of 6 integers representing the extent + + Returns: + The validated extent list + + Raises: + argparse.ArgumentTypeError: If validation fails + """ + if len(extent_list) != 6: + raise argparse.ArgumentTypeError('set_extent requires exactly 6 integers') + + xmin, xmax, ymin, ymax, zmin, zmax = extent_list + + if xmin >= xmax: + raise argparse.ArgumentTypeError(f'xmin ({xmin}) must be less than xmax ({xmax})') + if ymin >= ymax: + raise argparse.ArgumentTypeError(f'ymin ({ymin}) must be less than ymax ({ymax})') + if zmin >= zmax: + raise argparse.ArgumentTypeError(f'zmin ({zmin}) must be less than zmax ({zmax})') + + return extent_list + +def load_config_from_json(json_file): + """ + Load configuration from JSON file. + + Args: + json_file: Path to JSON configuration file + + Returns: + Dictionary of configuration parameters + """ + with open(json_file, 'r') as f: + return json.load(f) + +def save_config_to_json(args, json_file): + """ + Save configuration to JSON file or output to stdout. + + Args: + args: Parsed arguments namespace + json_file: Path to output JSON file, or None to output to stdout + """ + config = vars(args).copy() + # Remove json-related arguments + config.pop('json', None) + + if json_file is None: + # Output to stdout + print(json.dumps(config, indent=2)) + else: + # Validate file extension + if not json_file.endswith('.json'): + raise ValueError(f'Configuration file must have .json extension: {json_file}') + + # Save to file + with open(json_file, 'w') as f: + json.dump(config, f, indent=2) + print(f'Configuration saved to: {json_file}') + +def print_image_info(infile, image, brief=False): + """ + Print formatted image information including dimensions, spacing, orientation, and statistics. + Supports 2D, 3D, and 4D images. + + Args: + infile: Path to the image file + image: SimpleITK image object + brief: If True, print only brief summary (File, dim, pos, element size, Type of data). + If False (default), print full detailed information. + """ + guard = '!-------------------------------------------------------------------------------' + + # Get basic image properties + size = list(image.GetSize()) + spacing = list(image.GetSpacing()) + origin = list(image.GetOrigin()) + + # Determine actual dimensionality + actual_ndim = len(size) + + # Calculate position + position = [math.floor(o / s) if s != 0 else 0 for o, s in zip(origin, spacing)] + + # Determine effective dimensionality (count non-singleton dimensions) + effective_ndim = sum(1 for s in size[:actual_ndim] if s > 1) + dim_str = f'{effective_ndim}D' + + # Print formatted information + print(guard) + print('!>') + print(f'!> File {os.path.basename(infile)}') + + # Format dimension-dependent lines based on actual dimensionality + if actual_ndim == 2: + print(f'!> dim {size[0]:>8} {size[1]:>8}') + elif actual_ndim == 3: + print(f'!> dim {size[0]:>8} {size[1]:>8} {size[2]:>8}') + else: # 4D + print(f'!> dim {size[0]:>8} {size[1]:>8} {size[2]:>8} {size[3]:>8}') + + if brief: + # Brief mode: only essential fields + if actual_ndim == 2: + print(f'!> pos {position[0]:>8} {position[1]:>8}') + print(f'!> element size in mm {spacing[0]:>8.4f} {spacing[1]:>8.4f}') + elif actual_ndim == 3: + print(f'!> pos {position[0]:>8} {position[1]:>8} {position[2]:>8}') + print(f'!> element size in mm {spacing[0]:>8.4f} {spacing[1]:>8.4f} {spacing[2]:>8.4f}') + else: # 4D + print(f'!> pos {position[0]:>8} {position[1]:>8} {position[2]:>8} {position[3]:>8}') + print(f'!> element size in mm {spacing[0]:>8.4f} {spacing[1]:>8.4f} {spacing[2]:>8.4f} {spacing[3]:>8.4f}') + print(f'!> Type of data {image.GetPixelIDTypeAsString()}') + print('!>') + else: + # Full mode: all information + direction = list(image.GetDirection()) + + # Calculate derived properties (only for actual dimensions) + phys_dim = [s * sp for s, sp in zip(size[:actual_ndim], spacing[:actual_ndim])] + + n_image_voxels = np.prod(size[:actual_ndim]).astype(int) + voxel_volume = np.prod(spacing[:actual_ndim]) + total_physical_volume = np.prod(phys_dim[:actual_ndim]) + + # Calculate physical extent (bounding box) + phys_min = origin[:actual_ndim] + + phys_max_coords = [float(s) for s in size[:actual_ndim]] + phys_max = list(image.TransformContinuousIndexToPhysicalPoint(phys_max_coords)) + + # Handle direction matrix for different dimensions + if actual_ndim == 2: + # 2D: 2x2 matrix + direction_matrix = np.array(direction).reshape(2, 2) + # For 2D, we don't calculate handedness + handedness = "N/A (2D)" + elif actual_ndim == 3: + # 3D: 3x3 matrix + direction_matrix = np.array(direction).reshape(3, 3) + determinant = np.linalg.det(direction_matrix) + handedness = "right-handed" if determinant > 0 else "left-handed" + elif actual_ndim == 4: + # 4D: 4x4 matrix + direction_matrix = np.array(direction).reshape(4, 4) + # For 4D, calculate handedness from 3x3 submatrix + determinant = np.linalg.det(direction_matrix[:3, :3]) + handedness = "right-handed" if determinant > 0 else "left-handed" + else: + direction_matrix = np.eye(max(actual_ndim, 2)) + handedness = "unknown" + + # Get file size with proper unit + try: + file_size = os.path.getsize(infile) + size_names = ['Bytes', 'KBytes', 'MBytes', 'GBytes'] + size_unit_index = 0 + while file_size > 1024 and size_unit_index < len(size_names) - 1: + size_unit_index += 1 + file_size = file_size / 1024.0 + except OSError: + file_size = 0 + size_unit_index = 0 + size_names = ['Bytes', 'KBytes', 'MBytes', 'GBytes'] + + # Get image statistics + # Note: StatisticsImageFilter only works with 2D and 3D images + # For 4D images, extract the first 3D volume + stats = sitk.StatisticsImageFilter() + if actual_ndim == 4: + # Extract first 3D volume from 4D image for statistics + stats_image = image[:, :, :, 0] + stats.Execute(stats_image) + stats_note = " (from first volume)" + else: + stats.Execute(image) + stats_note = "" + + # Print dimension-dependent lines based on actual dimensionality + if actual_ndim == 2: + print(f'!> off {"-":>8} {"-":>8}') + print(f'!> pos {position[0]:>8} {position[1]:>8}') + print(f'!> element size in mm {spacing[0]:>8.4f} {spacing[1]:>8.4f}') + print(f'!> phys dim in mm {phys_dim[0]:>8.4f} {phys_dim[1]:>8.4f}') + print(f'!> origin in mm {origin[0]:>8.4f} {origin[1]:>8.4f}') + print('!>') + print(f'!> Physical extent (mm):') + print(f'!> min {phys_min[0]:>8.4f} {phys_min[1]:>8.4f}') + print(f'!> max {phys_max[0]:>8.4f} {phys_max[1]:>8.4f}') + print('!>') + print(f'!> Direction matrix ({handedness}):') + print(f'!> [{direction_matrix[0,0]:>7.4f} {direction_matrix[0,1]:>7.4f}]') + print(f'!> [{direction_matrix[1,0]:>7.4f} {direction_matrix[1,1]:>7.4f}]') + elif actual_ndim == 3: + print(f'!> off {"-":>8} {"-":>8} {"-":>8}') + print(f'!> pos {position[0]:>8} {position[1]:>8} {position[2]:>8}') + print(f'!> element size in mm {spacing[0]:>8.4f} {spacing[1]:>8.4f} {spacing[2]:>8.4f}') + print(f'!> phys dim in mm {phys_dim[0]:>8.4f} {phys_dim[1]:>8.4f} {phys_dim[2]:>8.4f}') + print(f'!> origin in mm {origin[0]:>8.4f} {origin[1]:>8.4f} {origin[2]:>8.4f}') + print('!>') + print(f'!> Physical extent (mm):') + print(f'!> min {phys_min[0]:>8.4f} {phys_min[1]:>8.4f} {phys_min[2]:>8.4f}') + print(f'!> max {phys_max[0]:>8.4f} {phys_max[1]:>8.4f} {phys_max[2]:>8.4f}') + print('!>') + print(f'!> Direction matrix ({handedness}):') + print(f'!> [{direction_matrix[0,0]:>7.4f} {direction_matrix[0,1]:>7.4f} {direction_matrix[0,2]:>7.4f}]') + print(f'!> [{direction_matrix[1,0]:>7.4f} {direction_matrix[1,1]:>7.4f} {direction_matrix[1,2]:>7.4f}]') + print(f'!> [{direction_matrix[2,0]:>7.4f} {direction_matrix[2,1]:>7.4f} {direction_matrix[2,2]:>7.4f}]') + else: # 4D + print(f'!> off {"-":>8} {"-":>8} {"-":>8} {"-":>8}') + print(f'!> pos {position[0]:>8} {position[1]:>8} {position[2]:>8} {position[3]:>8}') + print(f'!> element size in mm {spacing[0]:>8.4f} {spacing[1]:>8.4f} {spacing[2]:>8.4f} {spacing[3]:>8.4f}') + print(f'!> phys dim in mm {phys_dim[0]:>8.4f} {phys_dim[1]:>8.4f} {phys_dim[2]:>8.4f} {phys_dim[3]:>8.4f}') + print(f'!> origin in mm {origin[0]:>8.4f} {origin[1]:>8.4f} {origin[2]:>8.4f} {origin[3]:>8.4f}') + print('!>') + print(f'!> Physical extent (mm):') + print(f'!> min {phys_min[0]:>8.4f} {phys_min[1]:>8.4f} {phys_min[2]:>8.4f} {phys_min[3]:>8.4f}') + print(f'!> max {phys_max[0]:>8.4f} {phys_max[1]:>8.4f} {phys_max[2]:>8.4f} {phys_max[3]:>8.4f}') + print('!>') + print(f'!> Direction matrix ({handedness}):') + print(f'!> [{direction_matrix[0,0]:>7.4f} {direction_matrix[0,1]:>7.4f} {direction_matrix[0,2]:>7.4f} {direction_matrix[0,3]:>7.4f}]') + print(f'!> [{direction_matrix[1,0]:>7.4f} {direction_matrix[1,1]:>7.4f} {direction_matrix[1,2]:>7.4f} {direction_matrix[1,3]:>7.4f}]') + print(f'!> [{direction_matrix[2,0]:>7.4f} {direction_matrix[2,1]:>7.4f} {direction_matrix[2,2]:>7.4f} {direction_matrix[2,3]:>7.4f}]') + print(f'!> [{direction_matrix[3,0]:>7.4f} {direction_matrix[3,1]:>7.4f} {direction_matrix[3,2]:>7.4f} {direction_matrix[3,3]:>7.4f}]') + print('!>') + print(f'!> Dimensionality {dim_str} (actual: {actual_ndim}D)') + print(f'!> Number of voxels {n_image_voxels:,}') + + # Pre-format dimension-dependent labels for proper alignment + voxel_vol_label = f'Voxel volume (mm^{actual_ndim})' + phys_vol_label = f'Total physical volume (mm^{actual_ndim})' + print(f'!> {voxel_vol_label:<41} {voxel_volume:.6f}') + print(f'!> {phys_vol_label:<41} {total_physical_volume:.2f}') + if actual_ndim == 3: + print(f'!> Total physical volume (cm^3) {total_physical_volume/1000.0:.4f}') + print(f'!> Type of data {image.GetPixelIDTypeAsString()}') + print(f'!> Image min/max{stats_note:<28} {stats.GetMinimum():.2f} / {stats.GetMaximum():.2f}') + print(f'!> Image mean{stats_note:<30} {stats.GetMean():.2f}') + print(f'!> Image std dev{stats_note:<27} {stats.GetSigma():.2f}') + print(f'!> Total memory size {file_size:.1f} {size_names[size_unit_index]}') + print('!>') + + print(guard) diff --git a/ogo/util/graphcuts/ITK b/ogo/util/graphcuts/ITK new file mode 160000 index 0000000..ac080c9 --- /dev/null +++ b/ogo/util/graphcuts/ITK @@ -0,0 +1 @@ +Subproject commit ac080c97b1179f370ab15c8c2e8172f4479ecebd diff --git a/ogo/util/graphcuts/bin/EndostealSegmentation b/ogo/util/graphcuts/bin/EndostealSegmentation new file mode 100755 index 0000000..cba2fa7 Binary files /dev/null and b/ogo/util/graphcuts/bin/EndostealSegmentation differ diff --git a/ogo/util/graphcuts/bin/HU b/ogo/util/graphcuts/bin/HU new file mode 100755 index 0000000..924b77d Binary files /dev/null and b/ogo/util/graphcuts/bin/HU differ diff --git a/ogo/util/graphcuts/bin/PeriostealSegmentation b/ogo/util/graphcuts/bin/PeriostealSegmentation new file mode 100755 index 0000000..afae7da Binary files /dev/null and b/ogo/util/graphcuts/bin/PeriostealSegmentation differ diff --git a/ogo/util/graphcuts/bin/Sheetness b/ogo/util/graphcuts/bin/Sheetness new file mode 100755 index 0000000..dec293f Binary files /dev/null and b/ogo/util/graphcuts/bin/Sheetness differ diff --git a/ogo/util/graphcuts/bin/Sheetness2 b/ogo/util/graphcuts/bin/Sheetness2 new file mode 100755 index 0000000..dd22c3f Binary files /dev/null and b/ogo/util/graphcuts/bin/Sheetness2 differ diff --git a/ogo/util/spectral_util.py b/ogo/util/spectral_util.py new file mode 100644 index 0000000..1b41efc --- /dev/null +++ b/ogo/util/spectral_util.py @@ -0,0 +1,240 @@ +# /------------------------------------------------------------------------------+ +# | 28-OCT-2025 | +# | Copyright (c) Bone Imaging Laboratory | +# | All rights reserved | +# | bonelab@ucalgary.ca | +# +------------------------------------------------------------------------------+ + +import os +import numpy as np +from ogo.dat.MassAttenuationTables import mass_attenuation_tables +from scipy.interpolate import CubicSpline +from scipy.interpolate import UnivariateSpline + + +def mass_density(): # mass_density, g/cm3 + rho = { + "adipose": 9.500E-01, # 9.3 to 9.7 in ICRU + "air": 1.205E-03, + "blood": 1.060E+00, + "bone": 1.920E+00, + "calcium": 1.550E+00, + "cha": 3.225E+00, + "k2hpo4": 2.440E+00, + "muscle": 1.050E+00, + "water": 1.000E+00, + "ctwater": 1.033E+00, + "softtissue": 1.000E+00, + "redmarrow": 1.030E+00, + "yellowmarrow": 0.980E+00, + "spongiosa": 1.180E+00, + "iodine": 4.930, + "liver": 1.060E+00, + "brain": 1.040E+00, + "iron": 7.86 + } + return rho + + +def ct_number_to_mu(hu=0.0, mu_water=1.0): + """Convert CT number in HU to linear attenuation coefficient.""" + return mu_water * hu / 1000.0 + mu_water + + +def mu_to_ct_number(mu=0.0, mu_water=1.0): + """Convert linear attenuation coefficient to CT number in HU.""" + return (mu - mu_water) / mu_water * 1000.0 + + +def mu_from_mass_attenuation_coefficient(mu_rho=1.0, mass_density=1.0): + """Convert mass attenuation coefficient [cm^2/g] to linear attenuation [cm^-1].""" + return mu_rho * mass_density + + +def phantom_linear_fit(x=[1, 2, 3], y=[6, 20, 32]): + """Determine linear fit of density to Hounsfield units. + + x: densities (dependent), y: CT numbers (independent). Requires at least two points. + """ + if len(x) != len(y): + raise ValueError( + "phantom_linear_fit: Number of dependent and independent variables must be equal." + ) + + if len(x) < 2: + raise ValueError("phantom_linear_fit: Number of rods must be two or more.") + + deg = 1 + m, b = np.polyfit(x, y, deg) + return [m, b] + + +def interpolate_mass_attenuation(material="water", keV=90, method="cubic", smoothing_factor=None): + """Interpolate mass attenuation [cm2/g] from NIST tables at given energy (keV). + + method: 'linear', 'cubic', or 'log-log' (uses a univariate spline in log-log space). + smoothing_factor: smoothing parameter 's' for UnivariateSpline when using 'log-log'. + If None, s=0 (interpolating spline). + + Default method is 'cubic'. + """ + method_norm = method.lower().replace("_", "").replace("-", "") + if method_norm not in {"linear", "cubic", "loglog"}: + raise ValueError("interpolate_mass_attenuation: method must be 'linear', 'cubic' or 'log-log'.") + + table_name = material + "_table" + if table_name not in mass_attenuation_tables.keys(): + raise ValueError( + f"interpolate_mass_attenuation: Material {material} is not available in current NIST tables." + ) + + table = mass_attenuation_tables[table_name] + + # truncate table at material K-edge for selected materials (we need to be conservative to avoid interpolation errors) + material_key = material.lower() + k_edges = { + "iodine": 34.5, # I K-edge ~33.17 keV + "calcium": 4.5, # Ca K-edge ~4.04 keV, but we use a conservative estimate + "cha": 4.5, # calcium hydroxyapatite ~ Ca K-edge (4.038 keV) + "bone": 4.5, # bone dominated by Ca K-edge + "air": 3.5, # air ~3.2 keV + "blood": 7.2, # blood ~7.112 keV + "muscle": 4.0, # muscle - we seem to have an issue with our NIST data not monotonically increasing at values below 4 keV + "adipose": 3.0, # adipose - we seem to have an issue with our NIST data not monotonically increasing at values below 3 keV + "softtissue": 5.0, # soft tissue ~1.0 keV + "k2hpo4": 3.7 # k2hpo4 ~3.61 keV + } + + if material_key in k_edges: + k = float(k_edges[material_key]) + xp_min = float(table["Energy [keV]"].min()) + xp_max = float(table["Energy [keV]"].max()) + if k > xp_max: + raise ValueError( + f"interpolate_mass_attenuation: K-edge {k} keV for material {material} is outside table energy range [{xp_min}, {xp_max}]." + ) + # only apply truncation if K-edge is above the current minimum + # note that a problem with this code is that there are two mass attenuations at a given energy, and so we discard both. This affects, for example, iodine at 33 keV. + if k > xp_min: + table = table[table["Energy [keV]"] >= k] + + fp = table["Mass Attenuation [cm2/g]"].to_numpy() + xp = table["Energy [keV]"].to_numpy() + + keV_arr = np.asarray(keV) + if np.any(keV_arr < xp.min()) or np.any(keV_arr > xp.max()): + raise ValueError( + f"interpolate_mass_attenuation: Input keV {keV} is below or close to k-edge and cannot be interpreted from NIST table for material {material}." + ) + + if method_norm == "linear": + result = np.interp(keV_arr, xp, fp) + return result.item() if keV_arr.ndim == 0 else result + + if method_norm == "cubic": + spl = CubicSpline(xp, fp) + result = spl(keV_arr) + return result.item() if np.ndim(keV_arr) == 0 else result + + # method_norm == "loglog" -> perform interpolation in log-log space using UnivariateSpline + if np.any(xp <= 0) or np.any(fp <= 0) or np.any(keV_arr <= 0): + raise ValueError("interpolate_mass_attenuation: log-log interpolation requires positive energies and values.") + + + log_xp = np.log(xp) + log_fp = np.log(fp) + log_keV = np.log(keV_arr) + + s = 0 if smoothing_factor is None else float(smoothing_factor) + spline = UnivariateSpline(log_xp, log_fp, s=s, k=3) + log_result = spline(log_keV) + result = np.exp(log_result) + return result.item() if keV_arr.ndim == 0 else result + + +def solve_system_equations(A=np.ones((2, 2)), B=np.ones(2), use_numpy=False): + """Solve AX = B. Supports numpy solver or Cramer's rule for 2x2 and 3x3 systems.""" + if use_numpy: + a = np.array(A) + b = np.array(B) + X = np.linalg.solve(a, b) + return X + + size = len(A) + if size < 2 or size > 3: + raise ValueError("Matrix size must be 2 or 3.") + + if size == 2: + a = A[0, 0] + b = A[0, 1] + c = A[1, 0] + d = A[1, 1] + e = B[0] + f = B[1] + D = a * d - b * c + Dx = d * e - f * b + Dy = a * f - c * e + + x = Dx / D + y = Dy / D + return np.array([x, y]) + + if size == 3: + a = A[0, 0] + b = A[0, 1] + c = A[0, 2] + d = A[1, 0] + e = A[1, 1] + f = A[1, 2] + g = A[2, 0] + h = A[2, 1] + i = A[2, 2] + j = B[0] + k = B[1] + l = B[2] + D = a * (e * i - f * h) - b * (d * i - f * g) + c * (d * h - e * g) + Dx = j * (e * i - f * h) - b * (k * i - f * l) + c * (k * h - e * l) + Dy = a * (k * i - f * l) - j * (d * i - f * g) + c * (d * l - g * k) + Dz = a * (e * l - h * k) - b * (d * l - g * k) + j * (d * h - e * g) + + x = Dx / D + y = Dy / D + z = Dz / D + return np.array([x, y, z]) + + +def parse_filename(filename): + """Parse filename into basename, dirname, name (without extension), and extension. + Handles single extensions (.jpg, .nii, .tar) and common double/extensions (.nii.gz, .tar.gz, .tar.bz2). + """ + filename = str(filename) + basename = os.path.basename(filename) + dirname = os.path.dirname(filename) + + root, ext = os.path.splitext(basename) + ext = ext.lower() + + # common compressed suffixes that indicate a double extension + compressed_suffixes = {".gz", ".bz2", ".xz", ".zip"} + if ext in compressed_suffixes: + root2, ext2 = os.path.splitext(root) + if ext2: + ext = ext2.lower() + ext + root = root2 + + return basename, dirname, root, ext + + +def print_matrix(matrix, max_cols=6, float_fmt="{:12.4e}"): + """Print matrix rows (truncate to max_cols). Handles ints and floats uniformly.""" + for row in matrix: + row_to_print = row[:max_cols] if len(row) > max_cols else row + for element in row_to_print: + if isinstance(element, (int, np.integer)): + print(f"{element:4d}", end=" ") + else: + try: + print(float_fmt.format(float(element)), end=" ") + except Exception: + print(str(element), end=" ") + print() diff --git a/setup.cfg b/setup.cfg index 11ce1ca..9e8ab3e 100644 --- a/setup.cfg +++ b/setup.cfg @@ -51,6 +51,14 @@ console_scripts = ogoSegmentCalibrationRods = ogo.cli.SegmentCalibrationRods:main ogoConvertToRas = ogo.cli.ConvertToRas:main ogoCheckLabels = ogo.cli.CheckLabels:main + ogoCalculateLinearAttenuation = ogo.cli.CalculateLinearAttenuation:main + ogoMaterialDecomposition = ogo.cli.MaterialDecomposition:main + ogoStyleGuide = ogo.cli.StyleGuide:main + ogoImageCombine = ogo.cli.ImageCombine:main + ogoImageMirror = ogo.cli.ImageMirror:main + ogoImageChangeExtent = ogo.cli.ImageChangeExtent:main + ogoImageFilterByMask = ogo.cli.ImageFilterByMask:main + ogoAssignPhantomLabels = ogo.cli.AssignPhantomLabels:main [pbr] skip_changelog = 1