From c6b3c6c92fead5b165c5fe7f1f55049e1012423c Mon Sep 17 00:00:00 2001 From: Anthony Gatt Date: Wed, 29 Oct 2025 15:43:26 -0700 Subject: [PATCH 01/16] Add function to calculate distance to other surface along a unit vector Introduced `get_distance_other_surface_at_points_along_unit_vector` to compute distances from a surface to another along a specified unit vector. The function includes parameters for ray casting length and handling cases with no intersection. Updated documentation for clarity. --- pymskt/mesh/meshTools.py | 63 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 63 insertions(+) diff --git a/pymskt/mesh/meshTools.py b/pymskt/mesh/meshTools.py index f8eb790..9f764f3 100644 --- a/pymskt/mesh/meshTools.py +++ b/pymskt/mesh/meshTools.py @@ -519,6 +519,69 @@ def get_distance_other_surface_at_points( return np.asarray(distance_data, dtype=float) +def get_distance_other_surface_at_points_along_unit_vector( + surface, + other_surface, + unit_vector, + ray_cast_length=20.0, + percent_ray_length_opposite_direction=0.25, + no_distance_filler=0.0, +): # Could be nan?? + """ + Get distance to other surface at points along a unit vector. If no intersection is found, + the filler value is returned (default is 0.0). + + Parameters + ---------- + surface : Mesh + Mesh containing vtk.vtkPolyData - get distance to other surface at points along a unit vector + other_surface : Mesh + Mesh containing vtk.vtkPolyData - the other surface to get distance to + unit_vector : np.ndarray + Unit vector along which to get distance to other surface + ray_cast_length : float, optional + Length (mm) of ray to cast from surface to other surface when trying to find distance, by default 20.0 + percent_ray_length_opposite_direction : float, optional + How far to project ray inside of the surface. This is done just in case the other surface + ends up slightly inside of (or coincident with) the surface, by default 0.25 + no_distance_filler : float, optional + Value to use if no intersection is found, by default 0.0 + + Returns + ------- + np.ndarray + n_points x 1 array of the distance to the other surface at each point on the surface + """ + + points = surface.GetPoints() + obb_other_surface = get_obb_surface(other_surface) + + assert np.isclose(np.linalg.norm(unit_vector), 1.0), "Unit vector must be a unit vector" + + distance_data = [] + # Loop through all points + for idx in range(points.GetNumberOfPoints()): + point = points.GetPoint(idx) + + end_point_ray = n2l(l2n(point) + ray_cast_length * unit_vector) + start_point_ray = n2l( + l2n(point) + ray_cast_length * percent_ray_length_opposite_direction * -unit_vector + ) + + # Check if there are any intersections for the given ray + if is_hit(obb_other_surface, start_point_ray, end_point_ray): # intersections were found + # Retrieve coordinates of intersection points and intersected cell ids + points_intersect, cell_ids_intersect = get_intersect( + obb_other_surface, start_point_ray, end_point_ray + ) + distance_data.append(np.sqrt(np.sum(np.square(l2n(point) - l2n(points_intersect[0]))))) + + + else: + distance_data.append(no_distance_filler) + + return np.asarray(distance_data, dtype=float) + def set_mesh_physical_point_coords(mesh, new_points): """ Convenience function to update the x/y/z point coords of a mesh From b387c869b8373488f6697441fd5b651b4797e672 Mon Sep 17 00:00:00 2001 From: Anthony Gatt Date: Wed, 29 Oct 2025 15:50:56 -0700 Subject: [PATCH 02/16] Refactor distance calculation functions for improved clarity and functionality Created new helper function: _get_distance_with_directions This iterates over points and a direction vector per point to project rays along the vector, detect intersection(s) and return the distances. This is then wrapped by two other functions that project vectors normal to the surface, and project vectors using a fixed direction. --- pymskt/mesh/meshTools.py | 178 +++++++++++++++++++++------------------ 1 file changed, 94 insertions(+), 84 deletions(-) diff --git a/pymskt/mesh/meshTools.py b/pymskt/mesh/meshTools.py index 9f764f3..db84690 100644 --- a/pymskt/mesh/meshTools.py +++ b/pymskt/mesh/meshTools.py @@ -439,148 +439,158 @@ def get_cartilage_properties_at_points( ) -def get_distance_other_surface_at_points( - surface, - other_surface, +def _get_distance_with_directions( + points, + obb_other_surface, + directions, ray_cast_length=20.0, percent_ray_length_opposite_direction=0.25, no_distance_filler=0.0, -): # Could be nan?? +): """ - Extract cartilage outcomes (T2 & thickness) at all points on bone surface. + Private helper function to compute distances by ray casting along specified directions. Parameters ---------- - surface_bone : BoneMesh - Bone mesh containing vtk.vtkPolyData - get outcomes for nodes (vertices) on - this mesh - surface_cartilage : CartilageMesh - Cartilage mesh containing vtk.vtkPolyData - for obtaining cartilage outcomes. - t2_vtk_image : vtk.vtkImageData, optional - vtk object that contains our Cartilage T2 data, by default None - seg_vtk_image : vtk.vtkImageData, optional - vtk object that contains the segmentation mask(s) to help assign - labels to bone surface (e.g., most common), by default None - ray_cast_length : float, optional - Length (mm) of ray to cast from bone surface when trying to find cartilage (inner & - outter shell), by default 20.0 - percent_ray_length_opposite_direction : float, optional - How far to project ray inside of the bone. This is done just in case the cartilage - surface ends up slightly inside of (or coincident with) the bone surface, by default 0.25 - no_thickness_filler : float, optional - Value to use instead of thickness (if no cartilage), by default 0. - no_t2_filler : float, optional - Value to use instead of T2 (if no cartilage), by default 0. - no_seg_filler : int, optional - Value to use if no segmentation label available (because no cartilage?), by default 0 - line_resolution : int, optional - Number of points to have along line, by default 100 + points : vtk.vtkPoints + Points from the surface mesh + obb_other_surface : vtk.vtkOBBTree + OBB tree for the other surface + directions : np.ndarray + Array of direction vectors (n_points x 3) for ray casting + ray_cast_length : float + Length of ray to cast + percent_ray_length_opposite_direction : float + How far to project ray in opposite direction + no_distance_filler : float + Value to use when no intersection is found Returns ------- - list - Will return list of data for: - Cartilage thickness - Mean T2 at each point on bone - Most common cartilage label at each point on bone (normal to surface). + np.ndarray + Array of distances for each point """ - - normals = get_surface_normals(surface) - points = surface.GetPoints() - obb_other_surface = get_obb_surface(other_surface) - point_normals = normals.GetOutput().GetPointData().GetNormals() - distance_data = [] - # Loop through all points + for idx in range(points.GetNumberOfPoints()): point = points.GetPoint(idx) - normal = point_normals.GetTuple(idx) + direction = directions[idx] - end_point_ray = n2l(l2n(point) + ray_cast_length * l2n(normal)) + end_point_ray = n2l(l2n(point) + ray_cast_length * direction) start_point_ray = n2l( - l2n(point) + ray_cast_length * percent_ray_length_opposite_direction * (-l2n(normal)) + l2n(point) - ray_cast_length * percent_ray_length_opposite_direction * direction ) # Check if there are any intersections for the given ray - if is_hit(obb_other_surface, start_point_ray, end_point_ray): # intersections were found + if is_hit(obb_other_surface, start_point_ray, end_point_ray): # Retrieve coordinates of intersection points and intersected cell ids points_intersect, cell_ids_intersect = get_intersect( obb_other_surface, start_point_ray, end_point_ray ) - # points - # if len(points_intersect) == 1: distance_data.append(np.sqrt(np.sum(np.square(l2n(point) - l2n(points_intersect[0]))))) - # else: - # distance_data.append(no_distance_filler) - else: distance_data.append(no_distance_filler) return np.asarray(distance_data, dtype=float) -def get_distance_other_surface_at_points_along_unit_vector( +def get_distance_other_surface_at_points( surface, other_surface, - unit_vector, ray_cast_length=20.0, percent_ray_length_opposite_direction=0.25, no_distance_filler=0.0, -): # Could be nan?? +): """ - Get distance to other surface at points along a unit vector. If no intersection is found, - the filler value is returned (default is 0.0). + Get distance to another surface by projecting along surface normals at each point. Parameters ---------- surface : Mesh - Mesh containing vtk.vtkPolyData - get distance to other surface at points along a unit vector + Mesh containing vtk.vtkPolyData - get distance from this surface other_surface : Mesh Mesh containing vtk.vtkPolyData - the other surface to get distance to - unit_vector : np.ndarray - Unit vector along which to get distance to other surface ray_cast_length : float, optional - Length (mm) of ray to cast from surface to other surface when trying to find distance, by default 20.0 + Length (mm) of ray to cast from surface when trying to find distance, by default 20.0 percent_ray_length_opposite_direction : float, optional - How far to project ray inside of the surface. This is done just in case the other surface - ends up slightly inside of (or coincident with) the surface, by default 0.25 + How far to project ray in opposite direction. This is done just in case the other + surface ends up slightly inside of (or coincident with) the surface, by default 0.25 no_distance_filler : float, optional - Value to use if no intersection is found, by default 0.0 + Value to use when no intersection is found, by default 0.0 Returns ------- np.ndarray - n_points x 1 array of the distance to the other surface at each point on the surface + Array of distances (n_points,) to the other surface at each point """ - + normals = get_surface_normals(surface) points = surface.GetPoints() obb_other_surface = get_obb_surface(other_surface) - - assert np.isclose(np.linalg.norm(unit_vector), 1.0), "Unit vector must be a unit vector" + point_normals = normals.GetOutput().GetPointData().GetNormals() - distance_data = [] - # Loop through all points - for idx in range(points.GetNumberOfPoints()): - point = points.GetPoint(idx) + # Extract normals as numpy array + directions = np.array([point_normals.GetTuple(idx) for idx in range(points.GetNumberOfPoints())]) - end_point_ray = n2l(l2n(point) + ray_cast_length * unit_vector) - start_point_ray = n2l( - l2n(point) + ray_cast_length * percent_ray_length_opposite_direction * -unit_vector - ) + return _get_distance_with_directions( + points, + obb_other_surface, + directions, + ray_cast_length, + percent_ray_length_opposite_direction, + no_distance_filler, + ) - # Check if there are any intersections for the given ray - if is_hit(obb_other_surface, start_point_ray, end_point_ray): # intersections were found - # Retrieve coordinates of intersection points and intersected cell ids - points_intersect, cell_ids_intersect = get_intersect( - obb_other_surface, start_point_ray, end_point_ray - ) - distance_data.append(np.sqrt(np.sum(np.square(l2n(point) - l2n(points_intersect[0]))))) +def get_distance_other_surface_at_points_along_unit_vector( + surface, + other_surface, + unit_vector, + ray_cast_length=20.0, + percent_ray_length_opposite_direction=0.25, + no_distance_filler=0.0, +): + """ + Get distance to another surface by projecting along a specified unit vector direction. - else: - distance_data.append(no_distance_filler) + Parameters + ---------- + surface : Mesh + Mesh containing vtk.vtkPolyData - get distance from this surface + other_surface : Mesh + Mesh containing vtk.vtkPolyData - the other surface to get distance to + unit_vector : np.ndarray or list + Unit vector (3D) along which to project rays for distance calculation + ray_cast_length : float, optional + Length (mm) of ray to cast from surface when trying to find distance, by default 20.0 + percent_ray_length_opposite_direction : float, optional + How far to project ray in the opposite direction. This is done just in case the other surface + ends up slightly inside of (or coincident with) the surface, by default 0.25 + no_distance_filler : float, optional + Value to use when no intersection is found, by default 0.0 - return np.asarray(distance_data, dtype=float) + Returns + ------- + np.ndarray + Array of distances (n_points,) to the other surface at each point + """ + points = surface.GetPoints() + obb_other_surface = get_obb_surface(other_surface) + + unit_vector = np.asarray(unit_vector) + assert np.isclose(np.linalg.norm(unit_vector), 1.0), "unit_vector must have magnitude 1.0" + + # Create array of identical direction vectors for all points + n_points = points.GetNumberOfPoints() + directions = np.tile(unit_vector, (n_points, 1)) + + return _get_distance_with_directions( + points, + obb_other_surface, + directions, + ray_cast_length, + percent_ray_length_opposite_direction, + no_distance_filler, + ) def set_mesh_physical_point_coords(mesh, new_points): """ From bae3dddac6ad4039d2327de4381ae3a6bc432b72 Mon Sep 17 00:00:00 2001 From: Anthony Gatt Date: Wed, 29 Oct 2025 23:38:41 -0700 Subject: [PATCH 03/16] Add MeniscusMesh class and related functionality for meniscal analysis - Introduced MeniscusMesh class for creating and analyzing meniscus meshes, including extrusion and coverage metrics.... this MeniscusMesh class isnt really used. Maybe remove? or figure out if/how should update it individually? - Updated the mesh module to import MeniscusMesh and integrate it with existing functionality. - Added new functions for computing meniscal extrusion and meniscal coverage - Enhanced BoneMesh class to support meniscal outcomes computation and caching. - Directly call the new meniscal analysis functions to get extrusion/coverage. - Added tests for meniscal coverage and extrusion calculations, both individually as well as in the BoneMesh class. --- README.md | 34 + pymskt/mesh/__init__.py | 3 +- pymskt/mesh/mesh_meniscus.py | 685 +++++++++++++ pymskt/mesh/meshes.py | 472 ++++++++- .../meshMeniscus/MeniscusMesh_create_test.py | 33 + .../compute_meniscal_coverage_test.py | 52 + .../compute_meniscal_extrusion_test.py | 491 ++++++++++ .../meniscal_extrusion_Oct.29.2025.ipynb | 910 ++++++++++++++++++ ..._extrusion_function_test_Oct.20.2025.ipynb | 113 +++ 9 files changed, 2782 insertions(+), 11 deletions(-) create mode 100644 pymskt/mesh/mesh_meniscus.py create mode 100644 testing/mesh/meshMeniscus/MeniscusMesh_create_test.py create mode 100644 testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py create mode 100644 testing/mesh/meshMeniscus/compute_meniscal_extrusion_test.py create mode 100644 testing/scratch/meniscal_extrusion_Oct.29.2025.ipynb create mode 100644 testing/scratch/meniscal_extrusion_function_test_Oct.20.2025.ipynb diff --git a/README.md b/README.md index 5d1d2d7..ef42237 100644 --- a/README.md +++ b/README.md @@ -150,6 +150,40 @@ An example of how the cartilage thickness values are computed: ![](/images/cartilage_thickness_analysis.png) +### Meniscal Analysis + +Compute meniscal extrusion and coverage metrics: + +```python +import pymskt as mskt + +# Create tibia with cartilage labels +tibia = mskt.mesh.BoneMesh( + path_seg_image='tibia.nrrd', + label_idx=6, + dict_cartilage_labels={'medial': 2, 'lateral': 3} +) +tibia.create_mesh() +tibia.calc_cartilage_thickness() +tibia.assign_cartilage_regions() + +# Create meniscus meshes +med_meniscus = mskt.mesh.Mesh(path_seg_image='tibia.nrrd', label_idx=10) +med_meniscus.create_mesh() + +lat_meniscus = mskt.mesh.Mesh(path_seg_image='tibia.nrrd', label_idx=9) +lat_meniscus.create_mesh() + +# Set menisci (labels auto-inferred from dict_cartilage_labels) +tibia.set_menisci(medial_meniscus=med_meniscus, lateral_meniscus=lat_meniscus) + +# Access metrics (auto-computes on first access) +print(f"Medial extrusion: {tibia.med_men_extrusion:.2f} mm") +print(f"Medial coverage: {tibia.med_men_coverage:.1f}%") +print(f"Lateral extrusion: {tibia.lat_men_extrusion:.2f} mm") +print(f"Lateral coverage: {tibia.lat_men_coverage:.1f}%") +``` + # Development / Contributing General information for contributing can be found [here](https://github.com/gattia/pymskt/blob/main/CONTRIBUTING.md) diff --git a/pymskt/mesh/__init__.py b/pymskt/mesh/__init__.py index 901095f..b3c12a2 100644 --- a/pymskt/mesh/__init__.py +++ b/pymskt/mesh/__init__.py @@ -1,2 +1,3 @@ -from . import createMesh, io, meshCartilage, meshRegistration, meshTools, meshTransform, utils +from . import createMesh, io, meshCartilage, meshRegistration, meshTools, meshTransform, utils, mesh_meniscus from .meshes import * +from .mesh_meniscus import MeniscusMesh diff --git a/pymskt/mesh/mesh_meniscus.py b/pymskt/mesh/mesh_meniscus.py new file mode 100644 index 0000000..5024328 --- /dev/null +++ b/pymskt/mesh/mesh_meniscus.py @@ -0,0 +1,685 @@ +""" +Meniscus mesh class and analysis functions for computing meniscal outcomes, +including extrusion and coverage. + +This module provides functionality to analyze meniscal function using healthy cartilage +reference masks. Key metrics include: +- Meniscal extrusion: how far the meniscus extends beyond the cartilage rim +- Meniscal coverage: percentage of cartilage area covered by meniscus + +All distances are in mm, areas in mm², and coverage in mm² and percentage. +""" + +import numpy as np +from pymskt.mesh.meshes import Mesh + + +class MeniscusMesh(Mesh): + """ + Class to create, store, and process meniscus meshes with specialized + analysis functions for meniscal extrusion and coverage calculations. + + Parameters + ---------- + mesh : vtk.vtkPolyData, optional + vtkPolyData object that is basis of surface mesh, by default None + seg_image : SimpleITK.Image, optional + Segmentation image that can be used to create surface mesh, by default None + path_seg_image : str, optional + Path to a medical image (.nrrd) to load and create mesh from, by default None + label_idx : int, optional + Label of anatomy of interest, by default None + min_n_pixels : int, optional + All islands smaller than this size are dropped, by default 1000 + meniscus_type : str, optional + Type of meniscus ('medial' or 'lateral'), by default None + + Attributes + ---------- + meniscus_type : str + Type of meniscus ('medial' or 'lateral') + + Examples + -------- + >>> med_meniscus = MeniscusMesh( + ... path_seg_image='meniscus_seg.nrrd', + ... label_idx=1, + ... meniscus_type='medial' + ... ) + """ + + def __init__( + self, + mesh=None, + seg_image=None, + path_seg_image=None, + label_idx=None, + min_n_pixels=1000, + meniscus_type=None, + ): + super().__init__( + mesh=mesh, + seg_image=seg_image, + path_seg_image=path_seg_image, + label_idx=label_idx, + min_n_pixels=min_n_pixels, + ) + self._meniscus_type = meniscus_type + + @property + def meniscus_type(self): + """Get the meniscus type.""" + return self._meniscus_type + + @meniscus_type.setter + def meniscus_type(self, new_meniscus_type): + """Set the meniscus type with validation.""" + if new_meniscus_type not in [None, 'medial', 'lateral']: + raise ValueError("meniscus_type must be None, 'medial', or 'lateral'") + self._meniscus_type = new_meniscus_type + + +# ============================================================================ +# Helper Functions +# ============================================================================ + +def compute_tibia_axes( + tibia_mesh, + medial_cart_label, + lateral_cart_label, + scalar_array_name='labels', +): + """ + Compute anatomical axes (ML, IS, AP) from tibial cartilage regions. + + Uses PCA on combined cartilage points to find the tibial plateau normal (IS axis). + The superior direction is determined by checking which side the bone is on + relative to the cartilage. ML axis is from medial to lateral cartilage centers. + AP axis is the cross product of ML and IS. + + Parameters + ---------- + tibia_mesh : BoneMesh or Mesh + Tibia mesh with scalar values indicating cartilage regions + medial_cart_label : int or float + Scalar value indicating medial tibial cartilage region + lateral_cart_label : int or float + Scalar value indicating lateral tibial cartilage region + scalar_array_name : str, optional + Name of scalar array containing region labels, by default 'labels' + + Returns + ------- + dict + Dictionary containing: + - 'ml_axis': medial-lateral axis vector (medial to lateral, unit vector) + - 'is_axis': inferior-superior axis vector (unit vector pointing superior) + - 'ap_axis': anterior-posterior axis vector (unit vector) + - 'medial_center': medial cartilage center point + - 'lateral_center': lateral cartilage center point + + Examples + -------- + >>> axes = compute_tibia_axes(tibia, med_cart_label=2, lat_cart_label=3) + >>> ml_axis = axes['ml_axis'] + >>> is_axis = axes['is_axis'] + """ + # Get scalar array + region_array = tibia_mesh[scalar_array_name] + + # Extract cartilage points + med_tib_cart_mask = (region_array == medial_cart_label) + lat_tib_cart_mask = (region_array == lateral_cart_label) + + med_tib_cart_points = tibia_mesh.points[med_tib_cart_mask] + lat_tib_cart_points = tibia_mesh.points[lat_tib_cart_mask] + tib_cart_points = np.concatenate([med_tib_cart_points, lat_tib_cart_points], axis=0) + + # Do PCA to get the three axes of the tib_cart_points and take the last + # one as the inf/sup (normal to plateau) + X = tib_cart_points - tib_cart_points.mean(axis=0, keepdims=True) # (N,3) + # PCA via SVD: X = U S Vt, rows of Vt are PCs + U, S, Vt = np.linalg.svd(X, full_matrices=False) + pc1, pc2, pc3 = Vt # already orthonormal + + is_axis = pc3 + + # From the PCA we can't know what is up. Check which side the bone is on + # relative to the cartilage. The opposite direction from bone to cartilage is IS. + mean_tib = np.mean(tibia_mesh.points, axis=0) + mean_cart = np.mean(tib_cart_points, axis=0) + + # Update is_axis direction based on where mean_tib is relative to mean_cart + if np.dot(mean_tib - mean_cart, is_axis) > 0: + is_axis = -is_axis + + # Compute ML axis from cartilage centers + med_tib_center = np.mean(med_tib_cart_points, axis=0) + lat_tib_center = np.mean(lat_tib_cart_points, axis=0) + + ml_axis = lat_tib_center - med_tib_center + ml_axis = ml_axis / np.linalg.norm(ml_axis) + + # Compute AP axis as cross product + # NOTE: AP axis direction is not always same (front vs back) + # without inputting side (right.left). So, left it just as a general axis. + ap_axis = np.cross(ml_axis, is_axis) + ap_axis = ap_axis / np.linalg.norm(ap_axis) + + return { + 'ml_axis': ml_axis, + 'is_axis': is_axis, + 'ap_axis': ap_axis, + 'medial_center': med_tib_center, + 'lateral_center': lat_tib_center, + } + + +def _compute_extrusion_from_points( + cart_points, + men_points, + ml_axis, + side, +): + """ + Compute extrusion by comparing ML extremes of cartilage and meniscus. + + Helper function that projects points onto ML axis and computes the + signed extrusion distance. + + Parameters + ---------- + cart_points : np.ndarray + Cartilage points (N x 3) + men_points : np.ndarray + Meniscus points (M x 3) + ml_axis : np.ndarray + Medial-lateral axis vector + side : str + 'med', 'medial', 'lat', or 'lateral' + + Returns + ------- + float + Extrusion distance in mm (positive = extruded beyond cartilage) + """ + cart_points_ml = np.dot(cart_points, ml_axis) + men_points_ml = np.dot(men_points, ml_axis) + + if side in ['med', 'medial']: + cart_edge = np.min(cart_points_ml) + men_edge = np.min(men_points_ml) + extrusion = cart_edge - men_edge + elif side in ['lat', 'lateral']: + cart_edge = np.max(cart_points_ml) + men_edge = np.max(men_points_ml) + extrusion = men_edge - cart_edge + else: + raise ValueError(f'Invalid side: {side}, must be one of: med, medial, lat, lateral') + + return extrusion + + +def _compute_middle_region_extrusion( + cart_points, + men_points, + ap_axis, + ml_axis, + side, + middle_percentile_range, +): + """ + Compute extrusion using middle percentile range along AP axis. + + This helper function focuses on the central portion of the AP range + to avoid edge effects at the anterior and posterior extremes. + + Parameters + ---------- + cart_points : np.ndarray + Cartilage points (N x 3) + men_points : np.ndarray + Meniscus points (M x 3) + ap_axis : np.ndarray + Anterior-posterior axis vector + ml_axis : np.ndarray + Medial-lateral axis vector + side : str + 'med', 'medial', 'lat', or 'lateral' + middle_percentile_range : float + Fraction of AP range to use (centered on middle) + + Returns + ------- + float + Extrusion distance in mm (positive = extruded beyond cartilage) + """ + # Project cartilage points onto AP axis + cart_points_ap = np.dot(cart_points, ap_axis) + min_cart_ap = np.min(cart_points_ap) + max_cart_ap = np.max(cart_points_ap) + + # Get the middle +/- middle_percentile_range/2 of the cartilage along AP axis + middle_ap_cartilage = (min_cart_ap + max_cart_ap) / 2 + min_max_ap_cartilage_range = max_cart_ap - min_cart_ap + plus_minus_ap_cartilage_range = min_max_ap_cartilage_range * middle_percentile_range / 2 + lower_ap_cartilage = middle_ap_cartilage - plus_minus_ap_cartilage_range + upper_ap_cartilage = middle_ap_cartilage + plus_minus_ap_cartilage_range + + # Get points within the middle AP range for cartilage + ap_cart_indices = (cart_points_ap >= lower_ap_cartilage) & (cart_points_ap <= upper_ap_cartilage) + ml_cart_points = cart_points[ap_cart_indices] + + # Project meniscus points onto AP axis + men_points_ap = np.dot(men_points, ap_axis) + + # Get points within the middle AP range for meniscus + ap_men_indices = (men_points_ap >= lower_ap_cartilage) & (men_points_ap <= upper_ap_cartilage) + ml_men_points = men_points[ap_men_indices] + + # Compute extrusion + extrusion = _compute_extrusion_from_points( + cart_points=ml_cart_points, + men_points=ml_men_points, + ml_axis=ml_axis, + side=side, + ) + + return extrusion + + +def _get_single_compartment_coverage( + tibia_mesh, + meniscus_mesh, + cart_label, + is_direction, + side_name, + scalar_array_name, + ray_cast_length=10.0, +): + """ + Compute meniscal coverage for a single compartment. + + Helper function that performs ray casting from tibia to meniscus and + computes the area of cartilage covered by meniscus. + + Parameters + ---------- + tibia_mesh : BoneMesh or Mesh + Tibia mesh with cartilage region labels + meniscus_mesh : MeniscusMesh or Mesh + Meniscus mesh for this compartment + cart_label : int or float + Label value for this compartment's cartilage + is_direction : np.ndarray + Inferior-superior direction vector for ray casting + side_name : str + Name for this side ('med' or 'lat') used in output keys + scalar_array_name : str + Name of scalar array containing region labels + ray_cast_length : float, optional + Length of rays to cast, by default 20.0 mm + + Returns + ------- + dict + Dictionary containing: + - '{side_name}_cart_men_coverage': coverage percentage + - '{side_name}_cart_men_area': covered area (mm²) + - '{side_name}_cart_area': total cartilage area (mm²) + """ + # Calculate distance from tibia to meniscus along IS direction + tibia_mesh.calc_distance_to_other_mesh( + list_other_meshes=[meniscus_mesh], + ray_cast_length=ray_cast_length, + name=f'{side_name}_men_dist_mm', + direction=is_direction, + ) + + # Create binary masks + binary_mask_men_above = tibia_mesh[f'{side_name}_men_dist_mm'] > 0 + binary_mask_cart = tibia_mesh[scalar_array_name] == cart_label + + tibia_mesh[f'{side_name}_men_above'] = binary_mask_men_above.astype(float) + tibia_mesh[f'{side_name}_cart'] = binary_mask_cart.astype(float) + + # Extract cartilage submesh + tibia_cart = tibia_mesh.copy() + tibia_cart.remove_points(~binary_mask_cart, inplace=True) + tibia_cart.clean(inplace=True) + area_cart = tibia_cart.area + + # Extract covered cartilage submesh + tibia_cart_men = tibia_cart.copy() + tibia_cart_men.remove_points(tibia_cart_men[f'{side_name}_men_above'] == 0, inplace=True) + tibia_cart_men.clean(inplace=True) + area_cart_men = tibia_cart_men.area + + # Calculate coverage percentage + percent_cart_men_coverage = (area_cart_men / area_cart) * 100 if area_cart > 0 else 0.0 + + return { + f'{side_name}_cart_men_coverage': percent_cart_men_coverage, + f'{side_name}_cart_men_area': area_cart_men, + f'{side_name}_cart_area': area_cart, + } + + +# ============================================================================ +# Main Analysis Functions +# ============================================================================ + +def compute_meniscal_extrusion( + tibia_mesh, + medial_meniscus_mesh, + lateral_meniscus_mesh, + medial_cart_label, + lateral_cart_label, + scalar_array_name='labels', + middle_percentile_range=0.1, +): + """ + Compute meniscal extrusion for both medial and lateral menisci. + + Extrusion is computed by comparing the ML extremes of cartilage and meniscus + within the middle portion of the AP range. This avoids edge effects at the + anterior and posterior extremes. + + Parameters + ---------- + tibia_mesh : BoneMesh or Mesh + Tibia mesh with scalar values indicating cartilage regions from reference + medial_meniscus_mesh : MeniscusMesh or Mesh + Medial meniscus mesh + lateral_meniscus_mesh : MeniscusMesh or Mesh + Lateral meniscus mesh + medial_cart_label : int or float + Scalar value indicating medial cartilage region + lateral_cart_label : int or float + Scalar value indicating lateral cartilage region + scalar_array_name : str, optional + Name of scalar array containing region labels, by default 'labels' + middle_percentile_range : float, optional + Fraction of AP range to use for extrusion measurement (centered), by default 0.1 + + Returns + ------- + dict + Dictionary containing extrusion metrics (all distances in mm, positive = extruded): + - 'medial_extrusion_mm': medial extrusion distance + - 'lateral_extrusion_mm': lateral extrusion distance + - 'ml_axis': ML axis vector + - 'ap_axis': AP axis vector + - 'is_axis': IS axis vector + + Notes + ----- + Extrusion sign convention: positive values indicate meniscus extends + beyond the cartilage rim. Negative values indicate the meniscus is contained + within the cartilage boundaries. + + Examples + -------- + >>> results = compute_meniscal_extrusion( + ... tibia, med_meniscus, lat_meniscus, + ... medial_cart_label=2, lateral_cart_label=3 + ... ) + >>> print(f"Medial extrusion: {results['medial_extrusion_mm']:.2f} mm") + """ + # Compute anatomical axes + axes = compute_tibia_axes( + tibia_mesh, + medial_cart_label, + lateral_cart_label, + scalar_array_name + ) + + ml_axis = axes['ml_axis'] + ap_axis = axes['ap_axis'] + is_axis = axes['is_axis'] + + # Get cartilage points + region_array = tibia_mesh[scalar_array_name] + med_cart_indices = region_array == medial_cart_label + lat_cart_indices = region_array == lateral_cart_label + + med_cart_points = tibia_mesh.points[med_cart_indices] + lat_cart_points = tibia_mesh.points[lat_cart_indices] + + # Initialize results + results = { + 'ml_axis': ml_axis, + 'ap_axis': ap_axis, + 'is_axis': is_axis, + } + + # Compute medial extrusion (only if medial meniscus provided) + if medial_meniscus_mesh is not None: + med_men_points = medial_meniscus_mesh.points + med_men_extrusion = _compute_middle_region_extrusion( + cart_points=med_cart_points, + men_points=med_men_points, + ap_axis=ap_axis, + ml_axis=ml_axis, + side='med', + middle_percentile_range=middle_percentile_range, + ) + results['medial_extrusion_mm'] = med_men_extrusion + + # Compute lateral extrusion (only if lateral meniscus provided) + if lateral_meniscus_mesh is not None: + lat_men_points = lateral_meniscus_mesh.points + lat_men_extrusion = _compute_middle_region_extrusion( + cart_points=lat_cart_points, + men_points=lat_men_points, + ap_axis=ap_axis, + ml_axis=ml_axis, + side='lat', + middle_percentile_range=middle_percentile_range, + ) + results['lateral_extrusion_mm'] = lat_men_extrusion + + return results + + +def compute_meniscal_coverage( + tibia_mesh, + medial_meniscus_mesh, + lateral_meniscus_mesh, + medial_cart_label, + lateral_cart_label, + scalar_array_name='labels', + ray_cast_length=10.0, +): + """ + Compute meniscal coverage using superior-inferior ray casting. + + Coverage is computed by casting rays in the IS direction from tibial cartilage + reference points and checking for meniscus intersections. Areas are computed + using PyVista's mesh area calculations. + + Parameters + ---------- + tibia_mesh : BoneMesh or Mesh + Tibia mesh with scalar values indicating cartilage regions from reference + medial_meniscus_mesh : MeniscusMesh or Mesh + Medial meniscus mesh + lateral_meniscus_mesh : MeniscusMesh or Mesh + Lateral meniscus mesh + medial_cart_label : int or float + Scalar value indicating medial cartilage region + lateral_cart_label : int or float + Scalar value indicating lateral cartilage region + scalar_array_name : str, optional + Name of scalar array containing region labels, by default 'labels' + ray_cast_length : float, optional + Length of rays to cast in IS direction, by default 20.0 mm + + Returns + ------- + dict + Dictionary containing coverage metrics: + - 'medial_coverage_percent': percentage of medial cartilage covered by meniscus + - 'lateral_coverage_percent': percentage of lateral cartilage covered by meniscus + - 'medial_covered_area_mm2': area of medial cartilage covered (mm²) + - 'lateral_covered_area_mm2': area of lateral cartilage covered (mm²) + - 'medial_total_area_mm2': total medial cartilage area (mm²) + - 'lateral_total_area_mm2': total lateral cartilage area (mm²) + + Examples + -------- + >>> results = compute_meniscal_coverage( + ... tibia, med_meniscus, lat_meniscus, + ... medial_cart_label=2, lateral_cart_label=3 + ... ) + >>> print(f"Medial coverage: {results['medial_coverage_percent']:.1f}%") + """ + # Compute IS axis + axes = compute_tibia_axes( + tibia_mesh, + medial_cart_label, + lateral_cart_label, + scalar_array_name + ) + is_direction = axes['is_axis'] + + # Initialize results + results = {} + + # Compute medial coverage (only if medial meniscus provided) + if medial_meniscus_mesh is not None: + med_coverage = _get_single_compartment_coverage( + tibia_mesh=tibia_mesh, + meniscus_mesh=medial_meniscus_mesh, + cart_label=medial_cart_label, + is_direction=is_direction, + side_name='med', + scalar_array_name=scalar_array_name, + ray_cast_length=ray_cast_length, + ) + results['medial_coverage_percent'] = med_coverage['med_cart_men_coverage'] + results['medial_covered_area_mm2'] = med_coverage['med_cart_men_area'] + results['medial_total_area_mm2'] = med_coverage['med_cart_area'] + + # Compute lateral coverage (only if lateral meniscus provided) + if lateral_meniscus_mesh is not None: + lat_coverage = _get_single_compartment_coverage( + tibia_mesh=tibia_mesh, + meniscus_mesh=lateral_meniscus_mesh, + cart_label=lateral_cart_label, + is_direction=is_direction, + side_name='lat', + scalar_array_name=scalar_array_name, + ray_cast_length=ray_cast_length, + ) + results['lateral_coverage_percent'] = lat_coverage['lat_cart_men_coverage'] + results['lateral_covered_area_mm2'] = lat_coverage['lat_cart_men_area'] + results['lateral_total_area_mm2'] = lat_coverage['lat_cart_area'] + + return results + + +def analyze_meniscal_metrics( + tibia_mesh, + medial_meniscus_mesh, + lateral_meniscus_mesh, + medial_cart_label, + lateral_cart_label, + scalar_array_name='labels', + middle_percentile_range=0.1, + ray_cast_length=10.0, +): + """ + Comprehensive meniscal analysis computing both extrusion and coverage metrics. + + This is the main function for complete meniscal analysis. It computes + meniscal extrusion using the middle AP region and meniscal coverage + using IS-direction ray casting. + + Parameters + ---------- + tibia_mesh : BoneMesh or Mesh + Tibia mesh with scalar values indicating cartilage regions from reference + medial_meniscus_mesh : MeniscusMesh or Mesh + Medial meniscus mesh + lateral_meniscus_mesh : MeniscusMesh or Mesh + Lateral meniscus mesh + medial_cart_label : int or float + Scalar value indicating medial cartilage region + lateral_cart_label : int or float + Scalar value indicating lateral cartilage region + scalar_array_name : str, optional + Name of scalar array containing region labels, by default 'labels' + middle_percentile_range : float, optional + Fraction of AP range to use for extrusion measurement, by default 0.1 + ray_cast_length : float, optional + Length of rays to cast for coverage analysis, by default 20.0 mm + + Returns + ------- + dict + Dictionary containing all extrusion and coverage metrics: + + Extrusion metrics (mm, positive = extruded beyond cartilage rim): + - 'medial_extrusion_mm': medial extrusion distance + - 'lateral_extrusion_mm': lateral extrusion distance + + Coverage metrics: + - 'medial_coverage_percent': percentage of medial cartilage covered + - 'lateral_coverage_percent': percentage of lateral cartilage covered + - 'medial_covered_area_mm2': medial cartilage covered area (mm²) + - 'lateral_covered_area_mm2': lateral cartilage covered area (mm²) + - 'medial_total_area_mm2': total medial cartilage area (mm²) + - 'lateral_total_area_mm2': total lateral cartilage area (mm²) + + Reference frame: + - 'ml_axis': medial-lateral axis vector + - 'ap_axis': anterior-posterior axis vector + - 'is_axis': inferior-superior axis vector + + Notes + ----- + All meshes are automatically oriented with consistent normals before analysis. + + Examples + -------- + >>> results = analyze_meniscal_metrics( + ... tibia, med_meniscus, lat_meniscus, + ... medial_cart_label=2, lateral_cart_label=3 + ... ) + >>> print(f"Medial extrusion: {results['medial_extrusion_mm']:.2f} mm") + >>> print(f"Medial coverage: {results['medial_coverage_percent']:.1f}%") + """ + # Ensure tibia mesh is properly prepared + tibia_mesh.compute_normals(auto_orient_normals=True, inplace=True) + + # Ensure meniscus meshes are properly prepared (only if not None) + if medial_meniscus_mesh is not None: + medial_meniscus_mesh.compute_normals(auto_orient_normals=True, inplace=True) + if lateral_meniscus_mesh is not None: + lateral_meniscus_mesh.compute_normals(auto_orient_normals=True, inplace=True) + + # Check that at least one meniscus is provided + if medial_meniscus_mesh is None and lateral_meniscus_mesh is None: + raise ValueError("At least one meniscus mesh must be provided") + + # Compute extrusion metrics (only for menisci that are present) + extrusion_results = compute_meniscal_extrusion( + tibia_mesh, medial_meniscus_mesh, lateral_meniscus_mesh, + medial_cart_label, lateral_cart_label, scalar_array_name, middle_percentile_range + ) + + # Compute coverage metrics (only for menisci that are present) + coverage_results = compute_meniscal_coverage( + tibia_mesh, medial_meniscus_mesh, lateral_meniscus_mesh, + medial_cart_label, lateral_cart_label, scalar_array_name, ray_cast_length + ) + + # Combine results + results = {**extrusion_results, **coverage_results} + + return results + + + + diff --git a/pymskt/mesh/meshes.py b/pymskt/mesh/meshes.py index a84c689..44d9af4 100644 --- a/pymskt/mesh/meshes.py +++ b/pymskt/mesh/meshes.py @@ -32,6 +32,7 @@ gaussian_smooth_surface_scalars, get_cartilage_properties_at_points, get_distance_other_surface_at_points, + get_distance_other_surface_at_points_along_unit_vector, get_largest_connected_component, get_mesh_edge_lengths, get_mesh_physical_point_coords, @@ -802,6 +803,8 @@ def calc_distance_to_other_mesh( ray_cast_length=10.0, percent_ray_length_opposite_direction=0.25, name="thickness (mm)", + direction=None, + ): """ Using bone mesh (`_mesh`) and the list of cartilage meshes (`list_cartilage_meshes`) @@ -841,12 +844,26 @@ def calc_distance_to_other_mesh( # iterate over meshes and add their thicknesses to the thicknesses list. for other_mesh in list_other_meshes: - node_data = get_distance_other_surface_at_points( - self, - other_mesh, - ray_cast_length=ray_cast_length, - percent_ray_length_opposite_direction=percent_ray_length_opposite_direction, - ) + if direction is None: + node_data = get_distance_other_surface_at_points( + self, + other_mesh, + ray_cast_length=ray_cast_length, + percent_ray_length_opposite_direction=percent_ray_length_opposite_direction, + ) + + elif isinstance(direction, (np.ndarray, list, tuple)): + direction = np.array(direction) + direction = direction / np.linalg.norm(direction) + node_data = get_distance_other_surface_at_points_along_unit_vector( + self, + other_mesh, + unit_vector=direction, + ray_cast_length=ray_cast_length, + percent_ray_length_opposite_direction=percent_ray_length_opposite_direction, + ) + else: + raise ValueError(f"direction must be a numpy array, list, or tuple and received: {type(direction)}") distances += node_data @@ -1297,6 +1314,7 @@ def __init__( min_n_pixels=5000, list_cartilage_meshes=None, list_cartilage_labels=None, + dict_cartilage_labels=None, list_articular_surfaces=None, crop_percent=None, bone="femur", @@ -1325,6 +1343,11 @@ def __init__( list_cartilage_labels : list, optional List of `int` values that represent the different cartilage regions of interest appropriate for a single bone, by default None + dict_cartilage_labels : dict, optional + Dictionary mapping cartilage region names to label values. + For tibia: {'medial': 2, 'lateral': 3}. + This enables cleaner API for meniscal analysis without repeatedly + specifying labels, by default None crop_percent : float, optional Proportion value to crop long-axis of bone so it is proportional to the width of the bone for standardization purposes, by default 1.0 @@ -1332,13 +1355,19 @@ def __init__( String indicating what bone is being analyzed so that cropping can be applied appropriatey. {'femur', 'tibia'}, by default 'femur'. Patella is not an option because we do not need to crop for the patella. + tibia_idx : int, optional + Label index for tibia in segmentation (for registrations), by default None """ self._crop_percent = crop_percent self._bone = bone self._list_cartilage_meshes = list_cartilage_meshes self._list_cartilage_labels = list_cartilage_labels + self._dict_cartilage_labels = dict_cartilage_labels self._list_articular_surfaces = list_articular_surfaces self._tibia_idx = tibia_idx + self._meniscus_meshes = {} # Dictionary to store medial/lateral menisci + self._meniscal_outcomes = None # Cache for computed meniscal metrics + self._meniscal_cart_labels = None # Cache cartilage labels used for computation super().__init__( mesh=mesh, @@ -1362,7 +1391,11 @@ def copy(self, deep=True): copy_.bone = self.bone copy_.list_cartilage_meshes = self.list_cartilage_meshes copy_.list_cartilage_labels = self.list_cartilage_labels + copy_.dict_cartilage_labels = self.dict_cartilage_labels.copy() if self.dict_cartilage_labels else None copy_.list_articular_surfaces = self.list_articular_surfaces + copy_._meniscus_meshes = self._meniscus_meshes.copy() + copy_._meniscal_outcomes = self._meniscal_outcomes + copy_._meniscal_cart_labels = self._meniscal_cart_labels return copy_ def create_mesh( @@ -1477,7 +1510,8 @@ def create_cartilage_meshes( """ self._list_cartilage_meshes = [] - for cart_label_idx in self._list_cartilage_labels: + # Use property to handle both list_cartilage_labels and dict_cartilage_labels + for cart_label_idx in self.list_cartilage_labels: seg_array_view = sitk.GetArrayViewFromImage(self._seg_image) n_pixels_with_cart = np.sum(seg_array_view == cart_label_idx) if n_pixels_with_cart == 0: @@ -1556,9 +1590,10 @@ def calc_cartilage_thickness( self._list_cartilage_labels = list_cartilage_labels # If no cartilage stuff provided, then cant do this function - raise exception. - if (self._list_cartilage_meshes is None) & (self._list_cartilage_labels is None): + # Check using property to handle both list_cartilage_labels and dict_cartilage_labels + if (self._list_cartilage_meshes is None) & (self.list_cartilage_labels is None): raise Exception( - "No cartilage meshes or list of cartilage labels are provided! - These can be provided either to the class function `calc_cartilage_thickness` directly, or can be specified at the time of instantiating the `BoneMesh` class." + "No cartilage meshes or list of cartilage labels are provided! - These can be provided either to the class function `calc_cartilage_thickness` directly, or can be specified at the time of instantiating the `BoneMesh` class via list_cartilage_labels or dict_cartilage_labels." ) # if cartilage meshes don't exist yet, then make them. @@ -1950,13 +1985,25 @@ def list_cartilage_labels(self): """ Convenience function to get the list of labels for cartilage tissues associated with this bone. + + If list_cartilage_labels was not explicitly set but dict_cartilage_labels was, + this will return the values from dict_cartilage_labels in order. Returns ------- list list of `int`s for the cartilage tissues associated with this bone. """ - return self._list_cartilage_labels + # If explicit list provided, use it + if self._list_cartilage_labels is not None: + return self._list_cartilage_labels + + # Fall back to values from dict if available + if self._dict_cartilage_labels is not None: + return list(self._dict_cartilage_labels.values()) + + # Neither provided + return None @list_cartilage_labels.setter def list_cartilage_labels(self, new_list_cartilage_labels): @@ -1981,6 +2028,36 @@ def list_cartilage_labels(self, new_list_cartilage_labels): ] self._list_cartilage_labels = new_list_cartilage_labels + @property + def dict_cartilage_labels(self): + """ + Get the dictionary mapping cartilage region names to label values. + + Returns + ------- + dict or None + Dictionary mapping region names (e.g., 'medial', 'lateral') to label values. + Returns None if not set. + """ + return self._dict_cartilage_labels + + @dict_cartilage_labels.setter + def dict_cartilage_labels(self, new_dict_cartilage_labels): + """ + Set the dictionary mapping cartilage region names to label values. + + Parameters + ---------- + new_dict_cartilage_labels : dict or None + Dictionary mapping region names to label values. + For tibia: {'medial': 2, 'lateral': 3} + """ + if new_dict_cartilage_labels is not None and not isinstance(new_dict_cartilage_labels, dict): + raise TypeError( + f"dict_cartilage_labels must be a dict or None, got {type(new_dict_cartilage_labels)}" + ) + self._dict_cartilage_labels = new_dict_cartilage_labels + @property def crop_percent(self): """ @@ -2040,3 +2117,378 @@ def bone(self, new_bone): if not isinstance(new_bone, str): raise TypeError(f"New bone provided is type {type(new_bone)} - expected `str`") self._bone = new_bone + + # ============================================================================ + # Meniscus Analysis Methods (Tibia-specific) + # NOTE: Could be refactored into TibiaMesh class inheriting from BoneMesh + # ============================================================================ + + def set_menisci( + self, + medial_meniscus=None, + medial_cart_label=None, + lateral_meniscus=None, + lateral_cart_label=None, + scalar_array_name='labels', + ): + """ + Associate meniscus meshes and cartilage labels for meniscal analysis. + + This method stores references to meniscus meshes and their corresponding + cartilage labels. You can set one or both menisci, but BOTH cartilage labels + must be provided because tibial axes computation requires both cartilage regions. + + If dict_cartilage_labels was set during initialization, labels can be + automatically inferred and don't need to be explicitly provided. + + Parameters + ---------- + medial_meniscus : MeniscusMesh or Mesh, optional + Medial meniscus mesh, by default None + medial_cart_label : int or float, optional + Label value for medial tibial cartilage region. If None, uses value + from dict_cartilage_labels['medial'] if available. + lateral_meniscus : MeniscusMesh or Mesh, optional + Lateral meniscus mesh, by default None + lateral_cart_label : int or float, optional + Label value for lateral tibial cartilage region. If None, uses value + from dict_cartilage_labels['lateral'] if available. + scalar_array_name : str, optional + Name of scalar array containing region labels, by default 'labels' + + Raises + ------ + ValueError + If no menisci are provided or if cartilage labels cannot be determined + + Examples + -------- + >>> # With dict_cartilage_labels set at initialization + >>> tibia = BoneMesh( + ... path_seg_image='tibia.nrrd', label_idx=6, + ... dict_cartilage_labels={'medial': 2, 'lateral': 3} + ... ) + >>> tibia.set_menisci( + ... medial_meniscus=med_men, + ... lateral_meniscus=lat_men + ... ) # Labels auto-inferred! + + >>> # Or provide labels explicitly (overrides dict_cartilage_labels) + >>> tibia.set_menisci( + ... medial_meniscus=med_men, medial_cart_label=2, + ... lateral_meniscus=lat_men, lateral_cart_label=3 + ... ) + """ + # Must provide at least one meniscus + if medial_meniscus is None and lateral_meniscus is None: + raise ValueError("At least one meniscus must be provided") + + # Try to get labels from dict_cartilage_labels if not explicitly provided + if medial_cart_label is None and self._dict_cartilage_labels: + medial_cart_label = self._dict_cartilage_labels.get('medial') + if lateral_cart_label is None and self._dict_cartilage_labels: + lateral_cart_label = self._dict_cartilage_labels.get('lateral') + + # Both cartilage labels are required for axes computation + if medial_cart_label is None or lateral_cart_label is None: + raise ValueError( + "Both medial_cart_label and lateral_cart_label must be provided. " + "Tibial axes computation requires both cartilage regions, even if only " + "one meniscus is being analyzed. Either provide them explicitly or set " + "dict_cartilage_labels={'medial': X, 'lateral': Y} during initialization." + ) + + # Store menisci + if medial_meniscus is not None: + self._meniscus_meshes['medial'] = medial_meniscus + if lateral_meniscus is not None: + self._meniscus_meshes['lateral'] = lateral_meniscus + + # Store labels (always store both since both are required) + self._meniscal_cart_labels = { + 'medial': medial_cart_label, + 'lateral': lateral_cart_label, + 'scalar_array_name': scalar_array_name, + } + + # Clear cached outcomes when menisci/labels are updated + self._meniscal_outcomes = None + + def compute_meniscal_outcomes( + self, + medial_cart_label=None, + lateral_cart_label=None, + scalar_array_name=None, + middle_percentile_range=0.1, + ray_cast_length=20.0, + force_recompute=False, + ): + """ + Compute meniscal extrusion and coverage metrics. + + This method computes extrusion (how far meniscus extends beyond + cartilage rim) and coverage (percentage of cartilage covered by meniscus) + for menisci that have been set via set_menisci(). Can compute for one + or both compartments depending on what was set. + + Parameters + ---------- + medial_cart_label : int or float, optional + Label value for medial tibial cartilage region. + If None, uses label from set_menisci() call. + lateral_cart_label : int or float, optional + Label value for lateral tibial cartilage region. + If None, uses label from set_menisci() call. + scalar_array_name : str, optional + Name of scalar array containing region labels. + If None, uses value from set_menisci() call, by default 'labels' + middle_percentile_range : float, optional + Fraction of AP range to use for extrusion measurement, by default 0.1 + ray_cast_length : float, optional + Length of rays to cast for coverage analysis, by default 20.0 mm + force_recompute : bool, optional + Force recomputation even if cached results exist, by default False + + Returns + ------- + dict + Dictionary containing meniscal metrics for available compartments: + - 'medial_extrusion_mm': medial extrusion distance (mm) [if medial set] + - 'lateral_extrusion_mm': lateral extrusion distance (mm) [if lateral set] + - 'medial_coverage_percent': medial coverage percentage [if medial set] + - 'lateral_coverage_percent': lateral coverage percentage [if lateral set] + - 'medial_covered_area_mm2': medial covered area (mm²) [if medial set] + - 'lateral_covered_area_mm2': lateral covered area (mm²) [if lateral set] + - 'medial_total_area_mm2': total medial cartilage area (mm²) [if medial set] + - 'lateral_total_area_mm2': total lateral cartilage area (mm²) [if lateral set] + - 'ml_axis': medial-lateral axis vector + - 'ap_axis': anterior-posterior axis vector + - 'is_axis': inferior-superior axis vector + + Raises + ------ + ValueError + If no menisci are set or if required labels cannot be determined + + Examples + -------- + >>> # Set menisci with labels, then compute + >>> tibia.set_menisci( + ... medial_meniscus=med_men, medial_cart_label=2, + ... lateral_meniscus=lat_men, lateral_cart_label=3 + ... ) + >>> results = tibia.compute_meniscal_outcomes() + >>> print(f"Medial extrusion: {results['medial_extrusion_mm']:.2f} mm") + + >>> # Or provide labels explicitly + >>> results = tibia.compute_meniscal_outcomes( + ... medial_cart_label=2, lateral_cart_label=3 + ... ) + """ + # Return cached results if available and not forcing recompute + if self._meniscal_outcomes is not None and not force_recompute: + return self._meniscal_outcomes + + # Check that at least one meniscus is set + if not self._meniscus_meshes: + raise ValueError( + "No menisci have been set. Use set_menisci() to associate meniscus " + "meshes and cartilage labels before computing outcomes." + ) + + # Determine which menisci to compute for + has_medial = 'medial' in self._meniscus_meshes + has_lateral = 'lateral' in self._meniscus_meshes + + # Get labels (from parameters or cached values) + # Both labels are ALWAYS required for axes computation + if medial_cart_label is None: + if self._meniscal_cart_labels and 'medial' in self._meniscal_cart_labels: + medial_cart_label = self._meniscal_cart_labels['medial'] + else: + raise ValueError( + "medial_cart_label must be provided either in compute_meniscal_outcomes() " + "or previously in set_menisci(). Both cartilage labels are required for " + "tibial axes computation, even if only one meniscus is being analyzed." + ) + + if lateral_cart_label is None: + if self._meniscal_cart_labels and 'lateral' in self._meniscal_cart_labels: + lateral_cart_label = self._meniscal_cart_labels['lateral'] + else: + raise ValueError( + "lateral_cart_label must be provided either in compute_meniscal_outcomes() " + "or previously in set_menisci(). Both cartilage labels are required for " + "tibial axes computation, even if only one meniscus is being analyzed." + ) + + # Get scalar array name + if scalar_array_name is None: + if self._meniscal_cart_labels and 'scalar_array_name' in self._meniscal_cart_labels: + scalar_array_name = self._meniscal_cart_labels['scalar_array_name'] + else: + scalar_array_name = 'labels' + + # Import analysis functions + from pymskt.mesh.mesh_meniscus import analyze_meniscal_metrics + + # Always use the combined analysis function + # It will handle single meniscus cases by only computing metrics for present menisci + self._meniscal_outcomes = analyze_meniscal_metrics( + tibia_mesh=self, + medial_meniscus_mesh=self._meniscus_meshes.get('medial'), + lateral_meniscus_mesh=self._meniscus_meshes.get('lateral'), + medial_cart_label=medial_cart_label, + lateral_cart_label=lateral_cart_label, + scalar_array_name=scalar_array_name, + middle_percentile_range=middle_percentile_range, + ray_cast_length=ray_cast_length, + ) + + return self._meniscal_outcomes + + @property + def med_men_extrusion(self): + """ + Get medial meniscal extrusion value in mm. + + Automatically computes outcomes on first access if not already computed. + Positive values indicate meniscus extends beyond cartilage rim. + + Returns + ------- + float + Medial meniscal extrusion distance in mm + + Raises + ------ + ValueError + If menisci haven't been set or if medial meniscus was not included + + Examples + -------- + >>> tibia.set_menisci(medial_meniscus=med_men, lateral_meniscus=lat_men) + >>> print(f"Medial extrusion: {tibia.med_men_extrusion:.2f} mm") # Auto-computes! + """ + # Auto-compute on first access if not already computed + if self._meniscal_outcomes is None: + try: + self.compute_meniscal_outcomes() + except Exception as e: + raise ValueError( + f"Cannot compute meniscal outcomes automatically: {str(e)}\n" + "Ensure menisci are set via set_menisci() with appropriate labels." + ) + + if 'medial_extrusion_mm' not in self._meniscal_outcomes: + raise ValueError("Medial meniscus was not included in the analysis") + return self._meniscal_outcomes['medial_extrusion_mm'] + + @property + def lat_men_extrusion(self): + """ + Get lateral meniscal extrusion value in mm. + + Automatically computes outcomes on first access if not already computed. + Positive values indicate meniscus extends beyond cartilage rim. + + Returns + ------- + float + Lateral meniscal extrusion distance in mm + + Raises + ------ + ValueError + If menisci haven't been set or if lateral meniscus was not included + + Examples + -------- + >>> tibia.set_menisci(medial_meniscus=med_men, lateral_meniscus=lat_men) + >>> print(f"Lateral extrusion: {tibia.lat_men_extrusion:.2f} mm") # Auto-computes! + """ + # Auto-compute on first access if not already computed + if self._meniscal_outcomes is None: + try: + self.compute_meniscal_outcomes() + except Exception as e: + raise ValueError( + f"Cannot compute meniscal outcomes automatically: {str(e)}\n" + "Ensure menisci are set via set_menisci() with appropriate labels." + ) + + if 'lateral_extrusion_mm' not in self._meniscal_outcomes: + raise ValueError("Lateral meniscus was not included in the analysis") + return self._meniscal_outcomes['lateral_extrusion_mm'] + + @property + def med_men_coverage(self): + """ + Get medial meniscal coverage percentage. + + Automatically computes outcomes on first access if not already computed. + + Returns + ------- + float + Percentage of medial cartilage covered by meniscus + + Raises + ------ + ValueError + If menisci haven't been set or if medial meniscus was not included + + Examples + -------- + >>> tibia.set_menisci(medial_meniscus=med_men, lateral_meniscus=lat_men) + >>> print(f"Medial coverage: {tibia.med_men_coverage:.1f}%") # Auto-computes! + """ + # Auto-compute on first access if not already computed + if self._meniscal_outcomes is None: + try: + self.compute_meniscal_outcomes() + except Exception as e: + raise ValueError( + f"Cannot compute meniscal outcomes automatically: {str(e)}\n" + "Ensure menisci are set via set_menisci() with appropriate labels." + ) + + if 'medial_coverage_percent' not in self._meniscal_outcomes: + raise ValueError("Medial meniscus was not included in the analysis") + return self._meniscal_outcomes['medial_coverage_percent'] + + @property + def lat_men_coverage(self): + """ + Get lateral meniscal coverage percentage. + + Automatically computes outcomes on first access if not already computed. + + Returns + ------- + float + Percentage of lateral cartilage covered by meniscus + + Raises + ------ + ValueError + If menisci haven't been set or if lateral meniscus was not included + + Examples + -------- + >>> tibia.set_menisci(medial_meniscus=med_men, lateral_meniscus=lat_men) + >>> print(f"Lateral coverage: {tibia.lat_men_coverage:.1f}%") # Auto-computes! + """ + # Auto-compute on first access if not already computed + if self._meniscal_outcomes is None: + try: + self.compute_meniscal_outcomes() + except Exception as e: + raise ValueError( + f"Cannot compute meniscal outcomes automatically: {str(e)}\n" + "Ensure menisci are set via set_menisci() with appropriate labels." + ) + + if 'lateral_coverage_percent' not in self._meniscal_outcomes: + raise ValueError("Lateral meniscus was not included in the analysis") + return self._meniscal_outcomes['lateral_coverage_percent'] diff --git a/testing/mesh/meshMeniscus/MeniscusMesh_create_test.py b/testing/mesh/meshMeniscus/MeniscusMesh_create_test.py new file mode 100644 index 0000000..bbeb8a6 --- /dev/null +++ b/testing/mesh/meshMeniscus/MeniscusMesh_create_test.py @@ -0,0 +1,33 @@ +""" +Test file for MeniscusMesh class creation and basic functionality. + +TODO: Implement tests for: +- Creating MeniscusMesh from segmentation image +- Creating MeniscusMesh from existing mesh +- Setting and getting meniscus_type property +- Testing with different meniscus types ('medial', 'lateral') +- Validating that invalid meniscus_type raises appropriate error +""" + +import numpy as np +import pytest +from pymskt.mesh import MeniscusMesh + + +def test_meniscus_mesh_creation(): + """TODO: Test basic MeniscusMesh instantiation.""" + pass + + +def test_meniscus_type_property(): + """TODO: Test meniscus_type property setter and getter.""" + pass + + +def test_meniscus_type_validation(): + """TODO: Test that invalid meniscus_type values raise ValueError.""" + pass + + + + diff --git a/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py b/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py new file mode 100644 index 0000000..223e948 --- /dev/null +++ b/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py @@ -0,0 +1,52 @@ +""" +Test file for meniscal coverage computation. + +TODO: Implement tests for: +- Testing coverage calculation with synthetic data (known coverage values) +- Testing with 100% coverage (meniscus covers all cartilage) +- Testing with 0% coverage (no meniscus above cartilage) +- Testing with partial coverage +- Testing SI ray casting functionality +- Testing area calculation accuracy +- Testing with different SI tolerance values +- Edge cases: empty compartments, missing meniscus data +""" + +import numpy as np +import pytest +from pymskt.mesh import Mesh, MeniscusMesh +from pymskt.mesh.mesh_meniscus import compute_meniscal_coverage_si_ray + + +def test_coverage_synthetic_data(): + """TODO: Test coverage calculation with synthetic tibia and meniscus meshes.""" + pass + + +def test_coverage_full_coverage(): + """TODO: Test case where meniscus fully covers cartilage (100%).""" + pass + + +def test_coverage_no_coverage(): + """TODO: Test case where meniscus does not cover cartilage (0%).""" + pass + + +def test_coverage_partial_coverage(): + """TODO: Test case with partial meniscal coverage.""" + pass + + +def test_coverage_area_calculation(): + """TODO: Test that area calculation using cell sizes is accurate.""" + pass + + +def test_coverage_si_tolerance(): + """TODO: Test coverage with different SI tolerance values.""" + pass + + + + diff --git a/testing/mesh/meshMeniscus/compute_meniscal_extrusion_test.py b/testing/mesh/meshMeniscus/compute_meniscal_extrusion_test.py new file mode 100644 index 0000000..0650f57 --- /dev/null +++ b/testing/mesh/meshMeniscus/compute_meniscal_extrusion_test.py @@ -0,0 +1,491 @@ +""" +Test file for meniscal extrusion computation. + +TODO: Implement tests for: +- Testing extrusion calculation with synthetic data (known extrusion values) +- Testing with no extrusion (meniscus within cartilage rim) +- Edge cases: empty compartments, missing meniscus data +""" + +import numpy as np +import pytest +import os +from pymskt.mesh import Mesh, BoneMesh, MeniscusMesh +from pymskt.mesh.mesh_meniscus import compute_meniscal_extrusion + + +# ============================================================================ +# Fixtures for Meniscus Shift Tests +# ============================================================================ + +@pytest.fixture +def tibia_with_menisci(): + """ + Fixture that provides a tibia mesh with cartilage regions and menisci. + + Returns a tuple of (tibia, medial_meniscus, lateral_meniscus, baseline_results) + """ + test_dir = os.path.dirname(os.path.abspath(__file__)) + data_dir = os.path.join(test_dir, '..', '..', '..', 'data') + path_segmentation = os.path.join(data_dir, 'SAG_3D_DESS_RIGHT_bones_cart_men_fib-label.nrrd') + + if not os.path.exists(path_segmentation): + pytest.skip(f"Test data not found: {path_segmentation}") + + # Create tibia mesh with cartilage regions + tibia = BoneMesh( + path_seg_image=path_segmentation, + label_idx=6, + list_cartilage_labels=[2, 3] + ) + tibia.create_mesh() + tibia.calc_cartilage_thickness() + tibia.assign_cartilage_regions() + + # Create meniscus meshes + med_meniscus = Mesh( + path_seg_image=path_segmentation, + label_idx=10, + ) + med_meniscus.create_mesh() + med_meniscus.consistent_faces() + + lat_meniscus = Mesh( + path_seg_image=path_segmentation, + label_idx=9, + ) + lat_meniscus.create_mesh() + lat_meniscus.consistent_faces() + + # Compute baseline extrusion + baseline_results = compute_meniscal_extrusion( + tibia_mesh=tibia, + medial_meniscus_mesh=med_meniscus, + lateral_meniscus_mesh=lat_meniscus, + medial_cart_label=2, + lateral_cart_label=3, + scalar_array_name='labels', + middle_percentile_range=0.1, + ) + + return tibia, med_meniscus, lat_meniscus, baseline_results + + +# ============================================================================ +# Meniscus Shift Tests +# ============================================================================ + +def test_medial_meniscus_shift_medially_increases_extrusion(tibia_with_menisci): + """ + Test that shifting medial meniscus medially increases extrusion. + + When the medial meniscus is shifted 5mm medially (away from midline), + the extrusion value should increase (more positive or less negative). + """ + tibia, med_meniscus, lat_meniscus, baseline_results = tibia_with_menisci + baseline_extrusion = baseline_results['medial_extrusion_mm'] + + # Shift medial meniscus medially (5mm in +X direction for right knee) + med_meniscus_shifted = med_meniscus.copy() + med_meniscus_shifted.points = med_meniscus_shifted.points + np.array([5.0, 0.0, 0.0]) + + # Compute extrusion with shifted meniscus + results = compute_meniscal_extrusion( + tibia_mesh=tibia, + medial_meniscus_mesh=med_meniscus_shifted, + lateral_meniscus_mesh=lat_meniscus, + medial_cart_label=2, + lateral_cart_label=3, + scalar_array_name='labels', + middle_percentile_range=0.1, + ) + + # Verify extrusion increased + assert results['medial_extrusion_mm'] > baseline_extrusion, \ + f"Medial shift should increase extrusion. " \ + f"Baseline: {baseline_extrusion:.2f}, Shifted: {results['medial_extrusion_mm']:.2f}" + + print(f"\n✓ Medial meniscus medial shift test passed!") + print(f" Baseline: {baseline_extrusion:.2f} mm") + print(f" After medial shift: {results['medial_extrusion_mm']:.2f} mm") + print(f" Increase: {results['medial_extrusion_mm'] - baseline_extrusion:.2f} mm") + + +def test_medial_meniscus_shift_laterally_decreases_extrusion(tibia_with_menisci): + """ + Test that shifting medial meniscus laterally decreases extrusion. + + When the medial meniscus is shifted 5mm laterally (toward midline), + the extrusion value should decrease (less positive or more negative). + """ + tibia, med_meniscus, lat_meniscus, baseline_results = tibia_with_menisci + baseline_extrusion = baseline_results['medial_extrusion_mm'] + + # Shift medial meniscus laterally (5mm in -X direction for right knee) + med_meniscus_shifted = med_meniscus.copy() + med_meniscus_shifted.points = med_meniscus_shifted.points - np.array([5.0, 0.0, 0.0]) + + # Compute extrusion with shifted meniscus + results = compute_meniscal_extrusion( + tibia_mesh=tibia, + medial_meniscus_mesh=med_meniscus_shifted, + lateral_meniscus_mesh=lat_meniscus, + medial_cart_label=2, + lateral_cart_label=3, + scalar_array_name='labels', + middle_percentile_range=0.1, + ) + + # Verify extrusion decreased + assert results['medial_extrusion_mm'] < baseline_extrusion, \ + f"Lateral shift should decrease extrusion. " \ + f"Baseline: {baseline_extrusion:.2f}, Shifted: {results['medial_extrusion_mm']:.2f}" + + print(f"\n✓ Medial meniscus lateral shift test passed!") + print(f" Baseline: {baseline_extrusion:.2f} mm") + print(f" After lateral shift: {results['medial_extrusion_mm']:.2f} mm") + print(f" Decrease: {baseline_extrusion - results['medial_extrusion_mm']:.2f} mm") + + +def test_lateral_meniscus_shift_laterally_increases_extrusion(tibia_with_menisci): + """ + Test that shifting lateral meniscus laterally increases extrusion. + + When the lateral meniscus is shifted 5mm laterally (away from midline), + the extrusion value should increase (more positive or less negative). + """ + tibia, med_meniscus, lat_meniscus, baseline_results = tibia_with_menisci + baseline_extrusion = baseline_results['lateral_extrusion_mm'] + + # Shift lateral meniscus laterally (5mm in -X direction for right knee) + lat_meniscus_shifted = lat_meniscus.copy() + lat_meniscus_shifted.points = lat_meniscus_shifted.points - np.array([5.0, 0.0, 0.0]) + + # Compute extrusion with shifted meniscus + results = compute_meniscal_extrusion( + tibia_mesh=tibia, + medial_meniscus_mesh=med_meniscus, + lateral_meniscus_mesh=lat_meniscus_shifted, + medial_cart_label=2, + lateral_cart_label=3, + scalar_array_name='labels', + middle_percentile_range=0.1, + ) + + # Verify extrusion increased + assert results['lateral_extrusion_mm'] > baseline_extrusion, \ + f"Lateral shift should increase extrusion. " \ + f"Baseline: {baseline_extrusion:.2f}, Shifted: {results['lateral_extrusion_mm']:.2f}" + + print(f"\n✓ Lateral meniscus lateral shift test passed!") + print(f" Baseline: {baseline_extrusion:.2f} mm") + print(f" After lateral shift: {results['lateral_extrusion_mm']:.2f} mm") + print(f" Increase: {results['lateral_extrusion_mm'] - baseline_extrusion:.2f} mm") + + +def test_lateral_meniscus_shift_medially_decreases_extrusion(tibia_with_menisci): + """ + Test that shifting lateral meniscus medially decreases extrusion. + + When the lateral meniscus is shifted 5mm medially (toward midline), + the extrusion value should decrease (less positive or more negative). + """ + tibia, med_meniscus, lat_meniscus, baseline_results = tibia_with_menisci + baseline_extrusion = baseline_results['lateral_extrusion_mm'] + + # Shift lateral meniscus medially (5mm in +X direction for right knee) + lat_meniscus_shifted = lat_meniscus.copy() + lat_meniscus_shifted.points = lat_meniscus_shifted.points + np.array([5.0, 0.0, 0.0]) + + # Compute extrusion with shifted meniscus + results = compute_meniscal_extrusion( + tibia_mesh=tibia, + medial_meniscus_mesh=med_meniscus, + lateral_meniscus_mesh=lat_meniscus_shifted, + medial_cart_label=2, + lateral_cart_label=3, + scalar_array_name='labels', + middle_percentile_range=0.1, + ) + + # Verify extrusion decreased + assert results['lateral_extrusion_mm'] < baseline_extrusion, \ + f"Medial shift should decrease extrusion. " \ + f"Baseline: {baseline_extrusion:.2f}, Shifted: {results['lateral_extrusion_mm']:.2f}" + + print(f"\n✓ Lateral meniscus medial shift test passed!") + print(f" Baseline: {baseline_extrusion:.2f} mm") + print(f" After medial shift: {results['lateral_extrusion_mm']:.2f} mm") + print(f" Decrease: {baseline_extrusion - results['lateral_extrusion_mm']:.2f} mm") + + +# ============================================================================ +# Convenience API Tests +# ============================================================================ + +def test_dict_cartilage_labels_replaces_list(): + """ + Test that dict_cartilage_labels can replace list_cartilage_labels. + + Verifies that cartilage thickness and region assignment work with only + dict_cartilage_labels (no list_cartilage_labels needed). + """ + # Get path to test data + test_dir = os.path.dirname(os.path.abspath(__file__)) + data_dir = os.path.join(test_dir, '..', '..', '..', 'data') + path_segmentation = os.path.join(data_dir, 'SAG_3D_DESS_RIGHT_bones_cart_men_fib-label.nrrd') + + if not os.path.exists(path_segmentation): + pytest.skip(f"Test data not found: {path_segmentation}") + + # Create tibia with ONLY dict_cartilage_labels + tibia = BoneMesh( + path_seg_image=path_segmentation, + label_idx=6, + dict_cartilage_labels={'medial': 2, 'lateral': 3} + ) + + # Verify list_cartilage_labels property auto-extracts from dict + assert tibia.list_cartilage_labels == [2, 3] + + # Verify standard operations work + tibia.create_mesh() + tibia.calc_cartilage_thickness() # Should work with dict values + tibia.assign_cartilage_regions() # Should work with dict values + + # Verify thickness and labels were assigned + assert 'thickness (mm)' in tibia.point_data + assert 'labels' in tibia.point_data + + print("\n✓ dict_cartilage_labels successfully replaces list_cartilage_labels!") + + +def test_set_menisci_auto_infers_labels(): + """ + Test that set_menisci() automatically infers labels from dict_cartilage_labels. + + Verifies that cartilage labels don't need to be specified explicitly + when dict_cartilage_labels is set. + """ + test_dir = os.path.dirname(os.path.abspath(__file__)) + data_dir = os.path.join(test_dir, '..', '..', '..', 'data') + path_segmentation = os.path.join(data_dir, 'SAG_3D_DESS_RIGHT_bones_cart_men_fib-label.nrrd') + + if not os.path.exists(path_segmentation): + pytest.skip(f"Test data not found: {path_segmentation}") + + # Setup tibia + tibia = BoneMesh( + path_seg_image=path_segmentation, + label_idx=6, + dict_cartilage_labels={'medial': 2, 'lateral': 3} + ) + tibia.create_mesh() + tibia.calc_cartilage_thickness() + tibia.assign_cartilage_regions() + + # Create menisci + med_meniscus = Mesh(path_seg_image=path_segmentation, label_idx=10) + med_meniscus.create_mesh() + med_meniscus.consistent_faces() + + lat_meniscus = Mesh(path_seg_image=path_segmentation, label_idx=9) + lat_meniscus.create_mesh() + lat_meniscus.consistent_faces() + + # Test: set_menisci WITHOUT explicit labels (should auto-infer from dict) + tibia.set_menisci( + medial_meniscus=med_meniscus, + lateral_meniscus=lat_meniscus + ) + + # Verify labels were cached correctly + assert tibia._meniscal_cart_labels is not None + assert tibia._meniscal_cart_labels['medial'] == 2 + assert tibia._meniscal_cart_labels['lateral'] == 3 + + print("\n✓ set_menisci() successfully auto-infers labels from dict!") + + +def test_meniscal_properties_lazy_evaluation(): + """ + Test that meniscal properties auto-compute on first access (lazy evaluation). + + Verifies that calling properties like med_men_extrusion automatically + triggers computation without explicit compute_meniscal_outcomes() call. + """ + test_dir = os.path.dirname(os.path.abspath(__file__)) + data_dir = os.path.join(test_dir, '..', '..', '..', 'data') + path_segmentation = os.path.join(data_dir, 'SAG_3D_DESS_RIGHT_bones_cart_men_fib-label.nrrd') + + if not os.path.exists(path_segmentation): + pytest.skip(f"Test data not found: {path_segmentation}") + + # Setup + tibia = BoneMesh( + path_seg_image=path_segmentation, + label_idx=6, + dict_cartilage_labels={'medial': 2, 'lateral': 3} + ) + tibia.create_mesh() + tibia.calc_cartilage_thickness() + tibia.assign_cartilage_regions() + + med_meniscus = Mesh(path_seg_image=path_segmentation, label_idx=10) + med_meniscus.create_mesh() + med_meniscus.consistent_faces() + + lat_meniscus = Mesh(path_seg_image=path_segmentation, label_idx=9) + lat_meniscus.create_mesh() + lat_meniscus.consistent_faces() + + tibia.set_menisci(medial_meniscus=med_meniscus, lateral_meniscus=lat_meniscus) + + # Verify outcomes NOT computed yet + assert tibia._meniscal_outcomes is None + + # Access property - should trigger auto-computation + med_extrusion = tibia.med_men_extrusion + + # Verify outcomes NOW computed + assert tibia._meniscal_outcomes is not None + assert isinstance(med_extrusion, (int, float, np.number)) + + # Access another property - should use cached results (no recomputation) + lat_extrusion = tibia.lat_men_extrusion + assert isinstance(lat_extrusion, (int, float, np.number)) + + print("\n✓ Properties successfully auto-compute on first access!") + print(f" Medial extrusion: {med_extrusion:.2f} mm") + print(f" Lateral extrusion: {lat_extrusion:.2f} mm") + + +def test_meniscal_outcomes_caching(): + """ + Test that meniscal outcomes are properly cached and reused. + + Verifies that: + - Results are cached after first computation + - Property values match cached dictionary values + - All expected metrics are present in cache + """ + test_dir = os.path.dirname(os.path.abspath(__file__)) + data_dir = os.path.join(test_dir, '..', '..', '..', 'data') + path_segmentation = os.path.join(data_dir, 'SAG_3D_DESS_RIGHT_bones_cart_men_fib-label.nrrd') + + if not os.path.exists(path_segmentation): + pytest.skip(f"Test data not found: {path_segmentation}") + + # Setup + tibia = BoneMesh( + path_seg_image=path_segmentation, + label_idx=6, + dict_cartilage_labels={'medial': 2, 'lateral': 3} + ) + tibia.create_mesh() + tibia.calc_cartilage_thickness() + tibia.assign_cartilage_regions() + + med_meniscus = Mesh(path_seg_image=path_segmentation, label_idx=10) + med_meniscus.create_mesh() + med_meniscus.consistent_faces() + + lat_meniscus = Mesh(path_seg_image=path_segmentation, label_idx=9) + lat_meniscus.create_mesh() + lat_meniscus.consistent_faces() + + tibia.set_menisci(medial_meniscus=med_meniscus, lateral_meniscus=lat_meniscus) + + # Trigger computation via property access + med_extrusion = tibia.med_men_extrusion + lat_extrusion = tibia.lat_men_extrusion + med_coverage = tibia.med_men_coverage + lat_coverage = tibia.lat_men_coverage + + # Verify all metrics are cached + assert 'medial_extrusion_mm' in tibia._meniscal_outcomes + assert 'lateral_extrusion_mm' in tibia._meniscal_outcomes + assert 'medial_coverage_percent' in tibia._meniscal_outcomes + assert 'lateral_coverage_percent' in tibia._meniscal_outcomes + + # Verify property values match cached values + assert med_extrusion == tibia._meniscal_outcomes['medial_extrusion_mm'] + assert lat_extrusion == tibia._meniscal_outcomes['lateral_extrusion_mm'] + assert med_coverage == tibia._meniscal_outcomes['medial_coverage_percent'] + assert lat_coverage == tibia._meniscal_outcomes['lateral_coverage_percent'] + + print("\n✓ Meniscal outcomes properly cached and accessible!") + + +def test_meniscal_values_reasonable(): + """ + Test that computed meniscal values are reasonable. + + Verifies: + - Extrusion values are numeric types + - Coverage values are percentages (0-100) + """ + test_dir = os.path.dirname(os.path.abspath(__file__)) + data_dir = os.path.join(test_dir, '..', '..', '..', 'data') + path_segmentation = os.path.join(data_dir, 'SAG_3D_DESS_RIGHT_bones_cart_men_fib-label.nrrd') + + if not os.path.exists(path_segmentation): + pytest.skip(f"Test data not found: {path_segmentation}") + + # Setup + tibia = BoneMesh( + path_seg_image=path_segmentation, + label_idx=6, + dict_cartilage_labels={'medial': 2, 'lateral': 3} + ) + tibia.create_mesh() + tibia.calc_cartilage_thickness() + tibia.assign_cartilage_regions() + + med_meniscus = Mesh(path_seg_image=path_segmentation, label_idx=10) + med_meniscus.create_mesh() + med_meniscus.consistent_faces() + + lat_meniscus = Mesh(path_seg_image=path_segmentation, label_idx=9) + lat_meniscus.create_mesh() + lat_meniscus.consistent_faces() + + tibia.set_menisci(medial_meniscus=med_meniscus, lateral_meniscus=lat_meniscus) + + # Get values + med_extrusion = tibia.med_men_extrusion + lat_extrusion = tibia.lat_men_extrusion + med_coverage = tibia.med_men_coverage + lat_coverage = tibia.lat_men_coverage + + # Verify types (accept Python and numpy numeric types) + assert isinstance(med_extrusion, (int, float, np.number)) + assert isinstance(lat_extrusion, (int, float, np.number)) + assert isinstance(med_coverage, (int, float, np.number)) + assert isinstance(lat_coverage, (int, float, np.number)) + + # Verify coverage percentages are in valid range + assert 0 <= med_coverage <= 100, f"Medial coverage {med_coverage}% outside valid range" + assert 0 <= lat_coverage <= 100, f"Lateral coverage {lat_coverage}% outside valid range" + + print("\n✓ All meniscal values are reasonable!") + print(f" Medial: {med_extrusion:.2f} mm extrusion, {med_coverage:.1f}% coverage") + print(f" Lateral: {lat_extrusion:.2f} mm extrusion, {lat_coverage:.1f}% coverage") + + +# ============================================================================ +# TODO: Additional Tests +# ============================================================================ + +def test_extrusion_synthetic_data(): + """TODO: Test extrusion calculation with synthetic tibia and meniscus meshes.""" + pass + + +def test_extrusion_no_extrusion(): + """TODO: Test case where meniscus is fully within cartilage rim.""" + pass + + + + diff --git a/testing/scratch/meniscal_extrusion_Oct.29.2025.ipynb b/testing/scratch/meniscal_extrusion_Oct.29.2025.ipynb new file mode 100644 index 0000000..2ad7a2c --- /dev/null +++ b/testing/scratch/meniscal_extrusion_Oct.29.2025.ipynb @@ -0,0 +1,910 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "2bf43a6a", + "metadata": {}, + "outputs": [], + "source": [ + "import pymskt as mskt\n", + "from itkwidgets import view\n", + "import numpy as np\n", + "import os\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "843e0dda", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Initiating tibia mesh\n", + "Creating tibia mesh\n", + "Calculating cartilage thickness\n", + "WARNING: Mesh is now synonymous with pyvista.PolyData and thus this property is redundant and the Mesh object can be used for anything that pyvista.PolyData or vtk.vtkPolyData can be used for.\n", + "WARNING: Mesh is now synonymous with pyvista.PolyData and thus this property is redundant and the Mesh object can be used for anything that pyvista.PolyData or vtk.vtkPolyData can be used for.\n", + "INTERSECTION IS: 2\n", + "INTERSECTION IS: 2\n", + "Assigning cartilage regions\n", + "INTERSECTION IS: 2\n", + "INTERSECTION IS: 2\n" + ] + } + ], + "source": [ + "path_segmentation = '../../data/SAG_3D_DESS_RIGHT_bones_cart_men_fib-label.nrrd'\n", + "\n", + "print('Initiating tibia mesh')\n", + "tibia = mskt.mesh.BoneMesh(\n", + " path_seg_image=path_segmentation,\n", + " label_idx=6,\n", + " list_cartilage_labels=[2, 3]\n", + ")\n", + "\n", + "print('Creating tibia mesh')\n", + "tibia.create_mesh()\n", + "\n", + "print('Calculating cartilage thickness')\n", + "tibia.calc_cartilage_thickness()\n", + "print('Assigning cartilage regions')\n", + "tibia.assign_cartilage_regions()" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "8f070b73", + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "9690c36889bd41aba731b0080e5f899d", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'name': '_points', 'numberO…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "view(geometries=[tibia])" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "2ea8d441", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "regions_label = 'labels'\n", + "med_tib_cart_label = 2\n", + "lat_tib_cart_label = 3\n", + "\n", + "region_array = tibia[regions_label]\n", + "med_tib_cart_mask = (region_array == med_tib_cart_label)\n", + "lat_tib_cart_mask = (region_array == lat_tib_cart_label)\n", + "\n", + "med_tib_cart_points = tibia.points[med_tib_cart_mask]\n", + "lat_tib_cart_points = tibia.points[lat_tib_cart_mask]\n", + "tib_cart_points = np.concatenate([med_tib_cart_points, lat_tib_cart_points], axis=0)\n", + "\n", + "# do PCA to get the three axes of the tib_cart_points and take the last\n", + "# one as the inf/sup\n", + "X = tib_cart_points - tib_cart_points.mean(axis=0, keepdims=True) # (N,3)\n", + "# PCA via SVD: X = U S Vt, rows of Vt are PCs\n", + "U, S, Vt = np.linalg.svd(X, full_matrices=False)\n", + "pc1, pc2, pc3 = Vt # already orthonormal\n", + "\n", + "is_axis = pc3\n", + "\n", + "# from the PCA we cant know what it up. We should check which side the meniscus is on...\n", + "# or... which side the bone is on... So, from the middle of the cartilage, \n", + "# the opposide of the direction of the middle of the bone is the IS axis. \n", + "mean_tib = np.mean(tibia.points, axis=0)\n", + "mean_cart = np.mean(tib_cart_points, axis=0)\n", + "\n", + "# update is_axis direction based on where mean_tib is relative to mean_cart\n", + "if np.dot(mean_tib - mean_cart, is_axis) > 0:\n", + " is_axis = -is_axis\n", + "\n", + "\n", + "med_tib_center = np.mean(med_tib_cart_points, axis=0)\n", + "lat_tib_center = np.mean(lat_tib_cart_points, axis=0)\n", + "\n", + "ml_axis = lat_tib_center - med_tib_center\n", + "ml_axis = ml_axis / np.linalg.norm(ml_axis)\n", + "\n", + "ap_axis = np.cross(ml_axis, is_axis)\n", + "ap_axis = ap_axis / np.linalg.norm(ap_axis)\n", + "\n", + "dict_tibia_axes = {\n", + " 'ml_axis': ml_axis,\n", + " 'is_axis': is_axis,\n", + " 'ap_axis': ap_axis,\n", + " 'medial_center': med_tib_center,\n", + " 'lateral_center': lat_tib_center,\n", + "}\n", + "\n", + "def get_tibia_axes_meniscal_extrusion(\n", + " tibia_mesh, \n", + " regions_label,\n", + " med_tib_cart_label,\n", + " lat_tib_cart_label,\n", + "):\n", + " region_array = tibia_mesh[regions_label]\n", + " med_tib_cart_mask = (region_array == med_tib_cart_label)\n", + " lat_tib_cart_mask = (region_array == lat_tib_cart_label)\n", + "\n", + " med_tib_cart_points = tibia_mesh.points[med_tib_cart_mask]\n", + " lat_tib_cart_points = tibia_mesh.points[lat_tib_cart_mask]\n", + " tib_cart_points = np.concatenate([med_tib_cart_points, lat_tib_cart_points], axis=0)\n", + "\n", + " # do PCA to get the three axes of the tib_cart_points and take the last\n", + " # one as the inf/sup\n", + " X = tib_cart_points - tib_cart_points.mean(axis=0, keepdims=True) # (N,3)\n", + " # PCA via SVD: X = U S Vt, rows of Vt are PCs\n", + " U, S, Vt = np.linalg.svd(X, full_matrices=False)\n", + " pc1, pc2, pc3 = Vt # already orthonormal\n", + "\n", + " is_axis = pc3\n", + " # from the PCA we cant know what it up. We should check which side the meniscus is on...\n", + " # or... which side the bone is on... So, from the middle of the cartilage, \n", + " # the opposide of the direction of the middle of the bone is the IS axis. \n", + " mean_tib = np.mean(tibia.points, axis=0)\n", + " mean_cart = np.mean(tib_cart_points, axis=0)\n", + "\n", + " # update is_axis direction based on where mean_tib is relative to mean_cart\n", + " if np.dot(mean_tib - mean_cart, is_axis) > 0:\n", + " is_axis = -is_axis\n", + " \n", + " med_tib_center = np.mean(med_tib_cart_points, axis=0)\n", + " lat_tib_center = np.mean(lat_tib_cart_points, axis=0)\n", + "\n", + " ml_axis = lat_tib_center - med_tib_center\n", + " ml_axis = ml_axis / np.linalg.norm(ml_axis)\n", + " \n", + " ap_axis = np.cross(ml_axis, is_axis)\n", + " ap_axis = ap_axis / np.linalg.norm(ap_axis)\n", + "\n", + " dict_tibia_axes = {\n", + " 'ml_axis': ml_axis,\n", + " 'is_axis': is_axis,\n", + " 'ap_axis': ap_axis,\n", + " 'medial_center': med_tib_center,\n", + " 'lateral_center': lat_tib_center,\n", + " }\n", + " \n", + " return dict_tibia_axes\n" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "98c2f353", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{'ml_axis': pyvista_ndarray([-0.9919523 , 0.12175144, 0.0347462 ], dtype=float32), 'is_axis': array([0.03815307, 0.04018929, 0.9984634 ], dtype=float32), 'ap_axis': array([ 0.12016812, 0.99175525, -0.04451112], dtype=float32), 'medial_center': pyvista_ndarray([-57.089283, -5.439862, -9.311588], dtype=float32), 'lateral_center': pyvista_ndarray([-92.64511 , -1.0757675, -8.066135 ], dtype=float32)}\n", + "{'ml_axis': pyvista_ndarray([-0.9919523 , 0.12175144, 0.0347462 ], dtype=float32), 'is_axis': array([0.03815307, 0.04018929, 0.9984634 ], dtype=float32), 'ap_axis': array([ 0.12016812, 0.99175525, -0.04451112], dtype=float32), 'medial_center': pyvista_ndarray([-57.089283, -5.439862, -9.311588], dtype=float32), 'lateral_center': pyvista_ndarray([-92.64511 , -1.0757675, -8.066135 ], dtype=float32)}\n" + ] + } + ], + "source": [ + "dict_tibia_axes_func = get_tibia_axes_meniscal_extrusion(\n", + " tibia_mesh=tibia, \n", + " regions_label='labels',\n", + " med_tib_cart_label=2,\n", + " lat_tib_cart_label=3,\n", + ")\n", + "\n", + "print(dict_tibia_axes)\n", + "print(dict_tibia_axes_func)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "1164c084", + "metadata": {}, + "outputs": [], + "source": [ + "# now label tibia points as having meniscus above them along the si_axis\n", + "\n", + "med_meniscus = mskt.mesh.Mesh(\n", + " path_seg_image=path_segmentation,\n", + " label_idx=10,\n", + ")\n", + "med_meniscus.create_mesh()\n", + "med_meniscus.consistent_faces()\n", + "\n", + "lat_meniscus = mskt.mesh.Mesh(\n", + " path_seg_image=path_segmentation,\n", + " label_idx=9,\n", + ")\n", + "lat_meniscus.create_mesh()\n", + "lat_meniscus.consistent_faces()\n" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "4fcee9e4", + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "22b292b9804f4707ad5fabbefcc672a1", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'name': '_points', 'numberO…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "view(geometries=[tibia, med_meniscus, lat_meniscus])" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "id": "c460b7b3", + "metadata": {}, + "outputs": [], + "source": [ + "tibia.calc_distance_to_other_mesh(\n", + " list_other_meshes=[med_meniscus],\n", + " ray_cast_length=10, \n", + " name='med_men_dist_mm',\n", + " direction=dict_tibia_axes['is_axis'],\n", + ")\n", + "\n", + "tibia.calc_distance_to_other_mesh(\n", + " list_other_meshes=[lat_meniscus],\n", + " ray_cast_length=10, \n", + " name='lat_men_dist_mm',\n", + " direction=dict_tibia_axes['is_axis'],\n", + ")\n", + "\n", + "binary_mask_med_men_above = tibia['med_men_dist_mm'] > 0\n", + "binary_mask_lat_men_above = tibia['lat_men_dist_mm'] > 0\n", + "\n", + "binary_mask_med_cart = tibia['labels'] == med_tib_cart_label\n", + "binary_mask_lat_cart = tibia['labels'] == lat_tib_cart_label\n", + "\n", + "tibia['med_men_above'] = binary_mask_med_men_above.astype(float)\n", + "tibia['lat_men_above'] = binary_mask_lat_men_above.astype(float)\n", + " \n" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "id": "d76982df", + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "06c9a515106e40df9f0b2ef1a970ff75", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'name': '_points', 'numberO…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "view(geometries=[tibia])" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "id": "65bdf405", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Med Cart Area: 868.46 mm^2\n", + "Med Cart Men Area: 346.47 mm^2\n", + "Percent Med Men Coverage: 39.89%\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "c6dcfa1970814b7aae31496352b9c6a8", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'name': '_points', 'numberO…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "tibia_med_cart = tibia.copy()\n", + "# delete points that are not medial cartilage\n", + "tibia_med_cart.remove_points(~binary_mask_med_cart, inplace=True)\n", + "tibia_med_cart.clean(inplace=True)\n", + "\n", + "area_med_cart = tibia_med_cart.area\n", + "\n", + "tibia_med_cart_men = tibia_med_cart.copy()\n", + "tibia_med_cart_men.remove_points(tibia_med_cart_men['med_men_above'] == 0, inplace=True)\n", + "tibia_med_cart_men.clean(inplace=True)\n", + "\n", + "area_med_cart_men = tibia_med_cart_men.area\n", + "\n", + "percent_med_men_coverage = (area_med_cart_men / area_med_cart) * 100\n", + "\n", + "print(f'Med Cart Area: {area_med_cart:.2f} mm^2')\n", + "print(f'Med Cart Men Area: {area_med_cart_men:.2f} mm^2')\n", + "print(f'Percent Med Men Coverage: {percent_med_men_coverage:.2f}%')\n", + "\n", + "view(geometries=[tibia, tibia_med_cart, tibia_med_cart_men])" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "id": "0b5987ed", + "metadata": {}, + "outputs": [], + "source": [ + "def get_meniscal_coverage(\n", + " tibia_mesh,\n", + " meniscal_mesh,\n", + " tibia_regions_label,\n", + " tibia_cart_label,\n", + " is_direction,\n", + " meniscus_side,\n", + " ray_cast_length=20,\n", + "):\n", + " tibia_mesh.calc_distance_to_other_mesh(\n", + " list_other_meshes=[meniscal_mesh],\n", + " ray_cast_length=ray_cast_length, \n", + " name=f'{meniscus_side}_men_dist_mm',\n", + " direction=is_direction,\n", + " )\n", + " binary_mask_men_above = tibia_mesh[f'{meniscus_side}_men_dist_mm'] > 0\n", + " binary_mask_cart = tibia_mesh[tibia_regions_label] == tibia_cart_label\n", + " tibia_mesh[f'{meniscus_side}_men_above'] = binary_mask_men_above.astype(float)\n", + " tibia_mesh[f'{meniscus_side}_cart'] = binary_mask_cart.astype(float)\n", + "\n", + " tibia_cart = tibia_mesh.copy()\n", + " # delete points that are not medial cartilage\n", + " tibia_cart.remove_points(~binary_mask_cart, inplace=True)\n", + " tibia_cart.clean(inplace=True)\n", + "\n", + " area_cart = tibia_cart.area\n", + "\n", + " tibia_cart_men = tibia_cart.copy()\n", + " tibia_cart_men.remove_points(tibia_cart_men[f'{meniscus_side}_men_above'] == 0, inplace=True)\n", + " tibia_cart_men.clean(inplace=True)\n", + "\n", + " area_cart_men = tibia_cart_men.area\n", + "\n", + " percent_cart_men_coverage = (area_cart_men / area_cart) * 100\n", + " \n", + " dict_meniscal_coverage = {\n", + " f'{meniscus_side}_cart_men_coverage': percent_cart_men_coverage,\n", + " f'{meniscus_side}_cart_men_area': area_cart_men,\n", + " f'{meniscus_side}_cart_area': area_cart,\n", + " }\n", + " \n", + " return dict_meniscal_coverage\n", + "\n", + "dict_med_men_coverage = get_meniscal_coverage(\n", + " tibia_mesh=tibia,\n", + " meniscal_mesh=med_meniscus,\n", + " tibia_regions_label='labels',\n", + " tibia_cart_label=2,\n", + " meniscus_side='med',\n", + " is_direction=dict_tibia_axes['is_axis'],\n", + ")\n", + "\n", + "dict_lat_men_coverage = get_meniscal_coverage(\n", + " tibia_mesh=tibia,\n", + " meniscal_mesh=lat_meniscus,\n", + " tibia_regions_label='labels',\n", + " tibia_cart_label=3,\n", + " meniscus_side='lat',\n", + " is_direction=dict_tibia_axes['is_axis'],\n", + ")\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "id": "e350fef2", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{'med_cart_men_coverage': 39.89438019552111, 'med_cart_men_area': 346.4668587995976, 'med_cart_area': 868.460312208322}\n", + "{'lat_cart_men_coverage': 59.36730303771416, 'lat_cart_men_area': 438.1155627617297, 'lat_cart_area': 737.9745084317016}\n" + ] + } + ], + "source": [ + "print(dict_med_men_coverage)\n", + "print(dict_lat_men_coverage)" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "id": "d7afcb99", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "'\\nNow I want to work on extrusion. \\n- We have the ML axis from the tibia.\\n- We have the medial/lateral meniscus segmentation on the bone to define the edge that\\nextrusion will be compared against. \\n- We will use the meniscus itself to define the extrusion depth.\\n\\n- We are just going to bin X number of bins in the AP direction (disregarding what is front/back)\\n- We will compute extrusion +/- for each bin\\n- We will then return the: mean, median, min, max, std of extrusion. \\n'" + ] + }, + "execution_count": 28, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "\"\"\"\n", + "Now I want to work on extrusion. \n", + "- We have the ML axis from the tibia.\n", + "- We have the medial/lateral meniscus segmentation on the bone to define the edge that\n", + "extrusion will be compared against. \n", + "- We will use the meniscus itself to define the extrusion depth.\n", + "\n", + "- We are just going to bin X number of bins in the AP direction (disregarding what is front/back)\n", + "- We will compute extrusion +/- for each bin\n", + "- We will then return the: mean, median, min, max, std of extrusion. \n", + "\"\"\"" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "id": "0b2c988d", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Extrusion: -3.97 mm\n" + ] + } + ], + "source": [ + "# project all med cart points onto the ML axis\n", + "tib_cart_label = 'labels'\n", + "med_tib_cart_label = 2\n", + "lat_tib_cart_label = 3\n", + "\n", + "med_cart_points = tibia[tib_cart_label] == med_tib_cart_label\n", + "\n", + "med_cart_points = tibia.points[med_cart_points]\n", + "\n", + "# project med_cart_points onto the ML axis\n", + "ml_axis = dict_tibia_axes['ml_axis']\n", + "med_cart_points_ml = np.dot(med_cart_points, ml_axis)\n", + "\n", + "# get medial meniscus points & project onto the ML axis\n", + "med_men_points = med_meniscus.points\n", + "med_men_points_ml = np.dot(med_men_points, ml_axis)\n", + "\n", + "# get the min of the medial cart and the min of the medial men\n", + "min_med_cart_ml = np.min(med_cart_points_ml)\n", + "min_med_men_ml = np.min(med_men_points_ml)\n", + "\n", + "extrusion = min_med_men_ml - min_med_cart_ml\n", + "\n", + "print(f'Extrusion: {extrusion:.2f} mm')\n" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "id": "15637c08", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGdCAYAAAAMm0nCAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAJSJJREFUeJzt3Ql0VOX9//FvICQBJIEQstUQQGUnbGrMURCEJgIHtVK3sFkRBAMWgoixiCEuwdCDYEuxtCB6DJX6P4gKll9C2CxEllgMBI0EwUDJ0iokLCUh5P7P85wz0xkI+6Thmft+nXOZufe5M3PnnpuZD882PpZlWQIAAGCQRg19AAAAAFeLAAMAAIxDgAEAAMYhwAAAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMI6veKna2lo5evSotGjRQnx8fBr6cAAAwBVQ8+ueOHFCIiMjpVGjRvYLMCq8REVFNfRhAACAa3D48GG5+eab7RdgVM2L4wQEBgY29OEAAIArUFlZqSsgHN/jtgswjmYjFV4IMAAAmOVy3T/oxAsAAIxDgAEAAMYhwAAAAOMQYAAAgPcHmC1btsjw4cP1+GzVwWb16tVu5WpbXcu8efOc+7Rr1+6C8rlz57o9T35+vvTr108CAgJ0b+SMjIzreZ8AAMDOAebUqVPSs2dPWbRoUZ3lJSUlbsuyZct0QBkxYoTbfmlpaW77TZkyxW0IVXx8vERHR0teXp4OP6mpqbJkyZJreY8AAMDLXPUw6iFDhujlYsLDw93WP/nkExk4cKB06NDBbbsa333+vg6ZmZlSXV2tw4+fn59069ZNdu/eLfPnz5cJEyZc7SEDAAAvU699YMrKymTt2rUybty4C8pUk1Hr1q2ld+/euoalpqbGWZabmyv9+/fX4cUhISFBCgsL5dixY3W+VlVVla65cV0AAIB3qteJ7N577z1d0/Lwww+7bX/uueekT58+EhwcLNu2bZOUlBTdjKRqWJTS0lJp376922PCwsKcZa1atbrgtdLT02XOnDn1+XYAAIAdAoxqAho5cqTuiOsqOTnZeT8mJkbXtDzzzDM6hPj7+1/Ta6kQ5Pq8jqmIAQCA96m3APPFF1/oJp+VK1dedt/Y2FjdhHTo0CHp1KmT7hujmp9cOdYv1m9GBZ9rDT8AAMAs9dYHZunSpdK3b189YulyVAdd9ZPZoaGhej0uLk4P1z579qxzn+zsbB1u6mo+AgAA9nLVAebkyZM6cKhFOXjwoL5fXFzs1nzz0UcfydNPP33B41UH3QULFsjXX38t33//vR5xNG3aNBk1apQznCQmJupmJdX5t6CgQNfiLFy40K2JCAAA2NdVNyHt2rVLD4t2cISKsWPHyvLly/X9Dz/8UCzLkieeeOKCx6tmHlWu5nVRI4dUZ10VYFzDSVBQkGRlZUlSUpKuxQkJCZHZs2czhBoAAGg+lkoaXkjVAqkgVFFRIYGBgWJrG9PFGANTGvoIAAAGfH/zW0gAAMA4BBgAAGAcAgwAADAOAQYAABiHAAMAAIxDgAEAAMYhwAAAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjEOAAQAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwjm9DHwDgZmO6WSdkYEpDHwEA2BI1MAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjEOAAQAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAAHh/gNmyZYsMHz5cIiMjxcfHR1avXu1W/uSTT+rtrsv999/vts9PP/0kI0eOlMDAQGnZsqWMGzdOTp486bZPfn6+9OvXTwICAiQqKkoyMjKu9T0CAAC7B5hTp05Jz549ZdGiRRfdRwWWkpIS5/KXv/zFrVyFl4KCAsnOzpY1a9boUDRhwgRneWVlpcTHx0t0dLTk5eXJvHnzJDU1VZYsWXK1hwsAALyQ79U+YMiQIXq5FH9/fwkPD6+z7JtvvpF169bJzp075fbbb9fbfve738nQoUPlt7/9ra7ZyczMlOrqalm2bJn4+flJt27dZPfu3TJ//ny3oAMAAOypXvrAbNq0SUJDQ6VTp04yadIk+fHHH51lubm5utnIEV6UwYMHS6NGjWT79u3Offr376/Di0NCQoIUFhbKsWPH6nzNqqoqXXPjugAAAO/k8QCjmo/ef/99ycnJkTfffFM2b96sa2zOnTuny0tLS3W4ceXr6yvBwcG6zLFPWFiY2z6Odcc+50tPT5egoCDnovrNAAAA73TVTUiX8/jjjzvv9+jRQ2JiYuSWW27RtTKDBg2S+pKSkiLJycnOdVUDQ4gBAMA7eTzAnK9Dhw4SEhIiRUVFOsCovjHl5eVu+9TU1OiRSY5+M+q2rKzMbR/H+sX61qh+N2r5n9iY/r95HQAA0DDzwBw5ckT3gYmIiNDrcXFxcvz4cT26yGHDhg1SW1srsbGxzn3UyKSzZ88691EjllSfmlatWtX3IQMAAG8LMGq+FjUiSC3KwYMH9f3i4mJdNmPGDPnyyy/l0KFDuh/Mgw8+KLfeeqvuhKt06dJF95MZP3687NixQ7Zu3SqTJ0/WTU9qBJKSmJioO/Cq+WHUcOuVK1fKwoUL3ZqIAACAfV11gNm1a5f07t1bL4oKFer+7NmzpXHjxnoCugceeEA6duyoA0jfvn3liy++cGveUcOkO3furJuU1PDpe+65x22OF9UJNysrS4cj9fjp06fr52cINQAAUHwsy7K88VSoTrwqCFVUVOgZfz2KPjBwGJjCuQCABvj+5reQAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAABiHAAMAAIxDgAEAAMYhwAAAAOP4NvQBAEbbmC7GGJjS0EcAAB5DDQwAADAOAQYAABiHAAMAAIxDgAEAAMYhwAAAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjEOAAQAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMD7A8yWLVtk+PDhEhkZKT4+PrJ69Wpn2dmzZ2XmzJnSo0cPad68ud5nzJgxcvToUbfnaNeunX6s6zJ37ly3ffLz86Vfv34SEBAgUVFRkpGRcT3vEwAA2DnAnDp1Snr27CmLFi26oOz06dPy1Vdfycsvv6xvV61aJYWFhfLAAw9csG9aWpqUlJQ4lylTpjjLKisrJT4+XqKjoyUvL0/mzZsnqampsmTJkmt5jwAAwMv4Xu0DhgwZope6BAUFSXZ2ttu23//+93LnnXdKcXGxtG3b1rm9RYsWEh4eXufzZGZmSnV1tSxbtkz8/PykW7dusnv3bpk/f75MmDDhag8ZAAB4mXrvA1NRUaGbiFq2bOm2XTUZtW7dWnr37q1rWGpqapxlubm50r9/fx1eHBISEnRtzrFjx+p8naqqKl1z47oAAADvdNU1MFfjzJkzuk/ME088IYGBgc7tzz33nPTp00eCg4Nl27ZtkpKSopuRVA2LUlpaKu3bt3d7rrCwMGdZq1atLnit9PR0mTNnTn2+HQAA4O0BRnXoffTRR8WyLFm8eLFbWXJysvN+TEyMrml55plndAjx9/e/ptdTIcj1eVUNjOr8CwAAvI9vfYaXH374QTZs2OBW+1KX2NhY3YR06NAh6dSpk+4bU1ZW5raPY/1i/WZU8LnW8AMAAGzeB8YRXvbv3y/r16/X/VwuR3XQbdSokYSGhur1uLg4PVxbPZeD6hyswk1dzUcAAMBerroG5uTJk1JUVORcP3jwoA4gqj9LRESE/PKXv9RDqNesWSPnzp3TfVYUVa6ailQH3e3bt8vAgQP1SCS1Pm3aNBk1apQznCQmJur+LOPGjdN9aPbu3SsLFy6Ut956y5PvHQAAGMrHUp1UrsKmTZt0+Djf2LFj9Vwt53e+ddi4caMMGDBAh5tnn31Wvv32Wz1ySO0/evRo3X/FtQlITWSXlJQkO3fulJCQED1PjAozV0r1gVHDutUoqMs1YV21jemefT7gf2FgCucZwA3vSr+/rzrAmIIAA5yHAAPAi76/+S0kAABgHAIMAAAwDgEGAAAYp15n4gVwAzGp8zn9dQBcBjUwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAABiHAAMAAIxDgAEAAMYhwAAAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjEOAAQAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAAHh/gNmyZYsMHz5cIiMjxcfHR1avXu1WblmWzJ49WyIiIqRp06YyePBg2b9/v9s+P/30k4wcOVICAwOlZcuWMm7cODl58qTbPvn5+dKvXz8JCAiQqKgoycjIuNb3CAAA7B5gTp06JT179pRFixbVWa6Cxttvvy3vvPOObN++XZo3by4JCQly5swZ5z4qvBQUFEh2drasWbNGh6IJEyY4yysrKyU+Pl6io6MlLy9P5s2bJ6mpqbJkyZJrfZ8AAMCL+FiqyuRaH+zjIx9//LE89NBDel09laqZmT59ujz//PN6W0VFhYSFhcny5cvl8ccfl2+++Ua6du0qO3fulNtvv13vs27dOhk6dKgcOXJEP37x4sXym9/8RkpLS8XPz0/v8+KLL+ranm+//faKjk2FoKCgIP36qqbHozame/b5AJhtYEpDHwHgNa70+9ujfWAOHjyoQ4dqNnJQBxEbGyu5ubl6Xd2qZiNHeFHU/o0aNdI1No59+vfv7wwviqrFKSwslGPHjtX52lVVVfpNuy4AAMA7eTTAqPCiqBoXV2rdUaZuQ0ND3cp9fX0lODjYbZ+6nsP1Nc6Xnp6uw5JjUf1mAACAd/KaUUgpKSm6usmxHD58uKEPCQAAmBBgwsPD9W1ZWZnbdrXuKFO35eXlbuU1NTV6ZJLrPnU9h+trnM/f31+3lbkuAADAO3k0wLRv314HjJycHOc21RdF9W2Ji4vT6+r2+PHjenSRw4YNG6S2tlb3lXHso0YmnT171rmPGrHUqVMnadWqlScPGQAA2CHAqPladu/erRdHx111v7i4WI9Kmjp1qrz22mvy6aefyp49e2TMmDF6ZJFjpFKXLl3k/vvvl/Hjx8uOHTtk69atMnnyZD1CSe2nJCYm6g68an4YNdx65cqVsnDhQklOTvb0+wcAAAbyvdoH7Nq1SwYOHOhcd4SKsWPH6qHSL7zwgp4rRs3rompa7rnnHj1MWk1I55CZmalDy6BBg/TooxEjRui5YxxUJ9ysrCxJSkqSvn37SkhIiJ4cz3WuGAAAYF/XNQ/MjYx5YAD8zzAPDGD2PDAAAAD/CwQYAABgHAIMAADw/k68AACDmfZbbvQvwkVQAwMAAIxDgAEAAMYhwAAAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjMNMvABgt9ltAS9ADQwAADAOAQYAABiHAAMAAIxDgAEAAMYhwAAAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjEOAAQAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwjm9DHwAAABe1Md2ckzMwpaGPwFaogQEAAMYhwAAAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjOPxANOuXTvx8fG5YElKStLlAwYMuKBs4sSJbs9RXFwsw4YNk2bNmkloaKjMmDFDampqPH2oAADAUB6fiXfnzp1y7tw55/revXvl5z//uTzyyCPObePHj5e0tDTnugoqDuqxKryEh4fLtm3bpKSkRMaMGSNNmjSRN954w9OHCwAADOTxANOmTRu39blz58ott9wi9957r1tgUQGlLllZWbJv3z5Zv369hIWFSa9eveTVV1+VmTNnSmpqqvj5+Xn6kAEAgGHqtQ9MdXW1fPDBB/LUU0/ppiKHzMxMCQkJke7du0tKSoqcPn3aWZabmys9evTQ4cUhISFBKisrpaCgoD4PFwAAGKJef8xx9erVcvz4cXnyySed2xITEyU6OloiIyMlPz9f16wUFhbKqlWrdHlpaalbeFEc66rsYqqqqvTioAIPAADwTvUaYJYuXSpDhgzRYcVhwoQJzvuqpiUiIkIGDRokBw4c0E1N1yo9PV3mzJlz3ccMALCPBTnfee65/m+teNqhucM8/pzeot6akH744Qfdj+Xpp5++5H6xsbH6tqioSN+qvjFlZWVu+zjWL9ZvRlFNURUVFc7l8OHDHngXAADAVgHm3Xff1UOg1YiiS9m9e7e+VTUxSlxcnOzZs0fKy8ud+2RnZ0tgYKB07dr1os/j7++v93FdAACAd6qXJqTa2lodYMaOHSu+vv99CdVMtGLFChk6dKi0bt1a94GZNm2a9O/fX2JiYvQ+8fHxOqiMHj1aMjIydL+XWbNm6XlkVEgBAAColwCjmo7UZHRq9JErNQRalS1YsEBOnTolUVFRMmLECB1QHBo3bixr1qyRSZMm6dqY5s2b6yDkOm8MAACwt3oJMKoWxbKsC7arwLJ58+bLPl6NUvr888/r49AAAIAX4LeQAACAcQgwAADAOPU6DwwAAJ6cawVwoAYGAAAYhwADAACMQxMSAAAeMNX3/3n+PG7Ml3ozMEVMRg0MAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAABiHAAMAAIxDgAEAAMYhwAAAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjOPb0AcAALh+C3K+4zTCVqiBAQAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgHEYhAbhh3MgjaaYO6tjQhwDABTUwAADAOAQYAABgHAIMAAAwDn1gABu5kfuYAMDVoAYGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxPB5gUlNTxcfHx23p3Lmzs/zMmTOSlJQkrVu3lptuuklGjBghZWVlbs9RXFwsw4YNk2bNmkloaKjMmDFDampqPH2oAADAUPUyjLpbt26yfv36/76I739fZtq0abJ27Vr56KOPJCgoSCZPniwPP/ywbN26VZefO3dOh5fw8HDZtm2blJSUyJgxY6RJkybyxhtv1MfhAgAAw9RLgFGBRQWQ81VUVMjSpUtlxYoVct999+lt7777rnTp0kW+/PJLueuuuyQrK0v27dunA1BYWJj06tVLXn31VZk5c6au3fHz86uPQwYAAHYPMPv375fIyEgJCAiQuLg4SU9Pl7Zt20peXp6cPXtWBg8e7NxXNS+pstzcXB1g1G2PHj10eHFISEiQSZMmSUFBgfTu3bvO16yqqtKLQ2VlZX28NQA2xSSAgJf3gYmNjZXly5fLunXrZPHixXLw4EHp16+fnDhxQkpLS3UNSsuWLd0eo8KKKlPUrWt4cZQ7yi5GhSTVJOVYoqKiPP3WAACAt9bADBkyxHk/JiZGB5ro6Gj561//Kk2bNpX6kpKSIsnJyW41MIQYAAC8U70Po1a1LR07dpSioiLdL6a6ulqOHz/uto8aheToM6Nuzx+V5Fivq1+Ng7+/vwQGBrotAADAO9V7gDl58qQcOHBAIiIipG/fvno0UU5OjrO8sLBQD5tWfWUUdbtnzx4pLy937pOdna0DSdeuXev7cAEAgB2bkJ5//nkZPny4bjY6evSovPLKK9K4cWN54okndN+UcePG6aae4OBgHUqmTJmiQ4vqwKvEx8froDJ69GjJyMjQ/V5mzZql545RtSzAjY7OngBM+DxZ8H9rr+vxh+YOE68KMEeOHNFh5ccff5Q2bdrIPffco4dIq/vKW2+9JY0aNdIT2KlRQ2qE0R/+8Afn41XYWbNmjR51pIJN8+bNZezYsZKWlubpQwUAAIbysSzLEi+kOvGqGh8194zH+8NsTPfs88GrUAMDwAQLan55Q9bAXOn3N7+FBAAAjEOAAQAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAABiHAAMAAIxDgAEAAMYhwAAAAOMQYAAAgHEIMAAAwDi+DX0AwNVakPMdJw0AbI4aGAAAYBwCDAAAMA4BBgAAGIcAAwAAjEOAAQAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADCOxwNMenq63HHHHdKiRQsJDQ2Vhx56SAoLC932GTBggPj4+LgtEydOdNunuLhYhg0bJs2aNdPPM2PGDKmpqfH04QIAAAP5evoJN2/eLElJSTrEqMDx0ksvSXx8vOzbt0+aN2/u3G/8+PGSlpbmXFdBxeHcuXM6vISHh8u2bdukpKRExowZI02aNJE33njD04cMAADsHmDWrVvntr58+XJdg5KXlyf9+/d3CywqoNQlKytLB57169dLWFiY9OrVS1599VWZOXOmpKamip+fn6cPGwAAGKTe+8BUVFTo2+DgYLftmZmZEhISIt27d5eUlBQ5ffq0syw3N1d69Oihw4tDQkKCVFZWSkFBQZ2vU1VVpctdFwAA4J08XgPjqra2VqZOnSp33323DioOiYmJEh0dLZGRkZKfn69rVlQ/mVWrVuny0tJSt/CiONZV2cX63syZM6c+3w4AALBDgFF9Yfbu3St///vf3bZPmDDBeV/VtERERMigQYPkwIEDcsstt1zTa6lanOTkZOe6qoGJioq6jqMHAAC2a0KaPHmyrFmzRjZu3Cg333zzJfeNjY3Vt0VFRfpW9Y0pKytz28exfrF+M/7+/hIYGOi2AAAA7+TxAGNZlg4vH3/8sWzYsEHat29/2cfs3r1b36qaGCUuLk727Nkj5eXlzn2ys7N1KOnataunDxkAANi9CUk1G61YsUI++eQTPReMo89KUFCQNG3aVDcTqfKhQ4dK69atdR+YadOm6RFKMTExel817FoFldGjR0tGRoZ+jlmzZunnVjUtAADA3jxeA7N48WI98khNVqdqVBzLypUrdbkaAq2GR6uQ0rlzZ5k+fbqMGDFCPvvsM+dzNG7cWDc/qVtVGzNq1Cg9D4zrvDEAAMC+fOujCelSVMdaNdnd5ahRSp9//rkHjwwAAHgLfgsJAAAYhwADAACMU6/zwMBcC3K+a+hDAADgoqiBAQAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAABiHAAMAAIxDgAEAAMYhwAAAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjEOAAQAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGCcGzrALFq0SNq1aycBAQESGxsrO3bsaOhDAgAAN4AbNsCsXLlSkpOT5ZVXXpGvvvpKevbsKQkJCVJeXt7QhwYAABrYDRtg5s+fL+PHj5df/epX0rVrV3nnnXekWbNmsmzZsoY+NAAA0MB85QZUXV0teXl5kpKS4tzWqFEjGTx4sOTm5tb5mKqqKr04VFRU6NvKykrPH+CpM+LtzlRVN/QhAADqUW3N6et6fL18v7o8r2VZ5gWYf//733Lu3DkJCwtz267Wv/322zofk56eLnPmzLlge1RUVL0dJwAA5lpxXY8OWiD16sSJExIUFGRWgLkWqrZG9ZlxqK2tlZ9++klat24tPj4+F6Q7FWwOHz4sgYGBDXC0Ny7ODeeG64a/KT5vGp6dP4sty9LhJTIy8pL73ZABJiQkRBo3bixlZWVu29V6eHh4nY/x9/fXi6uWLVte8nXURWG3C+NKcW44N1w3/E3xedPw7PpZHHSJmpcbuhOvn5+f9O3bV3JyctxqVNR6XFxcgx4bAABoeDdkDYyimoPGjh0rt99+u9x5552yYMECOXXqlB6VBAAA7O2GDTCPPfaY/Otf/5LZs2dLaWmp9OrVS9atW3dBx95roZqa1Pwy5zc5gXPDdcPflKfxecO54bqpHz7W5cYpAQAA3GBuyD4wAAAAl0KAAQAAxiHAAAAA4xBgAACAcbw2wCxevFhiYmKckwCp+WP+9re/OcsHDBigZ+h1XSZOnCh2NHfuXP3+p06d6tx25swZSUpK0jMZ33TTTTJixIgLJha067mx67WTmpp6wfvu3Lmzs9zu18zlzo9drxuHf/7znzJq1Ch9fTRt2lR69Oghu3btcpar8SRq1GlERIQuV799t3//frGDy52bJ5988oJr5/777xe7u2GHUV+vm2++WX/53HbbbfoP47333pMHH3xQ/vGPf0i3bt30PurXrtPS0pyPUb92bTc7d+6UP/7xjzrsuZo2bZqsXbtWPvroIz0j4uTJk+Xhhx+WrVu3it3PjZ2vHfW3s379eue6r+9/P0K4Zi59fux83Rw7dkzuvvtuGThwoP6PZJs2bXQ4adWqlXOfjIwMefvtt/Vndfv27eXll1+WhIQE2bdvnwQEBIidz42iAsu7777rXPdnGhDvDTDDhw93W3/99dd1rcyXX37pDDDqw+NiP01gBydPnpSRI0fKn/70J3nttdfcfsl76dKlsmLFCrnvvvv0NvWH06VLF33+7rrrLrHruXGw67WjvpDret9cM5c+P3a/bt588039uz6uX8AqpDio/2SqyUpnzZql/6OpvP/++3rer9WrV8vjjz8udj03roHFjteOLZuQXKlftv7www/1TL6uP0WQmZmpf3epe/fu+scgT5++vp8WN42q7h82bJiuqnWVl5cnZ8+edduuqsLbtm0rubm5YudzY/drR/3PUP3AWocOHXTAKy4u1tu5Zi59fux+3Xz66ad6VvVHHnlEQkNDpXfv3vo/Bw4HDx7UE5a6/r2pmt/Y2Fiv/8y53Llx2LRpky7v1KmTTJo0SX788UexO6+tgVH27NmjA4tqm1dt8h9//LF07dpVlyUmJkp0dLT+sMnPz5eZM2dKYWGhrFq1SuxABbqvvvpKN5OcT32QqN+jOv/HMNX/hlSZnc+Nna8d9WWyfPly/QFaUlIic+bMkX79+snevXttf81c7vy0aNHCtteN8v333+sacPUTMS+99JL+23ruuef054z6yRjH58r5M63b4TPncufG0XykmvBVzcyBAwf0fkOGDNHhTv3wsW1ZXqyqqsrav3+/tWvXLuvFF1+0QkJCrIKCgjr3zcnJUTMSW0VFRZa3Ky4utkJDQ62vv/7aue3ee++1fv3rX+v7mZmZlp+f3wWPu+OOO6wXXnjBsvO5sfu14+rYsWNWYGCg9ec//9nW18yVnB+7XzdNmjSx4uLi3LZNmTLFuuuuu/T9rVu36nNx9OhRt30eeeQR69FHH7XsfG7qcuDAAX2+1q9fb9mZVzchqQR766236l+2Tk9Pl549e8rChQsv+r8npaioSLydqu4vLy+XPn366DZ7tWzevFl3oFP31f96qqur5fjx426PUyNKvL0N9nLnRjVH2vnacaVq6Dp27Kjft7ou7HrNXMn5qYudrhs1sshR++2g+tQ5mtgc18j5o9bscP1c7tzURTVRhoSE2OLauRSvDjDnq62tlaqqqjrLdu/e7byYvN2gQYN085p6z45FtcGqNnvH/SZNmkhOTo7zMaqqW/1BufYhsuO5qau61k7XzvkdnVV1tnrf6j8Jdr1mruT81MVO140aZaOuB1ffffedblJTVNOICiqu109lZaVs377d66+fy52buhw5ckT3gbHDtXNJlpdSTUabN2+2Dh48aOXn5+t1Hx8fKysrS1fZpqWl6aYlVf7JJ59YHTp0sPr372/Z1fnNJBMnTrTatm1rbdiwQZ8nVcV5fjWnHc+Nna+d6dOnW5s2bdLvW1X5Dx48WDfLlpeX63K7XzOXOj92vm6UHTt2WL6+vtbrr7+um/VVk2OzZs2sDz74wLnP3LlzrZYtW+pzoz6zH3zwQat9+/bWf/7zH8vO5+bEiRPW888/b+Xm5uprRzUb9enTx7rtttusM2fOWHbmtQHmqaeesqKjo3W7fJs2baxBgwbp8OLo56A+OIKDgy1/f3/r1ltvtWbMmGFVVFRYdnV+gFEfGs8++6zVqlUr/cf0i1/8wiopKbHsfm7sfO089thjVkREhP6b+tnPfqbXXftv2P2audT5sfN14/DZZ59Z3bt31++/c+fO1pIlS9zKa2trrZdfftkKCwvT+6jP7MLCQsvu5+b06dNWfHy8/h5T/WXU99r48eOt0tJSy+581D8NXQsEAABwNWzVBwYAAHgHAgwAADAOAQYAABiHAAMAAIxDgAEAAMYhwAAAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAxDT/HxHr1KyoGwHtAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.hist(med_cart_points_ml)\n", + "plt.hist(med_men_points_ml, alpha=0.5)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "id": "8cbede8b", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 31, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAOpNJREFUeJzt3Ql4VdXZ9vE7cyAkISGBhCSgyBAgzCCDCggoRbRgrVCqxVq1rdUW269+LR3etmrF1td+tbVF6VulbylFUJGqoCIoiIAyKvMsBMjEkBkyf9daGUyUIfM+5+z/77r2lZOTnZzNOcC5s9aznuVXUVFRIQAAAIf4O/XAAAAABmEEAAA4ijACAAAcRRgBAACOIowAAABHEUYAAICjCCMAAMBRhBEAAOCoQHmB8vJynTx5UuHh4fLz83P6cgAAQD2Yvqp5eXnq3Lmz/P39vTuMmCCSlJTk9GUAAIBGSE1NVWJioneHETMiUv2HiYiIcPpyAABAPeTm5trBhOr3ca8OI9VTMyaIEEYAAPAulyuxoIAVAAA4ijACAAAcRRgBAACOIowAAABHEUYAAICjCCMAAMBRhBEAAOAowggAAHAUYQQAADiKMAIAABxFGAEAAI4ijAAAAEe5NoycLynTvz86pu/+c4vKyyucvhwAAFzLtWHEeHz5Hr25K13rD512+lIAAHAt14aR0KAATR2YYG+/uDnV6csBAMC1XBtGjOnDkuzHt3alK7uw2OnLAQDAlVwdRlISItUnPkLFpeVatv2k05cDAIAruTqMGNOGJtqPL25iqgYAACe4PoxMHZSg4EB/7U7L1c4TOY68CAAAuJnrw0j7tsGa2DfOPhmLKWQFAKDVuT6M1J6qeXXbCdt/BAAAtB7CiKRrropRQvs2yj1falfWAACA1kMYMU+Cv59urxodYaoGAIDWRRip8tUhifLzkz44eFqpZwpb+WUAAMC9CCNVEqPa6truMfb2EgpZAQBoNYSRWqYNrezIumTLcZWxeR4AAK2CMFLLjX07qX3bIKXlnNe6g6da5xUAAMDlCCO1hAR+tnneYjqyAgDQKggjF5mqeXt3us4UsHkeAAAtjTDyOX06R6hfQqRKyipsEzQAANCyCCMXMG1YUk3PkYqKihZ+CQAAcDfCyAV8eUBnhQT6a296nj45zuZ5AAB4TBj59a9/LT8/vzpHcnLyRc+fP3/+F84PDQ2Vp4tsE6RJKZWb571IzxEAADxrZKRv375KS0urOdatW3fJ8yMiIuqcf/ToUXnTVM1r20/qXDGb5wEA0FICG/wNgYGKi6scNagPMxrSkPM9xYgrO6hLdFsdO1OoFTvT9JXBlXvXAAAAh0dGDhw4oM6dO6tbt2664447dOzYsUuen5+fr65duyopKUlTpkzRrl27LvsYRUVFys3NrXM4snnekMoA8iI9RwAA8IwwMnz4cFsH8uabb2ru3Lk6cuSIrrvuOuXl5V3w/F69eun555/XsmXLtGDBApWXl2vUqFE6fvz4JR9nzpw5ioyMrDlMkHHCV4cmyt9P+vDIGX16qsCRawAAwNf5VTRh7Wp2drYd9fjDH/6ge+6557Lnl5SUqHfv3poxY4YeffTRS46MmKOaGRkxgSQnJ8fWoLSmb77wkd7bl6UHrr9KD0+8eLEuAACoy7x/m0GFy71/N2lpb/v27dWzZ08dPHiwXucHBQVp0KBBlz0/JCTEXnTtw+mOrC9tOa7SsnLHrgMAAF/VpDBi6kEOHTqk+Pj4ep1fVlamHTt21Pt8TzChdydFhwUrI7dIaw9kOX05AAC4O4z8+Mc/1po1a/Tpp59q/fr1uvXWWxUQEGCnXYyZM2dq9uzZNec/8sgjevvtt3X48GFt3bpVd955p13ae++998pbBAf669ZB1ZvnXbrWBQAAtPDSXlN4aoLH6dOnFRsbq2uvvVYbN260tw2zssbf/7N8c/bsWd13331KT09XVFSUhgwZYkNMnz595E3MVM3f1x3RO3sydCq/SDHtQpy+JAAAfEaTClg9rQCmJU35ywf6ODVbP7+pt+4b3c2RawAAwJu0SgGrm0yvKmQ17eG9IL8BAOA1CCP1dMuAeLUJCtDBzHxtS81u2VcFAAAXIYzUU3hokG7qV7kKaDEdWQEAaDaEkQaYNrSyPfxrH59UQVFp870KAAC4GGGkAa6+MlpXxoSpoLhMy3ektdyrAgCAixBGGsDsQHx71ejI4s2pLfWaAADgKoSRBrptcOXmeZs+PatDWfkt86oAAOAihJEG6hQRqut7dbS3GR0BAKDpCCONMG1YZc+Rl7ecUAmb5wEA0CSEkUYYl9xRMe2CbWv49/axeR4AAE1BGGmEoAB/fWVwZSHri/QcAQCgSQgjTdg8z3h3X6Yyc8837VUAAMDFCCON1L1jOw3pGqWy8gq9su1E874qAAC4CGGkGTqymvbwbJ4HAEDjEEaaYHL/zmobHKDDpwq0+ejZpvwoAABcizDSBO1CAnVzfzbPAwCgKQgjTTS9qufIGzvSlM/meQAANBhhpIkGd4lSt9gwFRaX6fWPTzb1xwEA4DqEkWbYPG961TLfF9k8DwCABiOMNAPTAC3A30/bjmXrQEZec/xIAABcgzDSDGLDQ2yLeIPN8wAAaBjCSDOpnqp5ZesJFZeWN9ePBQDA5xFGmsnYXrHqGB6i0wXFWr03s7l+LAAAPo8w0kwCA/x125CqjqwUsgIAUG+EkWZ0e1UYeW9fptJz2DwPAID6IIw0o26x7XT1FdEqr5Be3nq8OX80AAA+izDSzKZVdWQ1UzVsngcAwOURRprZTf3i7J41R08X6sMjZ5r7xwMA4HMII82sbXCgbhnA5nkAANQXYaQFTKvqObJ8Z5pyz5e0xEMAAOAzCCMtYGBSe/Xs1E7nS8r1GpvnAQBwSYSRFto8r3p0ZPGm1JZ4CAAAfAZhpIXcOihBQQF++vh4jvam57bUwwAA4PUIIy2kQ7sQTejdyd5evImeIwAAXAxhpBV6jizddlxFpWUt+VAAALgjjPz617+29RC1j+Tk5Et+z5IlS+w5oaGh6tevn5YvXy63GN0jVnERoTpbWKJ3drN5HgAAzTIy0rdvX6WlpdUc69atu+i569ev14wZM3TPPfdo27Ztmjp1qj127twpNwjw99NXq/areZHN8wAAaJ4wEhgYqLi4uJojJibmouc+/fTT+tKXvqSHH35YvXv31qOPPqrBgwfrmWeekVvcPrQyjLx/IEsns885fTkAAHh/GDlw4IA6d+6sbt266Y477tCxY8cueu6GDRs0YcKEOvdNnDjR3n8pRUVFys3NrXN4q64dwjSiW7QqKqSXtlDICgBAk8LI8OHDNX/+fL355puaO3eujhw5ouuuu055eXkXPD89PV2dOlWuKKlmPjf3X8qcOXMUGRlZcyQlVRaCeqvptTbPKzdb+gIAgMaFkUmTJun2229X//797QiHKUbNzs7W4sWL1Zxmz56tnJycmiM11bsbh01KiVd4aKCOnz2nDYdPO305AAD4ztLe9u3bq2fPnjp48OAFv25qSjIyMurcZz43919KSEiIIiIi6hzeLDQoQFMGdq4ZHQEAAM0URvLz83Xo0CHFx1fuUvt5I0eO1KpVq+rct3LlSnu/21S3h1+xM105hWyeBwBAo8LIj3/8Y61Zs0affvqpXbZ76623KiAgwC7fNWbOnGmnWKrNmjXL1pc89dRT2rt3r+1TsnnzZj344INym34JkUqOC1dxabmWfXzC6csBAMA7w8jx48dt8OjVq5emTZumDh06aOPGjYqNjbVfNytrTO+RaqNGjdLChQs1b948DRgwQC+99JJeffVVpaSkyG1Mg7jahawAAKCSX0WFWXTq2czSXrOqxhSzenP9yNmCYg1/fJWKy8r1xg+uVd/OkU5fEgAAjr9/szdNK4oKC9YNfas3z2N0BAAAwogDplcVsr66/aTOl7B5HgAAjIy0smu6xyihfRvlnCvR27vrLnsGAMCNCCMObJ53W9XmeUzVAABAGHHE7UMS5ecnrTt4SqlnCvl7CABwNUZGHJAU3VbXXFW52zGb5wEA3I4w4pBpVT1HTBgpY/M8AICLEUYccmOfTopsE6QT2ef0wcFTTl0GAACOI4w4uHne1KrN816kIysAwMUIIx4wVbNyV4btzgoAgBsRRhxk2sH37Rxh28O/up3N8wAA7kQYcVj15nkvbkqVF2wTBABAsyOMOGzKgAQFB/prb3qedpzIcfpyAABodYQRh0W2DdKklDh7ezGFrAAAFyKMeIBpVZvnLWPzPACACxFGPMDIbh2UFN1GeedLtWJnmtOXAwBAqyKMeAB/fz/dPqRydGTxpuNOXw4AAK2KMOIhbqvaPG/D4dM6errA6csBAKDVEEY8REL7NrquR6y9vWQzoyMAAPcgjHiQ6VWFrGyeBwBwE8KIB5nQp6Oi2gYpPfe81h7IcvpyAABoFYQRDxISGKCpgxLs7cWbUp2+HAAAWgVhxEPbw7+zJ0On84ucvhwAAFocYcTDJMdFaEBipErKKrR0G5vnAQB8H2HEA01j8zwAgIsQRjzQLQM6KzTIXwcy87U9NdvpywEAoEURRjxQRGiQbkqJt7fZPA8A4OsIIx4+VfPax2kqLC51+nIAAGgxhBEPNfzKaHXt0Fb5RaVaviPd6csBAKDFEEY8lJ+fn6ZVdWSl5wgAwJcRRjzYbYMT5e8nffTpGR3Oynf6cgAAaBGEEQ8WFxmqsb062ttLtrB5HgDANxFGPNy0oYn248tbjqu0rNzpywEAoNkRRjzcuORO6hAWrMy8Ir23j83zAAC+hzDi4YID/fWVwVWb521m8zwAgO8hjHiB6lU1q/dmKiuPzfMAAL6lSWHkiSeesEtQH3rooYueM3/+fHtO7SM0NLQpD+s6PTqFa1CX9iotr9ArWylkBQD4lkaHkU2bNum5555T//79L3tuRESE0tLSao6jR4829mFda3rV6MiLm1NVUVHh9OUAAOBsGMnPz9cdd9yhv/3tb4qKirrs+WY0JC4urubo1KlTYx7W1W4e0FltggJ0OKtAW4+ddfpyAABwNow88MADmjx5siZMmFDv8NK1a1clJSVpypQp2rVr1yXPLyoqUm5ubp3D7dqFBGpy/8rN817cRCErAMDFYWTRokXaunWr5syZU6/ze/Xqpeeff17Lli3TggULVF5erlGjRun48YvXPpifHRkZWXOYEANpetXmea9/kmb3rAEAwHVhJDU1VbNmzdK//vWvehehjhw5UjNnztTAgQM1ZswYvfLKK4qNjbX1Jhcze/Zs5eTk1BzmcSEN7RqlbjFhKiwu0xufnOQpAQC4L4xs2bJFmZmZGjx4sAIDA+2xZs0a/elPf7K3y8rKLvszgoKCNGjQIB08ePCi54SEhNii19oHKmtvbq/ePG8zq2oAAC4MI+PHj9eOHTu0ffv2mmPo0KG2mNXcDggIuOzPMIHF/Iz4+Mr6BzTMbUMSFODvpy1Hz+pgZh5PHwDA6wU25OTw8HClpKTUuS8sLEwdOnSoud9MySQkJNTUlDzyyCMaMWKEunfvruzsbD355JN2ae+9997bnH8O1+gYHqrre3XUO3sy7OjIz27q7fQlAQDgWR1Yjx07ZnuJVDt79qzuu+8+9e7dWzfddJNdGbN+/Xr16dOnuR/adYWspgFaCZvnAQC8nF+FF3TQMgHGrKoxxazUj8gGkJFzVutUfpGe+8YQTewb5/RLBABAo9+/2ZvGCwUF+NvaEWMxPUcAAF6OMOLlm+e9uy9TGbnnnb4cAAAajTDipa6KbadhV0SpvEJ6mc3zAABejDDixap7jizZfJzN8wAAXosw4sUm94tXWHCAjpwq0EdHzjh9OQAANAphxIuFhQTqlgGd7e0XN9MyHwDgnQgjPjJVs3xHmnLPlzh9OQAANBhhxMsN7tJePTu10/mScs1+ZQe1IwAAr0MY8YHN8x6/tZ8C/f30xidpem7tYacvCQCABiGM+IChV0TrV1/ua2///s29Wrs/y+lLAgCg3ggjPuLO4V00fWiS7Tvy/X9v07HThU5fEgAA9UIY8aHpmkem9tXApPbKOVeib/9zswqLS52+LAAALosw4kNCAgP07J1DFNMuRHvT8/TwS59Q0AoA8HiEER8TFxmquXcOpqAVAOA1CCM+aBgFrQAAL0IYcUlB69HTBU5fEgAAF0QYcUlB63f+uYWCVgCARyKMuKCgNTacglYAgOcijLihoPUOCloBAJ6LMOICdGgFAHgywohLUNAKAPBUhBGXoKAVAOCpCCMuQkErAMATEUZcWtAaFOCnNz5J07NrDjt9SQAAlyOMuLWg9Za+9vbv39qrNfuznL4kn2L6uizbfkInss85fSkA4BUCnb4AOOOO4V2080SOFm1K1fcXbtVr379WXTuE8XI0UeqZQt31/Ec6fKqy4+2ApPaa3C9Ok1LilRTdlucXAC7Ar6KiokIeLjc3V5GRkcrJyVFERITTl+MzikrLNP25jdqemq3kuHC98r1RahtMPm2sHcdzdPf8TTqVX6TwkEDlF5eq9r+ufgmRmtQvTjelxOuKGIIfAN+XW8/3b8KIy6XnnNctz6xTVl6RJveP1zMzBtmVN2iY9/Zl6nv/2qrC4jL1jo/Q/LuHyTyNb+1M1/Id6frwyGm7T1C1PvERusmMmPSL11Wx7Xi6AfgkwgjqbfOnZzTjbxtVUlahn3wpWfePvYpnrwEWb0rV7KU7VFZeoWu7x2junYMVHhpU5xwzWvL2rgwt35GmDYdP23Or9eoUrpv6xdtw0qNTOM89AJ9BGEGDLNh4VL94daf9bX7+3VdrTM9YnsHLMDOcT686oD++c8B+/pVBCXritv4KDrx0XfiZgmKt3F05YvLBwVMqrRVMundsVxNMTEhhlAqANyOMoMFvrLNf2WELWiNCAylovYzSsnIb3szzZXxv7FV6eGKvBoeHnMISvb07XSt2puv9A1l2dKpat5iwyhqTfvF2WodgAsDbEEbQqILWr83bqG3Hsu1v5aagNSyEgtbPKygq1YMLt+rdfVny95N+MyVF3xjRtVmWBK/em6E3PknX2gNZKi4tr/la1w5t7YocM2JiCmEJJgC8AWEEjZKRe143/7mqoLVfvJ75OgWttZnn5Z5/bNInx3MUGuSvP31tkG7sG9fsf9vyzptgkqkVO9L17r5MFdUKJolRbexoyaSUOA1Mak8wAeCxCCNotC1Hz9gREgpa6zqcla9vvrBJx84UKqptkP7+zWEa3CWqVUZiTCAxwcQElHMlZTVf6xwZalfkmBGTQUlR8jdDNQDgIQgjaJJ/fXhUP19KQWu1rcfO6p75m3S2sERdotvapbvdHFiSe664zC4jXr4zXav3ZKig+LNg0ikipGoqJ15DukYpgGACwEvCSJPawT/xxBN2iPihhx665HlLlixRcnKyQkND1a9fPy1fvrwpD4tW8PWru+hrw5Js0y7TofXo6cqOom60cneGvv63jTaI9E+M1Mv3j3IkiBhtggPsSMifZwzSll/eoHnfGKJbByXYJmsZuUWav/5TTXtug0bMWaVfvrpT6w+dssW2AODJGt30bNOmTZo2bZpNOtdff73++Mc/XvC89evXa/To0ZozZ45uvvlmLVy4UL/73e+0detWpaSk1Oux6MDqDApapX9uPKpfLdtpG5Zd3ytWf7ljsEd2qTWv1boDp+xyYbNsOPd8ac3XOoQFa2JKZefXEd2iFRjAllQAfGCaJj8/X4MHD9Zf//pXPfbYYxo4cOBFw8j06dNVUFCg119/vea+ESNG2O959tlnm/UPg5YtaDV1CX/5+mBXFEyafxa/f2uf5r53yH5uRokem5riFW/kZhXOB4dOacWONL29O0PZhSU1XzO1Ljf2idNN/eM16qoOCvKCPw8A79Wi0zQPPPCAJk+erAkTJlz23A0bNnzhvIkTJ9r7L6aoqMj+AWofcEaniFA9e+dgBQX42d+6566pfHP2ZebN/EeLP64JIj+c0FNzvtLPK4KIYZquXd+ro37/1QHa9PMJ+uc9V2vG1UmKDgu2U00vbk61m/kNfewd/XjJx3p3b2adZcQA0NoaPN68aNEiO8VipmnqIz09XZ06dapzn/nc3H8xZkrnN7/5TUMvDS1kSNdo/frLfW1B65Nv7bMNuMb26uiTz7dZUnv/gq1ad/CULQA1IWTa0CR5KzPycV2PWHs8OqVcHx05ozd2pOmtXek6lV+sl7Yct0d4aKBu6N3J1qNc1yNGoUEBTl86ABdpUBhJTU3VrFmztHLlSluM2lJmz56tH/3oRzWfm5GRpCTvfUPwBXcM76qdJ3L0749S9YN/b/PJDq1mSsqMGOxNz1Pb4AD99Y7BPhW6zMjOqO4x9nhkSoo2fXrGTuWY7q+ZeUV6ZdsJe7QLCdT43h3typyxvWIJJgA8K4xs2bJFmZmZtl6kWllZmdauXatnnnnGTq8EBNT9jSouLk4ZGRl17jOfm/svJiQkxB7wLGZ0xLxRmw6t3/7fLT7VofVARp4NIidzziumXYhe+OYw9UuMlK8yoz4junWwx69u6astx87aTfxML5P03PNatv2kPUwoM7Ulo3vGanSPWNsJ1g01QwBaV4MKWPPy8nT06NE6991999122e5PfvKTC66OMQWshYWFeu2112ruGzVqlPr3708BqxfyxYLWDw+f1n3/u9muQDH7wfzjW1crKbqt3Ki8vELbUrNrRkxOZJ+r8/Wk6DY2lJhpn1HdOyjic7sTA4AjTc/Gjh1bZzXNzJkzlZCQYOs+qpf2jhkzxvYkMUWvpubk8ccfZ2mvj3Ro/b9f6qXvje0ub/XGJ2n64YvbVVxWbhuF/c/MoYoKC3b6sjyC+a9h18lcu0/O2v1Z2nL0bJ2N/MzoyqCk9nbUxNSZ9E9sT6M1AI0KI80+xn7s2DH5+/vXGQUxvUV+8Ytf6Gc/+5l69OihV199td49RuB5fKWg9e/rjuixN3bbxm4T+3bS018bRH1ELWbEKyUh0h4mcJq29BsPn7bB5P0Dp3T4VIE2Hz1rjz+s3K/2bYN0TfcYje4RYwNKfGQb515coAWYVWcns88p9WyhUs+c04nsQg1MitINfeou0kDDNXlkpDXQZ8QzzX7lE1vQGhEaqP88eK2uiAnzmqmI3y7fY8OIcdfIrvqvW/ryW30DpZ4ptKMm7+8/Zfua5NVqtGZ079jOTumM7hmj4Vd2sN1jAU9muhWbmikTNEzgOH72nI6fKay5bb72+XdM0/Zg889vUGRbpiwvhL1p0CpdP2fM26itx7LVs1M7Lf3eNR5f0Hq+pEz/Z/HHdnmr8dNJyfrO6G5eX/fiCf+Jf3w8W2v2n7IjJ58cz7Zda2v3Prn6img7nWNGTZLjwnnO4cgvIln5RTZI24BRK3SYj2nZ51Va+y/uBZjdupOi2trds3eezLX1c3+aMUhfHtC51f4c3oQwglbhTQWtOYUluu+fm22vDfPbzH/fPkBTBiY4fVk+KbuwWB8cPK33q+pNzCql2mLDQ2wwGdMz1k7tmBVMQFOZgf4zBcVKNSMaVVMptUc4jmefu2yDv+AAfyVEtbFhIzGqrS3ath+rPo9pF1zzf9wTK/bq2TWHNHVgZ/3xa4N4AS+AMIJW4w0FrWZViFm6ezAz324q99w3hth+G2idN4hDWQU2lJhpHVN3cr6k7htCSkKEXaFjpnVMIbEZSQEuJOdciR3ZOF4TOD4b2TAfC2vtZH0hpvA6PjLUhg0zwmFWztnbVR87hYfKv547XptePbc/u0GRbYK05RcTvKZLc2sijKBVLfzwmH62dIfMLwymR4cnFbTuPpmrb77wkW3sFRcRqhfuHqbe8exx5OT03uZPz1at0jmlPWl1t3swvU1GdutQM6VzZUyYx462ofkVFpfaEY3PB43q+2pvAnkh5q+KCRTVAaN6RCMxujJ8xEWGNtueTGZ6cuhv37H7Py3+zkhdfWV0s/xcX0IYQavzxIJWs5PtdxdsUX5Rqa1rmX/31ercnlUeniQz77x9napX6ZwuKK7zdfOmYkZNxvSM0cirYuxvofC+0bGi0nK7IsuMXJijslC01qhG1e3Pv/4XYqZKEqqmTmpGNqpGOTq3D1VIYOsVS5vWAEu3ndB3xnTT7Em9W+1xvQVhBHJ7QevSbcf18JJPbEHaiG7Reu4bQ3kj84ICw91puTaUmHCy+eiZL/Q2GWh6m5jGaz1jNIDeJs3+/J8rKVNBcanOFZepoKhM50pK7cfKEFEdJj4LFfZ21ddrvq+4TOeKS6s+Vp5zmbrQOkzgtLUa7WvVbFSNbJh6jrbBnlMo/9rHJ/X9f2+zq8fe+dEYpy/H4xBG4NqCVvNbmNld+Pdv7rOf3zKgs/779v6t+tsSmof5TfrDI6a3ySk7rXM4q6DO180o3LVmOseGk1gluGTUy0wPFJZUvtHXHm34QkgwH6u+XnChMGFCRMlnYcIEkZbWJijATsWZIubqItHPplQqp1O8qbOvqWEZ8uhK+0vPmofH+tyeXU1FGIFHFLQ+PLGXHri+9Qpay8or9Kv/7NSCjcfs598e3U0//VJyvQvS4NnMsL7ZUdmMmpiPn+9tclVsWNWUTqyGd4tu1d+gTUA4X1pul4+bkGBGCk2hrnmDN/fVvl1kzqm673yt25/d/7nzS8vtzzxfWhkaLrcipKnM7w9tTWgICbTBwTyPlR8rj7DgQNs3xox8mnARFhKgNsGBCqt1rr0vKLDqa1XfExTgk/8WzYjwhsOn9atb+ujua650+nI8CmEEHlPQ+vw3h+n6VihoNf9Z/2DRNq3cnWEf979u5j8G3+9tklNVa5Kl7amf620S4K+hV0TZIthru8fY3YjNm7l9k696Y29wKKgdLkpN6Kg819yuPZ3UWsy01RcCQk1Q+GIY+GKwCFTbkNpfq7zP9NKgaLj+/uf9w3rsjT3279mCe4e30KvtnQgjcNzsV3bo3x8da5WCVtNb4J5/bLI7CptloU9PH6hJ/eJb7PHgmX1kTCfYyt4mp76wyV9rCgn0t1sLmJEA88Zublcete//3NcCTWio+jwwQKHB5mPV+fZ25fm1A4d5HEKD8w5n5WvcU2ts/6Ktv7xB4V40zeSze9MA1X795T7al55rC1q//c/NLVbQevR0gb75wiYdOVVgC9/+566hGnYFS+zcxrTjvqlfvD1M3ZDZO6d6hc6mI2dUXlHxhVBQOzCE1Lr9WViodb4NC5/dDv1cQKi8XRkQfHEqAhfXLbad3fHb/J0zK8P4Rajh2JsGLSqzqqA1s4UKWj9OzbYjIqfyi23x4j++NUzdO4Y3288HgPp47PXd+p91R3Tb4EQ9NW0AT1oDR0ZoF4cW1TEiVHPvHGyHL5fvSNdf3zvUbD/73b2ZtlDWBBGzc/DS740iiABwxLjelXVx7+3LtIX0aBjCCFrckK7R+s2XU+zt/357n97dl9nkn/nipmO6938322JC06nzxe+MsMEHAJxgpobDQwNt0zazaSQahjCCVvH14V004+oudvvtWf/epk9P1e0XUV+mFuD/rdyvn7y8w/72YYZEzWodCsYAOMm0mDdLyo1VezJ4MRqIMIJWLWgd3KW93VvCFLSaFu0NUVJWrp+8/ImeXnXAfv79cd1tM7Pm2mcCAJpifNVUzao9TR/9dRv+F0erMR1Qn71ziDqGh2h/Rr4eXvKxHemoD9Nl8t5/bNbizcdlFir89tYU/Z8be7GsEYDHGNuzo/3/aW96nqNLy70RYQQOFLQOsQWtK3bWr6DVbKQ2fd4GrdmfZZdRzvvGUN0xvGurXC8A1FdUWLCGdI2yt1czVdMghBG0OvOP9ZEp9StoPZSVr6/8db12nshVdFiwFn17pCb06dSKVwsA9TcuufL/p1V7mappCMIIHGGKWU1R66UKWs0eN7fNXW+3Fe/aoa1euX+U3bEVADzVhKq6kfWHTtsNCVE/hBE4xmwqdbGC1jd3puvrf/tQ2YUlGpAYqZfvH9Wi7eQBoDl079hOSdFt7GaGphsr6ocwAo8raP3fDZ/q/n9tsTuVjk/uqH9/e4Ri2oXwSgHweKbD9PiqqZrVTNXUG2EEHlXQevuzG/Rfy3bZ6RszlfPcN4a06jbwANBsS3z3Zqqcbqz1QhiBRxW0bj561n788Y099fitKQqkhwgAL3P1ldEKCw5QVl6Rdp7McfpyvAJhBB7BjIJ8Z0w3tW8bpP++fYAeHNeDHiIAvHYKenRNN1ZW1dQHYQQeY/ak3tr2yxv01SGJTl8KADTJuOTqqRpaw9cHYQQeV/wFAN7u+uSOMv+dmR5J6Tnnnb4cj0cYAQCgmZkVgNV9kZpjp3JfRxgBAKAFmNYEBrv4Xh5hBACAFmwNv+7gKZ0vKeM5vgTCCAAALaB3fLg6R4bqfEm5Nhw6zXN8CYQRAABaqCB/XFUDtHfYxfeSCCMAALSQ2q3hzXYXuDDCCAAALWTkVR3UJihAaTnntSctj+f5IggjAAC0kNCgAF3TPcbeZlVNM4WRuXPnqn///oqIiLDHyJEjtWLFioueP3/+fDtnVvsIDQ1tyEMCAOAzG+fhwhq0HWpiYqKeeOIJ9ejRw859/eMf/9CUKVO0bds29e3b94LfY0LLvn37aj6nwyYAwI2t4T8+nm03z4sND3H6krx7ZOSWW27RTTfdZMNIz5499dvf/lbt2rXTxo0bL/o9JnzExcXVHJ06VRbzAADgBp0iQtUvIVKmfpVurM1cM1JWVqZFixapoKDATtdcTH5+vrp27aqkpCQ7irJr167GPiQAAF49OrKaXXybJ4zs2LHDjoaEhITou9/9rpYuXao+ffpc8NxevXrp+eef17Jly7RgwQKVl5dr1KhROn78+CUfo6ioSLm5uXUOAAC81YTelbMC7x/IUlEp3VibHEZMwNi+fbs+/PBD3X///brrrru0e/fuC55rRkxmzpypgQMHasyYMXrllVcUGxur55577pKPMWfOHEVGRtYcZlQFAABv1bdzhDqGh6iguEwfHj7j9OV4fxgJDg5W9+7dNWTIEBsaBgwYoKeffrpe3xsUFKRBgwbp4MGDlzxv9uzZysnJqTlSU1MbepkAAHgMf3+/z6ZqWFXT/H1GzNSLmVapb52JmeaJj4+/5HlmCqh6+XD1AQCANxtfNVWzam8G3VibsrTXjFhMmjRJXbp0UV5enhYuXKj33ntPb731lv26mZJJSEiwIybGI488ohEjRtiRlOzsbD355JM6evSo7r333oY8LAAAXu+a7h0UHOiv1DPndCAzXz07hTt9Sd4ZRjIzM23gSEtLs7UcpgGaCSI33HCD/fqxY8fk7//ZYMvZs2d13333KT09XVFRUXZqZ/369RcteAUAwFe1DQ7UqKs66L19WVq1J5MwUotfhRfs3GNW05jwY+pHmLIBAHirf248ql++ulPDrojSku+Okq/Lref7N3vTAADQSqqLWLccPauzBcU871UIIwAAtJKE9m2UHBeu8grpvf3sVVONMAIAgAMN0EzdCCoRRgAAaEXjqnbxXbM/SyVl5Tz3hBEAAFrXgMT26hAWrLzzpdr0Kd1YCSMAALSyAH8/Xc/GeXUwTQMAQCsbXxVGVtEa3iKMAADQyq7tEaOgAD8dOVWgw1n5rn/+CSMAALSy8NAgjejWwd5ezegIYQQAACcboL2zJ8P1LwAjIwAAOGB8cmW/kU2fnlXOuRJXvwaEEQAAHNClQ1v16NhOZeUVWrs/y9WvAWEEAACHG6CtcvlUDWEEAACHp2re25+lUhd3YyWMAADgkMFd2qt92yBlF5ZoW2q2a18HwggAAA4JDPDX2J6xcvuqGsIIAAAOGle1i+9qF+/iSxgBAMBBY3rG2v1qDmTm69jpQle+FoQRAAAcFNkmSMOuiLK3V+1151QNYQQAAA9ZVbPapa3hCSMAADhsfFW/kY2HTyu/qFRuQxgBAMBh3WLb6cqYMJWUVeh9F3ZjJYwAAOBBG+etcuFUDWEEAAAPmqp5d2+myssr5CaEEQAAPMCwK6IVHhqo0wXF2n7cXd1YCSMAAHiAoAB/23PEjQ3QCCMAAHjYVM07LmsNTxgBAMBDjO3ZUf5+0t70PJ3IPie3IIwAAOAhosKCNaRrlOsaoBFGAADwIOOqurGuctFUDWEEAAAPrBtZf+i0Covd0Y2VMAIAgAfp0bGdkqLbqLi0XB8cPC03IIwAAOBB/Pz8ajbOc8tUDWEEAAAPbQ2/2iXdWAkjAAB4mOHdohUWHKDMvCLtOpkrX0cYAQDAw4QEBui6HrGuaYDWoDAyd+5c9e/fXxEREfYYOXKkVqxYccnvWbJkiZKTkxUaGqp+/fpp+fLlTb1mAAB83riqVTVu6DfSoDCSmJioJ554Qlu2bNHmzZs1btw4TZkyRbt27brg+evXr9eMGTN0zz33aNu2bZo6dao9du7c2VzXDwCAT7q+V0f5+Uk7TuQoI/e8fJlfRUVFkypjoqOj9eSTT9rA8XnTp09XQUGBXn/99Zr7RowYoYEDB+rZZ5+t92Pk5uYqMjJSOTk5dkQGAAA3mPqXD7Q9NVtzvtJPM67uIm9T3/fvRteMlJWVadGiRTZsmOmaC9mwYYMmTJhQ576JEyfa+y+lqKjI/gFqHwAAuM34qlU1q3x8F98Gh5EdO3aoXbt2CgkJ0Xe/+10tXbpUffr0ueC56enp6tSpcq10NfO5uf9S5syZY5NU9ZGUlNTQywQAwOuN7135HvrBwVM6X1ImX9XgMNKrVy9t375dH374oe6//37ddddd2r17d7Ne1OzZs+2QTvWRmprarD8fAABv0Ds+XPGRoTpXUqYNh3y3G2uDw0hwcLC6d++uIUOG2BGMAQMG6Omnn77guXFxccrIqLskyXxu7r8UM+pSvWKn+gAAwI3dWMdVT9Xs9d0lvk3uM1JeXm5rPC7E1JKsWrWqzn0rV668aI0JAACoa0LVVM3qPZlq4poTjxXY0OmTSZMmqUuXLsrLy9PChQv13nvv6a233rJfnzlzphISEuyIiTFr1iyNGTNGTz31lCZPnmwLXs2S4Hnz5rXMnwYAAB8z8qoOCg3y18mc89qTlqc+nSPcPTKSmZlpA4epGxk/frw2bdpkg8gNN9xgv37s2DGlpaXVnD9q1CgbWEz4MNM5L730kl599VWlpKQ0/58EAAAfFBoUoGu7x9jbq310qqbJfUZaA31GAABu9u+Pjmn2Kzs0qEt7Lf3eNfIWLd5nBAAAtF43VsM0QDuVf+E6TW9GGAEAwMPFRYYqJSFCZi7jXR/cq4YwAgCAFxifXLWqhjACAACcML5qF9+1+7NUVOpb3VgZGQEAwAukdI5UbHiICorL9NGRM/IlhBEAALyAv7+fz26cRxgBAMBLjKvVGt4LOnPUG2EEAAAvcW2PGAUH+iv1zDkdzMyXryCMAADgJdoGB2rUVR3s7VU+tKqGMAIAgBcZX1M34jut4QkjAAB4keurwsiWo2d1tqBYvoAwAgCAF0mMaqvkuHCVV0hr9mfJFxBGAADw0gZo7/jIVA1hBAAALzOuqjW8GRkpKSuXtyOMAADgZQYmtVeHsGDlnS/V5k/PytsRRgAA8DIB/n4a28t3VtUQRgAA8OK6kdU+0G+EMAIAgBe6rkeMggL8dPhUgQ5neXc3VsIIAABeKDw0SMOv7OAToyOEEQAAvH3jvD2EEQAA4GDdyKZPzyjnXInXvgaMjAAA4KW6dghT947tVFpeobVe3I2VMAIAgA9snLfai+tGCCMAAHix8b0ru7G+uy9TZWbDGi9EGAEAwIsN7tJekW2ClF1Yoq3HvLMbK2EEAAAvFhjgr7G9Yr16VQ1hBAAAH5mqWb3XO1vDE0YAAPByY3rE2v1q9mfkK/VMobwNYQQAAC8X2TZIQ7tGee3GeYQRAAB8wISqqZpVXrjElzACAIAPGFfVjfXDw2eUX1Qqb0IYAQDAB1wV205XxoSpuKxc6w54VzdWwggAAD62cd47XrbElzACAICPtYZ/d2+myr2oGythBAAAHzHsymiFhwTqdEGxPj6eLZ8MI3PmzNGwYcMUHh6ujh07aurUqdq3b98lv2f+/Pny8/Orc4SGhjb1ugEAwOcEBfhrtBd2Y21QGFmzZo0eeOABbdy4UStXrlRJSYluvPFGFRQUXPL7IiIilJaWVnMcPXq0qdcNAAAuMVXjTUt8Axty8ptvvvmFUQ8zQrJlyxaNHj36ot9nRkPi4uIaf5UAAKBexvbqKH8/aU9ark5mn1Pn9m3k0zUjOTk59mN0dPQlz8vPz1fXrl2VlJSkKVOmaNeuXZc8v6ioSLm5uXUOAABwedFhwRrcJcqrRkcaHUbKy8v10EMP6ZprrlFKSspFz+vVq5eef/55LVu2TAsWLLDfN2rUKB0/fvyStSmRkZE1hwkxAACgYQ3QVntJa3i/ioqKRq39uf/++7VixQqtW7dOiYmJ9f4+U2fSu3dvzZgxQ48++uhFR0bMUc2MjJhAYkZiTP0JAAC4uP0Zebrx/61VcKC/tv/XDWob3KCqjGZj3r/NoMLl3r8bNTLy4IMP6vXXX9e7777boCBiBAUFadCgQTp48OBFzwkJCbEXXfsAAAD106NjOyVGtVFxabk+OHhanq5BYcQMopggsnTpUq1evVpXXnllgx+wrKxMO3bsUHx8fIO/FwAAXJ5ZOFK9qmb13gzfCiNmWa+p+1i4cKHtNZKenm6Pc+fO1Zwzc+ZMzZ49u+bzRx55RG+//bYOHz6srVu36s4777RLe++9997m/ZMAAIAa46t38d2TaQcTPFmDJpHmzp1rP44dO7bO/S+88IK++c1v2tvHjh2Tv/9nGefs2bO67777bGiJiorSkCFDtH79evXp06d5/gQAAOALhneLVtvgAGXmFWnniVz1S4yUzxWwemIBDAAA+Mx3/rlZb+3K0EMTeuihCT3lUwWsAADAe6ZqVnt4vxHCCAAAPur6XpVFrJ8cz1FG7nl5KsIIAAA+KjY8RAOS2tvb73rw6AhhBAAAHzbBCzbOI4wAAOCC1vDrDpzS+ZIyeSLCCAAAPqxPfITiI0N1rqRMGw57ZjdWwggAAD7ejXVcdTfWPZ45VUMYAQDAx42vmqpZtSfDI7uxEkYAAPBxo66KUWiQv07mnNfe9Dx5GsIIAAA+LjQoQNd2j/HYBmiEEQAAXGBccmU31nf2eN4uvoQRAABcYFxVEev21Gydyi+SJyGMAADgAnGRoUpJiJCpX31vX5Y8CWEEAACXTdWs8rCpGsIIAAAuMb5qqmbt/iwVl5bLUxBGAABwiX4JkXbzvILiMn105Iw8BWEEAACX8Pf307heHT1uVQ1hBAAAF26ct2qv53RjJYwAAOAi13aPUXCgv1LPnNOhrHx5AsIIAAAuEhYSqJHdOtjb73jIxnmEEQAAXLpx3mrCCAAAcLIb6+ajZ5RdWCynMTICAIDLJEa1VXJcuMo9pBsrYQQAABePjqzygF18CSMAALjQ+N6VreHX7MtUSZmz3VgJIwAAuNDApPaKDgtW7vlSbf70rKPXQhgBAMCFAvz9NLZXrL29eq+z3VgJIwAAuNSEqqkap+tGCCMAALjUdT1iFOjvp8NZBTpyqsCx6wh07JEBAICjwkODdO913dQpIkSRbYIcuw7CCAAALvbTSclOXwLTNAAAwFnUjAAAAEcRRgAAgKMIIwAAwFGEEQAA4D1hZM6cORo2bJjCw8PVsWNHTZ06Vfv27bvs9y1ZskTJyckKDQ1Vv379tHz58qZcMwAAcGsYWbNmjR544AFt3LhRK1euVElJiW688UYVFFy8Ucr69es1Y8YM3XPPPdq2bZsNMObYuXNnc1w/AADwcn4VFRUVjf3mrKwsO0JiQsro0aMveM706dNtWHn99ddr7hsxYoQGDhyoZ599tl6Pk5ubq8jISOXk5CgiIqKxlwsAAFpRfd+/m1QzYn64ER0dfdFzNmzYoAkTJtS5b+LEifb+iykqKrJ/gNoHAADwTY0OI+Xl5XrooYd0zTXXKCUl5aLnpaenq1Onyo14qpnPzf2Xqk0xSar6SEpKauxlAgAAXw0jpnbE1H0sWrSoea9I0uzZs+2oS/WRmpra7I8BAAA8Q6P2pnnwwQdtDcjatWuVmJh4yXPj4uKUkZFR5z7zubn/YkJCQuwBAAB8X4NGRkytqwkiS5cu1erVq3XllVde9ntGjhypVatW1bnPrMQx9wMAAAQ2dGpm4cKFWrZsme01Ul33Yeo62rRpY2/PnDlTCQkJtu7DmDVrlsaMGaOnnnpKkydPttM6mzdv1rx58+r9uNULfihkBQDAe1S/b1924W5FA5jTL3S88MILNeeMGTOm4q677qrzfYsXL67o2bNnRXBwcEXfvn0r3njjjYY8bEVqaupFH5uD54C/A/wd4O8Afwf4OyCPfg7M+/ilNKnPSGsxK3dOnjxpR2P8/PyaNbGZlTqmQJb+Jc7j9fA8vCaehdfDs/B6XJ6JGHl5eercubP8/f2bt4C1tZk/wOUKZZvCBBHCiOfg9fA8vCaehdfDs/B6XJop5bgcNsoDAACOIowAAABHuTqMmF4mv/rVr+hp4iF4PTwPr4ln4fXwLLwezccrClgBAIDvcvXICAAAcB5hBAAAOIowAgAAHEUYAQAAjnJ1GPnLX/6iK664QqGhoRo+fLg++ugjpy/Jlcw+RsOGDbMddjt27KipU6dq3759Tl8WqjzxxBO28/FDDz3Ec+KQEydO6M4771SHDh3sPmD9+vWze3zBGWVlZfrlL39pN4s1r8dVV12lRx999PL7r+CiXBtGXnzxRf3oRz+yS3u3bt2qAQMGaOLEicrMzHT60lxnzZo1dhPGjRs32h2dS0pKdOONN6qgoMDpS3O9TZs26bnnnlP//v1d/1w45ezZs7rmmmsUFBSkFStWaPfu3Xbj0aioKF4Th/zud7/T3Llz9cwzz2jPnj3289///vf685//zGvSSK5d2mtGQsxv4+YvU/X+N2afmu9///v66U9/6vTluVpWVpYdITEhZfTo0U5fjmvl5+dr8ODB+utf/6rHHntMAwcO1B//+EenL8t1zP9HH3zwgd5//32nLwVVbr75ZnXq1El///vfa56T2267zY6SLFiwgOepEVw5MlJcXKwtW7ZowoQJdfa/MZ9v2LDB0WuDlJOTY5+G6Ohong4HmdGqyZMn1/l3gtb3n//8R0OHDtXtt99uQ/qgQYP0t7/9jZfCQaNGjdKqVau0f/9++/nHH3+sdevWadKkSbwujeQVG+U1t1OnTtk5P5NsazOf792717HrQuUIlalNMMPSKSkpPCUOWbRokZ2+NNM0cNbhw4ftlICZVv7Zz35mX5Mf/OAHCg4O1l133cXL49BoldmxNzk5WQEBAfb95Le//a3uuOMOXo9GcmUYgWf/Nr5z5077WwackZqaqlmzZtn6HVPcDecDuhkZefzxx+3nZmTE/Bt59tlnCSMOWbx4sf71r39p4cKF6tu3r7Zv325/iercuTOvSSO5MozExMTYNJuRkVHnfvN5XFycY9fldg8++KBef/11rV27VomJiU5fjmuZKUxTyG3qRaqZ3/zM62JqrIqKiuy/H7SO+Ph49enTp859vXv31ssvv8xL4JCHH37Yjo587Wtfs5+b1U1Hjx61KwMZrWocV9aMmOHNIUOG2Dm/2r99mM9Hjhzp6LW5kamhNkFk6dKlWr16tV0uB+eMHz9eO3bssL/tVR/mN3MzBG1uE0Ral5my/PxSd1Or0LVr11a+ElQrLCy0dYa1mX8X5n0EjePKkRHDzL+aBGv+k7366qvtKgGzlPTuu+92+tJcOTVjhjuXLVtme42kp6fb+yMjI211OlqXeQ0+X68TFhZme1xQx9P6fvjDH9qCSTNNM23aNNsPad68efaAM2655RZbI9KlSxc7TbNt2zb94Q9/0Le+9S1eksaqcLE///nPFV26dKkIDg6uuPrqqys2btzo9CW5kvlreKHjhRdecPrSUGXMmDEVs2bN4vlwyGuvvVaRkpJSERISUpGcnFwxb948XgsH5ebm2n8P5v0jNDS0olu3bhU///nPK4qKinhdGsm1fUYAAIBncGXNCAAA8ByEEQAA4CjCCAAAcBRhBAAAOIowAgAAHEUYAQAAjiKMAAAARxFGAACAowgjAADAUYQRAADgKMIIAABwFGEEAADISf8fpZp7+lypWYUAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "def _get_extrusion_points(\n", + " cart_points,\n", + " men_points,\n", + " ml_axis,\n", + " side,\n", + "):\n", + " cart_points_ml = np.dot(cart_points, ml_axis)\n", + " men_points_ml = np.dot(men_points, ml_axis)\n", + "\n", + " if side in ['med', 'medial']:\n", + " cart_edge = np.min(cart_points_ml)\n", + " men_edge = np.min(men_points_ml)\n", + " extrusion = cart_edge - men_edge\n", + " elif side in ['lat', 'lateral']:\n", + " cart_edge = np.max(cart_points_ml)\n", + " men_edge = np.max(men_points_ml)\n", + " extrusion = men_edge - cart_edge\n", + " else:\n", + " raise ValueError(f'Invalid side: {side}, must be one of: med, medial, lat, lateral')\n", + "\n", + " return extrusion\n", + "\n", + "\n", + "# breakup points by bins in the AP direction. \n", + "n_bins = 10\n", + "\n", + "tib_cart_label = 'labels'\n", + "med_tib_cart_label = 2\n", + "lat_tib_cart_label = 3\n", + "\n", + "med_cart_points = tibia[tib_cart_label] == med_tib_cart_label\n", + "med_cart_points = tibia.points[med_cart_points]\n", + "\n", + "med_men_points = med_meniscus.points\n", + "\n", + "ap_axis = dict_tibia_axes['ap_axis']\n", + "\n", + "# project med_cart_points onto the AP axis\n", + "med_cart_points_ap = np.dot(med_cart_points, ap_axis)\n", + "min_med_cart_ap = np.min(med_cart_points_ap)\n", + "max_med_cart_ap = np.max(med_cart_points_ap)\n", + "bins = np.linspace(min_med_cart_ap, max_med_cart_ap, n_bins+1)\n", + "\n", + "med_men_points_ap = np.dot(med_men_points, ap_axis)\n", + "\n", + "list_extrusions = []\n", + "for i in range(n_bins):\n", + " bin_start = bins[i]\n", + " bin_end = bins[i+1]\n", + " bin_mask_cart = (med_cart_points_ap >= bin_start) & (med_cart_points_ap < bin_end)\n", + " bin_cart_points = med_cart_points[bin_mask_cart]\n", + " bin_mask_men = (med_men_points_ap >= bin_start) & (med_men_points_ap < bin_end)\n", + " bin_men_points = med_men_points[bin_mask_men]\n", + " \n", + " extrusion = _get_extrusion_points(\n", + " cart_points=bin_cart_points,\n", + " men_points=bin_men_points,\n", + " ml_axis=ml_axis,\n", + " side='med',\n", + " )\n", + " list_extrusions.append(extrusion)\n", + " \n", + "plt.plot(list_extrusions)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "id": "3c65964d", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Med Men Extrusion: 3.89 mm\n", + "Lat Men Extrusion: -0.47 mm\n", + "Lat Men Shift Extrusion: 4.33 mm\n", + "Med Men Shift Extrusion: -1.22 mm\n" + ] + } + ], + "source": [ + "def get_extrusion_stats_percentile(\n", + " tibia_mesh,\n", + " regions_label,\n", + " cart_label,\n", + " meniscus_mesh,\n", + " ap_axis,\n", + " ml_axis,\n", + " side,\n", + " middle_percentile_range=0.1,\n", + "):\n", + " cart_indices = tibia_mesh[regions_label] == cart_label\n", + " cart_points = tibia_mesh.points[cart_indices]\n", + " men_points = meniscus_mesh.points\n", + "\n", + " # project med_cart_points onto the AP axis\n", + " cart_points_ap = np.dot(cart_points, ap_axis)\n", + " min_cart_ap = np.min(cart_points_ap)\n", + " max_cart_ap = np.max(cart_points_ap)\n", + " \n", + " # get the middle +/- middle_percentile_range/2 of the med_cart_points_ap\n", + " # along the AP axis\n", + " middle_ap_cartilage = (min_cart_ap + max_cart_ap) / 2\n", + " min_max_ap_cartilage_range = max_cart_ap - min_cart_ap\n", + " plus_minus_ap_cartilage_range = min_max_ap_cartilage_range * middle_percentile_range / 2 \n", + " lower_ap_cartilage = middle_ap_cartilage - plus_minus_ap_cartilage_range\n", + " upper_ap_cartilage = middle_ap_cartilage + plus_minus_ap_cartilage_range\n", + " \n", + " # get the points along the AP axis that are within the lower_ap_cartilage and upper_ap_cartilage\n", + " ap_cart_indices = (cart_points_ap >= lower_ap_cartilage) & (cart_points_ap <= upper_ap_cartilage)\n", + " # ap_cart_points = med_cart_points[ap_cart_indices]\n", + " \n", + " # project meniscus points onto the AP axis\n", + " men_points_ap = np.dot(men_points, ap_axis)\n", + " \n", + " # get the points along the AP axis that are within the lower_ap_cartilage and upper_ap_cartilage\n", + " ap_men_indices = (men_points_ap >= lower_ap_cartilage) & (men_points_ap <= upper_ap_cartilage)\n", + " # ap_men_points = men_points[ap_men_indices]\n", + " \n", + " # we now have the ap_men_indices and ap_cart_indices\n", + " # we now need to extract the ml projected coordinates for these\n", + " # ap points\n", + " ml_cart_points = cart_points[ap_cart_indices]\n", + " ml_men_points = men_points[ap_men_indices]\n", + " \n", + " # get the extrusion for each point\n", + " extrusion = _get_extrusion_points(\n", + " cart_points=ml_cart_points,\n", + " men_points=ml_men_points,\n", + " ml_axis=ml_axis,\n", + " side=side,\n", + " )\n", + " \n", + " return extrusion\n", + "\n", + "\n", + "med_men_extrusion = get_extrusion_stats_percentile(\n", + " tibia_mesh=tibia,\n", + " regions_label='labels',\n", + " cart_label=2,\n", + " meniscus_mesh=med_meniscus,\n", + " ap_axis=dict_tibia_axes['ap_axis'],\n", + " ml_axis=dict_tibia_axes['ml_axis'],\n", + " middle_percentile_range=0.1,\n", + " side='med',\n", + ")\n", + "\n", + "lat_men_extrusion = get_extrusion_stats_percentile(\n", + " tibia_mesh=tibia,\n", + " regions_label='labels',\n", + " cart_label=3,\n", + " meniscus_mesh=lat_meniscus,\n", + " ap_axis=dict_tibia_axes['ap_axis'],\n", + " ml_axis=dict_tibia_axes['ml_axis'],\n", + " middle_percentile_range=0.1,\n", + " side='lat',\n", + ")\n", + "\n", + " \n", + " \n", + " \n", + "print(f'Med Men Extrusion: {med_men_extrusion:.2f} mm')\n", + "print(f'Lat Men Extrusion: {lat_men_extrusion:.2f} mm')\n", + "\n", + "late_men_shifted = lat_meniscus.copy()\n", + "late_men_shifted.points -= [5, 0, 0]\n", + "\n", + "lat_men_shift_extrusion = get_extrusion_stats_percentile(\n", + " tibia_mesh=tibia,\n", + " regions_label='labels',\n", + " cart_label=3,\n", + " meniscus_mesh=late_men_shifted,\n", + " ap_axis=dict_tibia_axes['ap_axis'],\n", + " ml_axis=dict_tibia_axes['ml_axis'],\n", + " middle_percentile_range=0.1,\n", + " side='lat',\n", + ")\n", + "\n", + "print(f'Lat Men Shift Extrusion: {lat_men_shift_extrusion:.2f} mm')\n", + "\n", + "med_men_shifted = med_meniscus.copy()\n", + "med_men_shifted.points -= [5, 0, 0]\n", + "\n", + "med_men_shift_extrusion = get_extrusion_stats_percentile(\n", + " tibia_mesh=tibia,\n", + " regions_label='labels',\n", + " cart_label=2,\n", + " meniscus_mesh=med_men_shifted,\n", + " ap_axis=dict_tibia_axes['ap_axis'],\n", + " ml_axis=dict_tibia_axes['ml_axis'],\n", + " middle_percentile_range=0.1,\n", + " side='med',\n", + ")\n", + "\n", + "print(f'Med Men Shift Extrusion: {med_men_shift_extrusion:.2f} mm')\n" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "id": "ca0e7dc9", + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "b5eebd1bcd784711b74e95a0fcce88ce", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'name': '_points', 'numberO…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "view(geometries=[tibia, late_men_shifted, med_men_shifted])" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "id": "1f35114b", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{'med_men_extrusion_middle_0.1p_mm': np.float32(3.8892097), 'lat_men_extrusion_middle_0.1p_mm': np.float32(-0.4709015)}\n" + ] + } + ], + "source": [ + "def get_med_lat_extrusion(\n", + " tibia_mesh,\n", + " regions_label,\n", + " med_cart_label,\n", + " lat_cart_label,\n", + " med_meniscus_mesh,\n", + " lat_meniscus_mesh,\n", + " ap_axis,\n", + " ml_axis,\n", + " middle_percentile_range=0.1,\n", + "):\n", + " med_men_extrusion = get_extrusion_stats_percentile(\n", + " tibia_mesh=tibia,\n", + " regions_label='labels',\n", + " cart_label=med_cart_label,\n", + " meniscus_mesh=med_meniscus_mesh,\n", + " ap_axis=ap_axis,\n", + " ml_axis=ml_axis,\n", + " middle_percentile_range=middle_percentile_range,\n", + " side='med',\n", + " )\n", + " \n", + " lat_men_extrusion = get_extrusion_stats_percentile(\n", + " tibia_mesh=tibia,\n", + " regions_label='labels',\n", + " cart_label=lat_cart_label,\n", + " meniscus_mesh=lat_meniscus_mesh,\n", + " ap_axis=ap_axis,\n", + " ml_axis=ml_axis,\n", + " middle_percentile_range=middle_percentile_range,\n", + " side='lat',\n", + " )\n", + " \n", + " dict_extrusions = {\n", + " f'med_men_extrusion_middle_{middle_percentile_range}p_mm': med_men_extrusion,\n", + " f'lat_men_extrusion_middle_{middle_percentile_range}p_mm': lat_men_extrusion,\n", + " }\n", + " return dict_extrusions\n", + " \n", + "\n", + "dict_extrusions = get_med_lat_extrusion(\n", + " tibia_mesh=tibia,\n", + " regions_label='labels',\n", + " med_cart_label=2,\n", + " lat_cart_label=3,\n", + " med_meniscus_mesh=med_meniscus,\n", + " lat_meniscus_mesh=lat_meniscus,\n", + " ap_axis=dict_tibia_axes['ap_axis'],\n", + " ml_axis=dict_tibia_axes['ml_axis'],\n", + " middle_percentile_range=0.1,\n", + ")\n", + "\n", + "\n", + "print(dict_extrusions)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "mskt", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.14" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/testing/scratch/meniscal_extrusion_function_test_Oct.20.2025.ipynb b/testing/scratch/meniscal_extrusion_function_test_Oct.20.2025.ipynb new file mode 100644 index 0000000..ad0e10e --- /dev/null +++ b/testing/scratch/meniscal_extrusion_function_test_Oct.20.2025.ipynb @@ -0,0 +1,113 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "e24e17aa", + "metadata": {}, + "outputs": [], + "source": [ + "import pymskt as mskt\n", + "from itkwidgets import view\n", + "import numpy as np\n", + "import os\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "e59a475b", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Initiating tibia mesh\n", + "Creating tibia mesh\n", + "Calculating cartilage thickness\n", + "WARNING: Mesh is now synonymous with pyvista.PolyData and thus this property is redundant and the Mesh object can be used for anything that pyvista.PolyData or vtk.vtkPolyData can be used for.\n", + "WARNING: Mesh is now synonymous with pyvista.PolyData and thus this property is redundant and the Mesh object can be used for anything that pyvista.PolyData or vtk.vtkPolyData can be used for.\n", + "INTERSECTION IS: 2\n", + "INTERSECTION IS: 2\n", + "Assigning cartilage regions\n", + "INTERSECTION IS: 2\n", + "INTERSECTION IS: 2\n", + "Medial extrusion: 3.89 mm\n", + "Medial coverage: 39.9%\n", + "Lateral extrusion: -0.47 mm\n", + "Lateral coverage: 59.4%\n" + ] + } + ], + "source": [ + "path_segmentation = '../../data/SAG_3D_DESS_RIGHT_bones_cart_men_fib-label.nrrd'\n", + "\n", + "print('Initiating tibia mesh')\n", + "tibia = mskt.mesh.BoneMesh(\n", + " path_seg_image=path_segmentation,\n", + " label_idx=6,\n", + " dict_cartilage_labels={'medial': 2, 'lateral': 3} # Only dict needed!\n", + ")\n", + "\n", + "print('Creating tibia mesh')\n", + "tibia.create_mesh()\n", + "\n", + "print('Calculating cartilage thickness')\n", + "tibia.calc_cartilage_thickness()\n", + "print('Assigning cartilage regions')\n", + "tibia.assign_cartilage_regions()\n", + "\n", + "med_meniscus = mskt.mesh.Mesh(\n", + " path_seg_image=path_segmentation,\n", + " label_idx=10,\n", + ")\n", + "med_meniscus.create_mesh()\n", + "med_meniscus.consistent_faces()\n", + "\n", + "lat_meniscus = mskt.mesh.Mesh(\n", + " path_seg_image=path_segmentation,\n", + " label_idx=9,\n", + ")\n", + "lat_meniscus.create_mesh()\n", + "lat_meniscus.consistent_faces()\n", + "\n", + "# Clean API - no need to specify labels again!\n", + "tibia.set_menisci(\n", + " medial_meniscus=med_meniscus,\n", + " lateral_meniscus=lat_meniscus\n", + ")\n", + "\n", + "# Access metrics via properties - auto-computes on first access!\n", + "print(f\"Medial extrusion: {tibia.med_men_extrusion:.2f} mm\")\n", + "print(f\"Medial coverage: {tibia.med_men_coverage:.1f}%\")\n", + "print(f\"Lateral extrusion: {tibia.lat_men_extrusion:.2f} mm\")\n", + "print(f\"Lateral coverage: {tibia.lat_men_coverage:.1f}%\")\n", + "\n", + "\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "mskt", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.14" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From ee1e251d5deff656b660e936194d3ad951380d90 Mon Sep 17 00:00:00 2001 From: Anthony Gatt Date: Wed, 29 Oct 2025 23:40:27 -0700 Subject: [PATCH 04/16] autoformat & bump version. --- pymskt/__init__.py | 2 +- pymskt/mesh/__init__.py | 13 +- pymskt/mesh/meshTools.py | 9 +- pymskt/mesh/mesh_meniscus.py | 291 +++++++++--------- pymskt/mesh/meshes.py | 181 +++++------ .../meshMeniscus/MeniscusMesh_create_test.py | 5 +- .../compute_meniscal_coverage_test.py | 7 +- .../compute_meniscal_extrusion_test.py | 270 ++++++++-------- 8 files changed, 394 insertions(+), 384 deletions(-) diff --git a/pymskt/__init__.py b/pymskt/__init__.py index 8aea147..e34cd29 100644 --- a/pymskt/__init__.py +++ b/pymskt/__init__.py @@ -12,4 +12,4 @@ RTOL = 1e-4 ATOL = 1e-5 -__version__ = "0.1.18" +__version__ = "0.1.19" diff --git a/pymskt/mesh/__init__.py b/pymskt/mesh/__init__.py index b3c12a2..2e58f76 100644 --- a/pymskt/mesh/__init__.py +++ b/pymskt/mesh/__init__.py @@ -1,3 +1,12 @@ -from . import createMesh, io, meshCartilage, meshRegistration, meshTools, meshTransform, utils, mesh_meniscus -from .meshes import * +from . import ( + createMesh, + io, + mesh_meniscus, + meshCartilage, + meshRegistration, + meshTools, + meshTransform, + utils, +) from .mesh_meniscus import MeniscusMesh +from .meshes import * diff --git a/pymskt/mesh/meshTools.py b/pymskt/mesh/meshTools.py index db84690..f2237f1 100644 --- a/pymskt/mesh/meshTools.py +++ b/pymskt/mesh/meshTools.py @@ -471,7 +471,7 @@ def _get_distance_with_directions( Array of distances for each point """ distance_data = [] - + for idx in range(points.GetNumberOfPoints()): point = points.GetPoint(idx) direction = directions[idx] @@ -529,7 +529,9 @@ def get_distance_other_surface_at_points( point_normals = normals.GetOutput().GetPointData().GetNormals() # Extract normals as numpy array - directions = np.array([point_normals.GetTuple(idx) for idx in range(points.GetNumberOfPoints())]) + directions = np.array( + [point_normals.GetTuple(idx) for idx in range(points.GetNumberOfPoints())] + ) return _get_distance_with_directions( points, @@ -575,7 +577,7 @@ def get_distance_other_surface_at_points_along_unit_vector( """ points = surface.GetPoints() obb_other_surface = get_obb_surface(other_surface) - + unit_vector = np.asarray(unit_vector) assert np.isclose(np.linalg.norm(unit_vector), 1.0), "unit_vector must have magnitude 1.0" @@ -592,6 +594,7 @@ def get_distance_other_surface_at_points_along_unit_vector( no_distance_filler, ) + def set_mesh_physical_point_coords(mesh, new_points): """ Convenience function to update the x/y/z point coords of a mesh diff --git a/pymskt/mesh/mesh_meniscus.py b/pymskt/mesh/mesh_meniscus.py index 5024328..5f29cd8 100644 --- a/pymskt/mesh/mesh_meniscus.py +++ b/pymskt/mesh/mesh_meniscus.py @@ -11,6 +11,7 @@ """ import numpy as np + from pymskt.mesh.meshes import Mesh @@ -18,7 +19,7 @@ class MeniscusMesh(Mesh): """ Class to create, store, and process meniscus meshes with specialized analysis functions for meniscal extrusion and coverage calculations. - + Parameters ---------- mesh : vtk.vtkPolyData, optional @@ -33,12 +34,12 @@ class MeniscusMesh(Mesh): All islands smaller than this size are dropped, by default 1000 meniscus_type : str, optional Type of meniscus ('medial' or 'lateral'), by default None - + Attributes ---------- meniscus_type : str Type of meniscus ('medial' or 'lateral') - + Examples -------- >>> med_meniscus = MeniscusMesh( @@ -47,7 +48,7 @@ class MeniscusMesh(Mesh): ... meniscus_type='medial' ... ) """ - + def __init__( self, mesh=None, @@ -65,16 +66,16 @@ def __init__( min_n_pixels=min_n_pixels, ) self._meniscus_type = meniscus_type - + @property def meniscus_type(self): """Get the meniscus type.""" return self._meniscus_type - + @meniscus_type.setter def meniscus_type(self, new_meniscus_type): """Set the meniscus type with validation.""" - if new_meniscus_type not in [None, 'medial', 'lateral']: + if new_meniscus_type not in [None, "medial", "lateral"]: raise ValueError("meniscus_type must be None, 'medial', or 'lateral'") self._meniscus_type = new_meniscus_type @@ -83,20 +84,21 @@ def meniscus_type(self, new_meniscus_type): # Helper Functions # ============================================================================ + def compute_tibia_axes( tibia_mesh, medial_cart_label, lateral_cart_label, - scalar_array_name='labels', + scalar_array_name="labels", ): """ Compute anatomical axes (ML, IS, AP) from tibial cartilage regions. - + Uses PCA on combined cartilage points to find the tibial plateau normal (IS axis). The superior direction is determined by checking which side the bone is on relative to the cartilage. ML axis is from medial to lateral cartilage centers. AP axis is the cross product of ML and IS. - + Parameters ---------- tibia_mesh : BoneMesh or Mesh @@ -107,7 +109,7 @@ def compute_tibia_axes( Scalar value indicating lateral tibial cartilage region scalar_array_name : str, optional Name of scalar array containing region labels, by default 'labels' - + Returns ------- dict @@ -117,7 +119,7 @@ def compute_tibia_axes( - 'ap_axis': anterior-posterior axis vector (unit vector) - 'medial_center': medial cartilage center point - 'lateral_center': lateral cartilage center point - + Examples -------- >>> axes = compute_tibia_axes(tibia, med_cart_label=2, lat_cart_label=3) @@ -126,52 +128,52 @@ def compute_tibia_axes( """ # Get scalar array region_array = tibia_mesh[scalar_array_name] - + # Extract cartilage points - med_tib_cart_mask = (region_array == medial_cart_label) - lat_tib_cart_mask = (region_array == lateral_cart_label) - + med_tib_cart_mask = region_array == medial_cart_label + lat_tib_cart_mask = region_array == lateral_cart_label + med_tib_cart_points = tibia_mesh.points[med_tib_cart_mask] lat_tib_cart_points = tibia_mesh.points[lat_tib_cart_mask] tib_cart_points = np.concatenate([med_tib_cart_points, lat_tib_cart_points], axis=0) - + # Do PCA to get the three axes of the tib_cart_points and take the last # one as the inf/sup (normal to plateau) X = tib_cart_points - tib_cart_points.mean(axis=0, keepdims=True) # (N,3) # PCA via SVD: X = U S Vt, rows of Vt are PCs U, S, Vt = np.linalg.svd(X, full_matrices=False) pc1, pc2, pc3 = Vt # already orthonormal - + is_axis = pc3 - + # From the PCA we can't know what is up. Check which side the bone is on # relative to the cartilage. The opposite direction from bone to cartilage is IS. mean_tib = np.mean(tibia_mesh.points, axis=0) mean_cart = np.mean(tib_cart_points, axis=0) - + # Update is_axis direction based on where mean_tib is relative to mean_cart if np.dot(mean_tib - mean_cart, is_axis) > 0: is_axis = -is_axis - + # Compute ML axis from cartilage centers med_tib_center = np.mean(med_tib_cart_points, axis=0) lat_tib_center = np.mean(lat_tib_cart_points, axis=0) - + ml_axis = lat_tib_center - med_tib_center ml_axis = ml_axis / np.linalg.norm(ml_axis) - + # Compute AP axis as cross product # NOTE: AP axis direction is not always same (front vs back) # without inputting side (right.left). So, left it just as a general axis. ap_axis = np.cross(ml_axis, is_axis) ap_axis = ap_axis / np.linalg.norm(ap_axis) - + return { - 'ml_axis': ml_axis, - 'is_axis': is_axis, - 'ap_axis': ap_axis, - 'medial_center': med_tib_center, - 'lateral_center': lat_tib_center, + "ml_axis": ml_axis, + "is_axis": is_axis, + "ap_axis": ap_axis, + "medial_center": med_tib_center, + "lateral_center": lat_tib_center, } @@ -183,10 +185,10 @@ def _compute_extrusion_from_points( ): """ Compute extrusion by comparing ML extremes of cartilage and meniscus. - + Helper function that projects points onto ML axis and computes the signed extrusion distance. - + Parameters ---------- cart_points : np.ndarray @@ -197,7 +199,7 @@ def _compute_extrusion_from_points( Medial-lateral axis vector side : str 'med', 'medial', 'lat', or 'lateral' - + Returns ------- float @@ -205,18 +207,18 @@ def _compute_extrusion_from_points( """ cart_points_ml = np.dot(cart_points, ml_axis) men_points_ml = np.dot(men_points, ml_axis) - - if side in ['med', 'medial']: + + if side in ["med", "medial"]: cart_edge = np.min(cart_points_ml) men_edge = np.min(men_points_ml) extrusion = cart_edge - men_edge - elif side in ['lat', 'lateral']: + elif side in ["lat", "lateral"]: cart_edge = np.max(cart_points_ml) men_edge = np.max(men_points_ml) extrusion = men_edge - cart_edge else: - raise ValueError(f'Invalid side: {side}, must be one of: med, medial, lat, lateral') - + raise ValueError(f"Invalid side: {side}, must be one of: med, medial, lat, lateral") + return extrusion @@ -230,10 +232,10 @@ def _compute_middle_region_extrusion( ): """ Compute extrusion using middle percentile range along AP axis. - + This helper function focuses on the central portion of the AP range to avoid edge effects at the anterior and posterior extremes. - + Parameters ---------- cart_points : np.ndarray @@ -248,7 +250,7 @@ def _compute_middle_region_extrusion( 'med', 'medial', 'lat', or 'lateral' middle_percentile_range : float Fraction of AP range to use (centered on middle) - + Returns ------- float @@ -258,25 +260,27 @@ def _compute_middle_region_extrusion( cart_points_ap = np.dot(cart_points, ap_axis) min_cart_ap = np.min(cart_points_ap) max_cart_ap = np.max(cart_points_ap) - + # Get the middle +/- middle_percentile_range/2 of the cartilage along AP axis middle_ap_cartilage = (min_cart_ap + max_cart_ap) / 2 min_max_ap_cartilage_range = max_cart_ap - min_cart_ap plus_minus_ap_cartilage_range = min_max_ap_cartilage_range * middle_percentile_range / 2 lower_ap_cartilage = middle_ap_cartilage - plus_minus_ap_cartilage_range upper_ap_cartilage = middle_ap_cartilage + plus_minus_ap_cartilage_range - + # Get points within the middle AP range for cartilage - ap_cart_indices = (cart_points_ap >= lower_ap_cartilage) & (cart_points_ap <= upper_ap_cartilage) + ap_cart_indices = (cart_points_ap >= lower_ap_cartilage) & ( + cart_points_ap <= upper_ap_cartilage + ) ml_cart_points = cart_points[ap_cart_indices] - + # Project meniscus points onto AP axis men_points_ap = np.dot(men_points, ap_axis) - + # Get points within the middle AP range for meniscus ap_men_indices = (men_points_ap >= lower_ap_cartilage) & (men_points_ap <= upper_ap_cartilage) ml_men_points = men_points[ap_men_indices] - + # Compute extrusion extrusion = _compute_extrusion_from_points( cart_points=ml_cart_points, @@ -284,7 +288,7 @@ def _compute_middle_region_extrusion( ml_axis=ml_axis, side=side, ) - + return extrusion @@ -299,10 +303,10 @@ def _get_single_compartment_coverage( ): """ Compute meniscal coverage for a single compartment. - + Helper function that performs ray casting from tibia to meniscus and computes the area of cartilage covered by meniscus. - + Parameters ---------- tibia_mesh : BoneMesh or Mesh @@ -319,7 +323,7 @@ def _get_single_compartment_coverage( Name of scalar array containing region labels ray_cast_length : float, optional Length of rays to cast, by default 20.0 mm - + Returns ------- dict @@ -332,36 +336,36 @@ def _get_single_compartment_coverage( tibia_mesh.calc_distance_to_other_mesh( list_other_meshes=[meniscus_mesh], ray_cast_length=ray_cast_length, - name=f'{side_name}_men_dist_mm', + name=f"{side_name}_men_dist_mm", direction=is_direction, ) - + # Create binary masks - binary_mask_men_above = tibia_mesh[f'{side_name}_men_dist_mm'] > 0 + binary_mask_men_above = tibia_mesh[f"{side_name}_men_dist_mm"] > 0 binary_mask_cart = tibia_mesh[scalar_array_name] == cart_label - - tibia_mesh[f'{side_name}_men_above'] = binary_mask_men_above.astype(float) - tibia_mesh[f'{side_name}_cart'] = binary_mask_cart.astype(float) - + + tibia_mesh[f"{side_name}_men_above"] = binary_mask_men_above.astype(float) + tibia_mesh[f"{side_name}_cart"] = binary_mask_cart.astype(float) + # Extract cartilage submesh tibia_cart = tibia_mesh.copy() tibia_cart.remove_points(~binary_mask_cart, inplace=True) tibia_cart.clean(inplace=True) area_cart = tibia_cart.area - + # Extract covered cartilage submesh tibia_cart_men = tibia_cart.copy() - tibia_cart_men.remove_points(tibia_cart_men[f'{side_name}_men_above'] == 0, inplace=True) + tibia_cart_men.remove_points(tibia_cart_men[f"{side_name}_men_above"] == 0, inplace=True) tibia_cart_men.clean(inplace=True) area_cart_men = tibia_cart_men.area - + # Calculate coverage percentage percent_cart_men_coverage = (area_cart_men / area_cart) * 100 if area_cart > 0 else 0.0 - + return { - f'{side_name}_cart_men_coverage': percent_cart_men_coverage, - f'{side_name}_cart_men_area': area_cart_men, - f'{side_name}_cart_area': area_cart, + f"{side_name}_cart_men_coverage": percent_cart_men_coverage, + f"{side_name}_cart_men_area": area_cart_men, + f"{side_name}_cart_area": area_cart, } @@ -369,22 +373,23 @@ def _get_single_compartment_coverage( # Main Analysis Functions # ============================================================================ + def compute_meniscal_extrusion( tibia_mesh, medial_meniscus_mesh, lateral_meniscus_mesh, medial_cart_label, lateral_cart_label, - scalar_array_name='labels', + scalar_array_name="labels", middle_percentile_range=0.1, ): """ Compute meniscal extrusion for both medial and lateral menisci. - + Extrusion is computed by comparing the ML extremes of cartilage and meniscus within the middle portion of the AP range. This avoids edge effects at the anterior and posterior extremes. - + Parameters ---------- tibia_mesh : BoneMesh or Mesh @@ -401,7 +406,7 @@ def compute_meniscal_extrusion( Name of scalar array containing region labels, by default 'labels' middle_percentile_range : float, optional Fraction of AP range to use for extrusion measurement (centered), by default 0.1 - + Returns ------- dict @@ -411,48 +416,43 @@ def compute_meniscal_extrusion( - 'ml_axis': ML axis vector - 'ap_axis': AP axis vector - 'is_axis': IS axis vector - + Notes ----- Extrusion sign convention: positive values indicate meniscus extends beyond the cartilage rim. Negative values indicate the meniscus is contained within the cartilage boundaries. - + Examples -------- >>> results = compute_meniscal_extrusion( - ... tibia, med_meniscus, lat_meniscus, + ... tibia, med_meniscus, lat_meniscus, ... medial_cart_label=2, lateral_cart_label=3 ... ) >>> print(f"Medial extrusion: {results['medial_extrusion_mm']:.2f} mm") """ # Compute anatomical axes - axes = compute_tibia_axes( - tibia_mesh, - medial_cart_label, - lateral_cart_label, - scalar_array_name - ) - - ml_axis = axes['ml_axis'] - ap_axis = axes['ap_axis'] - is_axis = axes['is_axis'] - + axes = compute_tibia_axes(tibia_mesh, medial_cart_label, lateral_cart_label, scalar_array_name) + + ml_axis = axes["ml_axis"] + ap_axis = axes["ap_axis"] + is_axis = axes["is_axis"] + # Get cartilage points region_array = tibia_mesh[scalar_array_name] med_cart_indices = region_array == medial_cart_label lat_cart_indices = region_array == lateral_cart_label - + med_cart_points = tibia_mesh.points[med_cart_indices] lat_cart_points = tibia_mesh.points[lat_cart_indices] - + # Initialize results results = { - 'ml_axis': ml_axis, - 'ap_axis': ap_axis, - 'is_axis': is_axis, + "ml_axis": ml_axis, + "ap_axis": ap_axis, + "is_axis": is_axis, } - + # Compute medial extrusion (only if medial meniscus provided) if medial_meniscus_mesh is not None: med_men_points = medial_meniscus_mesh.points @@ -461,11 +461,11 @@ def compute_meniscal_extrusion( men_points=med_men_points, ap_axis=ap_axis, ml_axis=ml_axis, - side='med', + side="med", middle_percentile_range=middle_percentile_range, ) - results['medial_extrusion_mm'] = med_men_extrusion - + results["medial_extrusion_mm"] = med_men_extrusion + # Compute lateral extrusion (only if lateral meniscus provided) if lateral_meniscus_mesh is not None: lat_men_points = lateral_meniscus_mesh.points @@ -474,11 +474,11 @@ def compute_meniscal_extrusion( men_points=lat_men_points, ap_axis=ap_axis, ml_axis=ml_axis, - side='lat', + side="lat", middle_percentile_range=middle_percentile_range, ) - results['lateral_extrusion_mm'] = lat_men_extrusion - + results["lateral_extrusion_mm"] = lat_men_extrusion + return results @@ -488,16 +488,16 @@ def compute_meniscal_coverage( lateral_meniscus_mesh, medial_cart_label, lateral_cart_label, - scalar_array_name='labels', + scalar_array_name="labels", ray_cast_length=10.0, ): """ Compute meniscal coverage using superior-inferior ray casting. - + Coverage is computed by casting rays in the IS direction from tibial cartilage - reference points and checking for meniscus intersections. Areas are computed + reference points and checking for meniscus intersections. Areas are computed using PyVista's mesh area calculations. - + Parameters ---------- tibia_mesh : BoneMesh or Mesh @@ -514,7 +514,7 @@ def compute_meniscal_coverage( Name of scalar array containing region labels, by default 'labels' ray_cast_length : float, optional Length of rays to cast in IS direction, by default 20.0 mm - + Returns ------- dict @@ -525,27 +525,22 @@ def compute_meniscal_coverage( - 'lateral_covered_area_mm2': area of lateral cartilage covered (mm²) - 'medial_total_area_mm2': total medial cartilage area (mm²) - 'lateral_total_area_mm2': total lateral cartilage area (mm²) - + Examples -------- >>> results = compute_meniscal_coverage( - ... tibia, med_meniscus, lat_meniscus, + ... tibia, med_meniscus, lat_meniscus, ... medial_cart_label=2, lateral_cart_label=3 ... ) >>> print(f"Medial coverage: {results['medial_coverage_percent']:.1f}%") """ # Compute IS axis - axes = compute_tibia_axes( - tibia_mesh, - medial_cart_label, - lateral_cart_label, - scalar_array_name - ) - is_direction = axes['is_axis'] - + axes = compute_tibia_axes(tibia_mesh, medial_cart_label, lateral_cart_label, scalar_array_name) + is_direction = axes["is_axis"] + # Initialize results results = {} - + # Compute medial coverage (only if medial meniscus provided) if medial_meniscus_mesh is not None: med_coverage = _get_single_compartment_coverage( @@ -553,14 +548,14 @@ def compute_meniscal_coverage( meniscus_mesh=medial_meniscus_mesh, cart_label=medial_cart_label, is_direction=is_direction, - side_name='med', + side_name="med", scalar_array_name=scalar_array_name, ray_cast_length=ray_cast_length, ) - results['medial_coverage_percent'] = med_coverage['med_cart_men_coverage'] - results['medial_covered_area_mm2'] = med_coverage['med_cart_men_area'] - results['medial_total_area_mm2'] = med_coverage['med_cart_area'] - + results["medial_coverage_percent"] = med_coverage["med_cart_men_coverage"] + results["medial_covered_area_mm2"] = med_coverage["med_cart_men_area"] + results["medial_total_area_mm2"] = med_coverage["med_cart_area"] + # Compute lateral coverage (only if lateral meniscus provided) if lateral_meniscus_mesh is not None: lat_coverage = _get_single_compartment_coverage( @@ -568,14 +563,14 @@ def compute_meniscal_coverage( meniscus_mesh=lateral_meniscus_mesh, cart_label=lateral_cart_label, is_direction=is_direction, - side_name='lat', + side_name="lat", scalar_array_name=scalar_array_name, ray_cast_length=ray_cast_length, ) - results['lateral_coverage_percent'] = lat_coverage['lat_cart_men_coverage'] - results['lateral_covered_area_mm2'] = lat_coverage['lat_cart_men_area'] - results['lateral_total_area_mm2'] = lat_coverage['lat_cart_area'] - + results["lateral_coverage_percent"] = lat_coverage["lat_cart_men_coverage"] + results["lateral_covered_area_mm2"] = lat_coverage["lat_cart_men_area"] + results["lateral_total_area_mm2"] = lat_coverage["lat_cart_area"] + return results @@ -585,17 +580,17 @@ def analyze_meniscal_metrics( lateral_meniscus_mesh, medial_cart_label, lateral_cart_label, - scalar_array_name='labels', + scalar_array_name="labels", middle_percentile_range=0.1, ray_cast_length=10.0, ): """ Comprehensive meniscal analysis computing both extrusion and coverage metrics. - + This is the main function for complete meniscal analysis. It computes meniscal extrusion using the middle AP region and meniscal coverage using IS-direction ray casting. - + Parameters ---------- tibia_mesh : BoneMesh or Mesh @@ -614,37 +609,37 @@ def analyze_meniscal_metrics( Fraction of AP range to use for extrusion measurement, by default 0.1 ray_cast_length : float, optional Length of rays to cast for coverage analysis, by default 20.0 mm - + Returns ------- dict Dictionary containing all extrusion and coverage metrics: - + Extrusion metrics (mm, positive = extruded beyond cartilage rim): - 'medial_extrusion_mm': medial extrusion distance - 'lateral_extrusion_mm': lateral extrusion distance - + Coverage metrics: - 'medial_coverage_percent': percentage of medial cartilage covered - - 'lateral_coverage_percent': percentage of lateral cartilage covered + - 'lateral_coverage_percent': percentage of lateral cartilage covered - 'medial_covered_area_mm2': medial cartilage covered area (mm²) - 'lateral_covered_area_mm2': lateral cartilage covered area (mm²) - 'medial_total_area_mm2': total medial cartilage area (mm²) - 'lateral_total_area_mm2': total lateral cartilage area (mm²) - + Reference frame: - 'ml_axis': medial-lateral axis vector - 'ap_axis': anterior-posterior axis vector - 'is_axis': inferior-superior axis vector - + Notes ----- All meshes are automatically oriented with consistent normals before analysis. - + Examples -------- >>> results = analyze_meniscal_metrics( - ... tibia, med_meniscus, lat_meniscus, + ... tibia, med_meniscus, lat_meniscus, ... medial_cart_label=2, lateral_cart_label=3 ... ) >>> print(f"Medial extrusion: {results['medial_extrusion_mm']:.2f} mm") @@ -652,34 +647,40 @@ def analyze_meniscal_metrics( """ # Ensure tibia mesh is properly prepared tibia_mesh.compute_normals(auto_orient_normals=True, inplace=True) - + # Ensure meniscus meshes are properly prepared (only if not None) if medial_meniscus_mesh is not None: medial_meniscus_mesh.compute_normals(auto_orient_normals=True, inplace=True) if lateral_meniscus_mesh is not None: lateral_meniscus_mesh.compute_normals(auto_orient_normals=True, inplace=True) - + # Check that at least one meniscus is provided if medial_meniscus_mesh is None and lateral_meniscus_mesh is None: raise ValueError("At least one meniscus mesh must be provided") - + # Compute extrusion metrics (only for menisci that are present) extrusion_results = compute_meniscal_extrusion( - tibia_mesh, medial_meniscus_mesh, lateral_meniscus_mesh, - medial_cart_label, lateral_cart_label, scalar_array_name, middle_percentile_range + tibia_mesh, + medial_meniscus_mesh, + lateral_meniscus_mesh, + medial_cart_label, + lateral_cart_label, + scalar_array_name, + middle_percentile_range, ) - + # Compute coverage metrics (only for menisci that are present) coverage_results = compute_meniscal_coverage( - tibia_mesh, medial_meniscus_mesh, lateral_meniscus_mesh, - medial_cart_label, lateral_cart_label, scalar_array_name, ray_cast_length + tibia_mesh, + medial_meniscus_mesh, + lateral_meniscus_mesh, + medial_cart_label, + lateral_cart_label, + scalar_array_name, + ray_cast_length, ) - + # Combine results results = {**extrusion_results, **coverage_results} - - return results - - - + return results diff --git a/pymskt/mesh/meshes.py b/pymskt/mesh/meshes.py index 44d9af4..1b3e43f 100644 --- a/pymskt/mesh/meshes.py +++ b/pymskt/mesh/meshes.py @@ -804,7 +804,6 @@ def calc_distance_to_other_mesh( percent_ray_length_opposite_direction=0.25, name="thickness (mm)", direction=None, - ): """ Using bone mesh (`_mesh`) and the list of cartilage meshes (`list_cartilage_meshes`) @@ -863,7 +862,9 @@ def calc_distance_to_other_mesh( percent_ray_length_opposite_direction=percent_ray_length_opposite_direction, ) else: - raise ValueError(f"direction must be a numpy array, list, or tuple and received: {type(direction)}") + raise ValueError( + f"direction must be a numpy array, list, or tuple and received: {type(direction)}" + ) distances += node_data @@ -1391,7 +1392,9 @@ def copy(self, deep=True): copy_.bone = self.bone copy_.list_cartilage_meshes = self.list_cartilage_meshes copy_.list_cartilage_labels = self.list_cartilage_labels - copy_.dict_cartilage_labels = self.dict_cartilage_labels.copy() if self.dict_cartilage_labels else None + copy_.dict_cartilage_labels = ( + self.dict_cartilage_labels.copy() if self.dict_cartilage_labels else None + ) copy_.list_articular_surfaces = self.list_articular_surfaces copy_._meniscus_meshes = self._meniscus_meshes.copy() copy_._meniscal_outcomes = self._meniscal_outcomes @@ -1985,7 +1988,7 @@ def list_cartilage_labels(self): """ Convenience function to get the list of labels for cartilage tissues associated with this bone. - + If list_cartilage_labels was not explicitly set but dict_cartilage_labels was, this will return the values from dict_cartilage_labels in order. @@ -1997,11 +2000,11 @@ def list_cartilage_labels(self): # If explicit list provided, use it if self._list_cartilage_labels is not None: return self._list_cartilage_labels - + # Fall back to values from dict if available if self._dict_cartilage_labels is not None: return list(self._dict_cartilage_labels.values()) - + # Neither provided return None @@ -2032,7 +2035,7 @@ def list_cartilage_labels(self, new_list_cartilage_labels): def dict_cartilage_labels(self): """ Get the dictionary mapping cartilage region names to label values. - + Returns ------- dict or None @@ -2045,14 +2048,16 @@ def dict_cartilage_labels(self): def dict_cartilage_labels(self, new_dict_cartilage_labels): """ Set the dictionary mapping cartilage region names to label values. - + Parameters ---------- new_dict_cartilage_labels : dict or None Dictionary mapping region names to label values. For tibia: {'medial': 2, 'lateral': 3} """ - if new_dict_cartilage_labels is not None and not isinstance(new_dict_cartilage_labels, dict): + if new_dict_cartilage_labels is not None and not isinstance( + new_dict_cartilage_labels, dict + ): raise TypeError( f"dict_cartilage_labels must be a dict or None, got {type(new_dict_cartilage_labels)}" ) @@ -2117,30 +2122,30 @@ def bone(self, new_bone): if not isinstance(new_bone, str): raise TypeError(f"New bone provided is type {type(new_bone)} - expected `str`") self._bone = new_bone - + # ============================================================================ # Meniscus Analysis Methods (Tibia-specific) # NOTE: Could be refactored into TibiaMesh class inheriting from BoneMesh # ============================================================================ - + def set_menisci( self, medial_meniscus=None, medial_cart_label=None, lateral_meniscus=None, lateral_cart_label=None, - scalar_array_name='labels', + scalar_array_name="labels", ): """ Associate meniscus meshes and cartilage labels for meniscal analysis. - + This method stores references to meniscus meshes and their corresponding cartilage labels. You can set one or both menisci, but BOTH cartilage labels must be provided because tibial axes computation requires both cartilage regions. - + If dict_cartilage_labels was set during initialization, labels can be automatically inferred and don't need to be explicitly provided. - + Parameters ---------- medial_meniscus : MeniscusMesh or Mesh, optional @@ -2155,12 +2160,12 @@ def set_menisci( from dict_cartilage_labels['lateral'] if available. scalar_array_name : str, optional Name of scalar array containing region labels, by default 'labels' - + Raises ------ ValueError If no menisci are provided or if cartilage labels cannot be determined - + Examples -------- >>> # With dict_cartilage_labels set at initialization @@ -2172,7 +2177,7 @@ def set_menisci( ... medial_meniscus=med_men, ... lateral_meniscus=lat_men ... ) # Labels auto-inferred! - + >>> # Or provide labels explicitly (overrides dict_cartilage_labels) >>> tibia.set_menisci( ... medial_meniscus=med_men, medial_cart_label=2, @@ -2182,13 +2187,13 @@ def set_menisci( # Must provide at least one meniscus if medial_meniscus is None and lateral_meniscus is None: raise ValueError("At least one meniscus must be provided") - + # Try to get labels from dict_cartilage_labels if not explicitly provided if medial_cart_label is None and self._dict_cartilage_labels: - medial_cart_label = self._dict_cartilage_labels.get('medial') + medial_cart_label = self._dict_cartilage_labels.get("medial") if lateral_cart_label is None and self._dict_cartilage_labels: - lateral_cart_label = self._dict_cartilage_labels.get('lateral') - + lateral_cart_label = self._dict_cartilage_labels.get("lateral") + # Both cartilage labels are required for axes computation if medial_cart_label is None or lateral_cart_label is None: raise ValueError( @@ -2197,23 +2202,23 @@ def set_menisci( "one meniscus is being analyzed. Either provide them explicitly or set " "dict_cartilage_labels={'medial': X, 'lateral': Y} during initialization." ) - + # Store menisci if medial_meniscus is not None: - self._meniscus_meshes['medial'] = medial_meniscus + self._meniscus_meshes["medial"] = medial_meniscus if lateral_meniscus is not None: - self._meniscus_meshes['lateral'] = lateral_meniscus - + self._meniscus_meshes["lateral"] = lateral_meniscus + # Store labels (always store both since both are required) self._meniscal_cart_labels = { - 'medial': medial_cart_label, - 'lateral': lateral_cart_label, - 'scalar_array_name': scalar_array_name, + "medial": medial_cart_label, + "lateral": lateral_cart_label, + "scalar_array_name": scalar_array_name, } - + # Clear cached outcomes when menisci/labels are updated self._meniscal_outcomes = None - + def compute_meniscal_outcomes( self, medial_cart_label=None, @@ -2225,12 +2230,12 @@ def compute_meniscal_outcomes( ): """ Compute meniscal extrusion and coverage metrics. - + This method computes extrusion (how far meniscus extends beyond cartilage rim) and coverage (percentage of cartilage covered by meniscus) for menisci that have been set via set_menisci(). Can compute for one or both compartments depending on what was set. - + Parameters ---------- medial_cart_label : int or float, optional @@ -2248,7 +2253,7 @@ def compute_meniscal_outcomes( Length of rays to cast for coverage analysis, by default 20.0 mm force_recompute : bool, optional Force recomputation even if cached results exist, by default False - + Returns ------- dict @@ -2264,12 +2269,12 @@ def compute_meniscal_outcomes( - 'ml_axis': medial-lateral axis vector - 'ap_axis': anterior-posterior axis vector - 'is_axis': inferior-superior axis vector - + Raises ------ ValueError If no menisci are set or if required labels cannot be determined - + Examples -------- >>> # Set menisci with labels, then compute @@ -2279,7 +2284,7 @@ def compute_meniscal_outcomes( ... ) >>> results = tibia.compute_meniscal_outcomes() >>> print(f"Medial extrusion: {results['medial_extrusion_mm']:.2f} mm") - + >>> # Or provide labels explicitly >>> results = tibia.compute_meniscal_outcomes( ... medial_cart_label=2, lateral_cart_label=3 @@ -2288,83 +2293,83 @@ def compute_meniscal_outcomes( # Return cached results if available and not forcing recompute if self._meniscal_outcomes is not None and not force_recompute: return self._meniscal_outcomes - + # Check that at least one meniscus is set if not self._meniscus_meshes: raise ValueError( "No menisci have been set. Use set_menisci() to associate meniscus " "meshes and cartilage labels before computing outcomes." ) - + # Determine which menisci to compute for - has_medial = 'medial' in self._meniscus_meshes - has_lateral = 'lateral' in self._meniscus_meshes - + has_medial = "medial" in self._meniscus_meshes + has_lateral = "lateral" in self._meniscus_meshes + # Get labels (from parameters or cached values) # Both labels are ALWAYS required for axes computation if medial_cart_label is None: - if self._meniscal_cart_labels and 'medial' in self._meniscal_cart_labels: - medial_cart_label = self._meniscal_cart_labels['medial'] + if self._meniscal_cart_labels and "medial" in self._meniscal_cart_labels: + medial_cart_label = self._meniscal_cart_labels["medial"] else: raise ValueError( "medial_cart_label must be provided either in compute_meniscal_outcomes() " "or previously in set_menisci(). Both cartilage labels are required for " "tibial axes computation, even if only one meniscus is being analyzed." ) - + if lateral_cart_label is None: - if self._meniscal_cart_labels and 'lateral' in self._meniscal_cart_labels: - lateral_cart_label = self._meniscal_cart_labels['lateral'] + if self._meniscal_cart_labels and "lateral" in self._meniscal_cart_labels: + lateral_cart_label = self._meniscal_cart_labels["lateral"] else: raise ValueError( "lateral_cart_label must be provided either in compute_meniscal_outcomes() " "or previously in set_menisci(). Both cartilage labels are required for " "tibial axes computation, even if only one meniscus is being analyzed." ) - + # Get scalar array name if scalar_array_name is None: - if self._meniscal_cart_labels and 'scalar_array_name' in self._meniscal_cart_labels: - scalar_array_name = self._meniscal_cart_labels['scalar_array_name'] + if self._meniscal_cart_labels and "scalar_array_name" in self._meniscal_cart_labels: + scalar_array_name = self._meniscal_cart_labels["scalar_array_name"] else: - scalar_array_name = 'labels' - + scalar_array_name = "labels" + # Import analysis functions from pymskt.mesh.mesh_meniscus import analyze_meniscal_metrics - + # Always use the combined analysis function # It will handle single meniscus cases by only computing metrics for present menisci self._meniscal_outcomes = analyze_meniscal_metrics( tibia_mesh=self, - medial_meniscus_mesh=self._meniscus_meshes.get('medial'), - lateral_meniscus_mesh=self._meniscus_meshes.get('lateral'), + medial_meniscus_mesh=self._meniscus_meshes.get("medial"), + lateral_meniscus_mesh=self._meniscus_meshes.get("lateral"), medial_cart_label=medial_cart_label, lateral_cart_label=lateral_cart_label, scalar_array_name=scalar_array_name, middle_percentile_range=middle_percentile_range, ray_cast_length=ray_cast_length, ) - + return self._meniscal_outcomes - + @property def med_men_extrusion(self): """ Get medial meniscal extrusion value in mm. - + Automatically computes outcomes on first access if not already computed. Positive values indicate meniscus extends beyond cartilage rim. - + Returns ------- float Medial meniscal extrusion distance in mm - + Raises ------ ValueError If menisci haven't been set or if medial meniscus was not included - + Examples -------- >>> tibia.set_menisci(medial_meniscus=med_men, lateral_meniscus=lat_men) @@ -2379,29 +2384,29 @@ def med_men_extrusion(self): f"Cannot compute meniscal outcomes automatically: {str(e)}\n" "Ensure menisci are set via set_menisci() with appropriate labels." ) - - if 'medial_extrusion_mm' not in self._meniscal_outcomes: + + if "medial_extrusion_mm" not in self._meniscal_outcomes: raise ValueError("Medial meniscus was not included in the analysis") - return self._meniscal_outcomes['medial_extrusion_mm'] - + return self._meniscal_outcomes["medial_extrusion_mm"] + @property def lat_men_extrusion(self): """ Get lateral meniscal extrusion value in mm. - + Automatically computes outcomes on first access if not already computed. Positive values indicate meniscus extends beyond cartilage rim. - + Returns ------- float Lateral meniscal extrusion distance in mm - + Raises ------ ValueError If menisci haven't been set or if lateral meniscus was not included - + Examples -------- >>> tibia.set_menisci(medial_meniscus=med_men, lateral_meniscus=lat_men) @@ -2416,28 +2421,28 @@ def lat_men_extrusion(self): f"Cannot compute meniscal outcomes automatically: {str(e)}\n" "Ensure menisci are set via set_menisci() with appropriate labels." ) - - if 'lateral_extrusion_mm' not in self._meniscal_outcomes: + + if "lateral_extrusion_mm" not in self._meniscal_outcomes: raise ValueError("Lateral meniscus was not included in the analysis") - return self._meniscal_outcomes['lateral_extrusion_mm'] - + return self._meniscal_outcomes["lateral_extrusion_mm"] + @property def med_men_coverage(self): """ Get medial meniscal coverage percentage. - + Automatically computes outcomes on first access if not already computed. - + Returns ------- float Percentage of medial cartilage covered by meniscus - + Raises ------ ValueError If menisci haven't been set or if medial meniscus was not included - + Examples -------- >>> tibia.set_menisci(medial_meniscus=med_men, lateral_meniscus=lat_men) @@ -2452,28 +2457,28 @@ def med_men_coverage(self): f"Cannot compute meniscal outcomes automatically: {str(e)}\n" "Ensure menisci are set via set_menisci() with appropriate labels." ) - - if 'medial_coverage_percent' not in self._meniscal_outcomes: + + if "medial_coverage_percent" not in self._meniscal_outcomes: raise ValueError("Medial meniscus was not included in the analysis") - return self._meniscal_outcomes['medial_coverage_percent'] - + return self._meniscal_outcomes["medial_coverage_percent"] + @property def lat_men_coverage(self): """ Get lateral meniscal coverage percentage. - + Automatically computes outcomes on first access if not already computed. - + Returns ------- float Percentage of lateral cartilage covered by meniscus - + Raises ------ ValueError If menisci haven't been set or if lateral meniscus was not included - + Examples -------- >>> tibia.set_menisci(medial_meniscus=med_men, lateral_meniscus=lat_men) @@ -2488,7 +2493,7 @@ def lat_men_coverage(self): f"Cannot compute meniscal outcomes automatically: {str(e)}\n" "Ensure menisci are set via set_menisci() with appropriate labels." ) - - if 'lateral_coverage_percent' not in self._meniscal_outcomes: + + if "lateral_coverage_percent" not in self._meniscal_outcomes: raise ValueError("Lateral meniscus was not included in the analysis") - return self._meniscal_outcomes['lateral_coverage_percent'] + return self._meniscal_outcomes["lateral_coverage_percent"] diff --git a/testing/mesh/meshMeniscus/MeniscusMesh_create_test.py b/testing/mesh/meshMeniscus/MeniscusMesh_create_test.py index bbeb8a6..7143fda 100644 --- a/testing/mesh/meshMeniscus/MeniscusMesh_create_test.py +++ b/testing/mesh/meshMeniscus/MeniscusMesh_create_test.py @@ -11,6 +11,7 @@ import numpy as np import pytest + from pymskt.mesh import MeniscusMesh @@ -27,7 +28,3 @@ def test_meniscus_type_property(): def test_meniscus_type_validation(): """TODO: Test that invalid meniscus_type values raise ValueError.""" pass - - - - diff --git a/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py b/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py index 223e948..b8fa957 100644 --- a/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py +++ b/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py @@ -14,7 +14,8 @@ import numpy as np import pytest -from pymskt.mesh import Mesh, MeniscusMesh + +from pymskt.mesh import MeniscusMesh, Mesh from pymskt.mesh.mesh_meniscus import compute_meniscal_coverage_si_ray @@ -46,7 +47,3 @@ def test_coverage_area_calculation(): def test_coverage_si_tolerance(): """TODO: Test coverage with different SI tolerance values.""" pass - - - - diff --git a/testing/mesh/meshMeniscus/compute_meniscal_extrusion_test.py b/testing/mesh/meshMeniscus/compute_meniscal_extrusion_test.py index 0650f57..4af0cee 100644 --- a/testing/mesh/meshMeniscus/compute_meniscal_extrusion_test.py +++ b/testing/mesh/meshMeniscus/compute_meniscal_extrusion_test.py @@ -7,41 +7,39 @@ - Edge cases: empty compartments, missing meniscus data """ +import os + import numpy as np import pytest -import os -from pymskt.mesh import Mesh, BoneMesh, MeniscusMesh -from pymskt.mesh.mesh_meniscus import compute_meniscal_extrusion +from pymskt.mesh import BoneMesh, MeniscusMesh, Mesh +from pymskt.mesh.mesh_meniscus import compute_meniscal_extrusion # ============================================================================ # Fixtures for Meniscus Shift Tests # ============================================================================ + @pytest.fixture def tibia_with_menisci(): """ Fixture that provides a tibia mesh with cartilage regions and menisci. - + Returns a tuple of (tibia, medial_meniscus, lateral_meniscus, baseline_results) """ test_dir = os.path.dirname(os.path.abspath(__file__)) - data_dir = os.path.join(test_dir, '..', '..', '..', 'data') - path_segmentation = os.path.join(data_dir, 'SAG_3D_DESS_RIGHT_bones_cart_men_fib-label.nrrd') - + data_dir = os.path.join(test_dir, "..", "..", "..", "data") + path_segmentation = os.path.join(data_dir, "SAG_3D_DESS_RIGHT_bones_cart_men_fib-label.nrrd") + if not os.path.exists(path_segmentation): pytest.skip(f"Test data not found: {path_segmentation}") - + # Create tibia mesh with cartilage regions - tibia = BoneMesh( - path_seg_image=path_segmentation, - label_idx=6, - list_cartilage_labels=[2, 3] - ) + tibia = BoneMesh(path_seg_image=path_segmentation, label_idx=6, list_cartilage_labels=[2, 3]) tibia.create_mesh() tibia.calc_cartilage_thickness() tibia.assign_cartilage_regions() - + # Create meniscus meshes med_meniscus = Mesh( path_seg_image=path_segmentation, @@ -49,14 +47,14 @@ def tibia_with_menisci(): ) med_meniscus.create_mesh() med_meniscus.consistent_faces() - + lat_meniscus = Mesh( path_seg_image=path_segmentation, label_idx=9, ) lat_meniscus.create_mesh() lat_meniscus.consistent_faces() - + # Compute baseline extrusion baseline_results = compute_meniscal_extrusion( tibia_mesh=tibia, @@ -64,10 +62,10 @@ def tibia_with_menisci(): lateral_meniscus_mesh=lat_meniscus, medial_cart_label=2, lateral_cart_label=3, - scalar_array_name='labels', + scalar_array_name="labels", middle_percentile_range=0.1, ) - + return tibia, med_meniscus, lat_meniscus, baseline_results @@ -75,20 +73,21 @@ def tibia_with_menisci(): # Meniscus Shift Tests # ============================================================================ + def test_medial_meniscus_shift_medially_increases_extrusion(tibia_with_menisci): """ Test that shifting medial meniscus medially increases extrusion. - + When the medial meniscus is shifted 5mm medially (away from midline), the extrusion value should increase (more positive or less negative). """ tibia, med_meniscus, lat_meniscus, baseline_results = tibia_with_menisci - baseline_extrusion = baseline_results['medial_extrusion_mm'] - + baseline_extrusion = baseline_results["medial_extrusion_mm"] + # Shift medial meniscus medially (5mm in +X direction for right knee) med_meniscus_shifted = med_meniscus.copy() med_meniscus_shifted.points = med_meniscus_shifted.points + np.array([5.0, 0.0, 0.0]) - + # Compute extrusion with shifted meniscus results = compute_meniscal_extrusion( tibia_mesh=tibia, @@ -96,15 +95,16 @@ def test_medial_meniscus_shift_medially_increases_extrusion(tibia_with_menisci): lateral_meniscus_mesh=lat_meniscus, medial_cart_label=2, lateral_cart_label=3, - scalar_array_name='labels', + scalar_array_name="labels", middle_percentile_range=0.1, ) - + # Verify extrusion increased - assert results['medial_extrusion_mm'] > baseline_extrusion, \ - f"Medial shift should increase extrusion. " \ + assert results["medial_extrusion_mm"] > baseline_extrusion, ( + f"Medial shift should increase extrusion. " f"Baseline: {baseline_extrusion:.2f}, Shifted: {results['medial_extrusion_mm']:.2f}" - + ) + print(f"\n✓ Medial meniscus medial shift test passed!") print(f" Baseline: {baseline_extrusion:.2f} mm") print(f" After medial shift: {results['medial_extrusion_mm']:.2f} mm") @@ -114,17 +114,17 @@ def test_medial_meniscus_shift_medially_increases_extrusion(tibia_with_menisci): def test_medial_meniscus_shift_laterally_decreases_extrusion(tibia_with_menisci): """ Test that shifting medial meniscus laterally decreases extrusion. - + When the medial meniscus is shifted 5mm laterally (toward midline), the extrusion value should decrease (less positive or more negative). """ tibia, med_meniscus, lat_meniscus, baseline_results = tibia_with_menisci - baseline_extrusion = baseline_results['medial_extrusion_mm'] - + baseline_extrusion = baseline_results["medial_extrusion_mm"] + # Shift medial meniscus laterally (5mm in -X direction for right knee) med_meniscus_shifted = med_meniscus.copy() med_meniscus_shifted.points = med_meniscus_shifted.points - np.array([5.0, 0.0, 0.0]) - + # Compute extrusion with shifted meniscus results = compute_meniscal_extrusion( tibia_mesh=tibia, @@ -132,15 +132,16 @@ def test_medial_meniscus_shift_laterally_decreases_extrusion(tibia_with_menisci) lateral_meniscus_mesh=lat_meniscus, medial_cart_label=2, lateral_cart_label=3, - scalar_array_name='labels', + scalar_array_name="labels", middle_percentile_range=0.1, ) - + # Verify extrusion decreased - assert results['medial_extrusion_mm'] < baseline_extrusion, \ - f"Lateral shift should decrease extrusion. " \ + assert results["medial_extrusion_mm"] < baseline_extrusion, ( + f"Lateral shift should decrease extrusion. " f"Baseline: {baseline_extrusion:.2f}, Shifted: {results['medial_extrusion_mm']:.2f}" - + ) + print(f"\n✓ Medial meniscus lateral shift test passed!") print(f" Baseline: {baseline_extrusion:.2f} mm") print(f" After lateral shift: {results['medial_extrusion_mm']:.2f} mm") @@ -150,17 +151,17 @@ def test_medial_meniscus_shift_laterally_decreases_extrusion(tibia_with_menisci) def test_lateral_meniscus_shift_laterally_increases_extrusion(tibia_with_menisci): """ Test that shifting lateral meniscus laterally increases extrusion. - + When the lateral meniscus is shifted 5mm laterally (away from midline), the extrusion value should increase (more positive or less negative). """ tibia, med_meniscus, lat_meniscus, baseline_results = tibia_with_menisci - baseline_extrusion = baseline_results['lateral_extrusion_mm'] - + baseline_extrusion = baseline_results["lateral_extrusion_mm"] + # Shift lateral meniscus laterally (5mm in -X direction for right knee) lat_meniscus_shifted = lat_meniscus.copy() lat_meniscus_shifted.points = lat_meniscus_shifted.points - np.array([5.0, 0.0, 0.0]) - + # Compute extrusion with shifted meniscus results = compute_meniscal_extrusion( tibia_mesh=tibia, @@ -168,15 +169,16 @@ def test_lateral_meniscus_shift_laterally_increases_extrusion(tibia_with_menisci lateral_meniscus_mesh=lat_meniscus_shifted, medial_cart_label=2, lateral_cart_label=3, - scalar_array_name='labels', + scalar_array_name="labels", middle_percentile_range=0.1, ) - + # Verify extrusion increased - assert results['lateral_extrusion_mm'] > baseline_extrusion, \ - f"Lateral shift should increase extrusion. " \ + assert results["lateral_extrusion_mm"] > baseline_extrusion, ( + f"Lateral shift should increase extrusion. " f"Baseline: {baseline_extrusion:.2f}, Shifted: {results['lateral_extrusion_mm']:.2f}" - + ) + print(f"\n✓ Lateral meniscus lateral shift test passed!") print(f" Baseline: {baseline_extrusion:.2f} mm") print(f" After lateral shift: {results['lateral_extrusion_mm']:.2f} mm") @@ -186,17 +188,17 @@ def test_lateral_meniscus_shift_laterally_increases_extrusion(tibia_with_menisci def test_lateral_meniscus_shift_medially_decreases_extrusion(tibia_with_menisci): """ Test that shifting lateral meniscus medially decreases extrusion. - + When the lateral meniscus is shifted 5mm medially (toward midline), the extrusion value should decrease (less positive or more negative). """ tibia, med_meniscus, lat_meniscus, baseline_results = tibia_with_menisci - baseline_extrusion = baseline_results['lateral_extrusion_mm'] - + baseline_extrusion = baseline_results["lateral_extrusion_mm"] + # Shift lateral meniscus medially (5mm in +X direction for right knee) lat_meniscus_shifted = lat_meniscus.copy() lat_meniscus_shifted.points = lat_meniscus_shifted.points + np.array([5.0, 0.0, 0.0]) - + # Compute extrusion with shifted meniscus results = compute_meniscal_extrusion( tibia_mesh=tibia, @@ -204,15 +206,16 @@ def test_lateral_meniscus_shift_medially_decreases_extrusion(tibia_with_menisci) lateral_meniscus_mesh=lat_meniscus_shifted, medial_cart_label=2, lateral_cart_label=3, - scalar_array_name='labels', + scalar_array_name="labels", middle_percentile_range=0.1, ) - + # Verify extrusion decreased - assert results['lateral_extrusion_mm'] < baseline_extrusion, \ - f"Medial shift should decrease extrusion. " \ + assert results["lateral_extrusion_mm"] < baseline_extrusion, ( + f"Medial shift should decrease extrusion. " f"Baseline: {baseline_extrusion:.2f}, Shifted: {results['lateral_extrusion_mm']:.2f}" - + ) + print(f"\n✓ Lateral meniscus medial shift test passed!") print(f" Baseline: {baseline_extrusion:.2f} mm") print(f" After medial shift: {results['lateral_extrusion_mm']:.2f} mm") @@ -223,138 +226,136 @@ def test_lateral_meniscus_shift_medially_decreases_extrusion(tibia_with_menisci) # Convenience API Tests # ============================================================================ + def test_dict_cartilage_labels_replaces_list(): """ Test that dict_cartilage_labels can replace list_cartilage_labels. - + Verifies that cartilage thickness and region assignment work with only dict_cartilage_labels (no list_cartilage_labels needed). """ # Get path to test data test_dir = os.path.dirname(os.path.abspath(__file__)) - data_dir = os.path.join(test_dir, '..', '..', '..', 'data') - path_segmentation = os.path.join(data_dir, 'SAG_3D_DESS_RIGHT_bones_cart_men_fib-label.nrrd') - + data_dir = os.path.join(test_dir, "..", "..", "..", "data") + path_segmentation = os.path.join(data_dir, "SAG_3D_DESS_RIGHT_bones_cart_men_fib-label.nrrd") + if not os.path.exists(path_segmentation): pytest.skip(f"Test data not found: {path_segmentation}") - + # Create tibia with ONLY dict_cartilage_labels tibia = BoneMesh( path_seg_image=path_segmentation, label_idx=6, - dict_cartilage_labels={'medial': 2, 'lateral': 3} + dict_cartilage_labels={"medial": 2, "lateral": 3}, ) - + # Verify list_cartilage_labels property auto-extracts from dict assert tibia.list_cartilage_labels == [2, 3] - + # Verify standard operations work tibia.create_mesh() tibia.calc_cartilage_thickness() # Should work with dict values - tibia.assign_cartilage_regions() # Should work with dict values - + tibia.assign_cartilage_regions() # Should work with dict values + # Verify thickness and labels were assigned - assert 'thickness (mm)' in tibia.point_data - assert 'labels' in tibia.point_data - + assert "thickness (mm)" in tibia.point_data + assert "labels" in tibia.point_data + print("\n✓ dict_cartilage_labels successfully replaces list_cartilage_labels!") def test_set_menisci_auto_infers_labels(): """ Test that set_menisci() automatically infers labels from dict_cartilage_labels. - + Verifies that cartilage labels don't need to be specified explicitly when dict_cartilage_labels is set. """ test_dir = os.path.dirname(os.path.abspath(__file__)) - data_dir = os.path.join(test_dir, '..', '..', '..', 'data') - path_segmentation = os.path.join(data_dir, 'SAG_3D_DESS_RIGHT_bones_cart_men_fib-label.nrrd') - + data_dir = os.path.join(test_dir, "..", "..", "..", "data") + path_segmentation = os.path.join(data_dir, "SAG_3D_DESS_RIGHT_bones_cart_men_fib-label.nrrd") + if not os.path.exists(path_segmentation): pytest.skip(f"Test data not found: {path_segmentation}") - + # Setup tibia tibia = BoneMesh( path_seg_image=path_segmentation, label_idx=6, - dict_cartilage_labels={'medial': 2, 'lateral': 3} + dict_cartilage_labels={"medial": 2, "lateral": 3}, ) tibia.create_mesh() tibia.calc_cartilage_thickness() tibia.assign_cartilage_regions() - + # Create menisci med_meniscus = Mesh(path_seg_image=path_segmentation, label_idx=10) med_meniscus.create_mesh() med_meniscus.consistent_faces() - + lat_meniscus = Mesh(path_seg_image=path_segmentation, label_idx=9) lat_meniscus.create_mesh() lat_meniscus.consistent_faces() - + # Test: set_menisci WITHOUT explicit labels (should auto-infer from dict) - tibia.set_menisci( - medial_meniscus=med_meniscus, - lateral_meniscus=lat_meniscus - ) - + tibia.set_menisci(medial_meniscus=med_meniscus, lateral_meniscus=lat_meniscus) + # Verify labels were cached correctly assert tibia._meniscal_cart_labels is not None - assert tibia._meniscal_cart_labels['medial'] == 2 - assert tibia._meniscal_cart_labels['lateral'] == 3 - + assert tibia._meniscal_cart_labels["medial"] == 2 + assert tibia._meniscal_cart_labels["lateral"] == 3 + print("\n✓ set_menisci() successfully auto-infers labels from dict!") def test_meniscal_properties_lazy_evaluation(): """ Test that meniscal properties auto-compute on first access (lazy evaluation). - + Verifies that calling properties like med_men_extrusion automatically triggers computation without explicit compute_meniscal_outcomes() call. """ test_dir = os.path.dirname(os.path.abspath(__file__)) - data_dir = os.path.join(test_dir, '..', '..', '..', 'data') - path_segmentation = os.path.join(data_dir, 'SAG_3D_DESS_RIGHT_bones_cart_men_fib-label.nrrd') - + data_dir = os.path.join(test_dir, "..", "..", "..", "data") + path_segmentation = os.path.join(data_dir, "SAG_3D_DESS_RIGHT_bones_cart_men_fib-label.nrrd") + if not os.path.exists(path_segmentation): pytest.skip(f"Test data not found: {path_segmentation}") - + # Setup tibia = BoneMesh( path_seg_image=path_segmentation, label_idx=6, - dict_cartilage_labels={'medial': 2, 'lateral': 3} + dict_cartilage_labels={"medial": 2, "lateral": 3}, ) tibia.create_mesh() tibia.calc_cartilage_thickness() tibia.assign_cartilage_regions() - + med_meniscus = Mesh(path_seg_image=path_segmentation, label_idx=10) med_meniscus.create_mesh() med_meniscus.consistent_faces() - + lat_meniscus = Mesh(path_seg_image=path_segmentation, label_idx=9) lat_meniscus.create_mesh() lat_meniscus.consistent_faces() - + tibia.set_menisci(medial_meniscus=med_meniscus, lateral_meniscus=lat_meniscus) - + # Verify outcomes NOT computed yet assert tibia._meniscal_outcomes is None - + # Access property - should trigger auto-computation med_extrusion = tibia.med_men_extrusion - + # Verify outcomes NOW computed assert tibia._meniscal_outcomes is not None assert isinstance(med_extrusion, (int, float, np.number)) - + # Access another property - should use cached results (no recomputation) lat_extrusion = tibia.lat_men_extrusion assert isinstance(lat_extrusion, (int, float, np.number)) - + print("\n✓ Properties successfully auto-compute on first access!") print(f" Medial extrusion: {med_extrusion:.2f} mm") print(f" Lateral extrusion: {lat_extrusion:.2f} mm") @@ -363,111 +364,111 @@ def test_meniscal_properties_lazy_evaluation(): def test_meniscal_outcomes_caching(): """ Test that meniscal outcomes are properly cached and reused. - + Verifies that: - Results are cached after first computation - Property values match cached dictionary values - All expected metrics are present in cache """ test_dir = os.path.dirname(os.path.abspath(__file__)) - data_dir = os.path.join(test_dir, '..', '..', '..', 'data') - path_segmentation = os.path.join(data_dir, 'SAG_3D_DESS_RIGHT_bones_cart_men_fib-label.nrrd') - + data_dir = os.path.join(test_dir, "..", "..", "..", "data") + path_segmentation = os.path.join(data_dir, "SAG_3D_DESS_RIGHT_bones_cart_men_fib-label.nrrd") + if not os.path.exists(path_segmentation): pytest.skip(f"Test data not found: {path_segmentation}") - + # Setup tibia = BoneMesh( path_seg_image=path_segmentation, label_idx=6, - dict_cartilage_labels={'medial': 2, 'lateral': 3} + dict_cartilage_labels={"medial": 2, "lateral": 3}, ) tibia.create_mesh() tibia.calc_cartilage_thickness() tibia.assign_cartilage_regions() - + med_meniscus = Mesh(path_seg_image=path_segmentation, label_idx=10) med_meniscus.create_mesh() med_meniscus.consistent_faces() - + lat_meniscus = Mesh(path_seg_image=path_segmentation, label_idx=9) lat_meniscus.create_mesh() lat_meniscus.consistent_faces() - + tibia.set_menisci(medial_meniscus=med_meniscus, lateral_meniscus=lat_meniscus) - + # Trigger computation via property access med_extrusion = tibia.med_men_extrusion lat_extrusion = tibia.lat_men_extrusion med_coverage = tibia.med_men_coverage lat_coverage = tibia.lat_men_coverage - + # Verify all metrics are cached - assert 'medial_extrusion_mm' in tibia._meniscal_outcomes - assert 'lateral_extrusion_mm' in tibia._meniscal_outcomes - assert 'medial_coverage_percent' in tibia._meniscal_outcomes - assert 'lateral_coverage_percent' in tibia._meniscal_outcomes - + assert "medial_extrusion_mm" in tibia._meniscal_outcomes + assert "lateral_extrusion_mm" in tibia._meniscal_outcomes + assert "medial_coverage_percent" in tibia._meniscal_outcomes + assert "lateral_coverage_percent" in tibia._meniscal_outcomes + # Verify property values match cached values - assert med_extrusion == tibia._meniscal_outcomes['medial_extrusion_mm'] - assert lat_extrusion == tibia._meniscal_outcomes['lateral_extrusion_mm'] - assert med_coverage == tibia._meniscal_outcomes['medial_coverage_percent'] - assert lat_coverage == tibia._meniscal_outcomes['lateral_coverage_percent'] - + assert med_extrusion == tibia._meniscal_outcomes["medial_extrusion_mm"] + assert lat_extrusion == tibia._meniscal_outcomes["lateral_extrusion_mm"] + assert med_coverage == tibia._meniscal_outcomes["medial_coverage_percent"] + assert lat_coverage == tibia._meniscal_outcomes["lateral_coverage_percent"] + print("\n✓ Meniscal outcomes properly cached and accessible!") def test_meniscal_values_reasonable(): """ Test that computed meniscal values are reasonable. - + Verifies: - Extrusion values are numeric types - Coverage values are percentages (0-100) """ test_dir = os.path.dirname(os.path.abspath(__file__)) - data_dir = os.path.join(test_dir, '..', '..', '..', 'data') - path_segmentation = os.path.join(data_dir, 'SAG_3D_DESS_RIGHT_bones_cart_men_fib-label.nrrd') - + data_dir = os.path.join(test_dir, "..", "..", "..", "data") + path_segmentation = os.path.join(data_dir, "SAG_3D_DESS_RIGHT_bones_cart_men_fib-label.nrrd") + if not os.path.exists(path_segmentation): pytest.skip(f"Test data not found: {path_segmentation}") - + # Setup tibia = BoneMesh( path_seg_image=path_segmentation, label_idx=6, - dict_cartilage_labels={'medial': 2, 'lateral': 3} + dict_cartilage_labels={"medial": 2, "lateral": 3}, ) tibia.create_mesh() tibia.calc_cartilage_thickness() tibia.assign_cartilage_regions() - + med_meniscus = Mesh(path_seg_image=path_segmentation, label_idx=10) med_meniscus.create_mesh() med_meniscus.consistent_faces() - + lat_meniscus = Mesh(path_seg_image=path_segmentation, label_idx=9) lat_meniscus.create_mesh() lat_meniscus.consistent_faces() - + tibia.set_menisci(medial_meniscus=med_meniscus, lateral_meniscus=lat_meniscus) - + # Get values med_extrusion = tibia.med_men_extrusion lat_extrusion = tibia.lat_men_extrusion med_coverage = tibia.med_men_coverage lat_coverage = tibia.lat_men_coverage - + # Verify types (accept Python and numpy numeric types) assert isinstance(med_extrusion, (int, float, np.number)) assert isinstance(lat_extrusion, (int, float, np.number)) assert isinstance(med_coverage, (int, float, np.number)) assert isinstance(lat_coverage, (int, float, np.number)) - + # Verify coverage percentages are in valid range assert 0 <= med_coverage <= 100, f"Medial coverage {med_coverage}% outside valid range" assert 0 <= lat_coverage <= 100, f"Lateral coverage {lat_coverage}% outside valid range" - + print("\n✓ All meniscal values are reasonable!") print(f" Medial: {med_extrusion:.2f} mm extrusion, {med_coverage:.1f}% coverage") print(f" Lateral: {lat_extrusion:.2f} mm extrusion, {lat_coverage:.1f}% coverage") @@ -477,6 +478,7 @@ def test_meniscal_values_reasonable(): # TODO: Additional Tests # ============================================================================ + def test_extrusion_synthetic_data(): """TODO: Test extrusion calculation with synthetic tibia and meniscus meshes.""" pass @@ -485,7 +487,3 @@ def test_extrusion_synthetic_data(): def test_extrusion_no_extrusion(): """TODO: Test case where meniscus is fully within cartilage rim.""" pass - - - - From c4d839b1ce22c1b2b7c9f3c8594b27b4dd60c4c6 Mon Sep 17 00:00:00 2001 From: Anthony Gatt Date: Thu, 30 Oct 2025 10:57:20 -0700 Subject: [PATCH 05/16] Update Python version requirements and CI configurations - Bumped minimum required Python version from 3.7 to 3.9 in pyproject.toml. - Dropped support for Python 3.7 and 3.8 in build configurations due to end-of-life status and compatibility issues. - Updated GitHub Actions workflows to use macOS 14 and removed deprecated macOS 13, ensuring compatibility with Apple Silicon and improved testing across Python versions 3.9 to 3.12. --- .github/workflows/build-test.yml | 7 ++++--- .github/workflows/building-pypi-linux.yml | 12 +++++------- pyproject.toml | 5 +++-- 3 files changed, 12 insertions(+), 12 deletions(-) diff --git a/.github/workflows/build-test.yml b/.github/workflows/build-test.yml index 8495bd7..b11c2bc 100644 --- a/.github/workflows/build-test.yml +++ b/.github/workflows/build-test.yml @@ -40,9 +40,10 @@ jobs: strategy: fail-fast: false matrix: - os: [macos-13, ubuntu-latest] - #TODO: Fix macos-12 to be macos-latest. Issue with using latest (arm?) and point cloud utils install. . - python-version: ["3.8", "3.9", "3.10", "3.11"] + os: [macos-14, ubuntu-latest] + # Updated from macos-13 to macos-14 (Apple Silicon) due to macos-13 deprecation (Dec 2025) + # Dropped Python 3.8 support (EOL Oct 2024) due to ITK/NumPy compatibility issues + python-version: ["3.9", "3.10", "3.11", "3.12"] steps: - uses: actions/checkout@v4 diff --git a/.github/workflows/building-pypi-linux.yml b/.github/workflows/building-pypi-linux.yml index 57cf813..111c64b 100644 --- a/.github/workflows/building-pypi-linux.yml +++ b/.github/workflows/building-pypi-linux.yml @@ -57,13 +57,11 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - # macos-13 = intel mac runner - # macos-14 = apple silicon mac runner - os: [ubuntu-latest, windows-2022, macos-13] - # add macos-14 when apple silicon mac runners are available - # for all of the python versions. Seems like only 3.10 + is - # available now, but is intended to be fixed: - # https://github.com/actions/setup-python/issues/808 + # macos-14 = Apple Silicon (arm64) runner + # macos-14-large = Intel (x86_64) runner + # Updated from macos-13 due to deprecation (Dec 2025) + # Using macos-14 for arm64 wheels and macos-14-large for Intel wheels + os: [ubuntu-latest, windows-2022, macos-14, macos-14-large] steps: - uses: actions/checkout@v4 diff --git a/pyproject.toml b/pyproject.toml index b5c8865..39eb072 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -15,7 +15,7 @@ build-backend = "setuptools.build_meta" name = "mskt" description = "vtk helper tools/functions for musculoskeletal analyses" readme = "README.md" -requires-python = ">=3.7" +requires-python = ">=3.9" keywords = ["python"] license = {text = "MIT"} authors = [ @@ -99,7 +99,8 @@ exclude = ''' # Information needed for cibuildwheel [tool.cibuildwheel] # build options: https://cibuildwheel.readthedocs.io/en/stable/options/#build-selection -build = ["cp37-*", "cp38-*", "cp39-*", "cp310-*"] +# Dropped Python 3.7 and 3.8 support (EOL) due to dependency compatibility issues +build = ["cp39-*", "cp310-*", "cp311-*", "cp312-*"] skip = ["*-win32", "*i686", "*aarch64", "*ppc64le", "*s390x", "*musllinux*"] # testing info: https://cibuildwheel.readthedocs.io/en/stable/options/#testing From d9026a1b80c186beb799590a5c325f5c18b7dd35 Mon Sep 17 00:00:00 2001 From: Anthony Gatt Date: Thu, 30 Oct 2025 11:03:12 -0700 Subject: [PATCH 06/16] Updated function filename --- testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py b/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py index b8fa957..7bba84e 100644 --- a/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py +++ b/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py @@ -16,7 +16,7 @@ import pytest from pymskt.mesh import MeniscusMesh, Mesh -from pymskt.mesh.mesh_meniscus import compute_meniscal_coverage_si_ray +from pymskt.mesh.mesh_meniscus import compute_meniscal_coverage def test_coverage_synthetic_data(): From 1212f03df060dc2e36f56c96e76d81f624cfd2f6 Mon Sep 17 00:00:00 2001 From: Anthony Date: Thu, 30 Oct 2025 15:03:50 -0700 Subject: [PATCH 07/16] Update testing/mesh/meshMeniscus/compute_meniscal_extrusion_test.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- testing/mesh/meshMeniscus/compute_meniscal_extrusion_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testing/mesh/meshMeniscus/compute_meniscal_extrusion_test.py b/testing/mesh/meshMeniscus/compute_meniscal_extrusion_test.py index 4af0cee..97bdc48 100644 --- a/testing/mesh/meshMeniscus/compute_meniscal_extrusion_test.py +++ b/testing/mesh/meshMeniscus/compute_meniscal_extrusion_test.py @@ -12,7 +12,7 @@ import numpy as np import pytest -from pymskt.mesh import BoneMesh, MeniscusMesh, Mesh +from pymskt.mesh import BoneMesh, Mesh from pymskt.mesh.mesh_meniscus import compute_meniscal_extrusion # ============================================================================ From 31f8871891c3112eb9112fd5515fdf002bf34b86 Mon Sep 17 00:00:00 2001 From: Anthony Date: Thu, 30 Oct 2025 15:04:13 -0700 Subject: [PATCH 08/16] Update testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py | 1 - 1 file changed, 1 deletion(-) diff --git a/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py b/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py index 7bba84e..7b82e84 100644 --- a/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py +++ b/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py @@ -16,7 +16,6 @@ import pytest from pymskt.mesh import MeniscusMesh, Mesh -from pymskt.mesh.mesh_meniscus import compute_meniscal_coverage def test_coverage_synthetic_data(): From d12f252c2c9f6dd0bf62a3c2e280ab49e1754e79 Mon Sep 17 00:00:00 2001 From: Anthony Date: Thu, 30 Oct 2025 15:04:24 -0700 Subject: [PATCH 09/16] Update testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py | 1 - 1 file changed, 1 deletion(-) diff --git a/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py b/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py index 7b82e84..8108c35 100644 --- a/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py +++ b/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py @@ -15,7 +15,6 @@ import numpy as np import pytest -from pymskt.mesh import MeniscusMesh, Mesh def test_coverage_synthetic_data(): From a9aac2f88facc6843708eba312b6e5d49307d9b5 Mon Sep 17 00:00:00 2001 From: Anthony Date: Thu, 30 Oct 2025 15:05:06 -0700 Subject: [PATCH 10/16] Update pymskt/mesh/meshes.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- pymskt/mesh/meshes.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/pymskt/mesh/meshes.py b/pymskt/mesh/meshes.py index 1b3e43f..9488021 100644 --- a/pymskt/mesh/meshes.py +++ b/pymskt/mesh/meshes.py @@ -853,7 +853,10 @@ def calc_distance_to_other_mesh( elif isinstance(direction, (np.ndarray, list, tuple)): direction = np.array(direction) - direction = direction / np.linalg.norm(direction) + norm = np.linalg.norm(direction) + if norm == 0: + raise ValueError("direction vector must have non-zero magnitude for normalization.") + direction = direction / norm node_data = get_distance_other_surface_at_points_along_unit_vector( self, other_mesh, From f4ac5324f2e219601a4b4a783504049913e5ecb6 Mon Sep 17 00:00:00 2001 From: Anthony Date: Thu, 30 Oct 2025 15:05:52 -0700 Subject: [PATCH 11/16] Update pymskt/mesh/mesh_meniscus.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- pymskt/mesh/mesh_meniscus.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/pymskt/mesh/mesh_meniscus.py b/pymskt/mesh/mesh_meniscus.py index 5f29cd8..b26afb4 100644 --- a/pymskt/mesh/mesh_meniscus.py +++ b/pymskt/mesh/mesh_meniscus.py @@ -360,7 +360,12 @@ def _get_single_compartment_coverage( area_cart_men = tibia_cart_men.area # Calculate coverage percentage - percent_cart_men_coverage = (area_cart_men / area_cart) * 100 if area_cart > 0 else 0.0 + if area_cart == 0: + raise ValueError( + f"Cartilage region is empty (area = 0) for compartment '{side_name}'. " + "Cannot compute meniscal coverage. This likely indicates invalid input data." + ) + percent_cart_men_coverage = (area_cart_men / area_cart) * 100 return { f"{side_name}_cart_men_coverage": percent_cart_men_coverage, From 37c66d367f146a1fc21610d6f3a0f97dad6ad1f3 Mon Sep 17 00:00:00 2001 From: Anthony Date: Thu, 30 Oct 2025 15:06:15 -0700 Subject: [PATCH 12/16] Update pymskt/mesh/mesh_meniscus.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- pymskt/mesh/mesh_meniscus.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pymskt/mesh/mesh_meniscus.py b/pymskt/mesh/mesh_meniscus.py index b26afb4..b97152d 100644 --- a/pymskt/mesh/mesh_meniscus.py +++ b/pymskt/mesh/mesh_meniscus.py @@ -164,7 +164,7 @@ def compute_tibia_axes( # Compute AP axis as cross product # NOTE: AP axis direction is not always same (front vs back) - # without inputting side (right.left). So, left it just as a general axis. + # without inputting side (right/left). So, left it just as a general axis. ap_axis = np.cross(ml_axis, is_axis) ap_axis = ap_axis / np.linalg.norm(ap_axis) From 1a57dbd610f943d65ca103c769625ab25e2493ff Mon Sep 17 00:00:00 2001 From: Anthony Date: Thu, 30 Oct 2025 15:07:15 -0700 Subject: [PATCH 13/16] Update pymskt/mesh/meshes.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- pymskt/mesh/meshes.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/pymskt/mesh/meshes.py b/pymskt/mesh/meshes.py index 9488021..108f662 100644 --- a/pymskt/mesh/meshes.py +++ b/pymskt/mesh/meshes.py @@ -2305,9 +2305,6 @@ def compute_meniscal_outcomes( ) # Determine which menisci to compute for - has_medial = "medial" in self._meniscus_meshes - has_lateral = "lateral" in self._meniscus_meshes - # Get labels (from parameters or cached values) # Both labels are ALWAYS required for axes computation if medial_cart_label is None: From 33916e66781e0b377d28b7b761ac54f1ae2810ee Mon Sep 17 00:00:00 2001 From: Anthony Date: Thu, 30 Oct 2025 15:07:26 -0700 Subject: [PATCH 14/16] Update testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py | 1 - 1 file changed, 1 deletion(-) diff --git a/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py b/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py index 8108c35..c0926a7 100644 --- a/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py +++ b/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py @@ -13,7 +13,6 @@ """ import numpy as np -import pytest From 963cbd927211ba5f4b077fbdfb1a0406a2c46fbb Mon Sep 17 00:00:00 2001 From: Anthony Date: Thu, 30 Oct 2025 15:07:35 -0700 Subject: [PATCH 15/16] Update testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py | 1 - 1 file changed, 1 deletion(-) diff --git a/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py b/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py index c0926a7..a45ca3d 100644 --- a/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py +++ b/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py @@ -12,7 +12,6 @@ - Edge cases: empty compartments, missing meniscus data """ -import numpy as np From 16dc752fd650c38d82e7e815baf41594f0a997c0 Mon Sep 17 00:00:00 2001 From: Anthony Gatt Date: Fri, 31 Oct 2025 13:23:49 -0700 Subject: [PATCH 16/16] lint/autoformat --- pymskt/mesh/meshes.py | 4 +++- testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py | 2 -- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/pymskt/mesh/meshes.py b/pymskt/mesh/meshes.py index 108f662..5e456fd 100644 --- a/pymskt/mesh/meshes.py +++ b/pymskt/mesh/meshes.py @@ -855,7 +855,9 @@ def calc_distance_to_other_mesh( direction = np.array(direction) norm = np.linalg.norm(direction) if norm == 0: - raise ValueError("direction vector must have non-zero magnitude for normalization.") + raise ValueError( + "direction vector must have non-zero magnitude for normalization." + ) direction = direction / norm node_data = get_distance_other_surface_at_points_along_unit_vector( self, diff --git a/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py b/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py index a45ca3d..6088576 100644 --- a/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py +++ b/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py @@ -13,8 +13,6 @@ """ - - def test_coverage_synthetic_data(): """TODO: Test coverage calculation with synthetic tibia and meniscus meshes.""" pass