retuve.hip_us.metrics.aca

  1# Copyright 2024 Adam McArthur
  2#
  3# Licensed under the Apache License, Version 2.0 (the "License");
  4# you may not use this file except in compliance with the License.
  5# You may obtain a copy of the License at
  6#
  7#     http://www.apache.org/licenses/LICENSE-2.0
  8#
  9# Unless required by applicable law or agreed to in writing, software
 10# distributed under the License is distributed on an "AS IS" BASIS,
 11# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 12# See the License for the specific language governing permissions and
 13# limitations under the License.
 14
 15"""
 16Metric: Acetabular Coverage Angle (ACA)
 17
 18Meant to be a more accurate version of alpha angle for 3D Ultrasounds.
 19
 20"""
 21
 22import copy
 23from enum import Enum
 24from typing import Dict, List, Tuple
 25
 26import numpy as np
 27import open3d as o3d
 28from numpy.typing import NDArray
 29
 30from retuve.hip_us.classes.enums import Side
 31from retuve.hip_us.classes.general import HipDataUS
 32from retuve.hip_us.typing import CoordinatesArray3D
 33from retuve.keyphrases.config import Config
 34from retuve.keyphrases.enums import ACASplit
 35
 36
 37class SideZ(Enum):
 38    """
 39    Enuum for side of the Illiac-Acetabular Bone view for each 2D Frame
 40    """
 41
 42    LEFT = "illiac wing"
 43    RIGHT = "acetabular roof"
 44
 45
 46"""
 47Matches the colors of the ACA to the side and third of the Illium.
 48"""
 49ACA_COLORS = {
 50    (SideZ.LEFT, Side.POST): "red",
 51    (SideZ.LEFT, Side.GRAF): "orange",
 52    (SideZ.LEFT, Side.ANT): "yellow",
 53    (SideZ.RIGHT, Side.POST): "blue",
 54    (SideZ.RIGHT, Side.GRAF): "cyan",
 55    (SideZ.RIGHT, Side.ANT): "lightgreen",
 56}
 57
 58
 59class Triangle:
 60    """
 61    Class to store the data of a triangle for the ACA calculation.
 62
 63    :attr apex_side: SideZ: The side of the Illium that the apex is on.
 64    :attr third: Side: The third of the Illium that the triangle is in.
 65    :attr normal: NDArray[np.float64]: The normal of the triangle.
 66    :attr centroid: NDArray[np.float64]: The centroid of the triangle.
 67    :attr color: str: The color of the triangle.
 68    """
 69
 70    def __init__(
 71        self,
 72        apex_side: SideZ,
 73        third: Side,
 74        normal: NDArray[np.float64],
 75        centroid: NDArray[np.float64],
 76    ):
 77        self.apex_side = apex_side
 78        self.third = third
 79        self.normal = normal
 80        self.centroid = centroid
 81        self.color = ACA_COLORS[(apex_side, third)]
 82
 83
 84def get_aca(
 85    illium_mesh: o3d.geometry.TriangleMesh,
 86    apex_points: NDArray[np.float64],
 87    graf_hip: HipDataUS,
 88    no_of_frames: int,
 89    config: Config,
 90) -> Tuple[
 91    Dict[Side, float],
 92    Tuple[CoordinatesArray3D, CoordinatesArray3D, str],
 93    List[Triangle],
 94    str,
 95]:
 96    """
 97    Calculate the acetabular coverage angle (ACA) for the hip.
 98
 99    :param illium_mesh (o3d.geometry.TriangleMesh): The Illium Mesh.
100    :param apex_points (NDArray[np.float64]): The apex points of the Illium.
101    :param graf_hip (HipDataUS): The Hip Data of the Graf Frame.
102    :param no_of_frames (int): The number of frames in the 3D Ultrasound.
103    :param config (Config): The Config object.
104
105    :return: Tuple of the ACA angles for each side, the average normals data,
106             the triangle data, and the recorded error.
107    """
108    z_gap = config.hip.z_gap * (200 / no_of_frames)
109
110    # add normals as arrows
111    illium_mesh.compute_vertex_normals()
112    normals = np.asarray(illium_mesh.vertex_normals)
113    vertices = np.asarray(illium_mesh.vertices)
114
115    # Compute triangle normals
116    triangle_normals = []
117    triangle_centroids = []
118    for tri in illium_mesh.triangles:
119        v0, v1, v2 = vertices[tri]
120        edge1 = v1 - v0
121        edge2 = v2 - v0
122        normal = np.cross(edge1, edge2)
123        centroid = (v0 + v1 + v2) / 3
124        triangle_normals.append(normal)
125        triangle_centroids.append(centroid)
126
127    triangle_normals = np.array(triangle_normals)
128    triangle_centroids = np.array(triangle_centroids)
129
130    # convert normals to unit vectors
131    triangle_normals /= np.linalg.norm(triangle_normals, axis=1)[:, None]
132
133    triangle_data = []
134
135    max_z = np.max(triangle_centroids[:, 2])
136    min_z = np.min(triangle_centroids[:, 2])
137
138    for normal, centroid in zip(triangle_normals, triangle_centroids):
139        if config.hip.aca_split == ACASplit.THIRDS:
140            third_z = (max_z - min_z) / 3
141
142            # find which third the centroid is in
143            if centroid[2] - min_z < third_z:
144                third = Side.POST
145            elif centroid[2] - min_z < third_z * 2:
146                third = Side.GRAF
147            else:
148                third = Side.ANT
149
150        elif config.hip.aca_split == ACASplit.GRAFS:
151            grafs_z = graf_hip.frame_no * z_gap
152
153            grafs_z = grafs_z
154
155            margin = (max_z - min_z) * 0.15
156
157            min_graf_z = grafs_z - margin
158
159            # if it is within 10% of the grafs_z, third = 2
160            if abs(centroid[2] - grafs_z) < margin:
161                third = Side.GRAF
162            elif centroid[2] < min_graf_z:
163                third = Side.POST
164            else:
165                third = Side.ANT
166
167        else:
168            raise ValueError("Invalid ACA Split Config Option")
169
170        apex_points = copy.deepcopy(np.array(apex_points))
171
172        # get the apex point with the closest z value
173        apex_point = apex_points[
174            np.argmin(np.abs(apex_points[:, 2] - centroid[2]))
175        ]
176
177        if centroid[0] < apex_point[0]:
178            apex_side = SideZ.LEFT
179        else:
180            apex_side = SideZ.RIGHT
181
182        triangle_data.append(
183            Triangle(
184                apex_side,
185                third,
186                normal,
187                centroid,
188            )
189        )
190
191    final_vectors = {
192        # (left, 1): vec
193    }
194
195    avg_normals_data = []
196
197    # for each color combination, find the average normal and add a cone for it
198    for apex_side in [SideZ.LEFT, SideZ.RIGHT]:
199        for third in Side.ALL():
200            normals = []
201            centroids = []
202
203            for triangle in triangle_data:
204                if triangle.apex_side == apex_side and triangle.third == third:
205                    normals.append(triangle.normal)
206                    centroids.append(triangle.centroid)
207
208            # This error is recorded in other places.
209            # And will happen when the grafs plane is detected
210            # At the edges of the measured region
211            if len(normals) == 0:
212                continue
213
214            normals = np.array(normals)
215            centroids = np.array(centroids)
216
217            avg_normal = np.mean(normals, axis=0)
218            avg_centroid = np.mean(centroids, axis=0)
219
220            avg_normals_data.append(
221                (
222                    avg_normal,
223                    avg_centroid,
224                    ACA_COLORS[(apex_side, third)],
225                )
226            )
227
228            final_vectors[apex_side, third] = avg_normal
229
230    aca_angles = {}
231
232    # for each third, find the angle between the apex left and right
233    for side in Side.ALL():
234        try:
235            apex_left = final_vectors[SideZ.LEFT, side]
236            apex_right = final_vectors[SideZ.RIGHT, side]
237            recorded_error = ""
238        except KeyError:
239            aca_angles[side] = 0
240            recorded_error = (
241                f"Not enough {Side.get_name(side)} to calculate ACA."
242            )
243            continue
244
245        # get the angle between the two apex points
246        angle = np.arccos(
247            np.dot(apex_left, apex_right)
248            / (np.linalg.norm(apex_left) * np.linalg.norm(apex_right))
249        )
250
251        # convert to degrees
252        angle = np.degrees(angle)
253
254        aca_angles[side] = round(angle, 2)
255
256    return aca_angles, avg_normals_data, triangle_data, recorded_error