retuve.hip_us.metrics.alpha

  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: Alpha Angle
 17
 18"""
 19
 20import numpy as np
 21
 22from retuve.classes.draw import Overlay
 23from retuve.classes.seg import SegObject
 24from retuve.hip_us.classes.general import HipDataUS, LandmarksUS
 25from retuve.keyphrases.config import Config
 26from retuve.keyphrases.enums import MetricUS
 27from retuve.utils import find_midline_extremes
 28
 29# To move the apex point left or right
 30APEX_RIGHT_FACTOR = 0
 31
 32
 33def find_alpha_landmarks(
 34    illium: SegObject, landmarks: LandmarksUS, config: Config
 35) -> LandmarksUS:
 36    """
 37    Algorithm to find the landmarks for the Alpha Angle.
 38
 39    :param illium: SegObject: The Illium SegObject.
 40    :param landmarks: LandmarksUS: The landmarks object.
 41    :param config: Config: The Config
 42
 43    :return: LandmarksUS: The updated landmarks object.
 44    """
 45    if not (
 46        illium
 47        and any(m in config.hip.measurements for m in MetricUS.ALL())
 48        and illium.mask is not None
 49    ):
 50        return landmarks
 51
 52    left_most, right_most = find_midline_extremes(illium.midline_moved)
 53
 54    if (
 55        right_most is None
 56        or left_most is None
 57        or (right_most[1] - left_most[1]) == 0
 58    ):
 59        return landmarks
 60
 61    # get the equation for the line between the two extreme points
 62    # y = mx + b
 63    m = (right_most[0] - left_most[0]) / (right_most[1] - left_most[1])
 64    b = left_most[0] - m * left_most[1]
 65
 66    # if m is 0, return landmarks
 67    if m == 0:
 68        return landmarks
 69
 70    # find the height
 71
 72    # find the orthogonal line to the line between the two extreme points
 73    # y = -1/m * x + b
 74    m_orth = -1 / m
 75    # for each point along the line, find the distance between that point and the
 76
 77    points_on_line = np.array(
 78        [
 79            [x, m * x + b]
 80            for x in range(int(left_most[1]), int(right_most[1]), 1)
 81        ]
 82    )
 83
 84    # Convert white_points to a convenient shape for vector operations
 85    midline_moved = np.array(illium.midline_moved)[
 86        :, ::-1
 87    ]  # Reverse each point only once
 88
 89    # Create an array for b_orth values for each point in points_on_line
 90    b_orth_array = points_on_line[:, 1] - m_orth * points_on_line[:, 0]
 91
 92    # Calculate y_values on the orthogonal line for each point in points_on_line
 93    y_values_orth_line = (
 94        m_orth * midline_moved[:, 0] + b_orth_array[:, np.newaxis]
 95    )
 96
 97    # Check which white points are close to each orthogonal line
 98    close_points = np.isclose(
 99        y_values_orth_line, midline_moved[:, 1], atol=0.8
100    )
101
102    # Calculate distances for all point pairs
103    distances = np.linalg.norm(
104        points_on_line[:, np.newaxis, :] - midline_moved, axis=2
105    )
106
107    # Mask distances with close_points to consider only relevant distances
108    masked_distances = np.where(close_points, distances, 0)
109
110    # Find the maximum distance and the corresponding points
111    max_distance = np.max(masked_distances)
112    max_distance = np.max(masked_distances)
113    if max_distance > 0:
114        point_index, apex_point_index = np.unravel_index(
115            np.argmax(masked_distances), masked_distances.shape
116        )
117        best_point = tuple(points_on_line[point_index])
118        # Check APEX_RIGHT_FACTOR is within bounds
119        if apex_point_index + APEX_RIGHT_FACTOR < len(midline_moved):
120            factor = APEX_RIGHT_FACTOR
121        else:
122            # max possible factor
123            factor = len(midline_moved) - apex_point_index - 1
124
125        best_apex_point = tuple(midline_moved[apex_point_index + factor])
126
127        left_most, right_most = find_midline_extremes(illium.midline_moved)
128        if right_most is None or left_most is None:
129            return landmarks
130
131        left_most, right_most = tuple(reversed(left_most)), tuple(
132            reversed(right_most)
133        )
134        mid_x = (left_most[0] + right_most[0]) / 2
135        if (
136            (best_apex_point[1] < best_point[1])  # apex above line
137            and right_most[0] > best_apex_point[0]  # apex left of right point
138            and best_apex_point[0] > mid_x  # mid_x left of apex
139        ):
140            # re-calculate the left and right points from original midline
141
142            landmarks.left = left_most
143            landmarks.right = right_most
144            landmarks.apex = best_apex_point
145
146    # reversed because cv2 uses (y, x)??
147    return landmarks
148
149
150def find_alpha_angle(points: LandmarksUS) -> float:
151    """
152    Calculate the Alpha Angle.
153
154    :param points: LandmarksUS: The landmarks object.
155
156    :return: float: The Alpha Angle.
157    """
158    if not (points and points.left and points.apex and points.right):
159        return 0
160
161    # find angle ABC of points
162    A, B, C = (
163        np.array(points.left),
164        np.array(points.apex),
165        np.array(points.right),
166    )
167    AB = np.linalg.norm(A - B)
168    BC = np.linalg.norm(B - C)
169    AC = np.linalg.norm(A - C)
170    angle = np.arccos((BC**2 + AB**2 - AC**2) / (2 * BC * AB))
171
172    angle = np.degrees(angle)
173    angle = round((180 - angle), 1)
174
175    return round(angle, 2)
176
177
178def draw_alpha(hip: HipDataUS, overlay: Overlay, config: Config) -> Overlay:
179    """
180    Draw the Alpha Angle on the Overlay.
181
182    :param hip: HipDataUS: The Hip Data.
183    :param overlay: Overlay: The Overlay.
184    :param config: Config: The Config.
185
186    :return: Overlay: The updated Overlay.
187    """
188    alpha = hip.get_metric(MetricUS.ALPHA)
189    if alpha != 0:
190        m1 = (hip.landmarks.apex[1] - hip.landmarks.left[1]) / (
191            hip.landmarks.apex[0] - hip.landmarks.left[0]
192        )
193        m2 = (hip.landmarks.right[1] - hip.landmarks.apex[1]) / (
194            hip.landmarks.right[0] - hip.landmarks.apex[0]
195        )
196
197        b1 = hip.landmarks.apex[1] - m1 * hip.landmarks.apex[0]
198        b2 = hip.landmarks.apex[1] - m2 * hip.landmarks.apex[0]
199
200        # x2 to account for difference in gradient
201        new_right_apex = (
202            hip.landmarks.apex[0] + 350,
203            m1 * (hip.landmarks.apex[0] + 350) + b1,
204        )
205        new_right_landmark = (
206            hip.landmarks.right[0],
207            m2 * (hip.landmarks.right[0]) + b2,
208        )
209
210        overlay.draw_lines(
211            (
212                (hip.landmarks.left, new_right_apex),
213                (hip.landmarks.apex, new_right_landmark),
214            ),
215        )
216
217        if config.visuals.display_full_metric_names:
218            title = "Alpha Angle"
219        else:
220            title = "a"
221
222        overlay.draw_text(
223            f"{title}: {alpha}",
224            int(hip.landmarks.left[0]),
225            int(hip.landmarks.left[1] - 40),
226            header="h2",
227        )
228
229    return overlay
230
231
232def bad_alpha(hip: HipDataUS) -> bool:
233    """
234    Check if the Alpha Angle is bad.
235
236    :param hip: HipDataUS: The Hip Data.
237
238    :return: bool: True if the Alpha Angle is bad.
239    """
240
241    if (
242        hip.get_metric(MetricUS.ALPHA) < 20
243        or hip.get_metric(MetricUS.ALPHA) > 100
244    ):
245        return True
246
247    return False