Slice thickness test xleyn

hazenlib/ACRObject.py

@@ -95,6 +95,40 @@
         # else:
         #     print("LR orientation swap not required.")
 
+    def calculate_rotation(img):
+        """Calculate the rotation angle of the phantom using a Gaussian blur, Canny edge detection and the Hough Transform.
+
+        Args: 
+            img (np.ndarray): pixel array of a DICOM object
+
+        Returns:
+            float: The rotation angle of the phantom wrt the positive x-axis in degrees.
+        """
+
+        img = (np.maximum(img,0) / img.max()) * 255
+        img = np.uint8(img)
+
+        imgBlur = cv.GaussianBlur(img, (13,13), 0)
+        canny_edge = cv.Canny(imgBlur, 5, 40)
+
+        test_angles = np.linspace(-pi/2, pi/2, 1801, endpoint = False)
+        hough_space, theta, rho = hough_line(canny_edge, test_angles)
+        _ , bestTheta, _ = hough_line_peaks(hough_space, theta, rho, min_angle = 10)
+        rotation = np.mean(bestTheta[:2]) * 180/pi 
+
+        # Transform to be angle wrt positive x axis.
+        rotation = 90 - (rotation + 180)
+
+        # Select positive angle that is less than 360 deg.
+        while rotation < 0:
+            rotation += 180
+
+        if rotation > 360:
+            rotation -= 360
+
+        # Angle in degrees from positive x-axis.
+        return rotation
+
     @staticmethod
     def determine_rotation(img):
         """Determine the rotation angle of the phantom using edge detection and the Hough transform.

hazenlib/tasks/acr_slice_thickness.py

@@ -17,6 +17,7 @@
 import sys
 import traceback
 import numpy as np
+import cv2
 
 import scipy
 import skimage.morphology
@@ -27,7 +28,6 @@
 from hazenlib.utils import get_image_orientation
 
 
-
 class ACRSliceThickness(HazenTask):
     """Slice width measurement class for DICOM images of the ACR phantom."""
 
@@ -340,3 +340,73 @@ def get_slice_thickness(self, dcm):
             self.report_files.append(img_path)
 
         return slice_thickness
+
+    @staticmethod
+    def offset_point(p1, p2, factor):
+        """Offsets a given point p1, by the vector between
+        points p2 and p1 divided by the factor parameter.
+
+        Args:
+            p1 (list): Point 1, [x, y]
+            p2 (list): Point 2, [x, y]
+            factor (int): Scaling factor for the vector offset.
+
+        Returns:
+            point (list): Point 1 after offset.
+        """
+        vector = [p2[0] - p1[0], p2[1] - p1[1]]
+        offset_vector = [x / factor for x in vector]
+        point = p1 + offset_vector
+
+        return point
+
+    def calculate_profiles(self, img):
+        """Calculates line profiles on image.
+        Works for a rotated phantom.
+
+        Args:
+            img (np.ndarray): Pixel array from DICOM image.
+
+        Returns:
+            profiles (list): A list of the two line profiles for ramps.
+        """
+        # Applying canny edge to uint8 representation of imgage and dilating.
+        img_uint8 = np.uint8(cv2.normalize(img, None, alpha=0, beta=255, norm_type=cv2.NORM_MINMAX))
+        img_uint8 = cv2.GaussianBlur(img_uint8, ksize=(15, 15), sigmaX=0, sigmaY=0)
+        canny = cv2.dilate(
+            cv2.Canny(img_uint8, threshold1=25, threshold2=50), cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (2, 2))
+        )
+
+        # Find Contours. Sort by horizontal span and select second in list (which will be central insert)
+        contours, _ = cv2.findContours(canny, mode=cv2.RETR_TREE, method=cv2.CHAIN_APPROX_NONE)
+        contours = sorted(contours, key=lambda cont: abs(np.max(cont[:, 0, 0]) - np.min(cont[:, 0, 0])), reverse=True)
+        rectCont = np.intp(cv2.boxPoints(cv2.minAreaRect(contours[1])))
+
+        # Offset points by 1/3 of distance to nearest point, towards that point
+        testPoint = rectCont[0]
+        _, closest, middle, furthest = sorted(rectCont, key=lambda x: np.linalg.norm(testPoint - x))
+        offset_points = [
+            self.offset_point(testPoint, closest, 3),
+            self.offset_point(closest, testPoint, 3),
+            self.offset_point(middle, furthest, 3),
+            self.offset_point(furthest, middle, 3),
+        ]
+
+        # Offset points by 1/8 of distance to line pair point, towards that point
+        testPoint = offset_points[0]
+        _, closest, middle, furthest = sorted(offset_points, key=lambda x: np.linalg.norm(testPoint - x))
+        offset_points = [
+            self.offset_point(testPoint, middle, 8),
+            self.offset_point(middle, testPoint, 8),
+            self.offset_point(closest, furthest, 8),
+            self.offset_point(furthest, closest, 8),
+        ]
+
+        # Determine which points to join to form the lines.
+        testPoint = offset_points[0]
+        _, closest, middle, furthest = sorted(offset_points, key=lambda x: np.linalg.norm(testPoint - x))
+
+        line1, line2 = [testPoint, middle], [closest, furthest]
+        profiles = [skimage.measure.profile_line(img, start[::-1], end[::-1]) for (start, end) in [line1, line2]]
+
+        return profiles