retuve.hip_us.metrics.coverage

  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: Coverage
 17"""
 18
 19import numpy as np
 20from networkx import diameter
 21from radstract.math import smart_find_intersection
 22
 23from retuve.classes.draw import Overlay
 24from retuve.classes.seg import SegObject
 25from retuve.hip_us.classes.general import HipDataUS, LandmarksUS
 26from retuve.keyphrases.config import Config
 27from retuve.keyphrases.enums import MetricUS
 28
 29FEM_HEAD_SCALE_FACTOR = 0.18
 30
 31
 32def find_cov_landmarks(
 33    femoral: SegObject, landmarks: LandmarksUS, config: Config
 34) -> LandmarksUS:
 35    """
 36    Algorithm to find the landmarks for the Coverage metric.
 37
 38    :param femoral: SegObject: The femoral SegObject.
 39    :param landmarks: LandmarksUS: The landmarks object.
 40    :param config: Config: The Config
 41
 42    :return: LandmarksUS: The updated landmarks object.
 43    """
 44    if not (
 45        femoral
 46        and any(m in config.hip.measurements for m in MetricUS.ALL())
 47        and femoral.points is not None
 48        and landmarks.apex is not None
 49    ):
 50        return landmarks
 51
 52    points = femoral.points
 53
 54    # Find the most right and left points
 55    most_right_point = max(points, key=lambda x: x[0])
 56    most_left_point = min(points, key=lambda x: x[0])
 57
 58    top_most_point = min(points, key=lambda x: x[1])
 59    bottom_most_point = max(points, key=lambda x: x[1])
 60
 61    # find diameter using right-left and top-bottom
 62    diameter_1 = abs(most_right_point[0] - most_left_point[0])
 63    diameter_2 = abs(bottom_most_point[1] - top_most_point[1])
 64
 65    if diameter_1 == 0 or diameter_2 == 0:
 66        return landmarks
 67
 68    # Reject frames with a non-circular femoral head
 69    if abs((diameter_1 - diameter_2) / diameter_1) > 0.35:
 70        return landmarks
 71
 72    diameter = diameter_2
 73
 74    # take half of the diameter
 75    radius = diameter / 2
 76
 77    center = (
 78        int(
 79            (abs(most_left_point[0] - most_right_point[0]) / 2)
 80            + most_left_point[0]
 81        ),
 82        int(top_most_point[1] + radius),
 83    )
 84
 85    # find the line from landmarks left to apex
 86    m = (landmarks.apex[1] - landmarks.left[1]) / (
 87        landmarks.apex[0] - landmarks.left[0]
 88    )
 89
 90    radius = radius + radius * config.hip.fem_extention
 91
 92    if m == 0:
 93        point_above = (center[0], center[1] + radius)
 94        point_below = (center[0], center[1] - radius)
 95    else:
 96        m_orth = -1 / m
 97        # equation of orthogonal line: y = m_orth * (x - center[0]) + center[1]
 98        # solve for x using the circle equation (x-center[0])**2 + (y-center[1])**2 = radius**2
 99
100        # break the radius into two parts, based on the
101        # sin and cos of the angle between the m_orth and the x-axis
102        angle = np.degrees(np.arctan(m_orth))
103
104        # break the radius into two parts, based on the
105        # sin and cos of the angle between the m_orth and the x-axis
106        radius_x = radius * np.cos(np.radians(angle))
107        radius_y = radius * np.sin(np.radians(angle))
108
109        point_above = (
110            int(center[0] + radius_x),
111            int(center[1] + radius_y),
112        )
113
114        point_below = (
115            int(center[0] - radius_x),
116            int(center[1] - radius_y),
117        )
118
119    # ensure that point_above and point_below are the right way round
120    # (Larger y represents down)
121    if point_above[1] > point_below[1]:
122        point_above, point_below = point_below, point_above
123
124    mid_cov_point = smart_find_intersection(
125        landmarks.apex,
126        landmarks.left,
127        point_above,
128        point_below,
129    )
130
131    landmarks.point_d = point_above
132    landmarks.point_D = point_below
133    landmarks.mid_cov_point = mid_cov_point
134
135    return landmarks
136
137
138def find_coverage(landmarks: LandmarksUS) -> float:
139    """
140    Calculate the Coverage metric.
141
142    :param landmarks: LandmarksUS: The landmarks object.
143
144    :return: float: The Coverage metric.
145    """
146    if not (
147        landmarks
148        and landmarks.mid_cov_point
149        and landmarks.point_D
150        and landmarks.point_d
151    ):
152        return 0
153
154    coverage = abs(landmarks.mid_cov_point[1] - landmarks.point_D[1]) / abs(
155        landmarks.point_D[1] - landmarks.point_d[1]
156    )
157
158    # if the mid_point is above the point_D, then the coverage is 0
159    if landmarks.mid_cov_point[1] > landmarks.point_D[1]:
160        coverage = 0
161
162    return round(coverage, 3)
163
164
165def draw_coverage(hip: HipDataUS, overlay: Overlay, config: Config) -> Overlay:
166    """
167    Draw the Coverage metric on the Overlay.
168
169    :param hip: HipDataUS: The Hip Data.
170    :param overlay: Overlay: The Overlay.
171    :param config: Config: The Config.
172
173    :return: Overlay: The updated Overlay.
174    """
175    coverage = hip.get_metric(MetricUS.COVERAGE)
176    if coverage != 0:
177        overlay.draw_lines(
178            ((hip.landmarks.point_D, hip.landmarks.point_d),),
179        )
180
181        xl, yl = hip.landmarks.mid_cov_point
182
183        if config.visuals.display_full_metric_names:
184            title = "Coverage"
185        else:
186            title = "cov"
187
188        overlay.draw_text(
189            f"{title}: {coverage:.2f}",
190            int(xl + 10),
191            int(yl - 30),
192            header="h2",
193        )
194
195    return overlay
196
197
198def bad_coverage(hip: HipDataUS) -> bool:
199    """
200    Check if the Coverage is bad.
201
202    :param hip: HipDataUS: The Hip Data.
203
204    :return: bool: True if the Coverage is bad.
205    """
206
207    if (
208        hip.get_metric(MetricUS.COVERAGE) < 0
209        or hip.get_metric(MetricUS.COVERAGE) > 1
210    ):
211        return True
212
213    return False