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/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/__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 901095f..2e58f76 100644 --- a/pymskt/mesh/__init__.py +++ b/pymskt/mesh/__init__.py @@ -1,2 +1,12 @@ -from . import createMesh, io, meshCartilage, meshRegistration, meshTools, meshTransform, utils +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 f8eb790..f2237f1 100644 --- a/pymskt/mesh/meshTools.py +++ b/pymskt/mesh/meshTools.py @@ -439,86 +439,162 @@ 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( + surface, + other_surface, + 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 surface normals at each point. + + 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 + 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 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 + + Returns + ------- + np.ndarray + 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) + point_normals = normals.GetOutput().GetPointData().GetNormals() + + # Extract normals as numpy array + directions = np.array( + [point_normals.GetTuple(idx) for idx in range(points.GetNumberOfPoints())] + ) + + return _get_distance_with_directions( + points, + obb_other_surface, + directions, + ray_cast_length, + percent_ray_length_opposite_direction, + no_distance_filler, + ) + + +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. + + 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 + + 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): """ 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 new file mode 100644 index 0000000..b97152d --- /dev/null +++ b/pymskt/mesh/mesh_meniscus.py @@ -0,0 +1,691 @@ +""" +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 + 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, + 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..5e456fd 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,7 @@ 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 +843,33 @@ 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) + 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, + 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 +1320,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 +1349,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 +1361,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 +1397,13 @@ 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 +1518,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 +1598,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. @@ -1951,12 +1994,24 @@ 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 +2036,38 @@ 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 +2127,375 @@ 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 + # 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/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 diff --git a/testing/mesh/meshMeniscus/MeniscusMesh_create_test.py b/testing/mesh/meshMeniscus/MeniscusMesh_create_test.py new file mode 100644 index 0000000..7143fda --- /dev/null +++ b/testing/mesh/meshMeniscus/MeniscusMesh_create_test.py @@ -0,0 +1,30 @@ +""" +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..6088576 --- /dev/null +++ b/testing/mesh/meshMeniscus/compute_meniscal_coverage_test.py @@ -0,0 +1,43 @@ +""" +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 +""" + + +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..97bdc48 --- /dev/null +++ b/testing/mesh/meshMeniscus/compute_meniscal_extrusion_test.py @@ -0,0 +1,489 @@ +""" +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 os + +import numpy as np +import pytest + +from pymskt.mesh import BoneMesh, 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") + + 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 +}