From 3f286b8b919ad87d5dc81be4b6f3772f2c810d78 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C5=BEenan=20Zuki=C4=87?= Date: Mon, 21 Feb 2022 16:34:26 -0500 Subject: [PATCH 1/7] ENH: Update required CMake and C++ standard to match ITK 5.3 --- CMakeLists.txt | 2 +- examples/CMakeLists.txt | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index ed18bd7..bc376cb 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,4 +1,4 @@ -cmake_minimum_required(VERSION 3.10.2) +cmake_minimum_required(VERSION 3.16.3) project(HASI) if(NOT ITK_SOURCE_DIR) diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 80f7356..c017a51 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,7 +1,7 @@ -cmake_minimum_required(VERSION 3.10.2) +cmake_minimum_required(VERSION 3.16.3) if(NOT CMAKE_CXX_STANDARD) - set(CMAKE_CXX_STANDARD 11) # Supported values are ``11``, ``14``, and ``17``. + set(CMAKE_CXX_STANDARD 14) # Supported values are ``11``, ``14``, and ``17``. endif() if(NOT CMAKE_CXX_STANDARD_REQUIRED) set(CMAKE_CXX_STANDARD_REQUIRED ON) From 6eccaed6d36ce9b7077442c44648d379e0c88fc0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C5=BEenan=20Zuki=C4=87?= Date: Wed, 30 Jun 2021 16:57:39 -0400 Subject: [PATCH 2/7] STYLE: remove unused code As ApplyToImageMetadata is not available in Python and might never be, remove the commented-out code which mentions it. See this commit for an attempt to make it available in Python: https://github.com/dzenanz/ITK/commit/db8157e9edcbad0d9d2a6fd4e98fab8ad1c04fb0 --- examples/ResampleAxisAligned.ipynb | 3 --- 1 file changed, 3 deletions(-) diff --git a/examples/ResampleAxisAligned.ipynb b/examples/ResampleAxisAligned.ipynb index 4bd8979..cc6261d 100644 --- a/examples/ResampleAxisAligned.ipynb +++ b/examples/ResampleAxisAligned.ipynb @@ -36,9 +36,6 @@ " transforms = itk.transformread(transform_filename)\n", " direct_transform = transforms[0]\n", " \n", - " # ApplyToImageMetadata method is not yet available in Python\n", - " # direct_transform.ApplyToImageMetadata(in_image)\n", - " # direct_transform.ApplyToImageMetadata(labels)\n", " in_image = itk.resample_in_place_image_filter(in_image, rigid_transform=direct_transform)\n", " labels = itk.resample_in_place_image_filter(labels, rigid_transform=direct_transform)\n", " \n", From e2555668a2c5b8bc88ffc6da75b63850b5b39d32 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C5=BEenan=20Zuki=C4=87?= Date: Thu, 1 Jul 2021 16:44:51 -0400 Subject: [PATCH 3/7] ENH: Add the overall segmentation pipeline in Python Handles laterality: Atlas and case under consideration can have different laterality (e.g. atlas is from right leg, case under consideration from left leg), so we appropriately mirror the atlas. Crop atlas to the bounding box of the labeled region (bone in question) Crop case to the bounding box of the bone in question --- src/hasi/mouse_femur_tibia_ct_morphometry.py | 307 +++++++++++++++++++ 1 file changed, 307 insertions(+) create mode 100644 src/hasi/mouse_femur_tibia_ct_morphometry.py diff --git a/src/hasi/mouse_femur_tibia_ct_morphometry.py b/src/hasi/mouse_femur_tibia_ct_morphometry.py new file mode 100644 index 0000000..d12c1a4 --- /dev/null +++ b/src/hasi/mouse_femur_tibia_ct_morphometry.py @@ -0,0 +1,307 @@ +#!/usr/bin/env python3 + +# Copyright NumFOCUS +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0.txt +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# Purpose: Overall segmentation pipeline + +import itk +import os +from pathlib import Path +import sys +import traceback + + +def sorted_file_list(folder, extension): + file_list = [] + for filename in os.listdir(folder): + if filename.endswith(extension): + filename = os.path.splitext(filename)[0] + filename = Path(filename).stem + file_list.append(filename) + + file_list = list(set(file_list)) # remove duplicates + file_list.sort() + return file_list + + +def read_slicer_fiducials(filename): + file = open(filename, 'r') + lines = file.readlines() + lines.pop(0) # Markups fiducial file version = 4.11 + + coordinate_system = lines[0][-4:-1] + if coordinate_system == 'RAS' or coordinate_system[-1:] == '0': + ras = True + elif coordinate_system == 'LPS' or coordinate_system[-1:] == '1': + ras = False + elif coordinate_system == 'IJK' or coordinate_system[-1:] == '2': + raise RuntimeError('Fiducials file with IJK coordinates is not supported') + else: + raise RuntimeError('Unrecognized coordinate system: ' + coordinate_system) + + lines.pop(0) # CoordinateSystem = 0 + lines.pop(0) # columns = id,x,y,z,ow,ox,oy,oz,vis,sel,lock,label,desc,associatedNodeID + + fiducials = [] + for line in lines: + e = line.split(',', 4) + p = itk.Point[itk.D, 3]() + for i in range(3): + p[i] = float(e[i + 1]) + if ras and i < 2: + p[i] = -p[i] + fiducials.append(p) + + return fiducials + + +rigid_transform_type = itk.VersorRigid3DTransform[itk.D] + +# create an atlas laterality changer transform +atlas_aa_laterality_inverter = itk.Rigid3DTransform.New() +invert_superior_inferior = atlas_aa_laterality_inverter.GetParameters() +# the canonical pose was chosen without regard for proper anatomical orientation +invert_superior_inferior[8] = -1 # so we mirror along SI axis +atlas_aa_laterality_inverter.SetParameters(invert_superior_inferior) + + +def register_landmarks(atlas_landmarks, input_landmarks): + transform_type = itk.Transform[itk.D, 3, 3] + landmark_transformer = itk.LandmarkBasedTransformInitializer[transform_type].New() + rigid_transform = rigid_transform_type.New() + landmark_transformer.SetFixedLandmarks(atlas_landmarks) + landmark_transformer.SetMovingLandmarks(input_landmarks) + landmark_transformer.SetTransform(rigid_transform) + landmark_transformer.InitializeTransform() + + # force rotation to be around center of femur head + rigid_transform.SetCenter(atlas_landmarks[0]) + # and make sure that the other corresponding point maps to it perfectly + rigid_transform.SetTranslation(input_landmarks[0] - atlas_landmarks[0]) + + return rigid_transform + + +# If label is non-zero, only the specified label participates +# in computation of the bounding box. +# Normally, all non-zero labels contribute to bounding box. +def label_bounding_box(segmentation, label=0): + if label != 0: + segmentation = itk.binary_threshold_image_filter( + segmentation, lower_threshold=label, upper_threshold=label) + image_mask_spatial_object = itk.ImageMaskSpatialObject[3].New(segmentation) + bounding_box = image_mask_spatial_object.ComputeMyBoundingBoxInIndexSpace() + return bounding_box + + +def process_case(root_dir, bone, case, bone_label, atlas): + case_base = root_dir + bone + '/' + case + '-' + atlas # prefix for case file names + + pose = read_slicer_fiducials(root_dir + bone + '/Pose.fcsv') + + case_landmarks = read_slicer_fiducials(root_dir + bone + '/' + case + '.fcsv') + pose_to_case = register_landmarks(case_landmarks, pose) + + if case[-1] != atlas[-1]: # last letter of file name is either L or R + print(f'Changing atlas laterality from {atlas[-1]} to {case[-1]}.') + # pose_to_case.Compose(atlas_aa_laterality_inverter, True) + pose_to_case.Compose(atlas_aa_laterality_inverter) + # we don't need to change laterality of atlas landmarks + # as they all lie in a plane with K coordinate of zero + + atlas_bone_label_filename = root_dir + bone + '/' + atlas + '-AA-' + bone + '-label.nrrd' + print(f'Reading {bone} variant of atlas labels from file: {atlas_bone_label_filename}') + atlas_aa_segmentation = itk.imread(atlas_bone_label_filename) + + atlas_bone_image_filename = root_dir + bone + '/' + atlas + '-AA-' + bone + '.nrrd' + print(f'Reading {bone} variant of atlas image from file: {atlas_bone_image_filename}') + atlas_aa_image = itk.imread(atlas_bone_image_filename) + + case_image_filename = root_dir + 'Data/' + case + '.nrrd' + print(f'Reading case image from file: {case_image_filename}') + case_image = itk.imread(case_image_filename) + + auto_segmentation_filename = root_dir + 'AutoSegmentations/' + case + '-label.nrrd' + print(f'Reading case bone segmentation from file: {auto_segmentation_filename}') + case_auto_segmentation = itk.imread(auto_segmentation_filename) + + print(f'Computing {bone} bounding box') + case_bounding_box = label_bounding_box(case_auto_segmentation, bone_label) + case_bone_image = itk.region_of_interest_image_filter( + case_image, + region_of_interest=case_bounding_box) + case_bone_image_filename = root_dir + 'Bones/' + case + '-' + bone + '.nrrd' + print(f'Writing case bone image to file: {case_bone_image_filename}') + itk.imwrite(case_bone_image, case_bone_image_filename) + + print('Writing atlas to case transform to file for initializing Elastix registration') + affine_pose_to_case = itk.AffineTransform[itk.D, 3].New() + affine_pose_to_case.SetCenter(pose_to_case.GetCenter()) + affine_pose_to_case.SetMatrix(pose_to_case.GetMatrix()) + affine_pose_to_case.SetOffset(pose_to_case.GetOffset()) + atlas_to_case_filename = case_base + '.tfm' + itk.transformwrite([affine_pose_to_case], atlas_to_case_filename) + out_elastix_transform = open(atlas_to_case_filename + '.txt', "w") + out_elastix_transform.writelines(['(Transform "File")\n', + '(TransformFileName "' + case + '-' + atlas + '.tfm")']) + out_elastix_transform.close() + + # Construct elastix parameter map + parameter_object = itk.ParameterObject.New() + resolutions = 4 + parameter_map_rigid = parameter_object.GetDefaultParameterMap('rigid', resolutions) + parameter_object.AddParameterMap(parameter_map_rigid) + parameter_map_bspline = parameter_object.GetDefaultParameterMap("bspline", resolutions, 1.0) + parameter_object.AddParameterMap(parameter_map_bspline) + parameter_object.SetParameter("DefaultPixelValue", "-1024") + parameter_object.SetParameter("Metric", "AdvancedMeanSquares") + # parameter_object.SetParameter("FixedInternalImagePixelType", "short") + # parameter_object.SetParameter("MovingInternalImagePixelType", "short") + # we still have to use float pixels + + print('Starting atlas registration') + registered, elastix_transform = itk.elastix_registration_method( + case_bone_image.astype(itk.F), # fixed image is used as primary input to the filter + moving_image=atlas_aa_image.astype(itk.F), + # moving_mask=atlas_aa_segmentation, + parameter_object=parameter_object, + initial_transform_parameter_file_name=atlas_to_case_filename + '.txt', + # log_to_console=True, + output_directory=root_dir + bone + '/', + log_file_name=case + '-' + atlas + '-elx.log' + ) + + # serialize each parameter map to a file. + for index in range(elastix_transform.GetNumberOfParameterMaps()): + parameter_map = elastix_transform.GetParameterMap(index) + elastix_transform.WriteParameterFile( + parameter_map, + case_base + f".{index}.txt") + + registered_filename = case_base + '-reg.nrrd' + print(f'Writing registered image to file {registered_filename}') + itk.imwrite(registered.astype(itk.SS), registered_filename) + + print('Running transformix') + elastix_transform.SetParameter('FinalBSplineInterpolationOrder', '0') + result_image_transformix = itk.transformix_filter( + atlas_aa_segmentation.astype(itk.F), # transformix only works with float pixels + elastix_transform, + # reference image? + ) + result_image = result_image_transformix.astype(itk.UC) + registered_label_file = case_base + '-label.nrrd' + print(f'Writing deformed atlas to {registered_label_file}') + itk.imwrite(result_image, registered_label_file, compression=True) + + + print('Computing morphometry features') + morphometry_filter = itk.BoneMorphometryFeaturesFilter[type(atlas_aa_image)].New(case_bone_image) + morphometry_filter.SetMaskImage(result_image) + morphometry_filter.Update() + print('BVTV', morphometry_filter.GetBVTV()) + print('TbN', morphometry_filter.GetTbN()) + print('TbTh', morphometry_filter.GetTbTh()) + print('TbSp', morphometry_filter.GetTbSp()) + print('BSBV', morphometry_filter.GetBSBV()) + + print('Generate the mesh from the segmented case') + padded_segmentation = itk.constant_pad_image_filter( + result_image, + PadUpperBound=1, + PadLowerBound=1, + Constant=0 + ) + + mesh = itk.cuberille_image_to_mesh_filter(padded_segmentation) + mesh_filename = case_base + '.vtk' + print(f'Writing the mesh to file {mesh_filename}') + itk.meshwrite(mesh, mesh_filename) + + canonical_pose_mesh = itk.transform_mesh_filter( + mesh, + transform=pose_to_case # TODO: we should use the result of Elastix registration here + ) + canonical_pose_filename = case_base + '.obj' + print(f'Writing canonical pose mesh to {canonical_pose_filename}') + itk.meshwrite(canonical_pose_mesh, canonical_pose_filename) + + +def main_processing(root_dir, bone, atlas, bone_label): + root_dir = os.path.abspath(root_dir) + '/' + data_list = sorted_file_list(root_dir + 'Data', '.nrrd') + if atlas not in data_list: + raise RuntimeError("Missing data file for the atlas") + data_list.remove(atlas) + + landmarks_list = sorted_file_list(root_dir + bone, '.fcsv') + if atlas not in landmarks_list: + raise RuntimeError("Missing landmarks file for the atlas") + landmarks_list.remove(atlas) + if 'Pose' not in landmarks_list: + raise RuntimeError("Missing Pose.fcsv file") + landmarks_list.remove('Pose') + + # check if there are any discrepancies + if data_list != landmarks_list: + print('There is a discrepancy between data_list and landmarks_list') + print('data_list:', data_list) + print('landmarks_list:', landmarks_list) + print(f'List of cases to process: {data_list}') + + atlas_image_filename = root_dir + bone + '/' + atlas + '-AA.nrrd' + print(f'Reading atlas image from file: {atlas_image_filename}') + atlas_aa_image = itk.imread(atlas_image_filename) + # atlas_aa_segmentation = itk.imread(root_dir + bone + '/' + atlas + '-AA.seg.nrrd', + # pixel_type=itk.VariableLengthVector[itk.UC]) + atlas_labels_filename = root_dir + bone + '/' + atlas + '-AA-label.nrrd' + print(f'Reading atlas labels from file: {atlas_labels_filename}') + atlas_aa_segmentation = itk.imread(atlas_labels_filename, pixel_type=itk.UC) + + print('Computing bounding box of the labels') + # reduce the image to a bounding box around the segmented bone + # the other content makes the registration more difficult + # because the knees will be bent to different degree etc + atlas_bounding_box = label_bounding_box(atlas_aa_segmentation) + + atlas_aa_segmentation = itk.region_of_interest_image_filter( + atlas_aa_segmentation, + region_of_interest=atlas_bounding_box) + atlas_bone_label_filename = root_dir + bone + '/' + atlas + '-AA-' + bone + '-label.nrrd' + print(f'Writing {bone} variant of atlas labels to file: {atlas_bone_label_filename}') + itk.imwrite(atlas_aa_segmentation, atlas_bone_label_filename) + + atlas_aa_image = itk.region_of_interest_image_filter( + atlas_aa_image, + region_of_interest=atlas_bounding_box) + atlas_bone_image_filename = root_dir + bone + '/' + atlas + '-AA-' + bone + '.nrrd' + print(f'Writing {bone} variant of atlas image to file: {atlas_bone_image_filename}') + itk.imwrite(atlas_aa_image.astype(itk.SS), atlas_bone_image_filename) + + + # now go through all the cases, doing main processing + for case in data_list: + print(u'\u2500' * 80) + print(f'Processing case {case}') + + process_case(root_dir, bone, case, bone_label, atlas) + + print(f'Done processing case {case}') + + +# main code +main_processing('../../', 'Tibia', '901-R', 2) +main_processing('../../', 'Tibia', '901-L', 2) +main_processing('../../', 'Femur', '907-L', 1) From f3bdd315a83931a26b47be4e6a9770b1001f29ab Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C5=BEenan=20Zuki=C4=87?= Date: Wed, 23 Feb 2022 16:37:58 -0500 Subject: [PATCH 4/7] BUG: Use subprocess invocation to avoid Elastix crash Without this, elastix crashes on the second iteration of the loop `for case in data_list`. --- src/hasi/mouse_femur_tibia_ct_morphometry.py | 29 ++++++++++++++------ 1 file changed, 20 insertions(+), 9 deletions(-) diff --git a/src/hasi/mouse_femur_tibia_ct_morphometry.py b/src/hasi/mouse_femur_tibia_ct_morphometry.py index d12c1a4..ea34a93 100644 --- a/src/hasi/mouse_femur_tibia_ct_morphometry.py +++ b/src/hasi/mouse_femur_tibia_ct_morphometry.py @@ -19,8 +19,8 @@ import itk import os from pathlib import Path +import subprocess import sys -import traceback def sorted_file_list(folder, extension): @@ -206,7 +206,6 @@ def process_case(root_dir, bone, case, bone_label, atlas): print(f'Writing deformed atlas to {registered_label_file}') itk.imwrite(result_image, registered_label_file, compression=True) - print('Computing morphometry features') morphometry_filter = itk.BoneMorphometryFeaturesFilter[type(atlas_aa_image)].New(case_bone_image) morphometry_filter.SetMaskImage(result_image) @@ -290,18 +289,30 @@ def main_processing(root_dir, bone, atlas, bone_label): print(f'Writing {bone} variant of atlas image to file: {atlas_bone_image_filename}') itk.imwrite(atlas_aa_image.astype(itk.SS), atlas_bone_image_filename) - # now go through all the cases, doing main processing for case in data_list: print(u'\u2500' * 80) print(f'Processing case {case}') - process_case(root_dir, bone, case, bone_label, atlas) + # process_case(root_dir, bone, case, bone_label, atlas) - print(f'Done processing case {case}') + # Elastix crashes on second iteration of this for loop, + # so we use subprocess as a work-around + status = subprocess.run(['python', __file__, root_dir, bone, case, str(bone_label), atlas]) + if status.returncode != 0: + print(f'Case {case} failed with error {status.returncode}') + else: + print(f'Success processing case {case}') -# main code -main_processing('../../', 'Tibia', '901-R', 2) -main_processing('../../', 'Tibia', '901-L', 2) -main_processing('../../', 'Femur', '907-L', 1) +if __name__ == '__main__': + if len(sys.argv) == 6: # this is the subprocess call + process_case(sys.argv[1], sys.argv[2], sys.argv[3], int(sys.argv[4]), sys.argv[5]) + elif len(sys.argv) == 1: # direct invocation + main_processing('../../', 'Tibia', '901-R', 2) + main_processing('../../', 'Tibia', '901-L', 2) + main_processing('../../', 'Femur', '907-L', 1) + else: + print(f'Invalid number of arguments: {len(sys.argv)}. Invoke the script with no arguments.') + sys.exit(len(sys.argv)) + From 3b4e7238bf96ae7e2e7d4aab401ef269e3477571 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C5=BEenan=20Zuki=C4=87?= Date: Thu, 24 Feb 2022 15:26:15 -0500 Subject: [PATCH 5/7] ENH: Further changes identified from a complete run Compress bone atlas As atlases are already axis-aligned, we don't need its landmarks --- src/hasi/mouse_femur_tibia_ct_morphometry.py | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/src/hasi/mouse_femur_tibia_ct_morphometry.py b/src/hasi/mouse_femur_tibia_ct_morphometry.py index ea34a93..2ea4fe6 100644 --- a/src/hasi/mouse_femur_tibia_ct_morphometry.py +++ b/src/hasi/mouse_femur_tibia_ct_morphometry.py @@ -242,15 +242,16 @@ def main_processing(root_dir, bone, atlas, bone_label): root_dir = os.path.abspath(root_dir) + '/' data_list = sorted_file_list(root_dir + 'Data', '.nrrd') if atlas not in data_list: - raise RuntimeError("Missing data file for the atlas") + raise RuntimeError('Missing data file for the atlas') data_list.remove(atlas) landmarks_list = sorted_file_list(root_dir + bone, '.fcsv') if atlas not in landmarks_list: - raise RuntimeError("Missing landmarks file for the atlas") - landmarks_list.remove(atlas) + print('Missing landmarks file for the atlas') + else: + landmarks_list.remove(atlas) if 'Pose' not in landmarks_list: - raise RuntimeError("Missing Pose.fcsv file") + raise RuntimeError('Missing Pose.fcsv file') landmarks_list.remove('Pose') # check if there are any discrepancies @@ -258,6 +259,7 @@ def main_processing(root_dir, bone, atlas, bone_label): print('There is a discrepancy between data_list and landmarks_list') print('data_list:', data_list) print('landmarks_list:', landmarks_list) + sys.exit(1) print(f'List of cases to process: {data_list}') atlas_image_filename = root_dir + bone + '/' + atlas + '-AA.nrrd' @@ -280,7 +282,7 @@ def main_processing(root_dir, bone, atlas, bone_label): region_of_interest=atlas_bounding_box) atlas_bone_label_filename = root_dir + bone + '/' + atlas + '-AA-' + bone + '-label.nrrd' print(f'Writing {bone} variant of atlas labels to file: {atlas_bone_label_filename}') - itk.imwrite(atlas_aa_segmentation, atlas_bone_label_filename) + itk.imwrite(atlas_aa_segmentation, atlas_bone_label_filename, compression=True) atlas_aa_image = itk.region_of_interest_image_filter( atlas_aa_image, @@ -304,14 +306,18 @@ def main_processing(root_dir, bone, atlas, bone_label): else: print(f'Success processing case {case}') + print(f'Processed {len(data_list)} cases for bone {bone} using {atlas} as atlas.\n\n\n') + if __name__ == '__main__': if len(sys.argv) == 6: # this is the subprocess call process_case(sys.argv[1], sys.argv[2], sys.argv[3], int(sys.argv[4]), sys.argv[5]) elif len(sys.argv) == 1: # direct invocation - main_processing('../../', 'Tibia', '901-R', 2) - main_processing('../../', 'Tibia', '901-L', 2) main_processing('../../', 'Femur', '907-L', 1) + main_processing('../../', 'Femur', '907-R', 1) + main_processing('../../', 'Tibia', '901-L', 2) + main_processing('../../', 'Tibia', '901-R', 2) + # TODO: add for loops here to do this for all available atlases else: print(f'Invalid number of arguments: {len(sys.argv)}. Invoke the script with no arguments.') sys.exit(len(sys.argv)) From fae2e7da986c2307d8308bc6646d1da5846a36ff Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C5=BEenan=20Zuki=C4=87?= Date: Thu, 24 Feb 2022 17:23:09 -0500 Subject: [PATCH 6/7] ENH: Record per-label bone morphometry into a CSV file --- src/hasi/mouse_femur_tibia_ct_morphometry.py | 36 +++++++++++++++----- 1 file changed, 28 insertions(+), 8 deletions(-) diff --git a/src/hasi/mouse_femur_tibia_ct_morphometry.py b/src/hasi/mouse_femur_tibia_ct_morphometry.py index 2ea4fe6..7cbb4b7 100644 --- a/src/hasi/mouse_femur_tibia_ct_morphometry.py +++ b/src/hasi/mouse_femur_tibia_ct_morphometry.py @@ -208,13 +208,30 @@ def process_case(root_dir, bone, case, bone_label, atlas): print('Computing morphometry features') morphometry_filter = itk.BoneMorphometryFeaturesFilter[type(atlas_aa_image)].New(case_bone_image) - morphometry_filter.SetMaskImage(result_image) - morphometry_filter.Update() - print('BVTV', morphometry_filter.GetBVTV()) - print('TbN', morphometry_filter.GetTbN()) - print('TbTh', morphometry_filter.GetTbTh()) - print('TbSp', morphometry_filter.GetTbSp()) - print('BSBV', morphometry_filter.GetBSBV()) + label_names = { + 1: 'Diaphysis', + 2: 'Metaphysis', + 3: 'Growth Plate', + 4: 'Epiphysis', + } + for label in label_names: + print(f'Label {label_names[label]}') + label_region = itk.binary_threshold_image_filter( + result_image, lower_threshold=label, upper_threshold=label) + morphometry_filter.SetMaskImage(label_region) + morphometry_filter.Update() + + csv_filename = root_dir + bone + '/' + atlas + '-BoneMorphometry.csv' + with open(csv_filename, 'a') as morphometry_csv: + # Bone,Atlas,Case,Label, BVTV,TbN,TbTh,TbSp,BSBV. For description of measurements see: + # https://github.com/InsightSoftwareConsortium/ITKBoneMorphometry/blob/v1.3.0/include/itkBoneMorphometryFeaturesFilter.h#L35-L36 + morphometry_csv.write(f'{bone},{atlas},{case},{label_names[label]}') + morphometry_csv.write(f',{morphometry_filter.GetBVTV()}') + morphometry_csv.write(f',{morphometry_filter.GetTbN()}') + morphometry_csv.write(f',{morphometry_filter.GetTbTh()}') + morphometry_csv.write(f',{morphometry_filter.GetTbSp()}') + morphometry_csv.write(f',{morphometry_filter.GetBSBV()}') + morphometry_csv.write('\n') print('Generate the mesh from the segmented case') padded_segmentation = itk.constant_pad_image_filter( @@ -291,6 +308,10 @@ def main_processing(root_dir, bone, atlas, bone_label): print(f'Writing {bone} variant of atlas image to file: {atlas_bone_image_filename}') itk.imwrite(atlas_aa_image.astype(itk.SS), atlas_bone_image_filename) + csv_filename = root_dir + bone + '/' + atlas + '-BoneMorphometry.csv' + with open(csv_filename, 'w') as morphometry_csv: + morphometry_csv.write('Bone,Atlas,Case,Label,BVTV,TbN,TbTh,TbSp,BSBV\n') + # now go through all the cases, doing main processing for case in data_list: print(u'\u2500' * 80) @@ -321,4 +342,3 @@ def main_processing(root_dir, bone, atlas, bone_label): else: print(f'Invalid number of arguments: {len(sys.argv)}. Invoke the script with no arguments.') sys.exit(len(sys.argv)) - From 0f50ab50d6fb532f4276ca88b244dd3c000a0155 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C5=BEenan=20Zuki=C4=87?= Date: Thu, 24 Feb 2022 17:41:11 -0500 Subject: [PATCH 7/7] ENH: Loop over all available atlases --- src/hasi/mouse_femur_tibia_ct_morphometry.py | 21 +++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/src/hasi/mouse_femur_tibia_ct_morphometry.py b/src/hasi/mouse_femur_tibia_ct_morphometry.py index 7cbb4b7..3fe67b7 100644 --- a/src/hasi/mouse_femur_tibia_ct_morphometry.py +++ b/src/hasi/mouse_femur_tibia_ct_morphometry.py @@ -214,6 +214,13 @@ def process_case(root_dir, bone, case, bone_label, atlas): 3: 'Growth Plate', 4: 'Epiphysis', } + # TODO: add proper reading of labels from .seg.nrrd files so this is not hardcoded + # label_names = { + # 1: 'Cortical', + # 2: 'Trabecular VOI', + # 3: 'New Trabecular VOI', + # 4: 'New Cortical VOI', + # } for label in label_names: print(f'Label {label_names[label]}') label_region = itk.binary_threshold_image_filter( @@ -331,14 +338,14 @@ def main_processing(root_dir, bone, atlas, bone_label): if __name__ == '__main__': - if len(sys.argv) == 6: # this is the subprocess call + if len(sys.argv) == 1: # direct invocation + atlas_list = ['901-L', '901-R', '907-L', '907-R', '917-L', '917-R', 'F9-3wk-02-L', 'F9-3wk-02-R'] + for atlas in atlas_list: + main_processing('../../', 'Femur', atlas, 1) + main_processing('../../', 'Tibia', atlas, 2) + + elif len(sys.argv) == 6: # this is the subprocess call process_case(sys.argv[1], sys.argv[2], sys.argv[3], int(sys.argv[4]), sys.argv[5]) - elif len(sys.argv) == 1: # direct invocation - main_processing('../../', 'Femur', '907-L', 1) - main_processing('../../', 'Femur', '907-R', 1) - main_processing('../../', 'Tibia', '901-L', 2) - main_processing('../../', 'Tibia', '901-R', 2) - # TODO: add for loops here to do this for all available atlases else: print(f'Invalid number of arguments: {len(sys.argv)}. Invoke the script with no arguments.') sys.exit(len(sys.argv))