diff --git a/02_CODE/cfg/hfe-nodaratis.yaml b/02_CODE/cfg/hfe-nodaratis.yaml index 5660b9c..fd627a3 100644 --- a/02_CODE/cfg/hfe-nodaratis.yaml +++ b/02_CODE/cfg/hfe-nodaratis.yaml @@ -9,7 +9,7 @@ hydra: mode: MULTIRUN sweeper: params: - simulations.grayscale_filenames: C0003093 + simulations.grayscale_filenames: C0003104 mesher: @@ -20,7 +20,7 @@ mesher: image_processing: origaim_separate: False # For Hosseini Dataset, Image parameters are read from original aim, not from processed BMD file, as there they were deleted by medtool pre-processing mask_separate: True # Standard evaluation gives two separate mask file (CORTMASK and TRABMASK), False if one maskfile with two different labels - imtype: BMD # NATIVE or BMD + imtype: NATIVE # NATIVE or BMD bvtv_scaling: 1 # 0: no scaling, 1: Scaling of BVTV 61um to BVTV 11.4um [Hosseini2017] bvtv_slope: 0.963 # bvtv intercept, added to each image voxel bvtv_intercept: 0.03814 # bvtv intercept, added to each image voxel diff --git a/02_CODE/cfg/mesh.yaml b/02_CODE/cfg/mesh.yaml index 4088d51..e331717 100644 --- a/02_CODE/cfg/mesh.yaml +++ b/02_CODE/cfg/mesh.yaml @@ -12,7 +12,7 @@ meshing_settings: outside_val: 1 # threshold value for the outside of the mask lower_thresh: 0 # lower threshold for the mask upper_thresh: 0.9 # upper threshold for the mask - s: 400 # smoothing factor of the spline + s: 100 # smoothing factor of the spline k: 3 # degree of the spline interp_points: 200 # number of points to interpolate the spline dp_simplification_outer: 3 # Ramer-Douglas-Peucker simplification factor for the periosteal contour @@ -26,9 +26,9 @@ meshing_settings: n_elms_longitudinal: 3 n_elms_transverse_trab: 15 n_elms_transverse_cort: 3 - n_elms_radial: 20 # ! Should be 10 if trab_refinement is True + n_elms_radial: 80 # ! Should be 40 if trab_refinement is True ellipsoid_fitting: True - + show_plots: False # show plots during construction show_gmsh: False # show gmsh GUI write_mesh: True # write mesh to file diff --git a/02_CODE/cfg/paths-nodaratis.yaml b/02_CODE/cfg/paths-nodaratis.yaml index 43d58db..4753f3e 100644 --- a/02_CODE/cfg/paths-nodaratis.yaml +++ b/02_CODE/cfg/paths-nodaratis.yaml @@ -10,11 +10,11 @@ paths: odb_python_script: 02_CODE/src/hfe_abq/readODB_acc.py filenames: - filename_postfix_cort_mask: _V3_CORT_MASK.mha - filename_postfix_trab_mask: _V3_TRAB_MASK.mha + filename_postfix_cort_mask: _CORT_MASK_UNCOMP.mha + filename_postfix_trab_mask: _TRAB_MASK_UNCOMP.mha filename_postfix_mask: _MASK.AIM # if mask separate == False, give this filename filename_postfix_bmd: _UNCOMP.mha - filename_postfix_seg: _SEG.mha + filename_postfix_seg: _UNCOMP_SEG.mha filename_postfix_common: _common_region_MASK.mhd # if registration is true -> file of common region filename_postfix_transform: _transformation.tfm # if registration is true -> transformation file diff --git a/02_CODE/cfg/simulations-nodaratis.yaml b/02_CODE/cfg/simulations-nodaratis.yaml index 67ecce1..6680692 100644 --- a/02_CODE/cfg/simulations-nodaratis.yaml +++ b/02_CODE/cfg/simulations-nodaratis.yaml @@ -1,6 +1,8 @@ simulations: grayscale_filenames: C0003093 + C0003104 folder_id: - C0003093: 442_R_75_F \ No newline at end of file + C0003093: 442_R_75_F + C0003104: 438_L_71_F \ No newline at end of file diff --git a/02_CODE/cfg/socket.yaml b/02_CODE/cfg/socket.yaml index 383fd36..0da6160 100644 --- a/02_CODE/cfg/socket.yaml +++ b/02_CODE/cfg/socket.yaml @@ -2,24 +2,24 @@ solver: site: remote # local or remote # * simone-kubuntu - # abaqus: /var/DassaultSystemes/SIMULIA/Commands/abq2021hf6 + abaqus: /var/DassaultSystemes/SIMULIA/Commands/abq2021hf6 # * monterosa # abaqus: /var/DassaultSystemes/SIMULIA/Commands/abq2021hf4 # * Ubelix - abaqus: /storage/workspaces/artorg_msb/hpc_abaqus/Software/DassaultSystemes/SIMULIA/Commands/abq2024 + # abaqus: /storage/workspaces/artorg_msb/hpc_abaqus/Software/DassaultSystemes/SIMULIA/Commands/abq2024 socket_paths: # paths that are socket specific (you might need to change them) # * simone-macbook # workdir: /Users/msb/Documents/01_PHD/03_Methods/CLEAN-HFE-ACCURATE # * simone-kubuntu - # workdir: /home/simoneponcioni/Documents/01_PHD/03_Methods/HFE - # scratchdir: /home/simoneponcioni/.SCRATCH - # odb2vtk: "/home/simoneponcioni/Documents/04_TOOLS/ODB2VTK/python/odb2vtk.py" + workdir: /home/simoneponcioni/Documents/01_PHD/03_Methods/HFE + scratchdir: /home/simoneponcioni/.SCRATCH + odb2vtk: "/home/simoneponcioni/Documents/04_TOOLS/ODB2VTK/python/odb2vtk.py" # * monterosa # workdir: /home/sp20q110/HFE # scratchdir: /home/sp20q110/.SCRATCH # odb2vtk: "/home/sp20q110/TOOLS/ODB2VTK/python/odb2vtk.py" # * Ubelix - workdir: /storage/workspaces/artorg_msb/hpc_abaqus/poncioni/HFE - scratchdir: /storage/workspaces/artorg_msb/hpc_abaqus/poncioni/.SCRATCH - odb2vtk: /storage/workspaces/artorg_msb/hpc_abaqus/poncioni/TOOLS/ODB2VTK/python/odb2vtk.py + # workdir: /storage/workspaces/artorg_msb/hpc_abaqus/poncioni/HFE + # scratchdir: /storage/workspaces/artorg_msb/hpc_abaqus/poncioni/.SCRATCH + # odb2vtk: /storage/workspaces/artorg_msb/hpc_abaqus/poncioni/TOOLS/ODB2VTK/python/odb2vtk.py diff --git a/02_CODE/src/hfe_accurate/material_mapping.py b/02_CODE/src/hfe_accurate/material_mapping.py index b426008..bd85e6b 100644 --- a/02_CODE/src/hfe_accurate/material_mapping.py +++ b/02_CODE/src/hfe_accurate/material_mapping.py @@ -907,7 +907,8 @@ def material_mapping_spline( botnodes = bone["bnds_bot"] topnodes = bone["bnds_top"] RP_coords_s = bone["reference_point_coord"] - inp_filename = filenames["INPname"] + # inp_filename = filenames["INPname"] + inp_filename = filenames.inp_name save_dir = Path(inp_filename).parent abq = AbaqusWriter( diff --git a/02_CODE/src/hfe_accurate/preprocessing.py b/02_CODE/src/hfe_accurate/preprocessing.py index e81cf64..08976cb 100644 --- a/02_CODE/src/hfe_accurate/preprocessing.py +++ b/02_CODE/src/hfe_accurate/preprocessing.py @@ -451,7 +451,7 @@ def msl_triangulation(cfg, SEG_array, cortmask, trabmask, spacing, tolerance): try: dimZ = (np.shape(SEG_array) * spacing)[2] except TypeError: - spacing = spacing[0] # assuming isotropic spacing + spacing = np.array([spacing[0], spacing[1], spacing[2]]) dimZ = (np.shape(SEG_array) * spacing)[2] SEG_vtk, trabmask, cortmask, spacing, tolerance, dimZ = input_sanity_check( diff --git a/02_CODE/src/hfe_utils/imutils.py b/02_CODE/src/hfe_utils/imutils.py index 3bc22cc..90089a6 100644 --- a/02_CODE/src/hfe_utils/imutils.py +++ b/02_CODE/src/hfe_utils/imutils.py @@ -279,48 +279,45 @@ def AIMreader(fileINname, spacing): imvtk = reader.GetOutput() return imvtk, spacing, scaling, slope, intercept, header -def mha_reader(filename): - imsitk = sitk.ReadImage(filename) - print(type(imsitk)) - spacing = imsitk.GetSpacing() - scaling = None - slope = None - intercept = None - header = None - return imsitk, spacing, scaling, slope, intercept, header def mha_header_reader(filename): reader = sitk.ImageFileReader() - reader.SetFileName(filename) + reader.SetFileName(str(filename.resolve())) reader.LoadPrivateTagsOn() reader.ReadImageInformation() # Extract metadata and replace _LINEBREAK_ with actual line breaks for SimpleITK compatibility metadata = {} for k in reader.GetMetaDataKeys(): - v = reader.GetMetaData(k).replace('_LINEBREAK_', '\n') + v = reader.GetMetaData(k).replace("_LINEBREAK_", "\n") metadata[k] = v - print("({0}) = = \"{1}\"".format(k, v)) + print('({0}) = = "{1}"'.format(k, v)) # Extract Density: slope from the metadata density_slope = None density_intercept = None mu_scaling = None for _, value in metadata.items(): - if 'Density: slope' in value or 'Density: intercept' in value or 'Mu_Scaling' in value: - lines = value.split('\n') + if ( + "Density: slope" in value + or "Density: intercept" in value + or "Mu_Scaling" in value + ): + lines = value.split("\n") for line in lines: - if 'Density: slope' in line: + if "Density: slope" in line: density_slope = float(line.split()[-1]) - elif 'Density: intercept' in line: + elif "Density: intercept" in line: density_intercept = float(line.split()[-1]) - elif 'Mu_Scaling' in line: + elif "Mu_Scaling" in line: mu_scaling = int(line.split()[-1]) spacing = reader.GetSpacing() - + spacing = np.array(spacing) + return spacing, mu_scaling, density_slope, density_intercept + def read_img_param(filenames: FileConfig): """ Reads image parameters from the AIM image header. @@ -340,20 +337,20 @@ def read_img_param(filenames: FileConfig): print("\n ... read images") raw_filename = filenames.raw_name - print(f'raw_filename: {raw_filename}') - if raw_filename.suffix == 'AIM*': + print(f"raw_filename: {raw_filename}") + if ".AIM" in raw_filename.suffix: try: _, spacing, scaling, slope, intercept, _ = AIMreader( - raw_filename, np.array([0.0606997, 0.0606997, 0.0606997]) + str(raw_filename.resolve()), np.array([0.0606997, 0.0606997, 0.0606997]) ) except Exception as e: logger.exception(f"An error occurred while using AIMreader: {e}") raise - if raw_filename.suffix in ('.mha', '.mhd'): + if raw_filename.suffix in (".mha", ".mhd"): try: spacing, scaling, slope, intercept = mha_header_reader(raw_filename) except Exception as e: - logger.exception(f"An error occurred while using mhd_reader: {e}") + logger.exception(f"An error occurred while using mha_header_reader: {e}") raise return spacing, scaling, slope, intercept @@ -383,38 +380,48 @@ def read_image(name, filenames, bone, lock): 4. Updates the bone dictionary with the processed image array. """ + if name == "BMD" or name == "NATIVE": + image = filenames.raw_name + elif name == "SEG": + image = filenames.seg_name + elif name == "TRABMASK": + image = filenames.trab_mask_name + elif name == "CORTMASK": + image = filenames.cort_mask_name + print("\n ... read file: " + name) spacing = bone["Spacing"] - if name.suffix == '*AIM*': - IMG_vtk = AIMreader(filenames[name + "name"], spacing)[0] + if "AIM" in filenames.filename_postfix_bmd: + image = str(image.resolve()) + IMG_vtk = AIMreader(image, spacing)[0] IMG_array = vtk2numpy(IMG_vtk) IMG_sitk = sitk.GetImageFromArray(IMG_array) else: - IMG_sitk = general_image_reader(filenames[name]) - + IMG_sitk = general_image_reader(image) + IMG_sitk = sitk.PermuteAxes(IMG_sitk, [2, 1, 0]) # pad image to avoid having non-zero values at the boundary IMG_pad = pad_image(IMG_sitk, iso_pad_size=10) IMG_pad = sitk.Flip(IMG_pad, [True, False, False]) #! ONLY FOR TIBIA VALIDATION DATASET - if "C0003114" in filenames[name + "name"]: + if "C0003114" in filenames.sample: print("Removing 35 slices") IMG_pad = IMG_pad[:-35, :, :] print(IMG_pad.GetSize()) - elif "C0003111" in filenames[name + "name"]: + elif "C0003111" in filenames.sample: print("Removing 35 slices") IMG_pad = IMG_pad[:-35, :, :] print(IMG_pad.GetSize()) - elif "C0003106" in filenames[name + "name"]: + elif "C0003106" in filenames.sample: print("Removing 35 slices") IMG_pad = IMG_pad[35:, :, :] print(IMG_pad.GetSize()) - elif "C0003096" in filenames[name + "name"]: + elif "C0003096" in filenames.sample: print("Removing 10 slices") IMG_pad = IMG_pad[:-10, :, :] print(IMG_pad.GetSize()) - elif "C0003094" in filenames[name + "name"]: + elif "C0003094" in filenames.sample: print("Removing 10 slices") IMG_pad = IMG_pad[5:-10, :, :] print(IMG_pad.GetSize()) diff --git a/02_CODE/src/hfe_utils/io_utils.py b/02_CODE/src/hfe_utils/io_utils.py index f52effa..e25717c 100644 --- a/02_CODE/src/hfe_utils/io_utils.py +++ b/02_CODE/src/hfe_utils/io_utils.py @@ -14,8 +14,10 @@ def ext(filename, new_ext): - """changes the file extension""" - return filename.replace("." + filename.split(".")[-1], new_ext) + """Changes the file extension.""" + if isinstance(filename, Path): + filename = str(filename.resolve()) + return filename.rsplit(".", 1)[0] + new_ext def print_mem_usage(): @@ -415,7 +417,9 @@ def log_summary(bone, config, filenames, var): class FileConfig: - def __init__(self, cfg, sample: str, pipeline: str = "fast", origaim_separate: bool = True): + def __init__( + self, cfg, sample: str, pipeline: str = "fast", origaim_separate: bool = True + ): self.cfg = cfg self.sample = sample self.pipeline = pipeline @@ -452,7 +456,9 @@ def boundary_conditions(self) -> str: @property def file_mask(self) -> Path: - if self.pipeline == "fast" or (self.pipeline == "accurate" and not self.cfg.image_processing.mask_separate): + if self.pipeline == "fast" or ( + self.pipeline == "accurate" and not self.cfg.image_processing.mask_separate + ): return Path(f"{self.sample}{self.cfg.filenames.filename_postfix_mask}") return None @@ -524,7 +530,7 @@ def set_filenames(self): "INPname": str(self.inp_name), "VTKname": str(self.vtk_name), "SUMname": str(self.sum_name), - "VER_BPVname": str(self.ver_bpv_name) + "VER_BPVname": str(self.ver_bpv_name), } if self.pipeline == "fast": @@ -544,10 +550,10 @@ def set_filenames(self): filename["MASKname"] = str(self.mask_name) return filename - + # Usage # cfg = ... # Your configuration object # sample = "sample_name" # file_config = FileConfig(cfg, sample) # print(file_config.bmd_name) - # print(file_config.inp_name) \ No newline at end of file + # print(file_config.inp_name)