Repository: GSTT-CSC/hazen

Issue:
acr ghosting module
problem centering ROIs. this datset is from bodycoil acquisition. @YassineAzma 
<img width="592" alt="Screenshot 2023-04-06 at 10 40 41" src="https://user-images.githubusercontent.com/93722461/231783117-9cce4c18-0783-4593-9596-5d989c6b728e.png">


Discussion:

Patch:
diff --git a/hazenlib/ACRObject.py b/hazenlib/ACRObject.py
new file mode 100644
index 0000000..6eb2828
--- /dev/null
+++ b/hazenlib/ACRObject.py
@@ -0,0 +1,248 @@
+import numpy as np
+import cv2
+import scipy
+import skimage
+
+
+class ACRObject:
+    def __init__(self, dcm):
+        self.dcm = dcm
+        self.images, self.dcm = self.sort_images()
+        self.rot_angle = self.determine_rotation()
+        self.centre, self.radius = self.find_phantom_center()
+
+    def sort_images(self):
+        """
+        Sort a stack of images based on slice position.
+
+        Returns
+        -------
+        img_stack : np.array
+            A sorted stack of images, where each image is represented as a 2D numpy array.
+        dcm_stack : pyd
+            A sorted stack of
+        """
+
+        z = np.array([dcm_file.ImagePositionPatient[2] for dcm_file in self.dcm])
+        dicom_stack = [self.dcm[i] for i in np.argsort(z)]
+        img_stack = [dicom.pixel_array for dicom in dicom_stack]
+
+        return img_stack, dicom_stack
+
+    def determine_rotation(self):
+        """
+        Determine the rotation angle of the phantom using edge detection and the Hough transform.
+
+        Returns
+        ------
+        rot_angle : float
+            The rotation angle in degrees.
+        """
+
+        thresh = cv2.threshold(self.images[0], 127, 255, cv2.THRESH_BINARY)[1]
+
+        kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (5, 5))
+        dilate = cv2.morphologyEx(thresh, cv2.MORPH_DILATE, kernel)
+        diff = cv2.absdiff(dilate, thresh)
+
+        h, theta, d = skimage.transform.hough_line(diff)
+        _, angles, _ = skimage.transform.hough_line_peaks(h, theta, d)
+
+        angle = np.rad2deg(scipy.stats.mode(angles)[0][0])
+        rot_angle = angle + 90 if angle < 0 else angle - 90
+
+        return rot_angle
+
+    def rotate_images(self):
+        """
+        Rotate the images by a specified angle. The value range and dimensions of the image are preserved.
+
+        Returns
+        -------
+        np.array:
+            The rotated images.
+        """
+
+        return skimage.transform.rotate(self.images, self.rot_angle, resize=False, preserve_range=True)
+
+    def mask_image(self, image, mag_threshold=0.05, open_threshold=500):
+        """
+        Mask an image by magnitude threshold before applying morphological opening to remove small unconnected
+        features. The convex hull is calculated in order to accommodate for potential air bubbles.
+
+        Returns
+        -------
+        np.array:
+            The masked image.
+        """
+        test_mask = self.circular_mask(self.centre, (80 // self.dcm[0].PixelSpacing[0]), image.shape)
+        test_image = image * test_mask
+        test_vals = test_image[np.nonzero(test_image)]
+        if np.percentile(test_vals, 80) - np.percentile(test_vals, 10) > 0.9 * np.max(image):
+            print('Large intensity variations detected in image. Using local thresholding!')
+            initial_mask = skimage.filters.threshold_sauvola(image, window_size=3, k=0.95)
+        else:
+            initial_mask = image > mag_threshold * np.max(image)
+
+        opened_mask = skimage.morphology.area_opening(initial_mask, area_threshold=open_threshold)
+        final_mask = skimage.morphology.convex_hull_image(opened_mask)
+
+        return final_mask
+
+    def find_phantom_center(self):
+        """
+        Find the center of the ACR phantom by filtering the uniformity slice and using the Hough circle detector.
+
+
+        Returns
+        -------
+        centre  : tuple
+            Tuple of ints representing the (x, y) center of the image.
+        """
+        img = self.images[6]
+        dx, dy = self.dcm[6].PixelSpacing
+        img_blur = cv2.GaussianBlur(img, (1, 1), 0)
+        img_grad = cv2.Sobel(img_blur, 0, dx=1, dy=1)
+
+        detected_circles = cv2.HoughCircles(img_grad, cv2.HOUGH_GRADIENT, 1,
+                                            minDist=int(180 / dy), param1=50, param2=30,
+                                            minRadius=int(180 / (2 * dy)), maxRadius=int(200 / (2 * dx))).flatten()
+
+        centre, radius = [int(i) for i in detected_circles[:2]], int(detected_circles[2])
+        return centre, radius
+
+    @staticmethod
+    def circular_mask(centre, radius, dims):
+        """
+        Sort a stack of images based on slice position.
+
+        Parameters
+        ----------
+        centre : tuple
+            The centre coordinates of the circular mask.
+        radius : int
+            The radius of the circular mask.
+        dims   : tuple
+            The dimensions of the circular mask.
+
+        Returns
+        -------
+        img_stack : np.array
+            A sorted stack of images, where each image is represented as a 2D numpy array.
+        """
+        # Define a circular logical mask
+        x = np.linspace(1, dims[0], dims[0])
+        y = np.linspace(1, dims[1], dims[1])
+
+        X, Y = np.meshgrid(x, y)
+        mask = (X - centre[0]) ** 2 + (Y - centre[1]) ** 2 <= radius ** 2
+
+        return mask
+
+    def measure_orthogonal_lengths(self, mask):
+        """
+        Compute the horizontal and vertical lengths of a mask, based on the centroid.
+
+        Parameters:
+        ----------
+        mask    : ndarray of bool
+            Boolean array of the image.
+
+        Returns:
+        ----------
+        length_dict : dict
+            A dictionary containing the following information for both horizontal and vertical line profiles:
+            'Horizontal Start'      | 'Vertical Start' : tuple of int
+                Horizontal/vertical starting point of the object.
+            'Horizontal End'        | 'Vertical End' : tuple of int
+                Horizontal/vertical ending point of the object.
+            'Horizontal Extent'     | 'Vertical Extent' : ndarray of int
+                Indices of the non-zero elements of the horizontal/vertical line profile.
+            'Horizontal Distance'   | 'Vertical Distance' : float
+                The horizontal/vertical length of the object.
+        """
+        dims = mask.shape
+        res = self.dcm[0].PixelSpacing
+
+        horizontal_start = (self.centre[1], 0)
+        horizontal_end = (self.centre[1], dims[0] - 1)
+        horizontal_line_profile = skimage.measure.profile_line(mask, horizontal_start, horizontal_end)
+        horizontal_extent = np.nonzero(horizontal_line_profile)[0]
+        horizontal_distance = (horizontal_extent[-1] - horizontal_extent[0]) * res[0]
+
+        vertical_start = (0, self.centre[0])
+        vertical_end = (dims[1] - 1, self.centre[0])
+        vertical_line_profile = skimage.measure.profile_line(mask, vertical_start, vertical_end)
+        vertical_extent = np.nonzero(vertical_line_profile)[0]
+        vertical_distance = (vertical_extent[-1] - vertical_extent[0]) * res[1]
+
+        length_dict = {
+            'Horizontal Start': horizontal_start,
+            'Horizontal End': horizontal_end,
+            'Horizontal Extent': horizontal_extent,
+            'Horizontal Distance': horizontal_distance,
+            'Vertical Start': vertical_start,
+            'Vertical End': vertical_end,
+            'Vertical Extent': vertical_extent,
+            'Vertical Distance': vertical_distance
+        }
+
+        return length_dict
+
+    @staticmethod
+    def rotate_point(origin, point, angle):
+        """
+        Compute the horizontal and vertical lengths of a mask, based on the centroid.
+
+        Parameters:
+        ----------
+        origin : tuple
+            The coordinates of the point around which the rotation is performed.
+        point  : tuple
+            The coordinates of the point to rotate.
+        angle  : int
+            Angle in degrees.
+
+        Returns:
+        ----------
+        x_prime : float
+            A float representing the x coordinate of the desired point after being rotated around an origin.
+        y_prime : float
+            A float representing the y coordinate of the desired point after being rotated around an origin.
+        """
+        theta = np.radians(angle)
+        c, s = np.cos(theta), np.sin(theta)
+
+        x_prime = origin[0] + c * (point[0] - origin[0]) - s * (point[1] - origin[1])
+        y_prime = origin[1] + s * (point[0] - origin[0]) + c * (point[1] - origin[1])
+        return x_prime, y_prime
+
+    @staticmethod
+    def find_n_highest_peaks(data, n, height=1):
+        """
+        Find the indices and amplitudes of the N highest peaks within a 1D array.
+
+        Parameters:
+        ----------
+        data    : np.array
+            The array containing the data to perform peak extraction on.
+        n       : int
+            The coordinates of the point to rotate.
+        height  : int or float
+            The amplitude threshold for peak identification.
+
+        Returns:
+        ----------
+        peak_locs       : np.array
+            A numpy array containing the indices of the N highest peaks identified.
+        peak_heights    : np.array
+            A numpy array containing the amplitudes of the N highest peaks identified.
+        """
+        peaks = scipy.signal.find_peaks(data, height)
+        pk_heights = peaks[1]['peak_heights']
+        pk_ind = peaks[0]
+
+        peak_heights = pk_heights[(-pk_heights).argsort()[:n]]  # find n highest peak amplitudes
+        peak_locs = pk_ind[(-pk_heights).argsort()[:n]]  # find n highest peak locations
+
+        return np.sort(peak_locs), np.sort(peak_heights)
diff --git a/hazenlib/tasks/acr_geometric_accuracy.py b/hazenlib/tasks/acr_geometric_accuracy.py
index 2f2ad0f..069200b 100644
--- a/hazenlib/tasks/acr_geometric_accuracy.py
+++ b/hazenlib/tasks/acr_geometric_accuracy.py
@@ -1,5 +1,6 @@
 """
 ACR Geometric Accuracy
+
 https://www.acraccreditation.org/-/media/acraccreditation/documents/mri/largephantomguidance.pdf
 
 Calculates geometric accuracy for slices 1 and 5 of the ACR phantom.
@@ -27,43 +28,40 @@ import skimage.measure
 import skimage.transform
 
 from hazenlib.HazenTask import HazenTask
+from hazenlib.ACRObject import ACRObject
 
 
 class ACRGeometricAccuracy(HazenTask):
 
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
+        self.ACR_obj = None
 
     def run(self) -> dict:
         results = {}
-        z = []
-        for dcm in self.data:
-            z.append(dcm.ImagePositionPatient[2])
-
-        idx_sort = np.argsort(z)
-
-        for dcm in self.data:
-            if dcm.ImagePositionPatient[2] == z[idx_sort[0]]:
-                try:
-                    result1 = self.get_geometric_accuracy_slice1(dcm)
-                except Exception as e:
-                    print(f"Could not calculate the geometric accuracy for {self.key(dcm)} because of : {e}")
-                    traceback.print_exc(file=sys.stdout)
-                    continue
-
-                results[self.key(dcm)] = result1
-            elif dcm.ImagePositionPatient[2] == z[idx_sort[4]]:
-                try:
-                    result5 = self.get_geometric_accuracy_slice5(dcm)
-                except Exception as e:
-                    print(f"Could not calculate the geometric accuracy for {self.key(dcm)} because of : {e}")
-                    traceback.print_exc(file=sys.stdout)
-                    continue
-
-                results[self.key(dcm)] = result5
-
-        L = result1 + result5
+        self.ACR_obj = ACRObject(self.data)
+        slice1_dcm = self.ACR_obj.dcm[0]
+        slice5_dcm = self.ACR_obj.dcm[4]
+
+        try:
+            lengths_1 = self.get_geometric_accuracy_slice1(slice1_dcm)
+            results[self.key(slice1_dcm)] = lengths_1
+        except Exception as e:
+            print(f"Could not calculate the geometric accuracy for {self.key(slice1_dcm)} because of : {e}")
+            traceback.print_exc(file=sys.stdout)
+
+        try:
+            lengths_5 = self.get_geometric_accuracy_slice5(slice5_dcm)
+            results[self.key(slice5_dcm)] = lengths_5
+        except Exception as e:
+            print(f"Could not calculate the geometric accuracy for {self.key(slice5_dcm)} because of : {e}")
+            traceback.print_exc(file=sys.stdout)
+
+
+        L = lengths_1 + lengths_5
+
         mean_err, max_err, cov_l = self.distortion_metric(L)
+
         print(f"Mean relative measurement error is equal to {np.round(mean_err, 2)}mm")
         print(f"Maximum absolute measurement error is equal to {np.round(max_err, 2)}mm")
         print(f"Coefficient of variation of measurements is equal to {np.round(cov_l, 2)}%")
@@ -74,72 +72,18 @@ class ACRGeometricAccuracy(HazenTask):
 
         return results
 
-    def centroid_com(self, dcm):
-        # Calculate centroid of object using a centre-of-mass calculation
-        thresh_img = dcm > 0.25 * np.max(dcm)
-        open_img = skimage.morphology.area_opening(thresh_img, area_threshold=500)
-        bhull = skimage.morphology.convex_hull_image(open_img)
-        coords = np.nonzero(bhull)  # row major - first array is columns
-
-        sum_x = np.sum(coords[1])
-        sum_y = np.sum(coords[0])
-        cx, cy = sum_x / coords[0].shape[0], sum_y / coords[1].shape[0]
-        cxy = (round(cx), round(cy))
-
-        return bhull, cxy
-
-    def horizontal_length(self, res, mask, cxy):
-        dims = mask.shape
-        start_h = (cxy[1], 0)
-        end_h = (cxy[1], dims[0] - 1)
-        line_profile_h = skimage.measure.profile_line(mask, start_h, end_h, mode='reflect')
-        extent_h = np.nonzero(line_profile_h)[0]
-        dist_h = (extent_h[-1] - extent_h[0]) * res[0]
-
-        h_dict = {
-            'Start': start_h,
-            'End': end_h,
-            'Extent': extent_h,
-            'Distance': dist_h
-        }
-        return h_dict
-
-    def vertical_length(self, res, mask, cxy):
-        dims = mask.shape
-        start_v = (0, cxy[0])
-        end_v = (dims[1] - 1, cxy[0])
-        line_profile_v = skimage.measure.profile_line(mask, start_v, end_v, mode='reflect')
-        extent_v = np.nonzero(line_profile_v)[0]
-        dist_v = (extent_v[-1] - extent_v[0]) * res[1]
-
-        v_dict = {
-            'Start': start_v,
-            'End': end_v,
-            'Extent': extent_v,
-            'Distance': dist_v
-        }
-        return v_dict
-
-    def rotate_point(self, origin, point, angle):
-        theta = np.radians(angle)
-        c, s = np.cos(theta), np.sin(theta)
-
-        x_prime = origin[0] + c * (point[0] - origin[0]) - s * (point[1] - origin[1])
-        y_prime = origin[1] + s * (point[0] - origin[0]) + c * (point[1] - origin[1])
-        return x_prime, y_prime
-
-    def diagonal_lengths(self, res, mask, cxy):
+    def diagonal_lengths(self, res, img, cxy):
         eff_res = np.sqrt(np.mean(np.square(res)))
-        mask_rotate = skimage.transform.rotate(mask, 45, center=(cxy[0], cxy[1]))
+        img_rotate = skimage.transform.rotate(img, 45, center=(cxy[0], cxy[1]))
 
-        h_dict = self.horizontal_length(res, mask_rotate, cxy)
-        extent_h = h_dict['Extent']
+        length_dict = self.ACR_obj.measure_orthogonal_lengths(img_rotate)
+        extent_h = length_dict['Horizontal Extent']
 
         origin = (cxy[0], cxy[1])
         start = (extent_h[0], cxy[1])
         end = (extent_h[-1], cxy[1])
-        se_x_start, se_y_start = self.rotate_point(origin, start, 45)
-        se_x_end, se_y_end = self.rotate_point(origin, end, 45)
+        se_x_start, se_y_start = ACRObject.rotate_point(origin, start, 45)
+        se_x_end, se_y_end = ACRObject.rotate_point(origin, end, 45)
 
         dist_se = np.sqrt(np.sum(np.square([se_x_end - se_x_start, se_y_end - se_y_start]))) * eff_res
         se_dict = {
@@ -149,13 +93,12 @@ class ACRGeometricAccuracy(HazenTask):
             'Distance': dist_se
         }
 
-        v_dict = self.vertical_length(res, mask_rotate, cxy)
-        extent_v = v_dict['Extent']
+        extent_v = length_dict['Vertical Extent']
 
         start = (cxy[0], extent_v[0])
         end = (cxy[0], extent_v[-1])
-        sw_x_start, sw_y_start = self.rotate_point(origin, start, 45)
-        sw_x_end, sw_y_end = self.rotate_point(origin, end, 45)
+        sw_x_start, sw_y_start = ACRObject.rotate_point(origin, start, 45)
+        sw_x_end, sw_y_end = ACRObject.rotate_point(origin, end, 45)
 
         dist_sw = np.sqrt(np.sum(np.square([sw_x_end - sw_x_start, sw_y_end - sw_y_start]))) * eff_res
         sw_dict = {
@@ -169,72 +112,92 @@ class ACRGeometricAccuracy(HazenTask):
 
     def get_geometric_accuracy_slice1(self, dcm):
         img = dcm.pixel_array
-        res = dcm.PixelSpacing
-        mask, cxy = self.centroid_com(img)
 
-        h_dict = self.horizontal_length(res, mask, cxy)
-        v_dict = self.vertical_length(res, mask, cxy)
+        mask = self.ACR_obj.mask_image(self.ACR_obj.images[6])
+        cxy = self.ACR_obj.centre
+        length_dict = self.ACR_obj.measure_orthogonal_lengths(mask)
 
         if self.report:
             import matplotlib.pyplot as plt
-            fig = plt.figure()
-            fig.set_size_inches(8, 8)
-            plt.imshow(img)
-
-            plt.arrow(h_dict['Extent'][0], cxy[1], h_dict['Extent'][-1] - h_dict['Extent'][0], 1, color='blue',
-                      length_includes_head=True, head_width=5)
-            plt.arrow(cxy[0], v_dict['Extent'][0], 1, v_dict['Extent'][-1] - v_dict['Extent'][0], color='orange',
-                      length_includes_head=True, head_width=5)
-            plt.legend([str(np.round(h_dict['Distance'], 2)) + 'mm',
-                        str(np.round(v_dict['Distance'], 2)) + 'mm'])
-            plt.axis('off')
-            plt.title('Geometric Accuracy for Slice 1')
+            fig, axes = plt.subplots(3, 1)
+            fig.set_size_inches(8, 24)
+            fig.tight_layout(pad=4)
+
+            axes[0].imshow(img)
+            axes[0].scatter(cxy[0], cxy[1], c='red')
+            axes[0].set_title('Centroid Location')
+
+            axes[1].imshow(mask)
+            axes[1].set_title('Thresholding Result')
+
+            axes[2].imshow(img)
+            axes[2].arrow(length_dict['Horizontal Extent'][0], cxy[1],
+                          length_dict['Horizontal Extent'][-1] - length_dict['Horizontal Extent'][0], 1, color='blue',
+                          length_includes_head=True, head_width=5)
+            axes[2].arrow(cxy[0], length_dict['Vertical Extent'][0], 1, length_dict['Vertical Extent'][-1] -
+                          length_dict['Vertical Extent'][0], color='orange', length_includes_head=True, head_width=5)
+            axes[2].legend([str(np.round(length_dict['Horizontal Distance'], 2)) + 'mm',
+                            str(np.round(length_dict['Vertical Distance'], 2)) + 'mm'])
+            axes[2].axis('off')
+            axes[2].set_title('Geometric Accuracy for Slice 1')
 
             img_path = os.path.realpath(os.path.join(self.report_path, f'{self.key(dcm)}.png'))
             fig.savefig(img_path)
             self.report_files.append(img_path)
 
-        return h_dict['Distance'], v_dict['Distance']
+        return length_dict['Horizontal Distance'], length_dict['Vertical Distance']
 
     def get_geometric_accuracy_slice5(self, dcm):
         img = dcm.pixel_array
         res = dcm.PixelSpacing
-        mask, cxy = self.centroid_com(img)
+        mask = self.ACR_obj.mask_image(self.ACR_obj.images[6])
+        cxy = self.ACR_obj.centre
 
-        h_dict = self.horizontal_length(res, mask, cxy)
-        v_dict = self.vertical_length(res, mask, cxy)
+        length_dict = self.ACR_obj.measure_orthogonal_lengths(mask)
         sw_dict, se_dict = self.diagonal_lengths(res, mask, cxy)
 
         if self.report:
             import matplotlib.pyplot as plt
-            fig = plt.figure()
-            fig.set_size_inches(8, 8)
-            plt.imshow(img)
-
-            plt.arrow(h_dict['Extent'][0], cxy[1], h_dict['Extent'][-1] - h_dict['Extent'][0], 1, color='blue',
-                      length_includes_head=True, head_width=5)
-            plt.arrow(cxy[0], v_dict['Extent'][0], 1, v_dict['Extent'][-1] - v_dict['Extent'][0], color='orange',
-                      length_includes_head=True, head_width=5)
-
-            plt.arrow(se_dict['Start'][0], se_dict['Start'][1], se_dict['Extent'][0], se_dict['Extent'][1],
-                      color='purple', length_includes_head=True, head_width=5)
-            plt.arrow(sw_dict['Start'][0], sw_dict['Start'][1], sw_dict['Extent'][0], sw_dict['Extent'][1],
-                      color='yellow', length_includes_head=True, head_width=5)
-
-            plt.legend([str(np.round(h_dict['Distance'], 2)) + 'mm',
-                        str(np.round(v_dict['Distance'], 2)) + 'mm',
-                        str(np.round(sw_dict['Distance'], 2)) + 'mm',
-                        str(np.round(se_dict['Distance'], 2)) + 'mm'])
-            plt.axis('off')
-            plt.title('Geometric Accuracy for Slice 5')
+            fig, axes = plt.subplots(3, 1)
+            fig.set_size_inches(8, 24)
+            fig.tight_layout(pad=4)
+
+            axes[0].imshow(img)
+            axes[0].scatter(cxy[0], cxy[1], c='red')
+            axes[0].axis('off')
+            axes[0].set_title('Centroid Location')
+
+            axes[1].imshow(mask)
+            axes[1].axis('off')
+            axes[1].set_title('Thresholding Result')
+
+            axes[2].imshow(img)
+            axes[2].arrow(length_dict['Horizontal Extent'][0], cxy[1], length_dict['Horizontal Extent'][-1]
+                          - length_dict['Horizontal Extent'][0], 1, color='blue', length_includes_head=True,
+                          head_width=5)
+            axes[2].arrow(cxy[0], length_dict['Vertical Extent'][0], 1, length_dict['Vertical Extent'][-1] -
+                          length_dict['Vertical Extent'][0], color='orange', length_includes_head=True, head_width=5)
+            axes[2].arrow(se_dict['Start'][0], se_dict['Start'][1], se_dict['Extent'][0], se_dict['Extent'][1],
+                          color='purple', length_includes_head=True, head_width=5)
+            axes[2].arrow(sw_dict['Start'][0], sw_dict['Start'][1], sw_dict['Extent'][0], sw_dict['Extent'][1],
+                          color='yellow', length_includes_head=True, head_width=5)
+
+            axes[2].legend([str(np.round(length_dict['Horizontal Distance'], 2)) + 'mm',
+                            str(np.round(length_dict['Vertical Distance'], 2)) + 'mm',
+                            str(np.round(sw_dict['Distance'], 2)) + 'mm',
+                            str(np.round(se_dict['Distance'], 2)) + 'mm'])
+            axes[2].axis('off')
+            axes[2].set_title('Geometric Accuracy for Slice 5')
 
             img_path = os.path.realpath(os.path.join(self.report_path, f'{self.key(dcm)}.png'))
             fig.savefig(img_path)
             self.report_files.append(img_path)
 
-        return h_dict['Distance'], v_dict['Distance'], sw_dict['Distance'], se_dict['Distance']
+        return length_dict['Horizontal Distance'], length_dict['Vertical Distance'], \
+            sw_dict['Distance'], se_dict['Distance']
 
-    def distortion_metric(self, L):
+    @staticmethod
+    def distortion_metric(L):
         err = [x - 190 for x in L]
         mean_err = np.mean(err)
 
diff --git a/hazenlib/tasks/acr_ghosting.py b/hazenlib/tasks/acr_ghosting.py
index 0433adb..8ea91ea 100644
--- a/hazenlib/tasks/acr_ghosting.py
+++ b/hazenlib/tasks/acr_ghosting.py
@@ -3,7 +3,7 @@ ACR Ghosting
 
 https://www.acraccreditation.org/-/media/acraccreditation/documents/mri/largephantomguidance.pdf
 
-Calculates percent-signal ghosting for slice 7 of the ACR phantom.
+Calculates the percent-signal ghosting for slice 7 of the ACR phantom.
 
 This script calculates the percentage signal ghosting in accordance with the ACR Guidance.
 This is done by first defining a large 200cm2 ROI before placing 10cm2 elliptical ROIs outside the phantom along the
@@ -19,34 +19,28 @@ import sys
 import traceback
 import os
 import numpy as np
-import skimage.morphology
 
 from hazenlib.HazenTask import HazenTask
+from hazenlib.ACRObject import ACRObject
 
 
 class ACRGhosting(HazenTask):
 
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
+        self.ACR_obj = None
 
     def run(self) -> dict:
         results = {}
-        z = []
-        for dcm in self.data:
-            z.append(dcm.ImagePositionPatient[2])
+        self.ACR_obj = ACRObject(self.data)
+        ghosting_dcm = self.ACR_obj.dcm[6]
 
-        idx_sort = np.argsort(z)
-
-        for dcm in self.data:
-            if dcm.ImagePositionPatient[2] == z[idx_sort[6]]:
-                try:
-                    result = self.get_signal_ghosting(dcm)
-                except Exception as e:
-                    print(f"Could not calculate the percent-signal ghosting for {self.key(dcm)} because of : {e}")
-                    traceback.print_exc(file=sys.stdout)
-                    continue
-
-                results[self.key(dcm)] = result
+        try:
+            result = self.get_signal_ghosting(ghosting_dcm)
+            results[self.key(ghosting_dcm)] = result
+        except Exception as e:
+            print(f"Could not calculate the percent-signal ghosting for {self.key(ghosting_dcm)} because of : {e}")
+            traceback.print_exc(file=sys.stdout)
 
         # only return reports if requested
         if self.report:
@@ -54,26 +48,14 @@ class ACRGhosting(HazenTask):
 
         return results
 
-    def centroid_com(self, dcm):
-        # Calculate centroid of object using a centre-of-mass calculation
-        thresh_img = dcm > 0.25 * np.max(dcm)
-        open_img = skimage.morphology.area_opening(thresh_img, area_threshold=500)
-        bhull = skimage.morphology.convex_hull_image(open_img)
-        coords = np.nonzero(bhull)  # row major - first array is columns
-
-        sum_x = np.sum(coords[1])
-        sum_y = np.sum(coords[0])
-        cxy = sum_x / coords[0].shape, sum_y / coords[1].shape
-
-        cxy = [cxy[0].astype(int), cxy[1].astype(int)]
-        return bhull, cxy
-
     def get_signal_ghosting(self, dcm):
         img = dcm.pixel_array
         res = dcm.PixelSpacing  # In-plane resolution from metadata
         r_large = np.ceil(80 / res[0]).astype(int)  # Required pixel radius to produce ~200cm2 ROI
         dims = img.shape
-        mask, cxy = self.centroid_com(img)
+
+        mask = self.ACR_obj.mask_image(img)
+        cxy = self.ACR_obj.centre
 
         nx = np.linspace(1, dims[0], dims[0])
         ny = np.linspace(1, dims[1], dims[1])
@@ -96,15 +78,15 @@ class ACRGhosting(HazenTask):
         else:
             w_factor = 1
 
-        w_ellipse = np.square((y - w_centre[0]) / (4 * w_factor)) + np.square(
-            (x - w_centre[1]) * w_factor) <= np.square(
-            10 / res[0])  # generate ellipse mask
+        w_ellipse = np.square((y - w_centre[0]) / (4 * w_factor)) + np.square((x - w_centre[1]) * w_factor) <= \
+                    np.square(10 / res[0])  # generate ellipse mask
 
         # EAST ELLIPSE
         e_point = np.argwhere(np.sum(mask, 0) > 0)[-1]  # find last column in mask
         e_centre = [cxy[1], e_point + np.ceil((dims[1] - e_point) / 2)]  # initialise centre of ellipse
         right_fov_to_centre = e_centre[1] + sad / 2 + 5  # edge of ellipse towards right FoV (+ tolerance)
-        centre_to_right_phantom = e_centre[1] - sad/2 - 5  # edge of ellipse towards right side of phantom (+ tolerance)
+        centre_to_right_phantom = e_centre[
+                                      1] - sad / 2 - 5  # edge of ellipse towards right side of phantom (+ tolerance)
         if right_fov_to_centre > dims[1] - 1 or centre_to_right_phantom < e_point:
             diffs = [dims[1] - 1 - right_fov_to_centre, centre_to_right_phantom - e_point]
             ind = diffs.index(max(diffs, key=abs))
@@ -157,40 +139,48 @@ class ACRGhosting(HazenTask):
         psg = 100 * np.absolute(
             ((n_ellipse_val + s_ellipse_val) - (w_ellipse_val + e_ellipse_val)) / (2 * large_roi_val))
 
-        psg = np.round(psg,3)
+        psg = np.round(psg, 3)
 
         if self.report:
             import matplotlib.pyplot as plt
+            fig, axes = plt.subplots(2, 1)
+            fig.set_size_inches(8, 16)
+            fig.tight_layout(pad=4)
+
             theta = np.linspace(0, 2 * np.pi, 360)
-            fig = plt.figure()
-            fig.set_size_inches(8, 8)
-            plt.imshow(img)
-
-            plt.plot(r_large * np.cos(theta) + cxy[0], r_large * np.sin(theta) + cxy[1] + 5 / res[1], c='black')
-            plt.text(cxy[0] - 3 * np.floor(10 / res[0]), cxy[1] + np.floor(10 / res[1]),
-                     "Mean = " + str(np.round(large_roi_val, 2)), c='white')
-
-            plt.plot(10. / res[0] * np.cos(theta) / w_factor + w_centre[1],
-                     10. / res[0] * np.sin(theta) * 4 * w_factor + w_centre[0], c='red')
-            plt.text(w_centre[1] - np.floor(10 / res[0]), w_centre[0], "Mean = " + str(np.round(w_ellipse_val, 2)),
-                     c='white')
-
-            plt.plot(10. / res[0] * np.cos(theta) / e_factor + e_centre[1],
-                     10. / res[0] * np.sin(theta) * 4 * e_factor + e_centre[0], c='red')
-            plt.text(e_centre[1] - np.floor(30 / res[0]), e_centre[0], "Mean = " + str(np.round(e_ellipse_val, 2)),
-                     c='white')
-
-            plt.plot(10. / res[0] * np.cos(theta) * 4 * n_factor + n_centre[1],
-                     10. / res[0] * np.sin(theta) / n_factor + n_centre[0], c='red')
-            plt.text(n_centre[1] - 5 * np.floor(10 / res[0]), n_centre[0], "Mean = " + str(np.round(n_ellipse_val, 2)),
-                     c='white')
-
-            plt.plot(10. / res[0] * np.cos(theta) * 4 * s_factor + s_centre[1],
-                     10. / res[0] * np.sin(theta) / s_factor + s_centre[0], c='red')
-            plt.text(s_centre[1], s_centre[0], "Mean = " + str(np.round(s_ellipse_val, 2)), c='white')
-
-            plt.axis('off')
-            plt.title('Percent Signal Ghosting = ' + str(np.round(psg, 3)) + '%')
+
+            axes[0].imshow(img)
+            axes[0].scatter(cxy[0], cxy[1], c='red')
+            axes[0].axis('off')
+            axes[0].set_title('Centroid Location')
+
+            axes[1].imshow(img)
+            axes[1].plot(r_large * np.cos(theta) + cxy[0], r_large * np.sin(theta) + cxy[1] + 5 / res[1], c='black')
+            axes[1].text(cxy[0] - 3 * np.floor(10 / res[0]), cxy[1] + np.floor(10 / res[1]),
+                         "Mean = " + str(np.round(large_roi_val, 2)), c='white')
+
+            axes[1].plot(10. / res[0] * np.cos(theta) / w_factor + w_centre[1],
+                         10. / res[0] * np.sin(theta) * 4 * w_factor + w_centre[0], c='red')
+            axes[1].text(w_centre[1] - np.floor(10 / res[0]), w_centre[0], "Mean = " + str(np.round(w_ellipse_val, 2)),
+                         c='white')
+
+            axes[1].plot(10. / res[0] * np.cos(theta) / e_factor + e_centre[1],
+                         10. / res[0] * np.sin(theta) * 4 * e_factor + e_centre[0], c='red')
+            axes[1].text(e_centre[1] - np.floor(30 / res[0]), e_centre[0], "Mean = " + str(np.round(e_ellipse_val, 2)),
+                         c='white')
+
+            axes[1].plot(10. / res[0] * np.cos(theta) * 4 * n_factor + n_centre[1],
+                         10. / res[0] * np.sin(theta) / n_factor + n_centre[0], c='red')
+            axes[1].text(n_centre[1] - 5 * np.floor(10 / res[0]), n_centre[0],
+                         "Mean = " + str(np.round(n_ellipse_val, 2)),
+                         c='white')
+
+            axes[1].plot(10. / res[0] * np.cos(theta) * 4 * s_factor + s_centre[1],
+                         10. / res[0] * np.sin(theta) / s_factor + s_centre[0], c='red')
+            axes[1].text(s_centre[1], s_centre[0], "Mean = " + str(np.round(s_ellipse_val, 2)), c='white')
+
+            axes[1].axis('off')
+            axes[1].set_title('Percent Signal Ghosting = ' + str(np.round(psg, 3)) + '%')
             img_path = os.path.realpath(os.path.join(self.report_path, f'{self.key(dcm)}.png'))
             fig.savefig(img_path)
             self.report_files.append(img_path)
diff --git a/hazenlib/tasks/acr_slice_position.py b/hazenlib/tasks/acr_slice_position.py
index 337916b..40fbd5e 100644
--- a/hazenlib/tasks/acr_slice_position.py
+++ b/hazenlib/tasks/acr_slice_position.py
@@ -33,62 +33,33 @@ import skimage.morphology
 import skimage.measure
 
 from hazenlib.HazenTask import HazenTask
-
-
-def find_n_peaks(data, n, height=1):
-    peaks = scipy.signal.find_peaks(data, height)
-    pk_heights = peaks[1]['peak_heights']
-    pk_ind = peaks[0]
-    highest_peaks = pk_ind[(-pk_heights).argsort()[:n]]  # find n highest peaks
-
-    return np.sort(highest_peaks)
+from hazenlib.ACRObject import ACRObject
 
 
 class ACRSlicePosition(HazenTask):
 
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
+        self.ACR_obj = None
 
     def run(self) -> dict:
         results = {}
-        z = []
-        for dcm in self.data:
-            z.append(dcm.ImagePositionPatient[2])
-
-        idx_sort = np.argsort(z)
-
-        for dcm in self.data:
-            curr_z = dcm.ImagePositionPatient[2]
-            if curr_z in (z[idx_sort[0]], z[idx_sort[10]]):
-                try:
-                    result = self.get_slice_position(dcm)
-                except Exception as e:
-                    print(f"Could not calculate the bar length difference for {self.key(dcm)} because of : {e}")
-                    traceback.print_exc(file=sys.stdout)
-                    continue
-
+        self.ACR_obj = ACRObject(self.data)
+        dcms = [self.ACR_obj.dcm[0], self.ACR_obj.dcm[-1]]
+        for dcm in dcms:
+            try:
+                result = self.get_slice_position(dcm)
                 results[self.key(dcm)] = result
-
+            except Exception as e:
+                print(f"Could not calculate the bar length difference for {self.key(dcm)} because of : {e}")
+                traceback.print_exc(file=sys.stdout)
+                continue
         # only return reports if requested
         if self.report:
             results['reports'] = {'images': self.report_files}
 
         return results
 
-    def centroid_com(self, dcm):
-        # Calculate centroid of object using a centre-of-mass calculation
-        thresh_img = dcm > 0.25 * np.max(dcm)
-        open_img = skimage.morphology.area_opening(thresh_img, area_threshold=500)
-        bhull = skimage.morphology.convex_hull_image(open_img)
-        coords = np.nonzero(bhull)  # row major - first array is columns
-
-        sum_x = np.sum(coords[1])
-        sum_y = np.sum(coords[0])
-        cxy = sum_x / coords[0].shape, sum_y / coords[1].shape
-
-        cxy = [cxy[0].astype(int), cxy[1].astype(int)]
-        return bhull, cxy
-
     def find_wedges(self, img, mask, res):
         # X COORDINATES
         x_investigate_region = np.ceil(35 / res[0]).astype(int)  # define width of region to test (comparable to wedges)
@@ -115,7 +86,7 @@ class ACRSlicePosition(HazenTask):
         mean_x_profile = np.mean(invest_x, 1)  # mean of horizontal projections of phantom
         abs_diff_x_profile = np.abs(np.diff(mean_x_profile))  # absolute first derivative of mean
 
-        x_peaks = find_n_peaks(abs_diff_x_profile, 2)  # find two highest peaks
+        x_peaks, _ = self.ACR_obj.find_n_highest_peaks(abs_diff_x_profile, 2)  # find two highest peaks
         x_locs = w_point + x_peaks  # x coordinates of these peaks in image coordinate system(before diff operation)
 
         width_pts = [x_locs[0], x_locs[1]]  # width of wedges
@@ -147,7 +118,7 @@ class ACRSlicePosition(HazenTask):
         mean_y_profile = np.mean(invest_y, 1)  # mean of vertical projections of phantom
         abs_diff_y_profile = np.abs(np.diff(mean_y_profile))  # absolute first derivative of mean
 
-        y_peaks = find_n_peaks(abs_diff_y_profile, 2)  # find two highest peaks
+        y_peaks, _ = self.ACR_obj.find_n_highest_peaks(abs_diff_y_profile, 2)  # find two highest peaks
         y_locs = w_point + y_peaks - 1  # y coordinates of these peaks in image coordinate system(before diff operation)
 
         if y_locs[1] - y_locs[0] < 5 / res[1]:
@@ -156,14 +127,15 @@ class ACRSlicePosition(HazenTask):
             y = np.round(np.min(y_locs) + 0.25 * np.abs(np.diff(y_locs)))  # define y coordinate
 
         dist_to_y = np.abs(n_point - y[0]) * res[1]  # distance to y from top of phantom
-        y_pts = np.append(y, np.round(y[0] + (47 - dist_to_y) / res[1])).astype(int)  # place 2nd y point 47mm from top of phantom
+        y_pts = np.append(y, np.round(y[0] + (47 - dist_to_y) / res[1])).astype(
+            int)  # place 2nd y point 47mm from top of phantom
 
         return x_pts, y_pts
 
     def get_slice_position(self, dcm):
         img = dcm.pixel_array
         res = dcm.PixelSpacing  # In-plane resolution from metadata
-        mask, cxy = self.centroid_com(img)
+        mask = self.ACR_obj.mask_image(self.ACR_obj.images[6])
         x_pts, y_pts = self.find_wedges(img, mask, res)
 
         line_prof_L = skimage.measure.profile_line(img, (y_pts[0], x_pts[0]), (y_pts[1], x_pts[0]),
@@ -179,7 +151,7 @@ class ACRSlicePosition(HazenTask):
         interp_line_prof_R = scipy.interpolate.interp1d(x, line_prof_R)(new_x)  # interpolate right line profile
 
         delta = interp_line_prof_L - interp_line_prof_R  # difference of line profiles
-        peaks = find_n_peaks(abs(delta), 2, 0.5 * np.max(abs(delta)))  # find two highest peaks
+        peaks, _ = ACRObject.find_n_highest_peaks(abs(delta), 2, 0.5 * np.max(abs(delta)))  # find two highest peaks
 
         if len(peaks) == 1:
             peaks = [peaks[0] - 50, peaks[0] + 50]  # if only one peak, set dummy range
@@ -213,30 +185,31 @@ class ACRSlicePosition(HazenTask):
 
         if self.report:
             import matplotlib.pyplot as plt
-            fig = plt.figure()
-            plt.suptitle('Bar Length Difference = ' + str(np.round(dL, 2)) + 'mm', x=0.5, ha='center')
-            fig.set_size_inches(8, 8)
-
-            plt.subplot(2, 2, (1, 3))
-            plt.imshow(img)
-            plt.plot([x_pts[0], x_pts[0]], [y_pts[0], y_pts[1]], 'b')
-            plt.plot([x_pts[1], x_pts[1]], [y_pts[0], y_pts[1]], 'r')
-            plt.axis('off')
-            plt.tight_layout()
-
-            plt.subplot(2, 2, 2)
-            plt.grid()
-            plt.plot((1 / interp_factor) * np.linspace(1, len(interp_line_prof_L), len(interp_line_prof_L)) * res[1],
-                     interp_line_prof_L, 'b')
-            plt.plot((1 / interp_factor) * np.linspace(1, len(interp_line_prof_R), len(interp_line_prof_R)) * res[1],
-                     interp_line_prof_R, 'r')
-            plt.title('Original Line Profiles')
-            plt.xlabel('Relative Pixel Position (mm)')
-            plt.tight_layout()
-
-            plt.subplot(2, 2, 4)
-            plt.grid()
-            plt.plot((1 / interp_factor) * np.linspace(1, len(interp_line_prof_L), len(interp_line_prof_L)) * res[1],
+            fig, axes = plt.subplots(4, 1)
+            fig.set_size_inches(8, 32)
+            fig.tight_layout(pad=4)
+
+            axes[0].imshow(mask)
+            axes[0].axis('off')
+            axes[0].set_title('Thresholding Result')
+
+            axes[1].imshow(img)
+            axes[1].plot([x_pts[0], x_pts[0]], [y_pts[0], y_pts[1]], 'b')
+            axes[1].plot([x_pts[1], x_pts[1]], [y_pts[0], y_pts[1]], 'r')
+            axes[1].axis('off')
+            axes[1].set_title('Line Profiles')
+
+            axes[2].grid()
+            axes[2].plot(
+                (1 / interp_factor) * np.linspace(1, len(interp_line_prof_L), len(interp_line_prof_L)) * res[1],
+                interp_line_prof_L, 'b')
+            axes[2].plot(
+                (1 / interp_factor) * np.linspace(1, len(interp_line_prof_R), len(interp_line_prof_R)) * res[1],
+                interp_line_prof_R, 'r')
+            axes[2].set_title('Original Line Profiles')
+            axes[2].set_xlabel('Relative Pixel Position (mm)')
+
+            axes[3].plot((1 / interp_factor) * np.linspace(1, len(interp_line_prof_L), len(interp_line_prof_L)) * res[1],
                      interp_line_prof_L, 'b')
 
             shift_line = np.roll(interp_line_prof_R, pos * shift)
@@ -249,11 +222,11 @@ class ACRSlicePosition(HazenTask):
             else:
                 shift_line[0:np.abs(pos) * shift] = np.nan
 
-            plt.plot((1 / interp_factor) * np.linspace(1, len(interp_line_prof_L), len(interp_line_prof_L)) * res[1],
+            axes[3].grid()
+            axes[3].plot((1 / interp_factor) * np.linspace(1, len(interp_line_prof_L), len(interp_line_prof_L)) * res[1],
                      shift_line, 'r')
-            plt.title('Shifted Line Profiles')
-            plt.xlabel('Relative Pixel Position (mm)')
-            plt.tight_layout()
+            axes[3].set_title('Shifted Line Profiles')
+            axes[3].set_xlabel('Relative Pixel Position (mm)')
 
             img_path = os.path.realpath(os.path.join(self.report_path, f'{self.key(dcm)}.png'))
             fig.savefig(img_path)
diff --git a/hazenlib/tasks/acr_slice_thickness.py b/hazenlib/tasks/acr_slice_thickness.py
index 8e9571a..709c5df 100644
--- a/hazenlib/tasks/acr_slice_thickness.py
+++ b/hazenlib/tasks/acr_slice_thickness.py
@@ -16,44 +16,30 @@ yassine.azma@rmh.nhs.uk
 import sys
 import traceback
 import os
-from hazenlib.HazenTask import HazenTask
 import numpy as np
 import skimage.morphology
 import skimage.measure
 import scipy
 
-
-def find_n_peaks(data, n, height=1):
-    peaks = scipy.signal.find_peaks(data, height)
-    pk_heights = peaks[1]['peak_heights']
-    pk_ind = peaks[0]
-    highest_peaks = pk_ind[(-pk_heights).argsort()[:n]]  # find n highest peaks
-
-    return np.sort(highest_peaks)
+from hazenlib.HazenTask import HazenTask
+from hazenlib.ACRObject import ACRObject
 
 
 class ACRSliceThickness(HazenTask):
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
+        self.ACR_obj = None
 
     def run(self) -> dict:
         results = {}
-        z = []
-        for dcm in self.data:
-            z.append(dcm.ImagePositionPatient[2])
-
-        idx_sort = np.argsort(z)
-        for dcm in self.data:
-            curr_z = dcm.ImagePositionPatient[2]
-            if curr_z == z[idx_sort[0]]:
-                try:
-                    result = self.get_slice_thickness(dcm)
-                except Exception as e:
-                    print(f"Could not calculate the slice thickness for {self.key(dcm)} because of : {e}")
-                    traceback.print_exc(file=sys.stdout)
-                    continue
-
-                results[self.key(dcm)] = result
+        self.ACR_obj = ACRObject(self.data)
+        slice_thickness_dcm = self.ACR_obj.dcm[0]
+        try:
+            result = self.get_slice_thickness(slice_thickness_dcm)
+            results[self.key(slice_thickness_dcm)] = result
+        except Exception as e:
+            print(f"Could not calculate the slice thickness for {self.key(slice_thickness_dcm)} because of : {e}")
+            traceback.print_exc(file=sys.stdout)
 
         # only return reports if requested
         if self.report:
@@ -61,20 +47,6 @@ class ACRSliceThickness(HazenTask):
 
         return results
 
-    def centroid_com(self, dcm):
-        # Calculate centroid of object using a centre-of-mass calculation
-        thresh_img = dcm > 0.25 * np.max(dcm)
-        open_img = skimage.morphology.area_opening(thresh_img, area_threshold=500)
-        bhull = skimage.morphology.convex_hull_image(open_img)
-        coords = np.nonzero(bhull)  # row major - first array is columns
-
-        sum_x = np.sum(coords[1])
-        sum_y = np.sum(coords[0])
-        cxy = sum_x / coords[0].shape, sum_y / coords[1].shape
-
-        cxy = [int(cxy[0]), int(cxy[1])]
-        return bhull, cxy
-
     def find_ramps(self, img, centre, res):
         # X
         investigate_region = int(np.ceil(5.5 / res[1]).item())
@@ -97,7 +69,7 @@ class ACRSliceThickness(HazenTask):
         # [2] - hypointense region with ramps and hyperintense phantom
         # [3] - hyperintense phantom and background
 
-        x_peaks = find_n_peaks(abs_diff_x_profile, 4)
+        x_peaks, _ = self.ACR_obj.find_n_highest_peaks(abs_diff_x_profile, 4)
         x_locs = np.sort(x_peaks) - 1
 
         width_pts = [x_locs[1], x_locs[2]]
@@ -105,14 +77,13 @@ class ACRSliceThickness(HazenTask):
 
         # take rough estimate of x points for later line profiles
         x = np.round([np.min(width_pts) + 0.2 * width, np.max(width_pts) - 0.2 * width])
-
         # Y
         c = skimage.measure.profile_line(img, (centre[1] - 2 * investigate_region, centre[0]),
                                          (centre[1] + 2 * investigate_region, centre[0]), mode='constant').flatten()
 
         abs_diff_y_profile = np.absolute(np.diff(c))
 
-        y_peaks = find_n_peaks(abs_diff_y_profile, 2)
+        y_peaks, _ = self.ACR_obj.find_n_highest_peaks(abs_diff_y_profile, 2)
         y_locs = centre[1] - 2 * investigate_region + 1 + y_peaks
         height = np.max(y_locs) - np.min(y_locs)
 
@@ -149,7 +120,7 @@ class ACRSliceThickness(HazenTask):
     def get_slice_thickness(self, dcm):
         img = dcm.pixel_array
         res = dcm.PixelSpacing  # In-plane resolution from metadata
-        mask, cxy = self.centroid_com(img)
+        cxy = self.ACR_obj.centre
         x_pts, y_pts = self.find_ramps(img, cxy, res)
 
         interp_factor = 5
@@ -162,11 +133,9 @@ class ACRSliceThickness(HazenTask):
         fwhm_store = []
         for i, offset in enumerate(offsets):
             lines = [skimage.measure.profile_line(img, (offset + y_pts[0], x_pts[0]), (offset + y_pts[0], x_pts[1]),
-                                                  linewidth=2,
-                                                  mode='constant').flatten(),
+                                                  linewidth=2, mode='constant').flatten(),
                      skimage.measure.profile_line(img, (offset + y_pts[1], x_pts[0]), (offset + y_pts[1], x_pts[1]),
-                                                  linewidth=2,
-                                                  mode='constant').flatten()]
+                                                  linewidth=2, mode='constant').flatten()]
 
             interp_lines = [scipy.interpolate.interp1d(sample, line)(new_sample) for line in lines]
             fwhm = [self.FWHM(interp_line) for interp_line in interp_lines]
@@ -186,38 +155,43 @@ class ACRSliceThickness(HazenTask):
 
         if self.report:
             import matplotlib.pyplot as plt
-            fig = plt.figure()
-
-            plt.subplot(2, 2, (1, 3))
-            plt.imshow(img)
-            plt.plot([x_pts[0], x_pts[1]], offsets[z_ind] + [y_pts[0], y_pts[0]], 'b-')
-            plt.plot([x_pts[0], x_pts[1]], offsets[z_ind] + [y_pts[1], y_pts[1]], 'r-')
-            plt.axis('off')
-            plt.title(f'Slice Thickness = {np.round(dz[z_ind], 2)}mm')
-            plt.tight_layout()
-
-            width = fwhm_store[z_ind][1][0] * (1 / interp_factor) * res[0], fwhm_store[z_ind][1][1] * (1 / interp_factor) * \
-                    res[0]
+            fig, axes = plt.subplots(4, 1)
+            fig.set_size_inches(8, 24)
+            fig.tight_layout(pad=4)
+
+            axes[0].imshow(img)
+            axes[0].scatter(cxy[0], cxy[1], c='red')
+            axes[0].axis('off')
+            axes[0].set_title('Centroid Location')
+
+            axes[1].imshow(img)
+            axes[1].plot([x_pts[0], x_pts[1]], offsets[z_ind] + [y_pts[0], y_pts[0]], 'b-')
+            axes[1].plot([x_pts[0], x_pts[1]], offsets[z_ind] + [y_pts[1], y_pts[1]], 'r-')
+            axes[1].axis('off')
+            axes[1].set_title('Line Profiles')
+
+            width = fwhm_store[z_ind][1][0] * (1 / interp_factor) * res[0], fwhm_store[z_ind][1][1] * (
+                    1 / interp_factor) * res[0]
             x_ramp = new_sample * res[0]
             x_extent = np.max(x_ramp)
             y_ramp = line_store[z_ind][1]
             y_extent = np.max(y_ramp)
             max_loc = np.argmax(y_ramp) * (1 / interp_factor) * res[0]
 
-            plt.subplot(2, 2, 2)
-            plt.plot(x_ramp, y_ramp, 'r', label=f'FWHM={np.round(ramp_length[1][z_ind], 2)}mm')
-            plt.axhline(0.5 * y_extent, xmin=width[0] / x_extent, xmax=width[1] / x_extent, linestyle='dashdot', color='k')
-            plt.axvline(max_loc, ymin=0, ymax=10 / 11, linestyle='dashdot', color='k')
+            axes[2].plot(x_ramp, y_ramp, 'r', label=f'FWHM={np.round(ramp_length[1][z_ind], 2)}mm')
+            axes[2].axhline(0.5 * y_extent, xmin=width[0] / x_extent, xmax=width[1] / x_extent, linestyle='dashdot',
+                        color='k')
+            axes[2].axvline(max_loc, ymin=0, ymax=10 / 11, linestyle='dashdot', color='k')
 
-            plt.xlabel('Relative Position (mm)')
-            plt.xlim([0, x_extent])
-            plt.ylim([0, y_extent * 1.1])
-            plt.title('Upper Ramp')
-            plt.grid()
-            plt.legend(loc='best')
-            plt.tight_layout()
+            axes[2].set_xlabel('Relative Position (mm)')
+            axes[2].set_xlim([0, x_extent])
+            axes[2].set_ylim([0, y_extent * 1.1])
+            axes[2].set_title('Upper Ramp')
+            axes[2].grid()
+            axes[2].legend(loc='best')
 
-            width = fwhm_store[z_ind][0][0] * (1 / interp_factor) * res[0], fwhm_store[z_ind][0][1] * (1 / interp_factor) * \
+            width = fwhm_store[z_ind][0][0] * (1 / interp_factor) * res[0], fwhm_store[z_ind][0][1] * (
+                    1 / interp_factor) * \
                     res[0]
             x_ramp = new_sample * res[0]
             x_extent = np.max(x_ramp)
@@ -225,23 +199,20 @@ class ACRSliceThickness(HazenTask):
             y_extent = np.max(y_ramp)
             max_loc = np.argmax(y_ramp) * (1 / interp_factor) * res[0]
 
-            plt.subplot(2, 2, 4)
-            plt.plot(x_ramp, y_ramp, 'b', label=f'FWHM={np.round(ramp_length[0][z_ind], 2)}mm')
-            plt.axhline(0.5 * y_extent, xmin=width[0] / x_extent, xmax=width[1] / x_extent, linestyle='dashdot', color='k')
-            plt.axvline(max_loc, ymin=0, ymax=10 / 11, linestyle='dashdot', color='k')
+            axes[3].plot(x_ramp, y_ramp, 'b', label=f'FWHM={np.round(ramp_length[0][z_ind], 2)}mm')
+            axes[3].axhline(0.5 * y_extent, xmin=width[0] / x_extent, xmax=width[1] / x_extent, linestyle='dashdot',
+                        color='k')
+            axes[3].axvline(max_loc, ymin=0, ymax=10 / 11, linestyle='dashdot', color='k')
 
-            plt.xlabel('Relative Position (mm)')
-            plt.xlim([0, x_extent])
-            plt.ylim([0, y_extent * 1.1])
-            plt.title('Lower Ramp')
-            plt.grid()
-            plt.legend(loc='best')
-            plt.tight_layout()
+            axes[3].set_xlabel('Relative Position (mm)')
+            axes[3].set_xlim([0, x_extent])
+            axes[3].set_ylim([0, y_extent * 1.1])
+            axes[3].set_title('Lower Ramp')
+            axes[3].grid()
+            axes[3].legend(loc='best')
 
             img_path = os.path.realpath(os.path.join(self.report_path, f'{self.key(dcm)}_slice_thickness.png'))
             fig.savefig(img_path)
             self.report_files.append(img_path)
 
         return slice_thickness
-
-
diff --git a/hazenlib/tasks/acr_snr.py b/hazenlib/tasks/acr_snr.py
index 1d45f85..087c3de 100644
--- a/hazenlib/tasks/acr_snr.py
+++ b/hazenlib/tasks/acr_snr.py
@@ -16,37 +16,45 @@ Created by Neil Heraghty (Adapted by Yassine Azma)
 import sys
 import traceback
 import os
+
 import hazenlib.utils
 from hazenlib.HazenTask import HazenTask
 from scipy import ndimage
+
 import numpy as np
 import skimage.morphology
 import pydicom
+from scipy import ndimage
+
+from hazenlib.HazenTask import HazenTask
+from hazenlib.ACRObject import ACRObject
 
 
 class ACRSNR(HazenTask):
 
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
+        self.ACR_obj = None
+
         self.data2 = []
 
     def run(self, measured_slice_width=None, subtract=None) -> dict:
+        
         if measured_slice_width is not None:
             measured_slice_width = float(measured_slice_width)
+        
         snr_results = {}
+        self.ACR_obj = [ACRObject(self.data)]
+        snr_dcm = self.ACR_obj[0].dcm[6]
 
         # SINGLE METHOD (SMOOTHING)
         if subtract is None:
-            z = [dcm.ImagePositionPatient[2] for dcm in self.data]
-
-            idx_sort = np.argsort(z)
-            ind = idx_sort[6]
             try:
-                snr, normalised_snr = self.snr_by_smoothing(self.data[ind], measured_slice_width)
-                snr_results[f"snr_smoothing_measured_{self.key(self.data[0])}"] = round(snr, 2)
-                snr_results[f"snr_smoothing_normalised_{self.key(self.data[0])}"] = round(normalised_snr, 2)
+                snr, normalised_snr = self.snr_by_smoothing(snr_dcm, measured_slice_width)
+                snr_results[f"snr_smoothing_measured_{self.key(snr_dcm)}"] = round(snr, 2)
+                snr_results[f"snr_smoothing_normalised_{self.key(snr_dcm)}"] = round(normalised_snr, 2)
             except Exception as e:
-                print(f"Could not calculate the SNR for {self.key(self.data)} because of : {e}")
+                print(f"Could not calculate the SNR for {self.key(snr_dcm)} because of : {e}")
                 traceback.print_exc(file=sys.stdout)
         # SUBTRACTION METHOD
         else:
@@ -54,22 +62,19 @@ class ACRSNR(HazenTask):
             filenames = [f'{subtract}/{file}' for file in temp]
 
             self.data2 = [pydicom.dcmread(dicom) for dicom in filenames]
-
-            z = [dcm.ImagePositionPatient[2] for dcm in self.data]
-            z2 = [dcm.ImagePositionPatient[2] for dcm in self.data2]
-
-            idx_sort, idx_sort2 = np.argsort(z), np.argsort(z2)
-            ind, ind2 = idx_sort[6], idx_sort2[6]
+            self.ACR_obj.append(ACRObject(self.data2))
+            snr_dcm2 = self.ACR_obj[1].dcm[6]
             try:
-                snr, normalised_snr = self.snr_by_subtraction(self.data[ind], self.data2[ind2])
-                snr_results[f"snr_subtraction_measured_{self.key(self.data[0])}"] = round(snr, 2)
-                snr_results[f"snr_subtraction_normalised_{self.key(self.data[0])}"] = round(normalised_snr, 2)
+                snr, normalised_snr = self.snr_by_subtraction(snr_dcm, snr_dcm2)
+                snr_results[f"snr_subtraction_measured_{self.key(snr_dcm)}"] = round(snr, 2)
+                snr_results[f"snr_subtraction_normalised_{self.key(snr_dcm)}"] = round(normalised_snr, 2)
             except Exception as e:
-                print(f"Could not calculate the SNR for {self.key(self.data)} and "
-                      f"{self.key(self.data2)} because of : {e}")
+                print(f"Could not calculate the SNR for {self.key(snr_dcm)} and "
+                      f"{self.key(snr_dcm2)} because of : {e}")
                 traceback.print_exc(file=sys.stdout)
 
-        results = {self.key(self.data[0]): snr_results}
+
+        results = {self.key(snr_dcm): snr_results, 'reports': {'images': self.report_files}}
 
         # only return reports if requested
         if self.report:
@@ -77,21 +82,6 @@ class ACRSNR(HazenTask):
 
         return results
 
-    def centroid(self, dcm):
-        img = dcm.pixel_array
-        mask = img > 0.25 * np.max(img)
-        open_img = skimage.morphology.area_opening(mask, area_threshold=500)
-        mask = skimage.morphology.convex_hull_image(open_img)
-
-        coords = np.nonzero(mask)  # row major - first array is columns
-
-        sum_x = np.sum(coords[1])
-        sum_y = np.sum(coords[0])
-        cxy = sum_x / coords[0].shape, sum_y / coords[1].shape
-
-        cxy = [cxy[0].astype(int), cxy[1].astype(int)]
-        return mask, cxy
-
     def get_normalised_snr_factor(self, dcm, measured_slice_width=None) -> float:
         dx, dy = hazenlib.utils.get_pixel_size(dcm)
         bandwidth = hazenlib.utils.get_bandwidth(dcm)
@@ -126,7 +116,8 @@ class ACRSNR(HazenTask):
         """
         a = dcm.pixel_array.astype('int')
 
-        # filter size = 9, following MATLAB code and McCann 2013 paper for head coil, although note McCann 2013 recommends 25x25 for body coil.
+        # filter size = 9, following MATLAB code and McCann 2013 paper for head coil, although note McCann 2013
+        # recommends 25x25 for body coil.
         filtered_array = ndimage.uniform_filter(a, 25, mode='constant')
         return filtered_array
 
@@ -196,7 +187,7 @@ class ACRSNR(HazenTask):
         normalised_snr: float
 
         """
-        _, centre = self.centroid(dcm)
+        centre = self.ACR_obj[0].centre
         col, row = centre
         noise_img = self.get_noise_image(dcm)
 
@@ -205,7 +196,8 @@ class ACRSNR(HazenTask):
 
         noise = [np.std(roi, ddof=1) for roi in
                  self.get_roi_samples(ax=None, dcm=noise_img, centre_col=int(col), centre_row=int(row))]
-        # note no root_2 factor in noise for smoothed subtraction (one image) method, replicating Matlab approach and McCann 2013
+        # note no root_2 factor in noise for smoothed subtraction (one image) method, replicating Matlab approach and
+        # McCann 2013
 
         snr = np.mean(np.divide(signal, noise))
 
@@ -213,15 +205,17 @@ class ACRSNR(HazenTask):
 
         if self.report:
             import matplotlib.pyplot as plt
-            fig, axes = plt.subplots(1, 1)
-            fig.set_size_inches(5, 5)
-            fig.tight_layout(pad=1)
+            fig, axes = plt.subplots(2, 1)
+            fig.set_size_inches(8, 16)
+            fig.tight_layout(pad=4)
+
+            axes[0].imshow(dcm.pixel_array)
+            axes[0].scatter(centre[0], centre[1], c='red')
+            axes[0].set_title('Centroid Location')
 
-            axes.set_title('smoothed noise image')
-            axes.imshow(noise_img, cmap='gray', label='smoothed noise image')
-            axes.scatter(col, row, 10, marker="+", label='centre')
-            self.get_roi_samples(axes, dcm, int(col), int(row))
-            axes.legend()
+            axes[1].set_title('Smoothed Noise Image')
+            axes[1].imshow(noise_img, cmap='gray')
+            self.get_roi_samples(axes[1], dcm, int(col), int(row))
 
             img_path = os.path.realpath(os.path.join(self.report_path, f'{self.key(dcm)}_smoothing.png'))
             fig.savefig(img_path)
@@ -242,7 +236,7 @@ class ACRSNR(HazenTask):
         -------
 
         """
-        _, centre = self.centroid(dcm1)
+        centre = self.ACR_obj[0].centre
         col, row = centre
 
         difference = np.subtract(dcm1.pixel_array.astype('int'), dcm2.pixel_array.astype('int'))
@@ -259,15 +253,19 @@ class ACRSNR(HazenTask):
 
         if self.report:
             import matplotlib.pyplot as plt
-            fig, axes = plt.subplots(1, 1)
-            fig.set_size_inches(5, 5)
-            fig.tight_layout(pad=1)
-
-            axes.set_title('difference image')
-            axes.imshow(difference, cmap='gray', label='difference image')
-            axes.scatter(col, row, 10, marker="+", label='centre')
-            self.get_roi_samples(axes, dcm1, int(col), int(row))
-            axes.legend()
+            fig, axes = plt.subplots(2, 1)
+            fig.set_size_inches(8, 16)
+            fig.tight_layout(pad=4)
+
+            axes[0].imshow(dcm1.pixel_array)
+            axes[0].scatter(centre[0], centre[1], c='red')
+            axes[0].axis('off')
+            axes[0].set_title('Centroid Location')
+
+            axes[1].set_title('Difference Image')
+            axes[1].imshow(difference, cmap='gray',)
+            self.get_roi_samples(axes[1], dcm1, int(col), int(row))
+            axes[1].axis('off')
 
             img_path = os.path.realpath(os.path.join(self.report_path, f'{self.key(dcm1)}_snr_subtraction.png'))
             fig.savefig(img_path)
diff --git a/hazenlib/tasks/acr_spatial_resolution.py b/hazenlib/tasks/acr_spatial_resolution.py
index f97d8a0..9b45572 100644
--- a/hazenlib/tasks/acr_spatial_resolution.py
+++ b/hazenlib/tasks/acr_spatial_resolution.py
@@ -44,42 +44,35 @@ import skimage.morphology
 import skimage.measure
 
 from hazenlib.HazenTask import HazenTask
-
-
-def find_n_peaks(data, n, height=1):
-    peaks = scipy.signal.find_peaks(data, height)
-    pk_heights = peaks[1]['peak_heights']
-    pk_ind = peaks[0]
-
-    peak_heights = pk_heights[(-pk_heights).argsort()[:n]]
-    peak_locs = pk_ind[(-pk_heights).argsort()[:n]]  # find n highest peaks
-
-    return np.sort(peak_locs), np.sort(peak_heights)
+from hazenlib.ACRObject import ACRObject
+from hazenlib.logger import logger
 
 
 class ACRSpatialResolution(HazenTask):
 
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
+        self.ACR_obj = None
 
     def run(self) -> dict:
         mtf_results = {}
-        z = []
-        for dcm in self.data:
-            z.append(dcm.ImagePositionPatient[2])
-
-        idx_sort = np.argsort(z)
-
-        for dcm in self.data:
-            if dcm.ImagePositionPatient[2] == z[idx_sort[0]]:
-                try:
-                    raw_res, fitted_res = self.get_mtf50(dcm)
-                    mtf_results[f"raw_mtf50_{self.key(self.data[0])}"] = raw_res
-                    mtf_results[f"fitted_mtf50_{self.key(self.data[0])}"] = fitted_res
-                except Exception as e:
-                    print(f"Could not calculate the spatial resolution for {self.key(dcm)} because of : {e}")
-                    traceback.print_exc(file=sys.stdout)
-                    continue
+        self.ACR_obj = ACRObject(self.data)
+        rot_ang = self.ACR_obj.rot_angle
+
+        if np.round(np.abs(rot_ang), 2) < 3:
+            logger.warning(f'The estimated rotation angle of the ACR phantom is {np.round(rot_ang, 3)} degrees, which '
+                           f'is less than the recommended 3 degrees. Results will be unreliable!')
+
+        mtf_dcm = self.ACR_obj.dcm[0]
+
+        try:
+            raw_res, fitted_res = self.get_mtf50(mtf_dcm)
+            mtf_results[f"estimated_rotation_angle_{self.key(mtf_dcm)}"] = rot_ang
+            mtf_results[f"raw_mtf50_{self.key(mtf_dcm)}"] = raw_res
+            mtf_results[f"fitted_mtf50_{self.key(mtf_dcm)}"] = fitted_res
+        except Exception as e:
+            print(f"Could not calculate the spatial resolution for {self.key(mtf_dcm)} because of : {e}")
+            traceback.print_exc(file=sys.stdout)
 
         results = {self.key(self.data[0]): mtf_results}
 
@@ -89,40 +82,6 @@ class ACRSpatialResolution(HazenTask):
 
         return results
 
-    def centroid_com(self, dcm):
-        # Calculate centroid of object using a centre-of-mass calculation
-        thresh_img = dcm > 0.25 * np.max(dcm)
-        open_img = skimage.morphology.area_opening(thresh_img, area_threshold=500)
-        bhull = skimage.morphology.convex_hull_image(open_img)
-        coords = np.nonzero(bhull)  # row major - first array is columns
-
-        sum_x = np.sum(coords[1])
-        sum_y = np.sum(coords[0])
-        cx, cy = sum_x / coords[0].shape[0], sum_y / coords[1].shape[0]
-        cxy = (round(cx), round(cy))
-
-        return bhull, cxy
-
-    def find_rotation(self, dcm):
-        thresh = dcm * (dcm > 0.2 * np.max(dcm))
-
-        kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (5, 5))
-        dilate = cv2.morphologyEx(thresh, cv2.MORPH_DILATE, kernel)
-        diff = cv2.absdiff(dilate, thresh)
-
-        edges = (diff >= 200) * 1
-
-        h, theta, d = skimage.transform.hough_line(edges)
-        accum, angles, dists = skimage.transform.hough_line_peaks(h, theta, d)
-
-        try:
-            angle = np.rad2deg(scipy.stats.mode(angles)[0][0])
-            rot_angle = angle + 90 if angle < 0 else angle - 90
-        except IndexError:
-            rot_angle = 0
-
-        return rot_angle
-
     def y_position_for_ramp(self, res, img, cxy):
         investigate_region = int(np.ceil(5.5 / res[1]).item())
 
@@ -156,8 +115,8 @@ class ACRSpatialResolution(HazenTask):
         edge_sum_rows = np.sum(crop_img, axis=1).astype(np.int_)
         edge_sum_cols = np.sum(crop_img, axis=0).astype(np.int_)
 
-        _, pk_rows_height = find_n_peaks(np.abs(np.diff(edge_sum_rows)), 1)
-        _, pk_cols_height = find_n_peaks(np.abs(np.diff(edge_sum_cols)), 1)
+        _, pk_rows_height = self.ACR_obj.find_n_highest_peaks(np.abs(np.diff(edge_sum_rows)), 1)
+        _, pk_cols_height = self.ACR_obj.find_n_highest_peaks(np.abs(np.diff(edge_sum_cols)), 1)
 
         edge_type = 'vertical' if pk_rows_height > pk_cols_height else 'horizontal'
 
@@ -262,7 +221,7 @@ class ACRSpatialResolution(HazenTask):
 
             return sigmoid
 
-        popt, pcov = scipy.optimize.curve_fit(func, np.arange(1, len(erf) + 1), erf, sigma=np.diag(1 / weights),
+        popt, pcov = scipy.optimize.curve_fit(func, np.arange(1, len(erf) + 1), erf, sigma=(1 / weights),
                                               p0=[np.min(erf), np.max(erf), 0, sum(turning_points) / 2, 1], maxfev=5000)
         erf_fit = func(np.arange(1, len(erf) + 1), popt[0], popt[1], popt[2], popt[3], popt[4])
 
@@ -296,14 +255,7 @@ class ACRSpatialResolution(HazenTask):
     def get_mtf50(self, dcm):
         img = dcm.pixel_array
         res = dcm.PixelSpacing
-        _, cxy = self.centroid_com(img)
-        rot_ang = self.find_rotation(img)
-
-        if np.round(np.abs(rot_ang), 2) < 3:
-            print(f'Rotation angle of the ACR phantom is {np.round(rot_ang, 3)}, which has an absolute is less than 3 '
-                  f'degree. Results will be unreliable!')
-        else:
-            print(f'Rotation angle of the ACR phantom is {np.round(rot_ang, 3)}')
+        cxy = self.ACR_obj.centre
 
         ramp_x, ramp_y = int(cxy[0]), self.y_position_for_ramp(res, img, cxy)
         width = int(13 * img.shape[0] / 256)
@@ -324,54 +276,50 @@ class ACRSpatialResolution(HazenTask):
             import matplotlib.pyplot as plt
             import matplotlib.patches as patches
 
-            fig = plt.figure(figsize=(8, 8))
-
-            gs = fig.add_gridspec(3, 2)
-            ax0 = fig.add_subplot(gs[:2, 0])
-            ax2 = fig.add_subplot(gs[2:, 0])
-            ax3 = fig.add_subplot(gs[0, 1])
-            ax4 = fig.add_subplot(gs[1, 1])
-            ax5 = fig.add_subplot(gs[2, 1])
+            fig, axes = plt.subplots(5, 1)
+            fig.set_size_inches(8, 40)
+            fig.tight_layout(pad=4)
 
-            ax0.imshow(img, interpolation='none')
+            axes[0].imshow(img, interpolation='none')
             rect = patches.Rectangle((ramp_x - width // 2 - 1, ramp_y - width // 2 - 1), width, width, linewidth=1,
                                      edgecolor='w', facecolor='none')
-            ax0.add_patch(rect)
-            ax0.axis('off')
+            axes[0].add_patch(rect)
+            axes[0].axis('off')
+            axes[0].set_title('Segmented Edge')
 
-            ax2.imshow(crop_img)
+            axes[1].imshow(crop_img)
             if edge_type == 'vertical':
-                ax2.plot(np.arange(0, width - 1), np.mean(edge_loc) - slope * np.arange(0, width - 1), color='r')
+                axes[1].plot(np.arange(0, width - 1), np.mean(edge_loc) - slope * np.arange(0, width - 1), color='r')
             else:
-                ax2.plot(np.mean(edge_loc) + slope * np.arange(0, width - 1), np.arange(0, width - 1), color='r')
-            ax2.axis('off')
-            ax2.set_title('Cropped Edge', fontsize=14)
-
-            ax3.plot(erf, 'rx', ms=5)
-            ax3.plot(erf_fit, 'k', lw=3)
-            ax3.set_ylabel('Signal Intensity')
-            ax3.set_xlabel('Pixel')
-            ax3.grid()
-            ax3.set_title('ERF', fontsize=14)
-
-            ax4.plot(lsf_raw, 'rx', ms=5)
-            ax4.plot(lsf_fit, 'k', lw=3)
-            ax4.set_ylabel(r'$\Delta$' + ' Signal Intensity')
-            ax4.set_xlabel('Pixel')
-            ax4.grid()
-            ax4.set_title('LSF', fontsize=14)
-
-            ax5.plot(freq, MTF_raw, 'rx', ms=8, label=f'Raw Data - {round(eff_raw_res, 2)}mm @ 50%')
-            ax5.plot(freq, MTF_fit, 'k', lw=3, label=f'Weighted Sigmoid Fit of ERF - {round(eff_fit_res, 2)}mm @ 50%')
-            ax5.set_xlabel('Spatial Frequency (lp/mm)')
-            ax5.set_ylabel('Modulation Transfer Ratio')
-            ax5.set_xlim([-0.05, 1])
-            ax5.set_ylim([0, 1.05])
-            ax5.grid()
-            ax5.legend(fancybox='true', bbox_to_anchor=[1.05, -0.25, 0, 0])
-            ax5.set_title('MTF', fontsize=14)
-
-            plt.tight_layout()
+                axes[1].plot(np.mean(edge_loc) + slope * np.arange(0, width - 1), np.arange(0, width - 1), color='r')
+            axes[1].axis('off')
+            axes[1].set_title('Cropped Edge', fontsize=14)
+
+            axes[2].plot(erf, 'rx', ms=5, label='Raw Data')
+            axes[2].plot(erf_fit, 'k', lw=3, label='Fitted Data')
+            axes[2].set_ylabel('Signal Intensity')
+            axes[2].set_xlabel('Pixel')
+            axes[2].grid()
+            axes[2].legend(fancybox='true')
+            axes[2].set_title('ERF', fontsize=14)
+
+            axes[3].plot(lsf_raw, 'rx', ms=5, label='Raw Data')
+            axes[3].plot(lsf_fit, 'k', lw=3, label='Fitted Data')
+            axes[3].set_ylabel(r'$\Delta$' + ' Signal Intensity')
+            axes[3].set_xlabel('Pixel')
+            axes[3].grid()
+            axes[3].legend(fancybox='true')
+            axes[3].set_title('LSF', fontsize=14)
+
+            axes[4].plot(freq, MTF_raw, 'rx', ms=8, label=f'Raw Data - {round(eff_raw_res, 2)}mm @ 50%')
+            axes[4].plot(freq, MTF_fit, 'k', lw=3, label=f'Weighted Sigmoid Fit of ERF - {round(eff_fit_res, 2)}mm @ 50%')
+            axes[4].set_xlabel('Spatial Frequency (lp/mm)')
+            axes[4].set_ylabel('Modulation Transfer Ratio')
+            axes[4].set_xlim([-0.05, 1])
+            axes[4].set_ylim([0, 1.05])
+            axes[4].grid()
+            axes[4].legend(fancybox='true')
+            axes[4].set_title('MTF', fontsize=14)
 
             img_path = os.path.realpath(os.path.join(self.report_path, f'{self.key(dcm)}.png'))
             fig.savefig(img_path)
diff --git a/hazenlib/tasks/acr_uniformity.py b/hazenlib/tasks/acr_uniformity.py
index 4d3250d..54354ed 100644
--- a/hazenlib/tasks/acr_uniformity.py
+++ b/hazenlib/tasks/acr_uniformity.py
@@ -3,9 +3,9 @@ ACR Uniformity
 
 https://www.acraccreditation.org/-/media/acraccreditation/documents/mri/largephantomguidance.pdf
 
-Calculates uniformity for slice 7 of the ACR phantom.
+Calculates the percentage integral uniformity for slice 7 of the ACR phantom.
 
-This script calculates the integral uniformity in accordance with the ACR Guidance.
+This script calculates the percentage integral uniformity in accordance with the ACR Guidance.
 This is done by first defining a large 200cm2 ROI before placing 1cm2 ROIs at every pixel within
 the large ROI. At each point, the mean of the 1cm2 ROI is calculated. The ROIs with the maximum and
 minimum mean value are used to calculate the integral uniformity. The results are also visualised.
@@ -20,34 +20,29 @@ import sys
 import traceback
 import os
 import numpy as np
-import skimage.morphology
 
 from hazenlib.HazenTask import HazenTask
+from hazenlib.ACRObject import ACRObject
 
 
 class ACRUniformity(HazenTask):
 
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
+        self.ACR_obj = None
 
     def run(self) -> dict:
         results = {}
-        z = []
-        for dcm in self.data:
-            z.append(dcm.ImagePositionPatient[2])
+        self.ACR_obj = ACRObject(self.data)
+        uniformity_dcm = self.ACR_obj.dcm[6]
 
-        idx_sort = np.argsort(z)
-
-        for dcm in self.data:
-            if dcm.ImagePositionPatient[2] == z[idx_sort[6]]:
-                try:
-                    result = self.get_integral_uniformity(dcm)
-                except Exception as e:
-                    print(f"Could not calculate the percent integral uniformity for {self.key(dcm)} because of : {e}")
-                    traceback.print_exc(file=sys.stdout)
-                    continue
-
-                results[self.key(dcm)] = result
+        try:
+            result = self.get_integral_uniformity(uniformity_dcm)
+            results[self.key(uniformity_dcm)] = result
+        except Exception as e:
+            print(
+                f"Could not calculate the percent integral uniformity for {self.key(uniformity_dcm)} because of : {e}")
+            traceback.print_exc(file=sys.stdout)
 
         # only return reports if requested
         if self.report:
@@ -55,29 +50,6 @@ class ACRUniformity(HazenTask):
 
         return results
 
-    def centroid_com(self, dcm):
-        # Calculate centroid of object using a centre-of-mass calculation
-        thresh_img = dcm > 0.25 * np.max(dcm)
-        open_img = skimage.morphology.area_opening(thresh_img, area_threshold=500)
-        bhull = skimage.morphology.convex_hull_image(open_img)
-        coords = np.nonzero(bhull)  # row major - first array is columns
-
-        sum_x = np.sum(coords[1])
-        sum_y = np.sum(coords[0])
-        cxy = sum_x / coords[0].shape, sum_y / coords[1].shape
-
-        cxy = [cxy[0].astype(int), cxy[1].astype(int)]
-        return cxy
-
-    def circular_mask(self, centre, radius, dims):
-        # Define a circular logical mask
-        nx = np.linspace(1, dims[0], dims[0])
-        ny = np.linspace(1, dims[1], dims[1])
-
-        x, y = np.meshgrid(nx, ny)
-        mask = np.square(x - centre[0]) + np.square(y - centre[1]) <= radius ** 2
-        return mask
-
     def get_integral_uniformity(self, dcm):
         # Calculate the integral uniformity in accordance with ACR guidance.
         img = dcm.pixel_array
@@ -87,54 +59,75 @@ class ACRUniformity(HazenTask):
         d_void = np.ceil(5 / res[0]).astype(int)  # Offset distance for rectangular void at top of phantom
         dims = img.shape  # Dimensions of image
 
-        cxy = self.centroid_com(img)
-        base_mask = self.circular_mask([cxy[0], cxy[1] + d_void], r_small, dims)  # Dummy circular mask at centroid
+        cxy = self.ACR_obj.centre
+        base_mask = ACRObject.circular_mask((cxy[0], cxy[1] + d_void), r_small, dims)  # Dummy circular mask at
+        # centroid
         coords = np.nonzero(base_mask)  # Coordinates of mask
 
-        lroi = self.circular_mask([cxy[0], cxy[1] + d_void], r_large, dims)
+        lroi = self.ACR_obj.circular_mask([cxy[0], cxy[1] + d_void], r_large, dims)
         img_masked = lroi * img
+        half_max = np.percentile(img_masked[np.nonzero(img_masked)], 50)
+
+        min_image = img_masked * (img_masked < half_max)
+        max_image = img_masked * (img_masked > half_max)
 
-        lroi_extent = np.nonzero(lroi)
+        min_rows, min_cols = np.nonzero(min_image)[0], np.nonzero(min_image)[1]
+        max_rows, max_cols = np.nonzero(max_image)[0], np.nonzero(max_image)[1]
 
-        mean_val = np.zeros(lroi_extent[0].shape)
         mean_array = np.zeros(img_masked.shape)
 
-        for ii in range(0, len(lroi_extent[0])):
-            centre = [lroi_extent[0][ii], lroi_extent[1][ii]]  # Extract coordinates of new mask centre within large ROI
-            trans_mask = [coords[0] + centre[0] - cxy[0] - d_void,
-                          coords[1] + centre[1] - cxy[1]]  # Translate mask within the limits of the large ROI
-            sroi_val = img_masked[trans_mask[0], trans_mask[1]]  # Extract values within translated mask
-            if np.count_nonzero(sroi_val) < np.count_nonzero(base_mask):
-                mean_val[ii] = 0
-            else:
-                mean_val[ii] = np.mean(sroi_val[np.nonzero(sroi_val)])
-            mean_array[lroi_extent[0][ii], lroi_extent[1][ii]] = mean_val[ii]
+        def uniformity_iterator(masked_image, sample_mask, rows, cols):
+            coords = np.nonzero(sample_mask)  # Coordinates of mask
+            for idx, (row, col) in enumerate(zip(rows, cols)):
+                centre = [row, col]
+                translate_mask = [coords[0] + centre[0] - cxy[0] - d_void,
+                                  coords[1] + centre[1] - cxy[1]]
+                values = masked_image[translate_mask[0], translate_mask[1]]
+                if np.count_nonzero(values) < np.count_nonzero(sample_mask):
+                    mean_val = 0
+                else:
+                    mean_val = np.mean(values[np.nonzero(values)])
 
-        sig_max = np.max(mean_val)
-        sig_min = np.min(mean_val[np.nonzero(mean_val)])
+                mean_array[row, col] = mean_val
 
-        max_loc = np.where(mean_array == sig_max)
-        min_loc = np.where(mean_array == sig_min)
+            return mean_array
+
+        min_data = uniformity_iterator(min_image, base_mask, min_rows, min_cols)
+        max_data = uniformity_iterator(max_image, base_mask, max_rows, max_cols)
+
+        sig_max = np.max(max_data)
+        sig_min = np.min(min_data[np.nonzero(min_data)])
+
+        max_loc = np.where(max_data == sig_max)
+        min_loc = np.where(min_data == sig_min)
 
         piu = 100 * (1 - (sig_max - sig_min) / (sig_max + sig_min))
 
-        piu = np.round(piu,2)
+        piu = np.round(piu, 2)
         if self.report:
             import matplotlib.pyplot as plt
+            fig, axes = plt.subplots(2, 1)
+            fig.set_size_inches(8, 16)
+            fig.tight_layout(pad=4)
+
             theta = np.linspace(0, 2 * np.pi, 360)
-            fig = plt.figure()
-            fig.set_size_inches(8, 8)
-            plt.imshow(img)
-
-            plt.scatter([max_loc[1], min_loc[1]], [max_loc[0], min_loc[0]], c='red', marker='x')
-            plt.plot(r_small * np.cos(theta) + max_loc[1], r_small * np.sin(theta) + max_loc[0], c='yellow')
-            plt.annotate('Min = ' + str(np.round(sig_min, 1)), [min_loc[1], min_loc[0] + 10 / res[0]], c='white')
-
-            plt.plot(r_small * np.cos(theta) + min_loc[1], r_small * np.sin(theta) + min_loc[0], c='yellow')
-            plt.annotate('Max = ' + str(np.round(sig_max, 1)), [max_loc[1], max_loc[0] + 10 / res[0]], c='white')
-            plt.plot(r_large * np.cos(theta) + cxy[1], r_large * np.sin(theta) + cxy[0] + 5 / res[1], c='black')
-            plt.axis('off')
-            plt.title('Percent Integral Uniformity = ' + str(np.round(piu, 2)) + '%')
+
+            axes[0].imshow(img)
+            axes[0].scatter(cxy[0], cxy[1], c='red')
+            axes[0].axis('off')
+            axes[0].set_title('Centroid Location')
+
+            axes[1].imshow(img)
+            axes[1].scatter([max_loc[1], min_loc[1]], [max_loc[0], min_loc[0]], c='red', marker='x')
+            axes[1].plot(r_small * np.cos(theta) + max_loc[1], r_small * np.sin(theta) + max_loc[0], c='yellow')
+            axes[1].annotate('Min = ' + str(np.round(sig_min, 1)), [min_loc[1], min_loc[0] + 10 / res[0]], c='white')
+
+            axes[1].plot(r_small * np.cos(theta) + min_loc[1], r_small * np.sin(theta) + min_loc[0], c='yellow')
+            axes[1].annotate('Max = ' + str(np.round(sig_max, 1)), [max_loc[1], max_loc[0] + 10 / res[0]], c='white')
+            axes[1].plot(r_large * np.cos(theta) + cxy[1], r_large * np.sin(theta) + cxy[0] + 5 / res[1], c='black')
+            axes[1].axis('off')
+            axes[1].set_title('Percent Integral Uniformity = ' + str(np.round(piu, 2)) + '%')
+
             img_path = os.path.realpath(os.path.join(self.report_path, f'{self.key(dcm)}.png'))
             fig.savefig(img_path)
             self.report_files.append(img_path)


Test patch:
diff --git a/tests/test_ACRObject.py b/tests/test_ACRObject.py
new file mode 100644
index 0000000..ecb8e35
--- /dev/null
+++ b/tests/test_ACRObject.py
@@ -0,0 +1,37 @@
+import os
+import unittest
+import pathlib
+import pydicom
+import numpy as np
+
+from hazenlib import HazenTask
+from hazenlib.ACRObject import ACRObject
+from tests import TEST_DATA_DIR, TEST_REPORT_DIR
+
+
+class TestACRTools(unittest.TestCase):
+    rotation = [-1.0, 0.0]
+    centre = [(129, 130), (253, 255)]
+    test_point = (-60.98, -45.62)
+
+    def setUp(self):
+        self.Siemens_data = [pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens', f'{i}')) for i in
+                             os.listdir(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens'))]
+        self.GE_data = [pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'GE', f'{i}')) for i in
+                        os.listdir(os.path.join(TEST_DATA_DIR, 'acr', 'GE'))]
+
+        self.Siemens_ACR_obj = ACRObject(self.Siemens_data)
+        self.GE_ACR_obj = ACRObject(self.GE_data)
+
+    def test_find_rotation(self):
+        assert self.rotation[0] == np.round(self.Siemens_ACR_obj.determine_rotation(), 1)
+        assert self.rotation[1] == np.round(self.GE_ACR_obj.determine_rotation(), 1)
+
+    def test_find_centre(self):
+        assert (self.centre[0] == np.round(self.Siemens_ACR_obj.centre, 1)).all() == True
+        assert (self.centre[1] == np.round(self.GE_ACR_obj.centre, 1)).all() == True
+
+    def test_rotate_point(self):
+        rotated_point = np.array(self.Siemens_ACR_obj.rotate_point((0, 0), (30, 70), 150))
+        rotated_point = np.round(rotated_point, 2)
+        assert (rotated_point == self.test_point).all() == True
diff --git a/tests/test_acr_geometric_accuracy.py b/tests/test_acr_geometric_accuracy.py
index b89aadc..ed8ff80 100644
--- a/tests/test_acr_geometric_accuracy.py
+++ b/tests/test_acr_geometric_accuracy.py
@@ -5,28 +5,27 @@ import pydicom
 import numpy as np
 
 from hazenlib.tasks.acr_geometric_accuracy import ACRGeometricAccuracy
+from hazenlib.ACRObject import ACRObject
 from tests import TEST_DATA_DIR, TEST_REPORT_DIR
 
 
 class TestACRGeometricAccuracySiemens(unittest.TestCase):
     ACR_GEOMETRIC_ACCURACY_DATA = pathlib.Path(TEST_DATA_DIR / 'acr')
-    centre = (128, 129)
-    L1 = 190.43, 186.52
-    L5 = 190.43, 186.52, 189.45, 191.41
-    test_point = (-60.98, -45.62)
+    L1 = 192.38, 188.48
+    L5 = 192.38, 188.48, 190.43, 192.38
 
     def setUp(self):
         self.acr_geometric_accuracy_task = ACRGeometricAccuracy(data_paths=[os.path.join(TEST_DATA_DIR, 'acr')],
                                                                 report_dir=pathlib.PurePath.joinpath(TEST_REPORT_DIR))
-        self.dcm = pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens', '0.dcm'))
-        self.dcm2 = pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens', '4.dcm'))
+        self.acr_geometric_accuracy_task.ACR_obj = ACRObject(
+            [pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens', f'{i}')) for i in
+             os.listdir(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens'))])
 
-    def test_object_centre(self):
-        data = self.dcm.pixel_array
-        assert self.acr_geometric_accuracy_task.centroid_com(data)[1] == self.centre
+        self.dcm_1 = self.acr_geometric_accuracy_task.ACR_obj.dcm[0]
+        self.dcm_5 = self.acr_geometric_accuracy_task.ACR_obj.dcm[4]
 
-    def test_geo_accuracy_slice1(self):
-        slice1_vals = np.array(self.acr_geometric_accuracy_task.get_geometric_accuracy_slice1(self.dcm))
+    def test_geometric_accuracy_slice_1(self):
+        slice1_vals = np.array(self.acr_geometric_accuracy_task.get_geometric_accuracy_slice1(self.dcm_1))
         slice1_vals = np.round(slice1_vals, 2)
 
         print("\ntest_geo_accuracy.py::TestGeoAccuracy::test_geo_accuracy_slice1")
@@ -35,8 +34,8 @@ class TestACRGeometricAccuracySiemens(unittest.TestCase):
 
         assert (slice1_vals == self.L1).all() == True
 
-    def test_geo_accuracy_slice5(self):
-        slice5_vals = np.array(self.acr_geometric_accuracy_task.get_geometric_accuracy_slice5(self.dcm2))
+    def test_geometric_accuracy_slice_5(self):
+        slice5_vals = np.array(self.acr_geometric_accuracy_task.get_geometric_accuracy_slice5(self.dcm_5))
         slice5_vals = np.round(slice5_vals, 2)
 
         print("\ntest_geo_accuracy.py::TestGeoAccuracy::test_geo_accuracy_slice1")
@@ -44,38 +43,36 @@ class TestACRGeometricAccuracySiemens(unittest.TestCase):
         print("fixed value:", self.L5)
         assert (slice5_vals == self.L5).all() == True
 
-    def test_rotate_point(self):
-        rotated_point = np.array(self.acr_geometric_accuracy_task.rotate_point((0, 0), (30, 70), 150))
-        rotated_point = np.round(rotated_point, 2)
-        print(rotated_point)
-        assert (rotated_point == self.test_point).all() == True
 
-
-# class TestACRUniformityPhilips(unittest.TestCase):
+# TODO: Add unit tests for Philips datasets.
 
 class TestACRGeometricAccuracyGE(unittest.TestCase):
     ACR_GEOMETRIC_ACCURACY_DATA = pathlib.Path(TEST_DATA_DIR / 'acr')
-    L1 = 189.92, 187.89
-    L5 = 189.92, 188.39, 190.43, 189.92
-    distortion_metrics = [-0.59, 2.11, 0.49]
+    L1 = 191.44, 191.44
+    L5 = 191.44, 191.44, 191.44, 189.41
+    distortion_metrics = [1.1, 1.44, 0.4]
 
     def setUp(self):
         self.acr_geometric_accuracy_task = ACRGeometricAccuracy(data_paths=[os.path.join(TEST_DATA_DIR, 'acr')],
                                                                 report_dir=pathlib.PurePath.joinpath(TEST_REPORT_DIR))
-        self.dcm = pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'GE', '10.dcm'))
-        self.dcm2 = pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'GE', '6.dcm'))
+        self.acr_geometric_accuracy_task.ACR_obj = ACRObject(
+            [pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'GE', f'{i}')) for i in
+             os.listdir(os.path.join(TEST_DATA_DIR, 'acr', 'GE'))])
+
+        self.dcm_1 = self.acr_geometric_accuracy_task.ACR_obj.dcm[0]
+        self.dcm_5 = self.acr_geometric_accuracy_task.ACR_obj.dcm[4]
 
     def test_geo_accuracy_slice1(self):
-        slice1_vals = np.array(self.acr_geometric_accuracy_task.get_geometric_accuracy_slice1(self.dcm))
+        slice1_vals = np.array(self.acr_geometric_accuracy_task.get_geometric_accuracy_slice1(self.dcm_1))
         slice1_vals = np.round(slice1_vals, 2)
         assert (slice1_vals == self.L1).all() == True
 
     def test_geo_accuracy_slice5(self):
-        slice5_vals = np.array(self.acr_geometric_accuracy_task.get_geometric_accuracy_slice5(self.dcm2))
+        slice5_vals = np.array(self.acr_geometric_accuracy_task.get_geometric_accuracy_slice5(self.dcm_5))
         slice5_vals = np.round(slice5_vals, 2)
         assert (slice5_vals == self.L5).all() == True
 
     def test_distortion_metrics(self):
-        metrics = np.array(self.acr_geometric_accuracy_task.distortion_metric(self.L1+self.L5))
+        metrics = np.array(self.acr_geometric_accuracy_task.distortion_metric(self.L1 + self.L5))
         metrics = np.round(metrics, 2)
         assert (metrics == self.distortion_metrics).all() == True
diff --git a/tests/test_acr_ghosting.py b/tests/test_acr_ghosting.py
index b011d28..96e4798 100644
--- a/tests/test_acr_ghosting.py
+++ b/tests/test_acr_ghosting.py
@@ -4,21 +4,24 @@ import pathlib
 import pydicom
 
 from hazenlib.tasks.acr_ghosting import ACRGhosting
+from hazenlib.ACRObject import ACRObject
 from tests import TEST_DATA_DIR, TEST_REPORT_DIR
 
 
 class TestACRGhostingSiemens(unittest.TestCase):
     ACR_GHOSTING_DATA = pathlib.Path(TEST_DATA_DIR / 'acr')
     centre = [129, 128]
-    psg = 0.056
+    psg = 0.035
 
     def setUp(self):
         self.acr_ghosting_task = ACRGhosting(data_paths=[os.path.join(TEST_DATA_DIR, 'acr')],
                                              report_dir=pathlib.PurePath.joinpath(TEST_REPORT_DIR))
-        self.dcm = pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens', '6.dcm'))
 
-    def test_object_centre(self):
-        assert self.acr_ghosting_task.centroid_com(self.dcm.pixel_array)[1] == self.centre
+        self.acr_ghosting_task.ACR_obj = ACRObject(
+            [pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens', f'{i}')) for i in
+             os.listdir(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens'))])
+
+        self.dcm = self.acr_ghosting_task.ACR_obj.dcm[6]
 
     def test_ghosting(self):
         ghosting_val = round(self.acr_ghosting_task.get_signal_ghosting(self.dcm), 3)
@@ -33,15 +36,16 @@ class TestACRGhostingSiemens(unittest.TestCase):
 class TestACRGhostingGE(unittest.TestCase):
     ACR_GHOSTING_DATA = pathlib.Path(TEST_DATA_DIR / 'acr')
     centre = [253, 256]
-    psg = 0.487
+    psg = 0.471
 
     def setUp(self):
         self.acr_ghosting_task = ACRGhosting(data_paths=[os.path.join(TEST_DATA_DIR, 'acr')],
                                              report_dir=pathlib.PurePath.joinpath(TEST_REPORT_DIR))
-        self.dcm = pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'GE', '4.dcm'))
+        self.acr_ghosting_task.ACR_obj = ACRObject(
+            [pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'GE', f'{i}')) for i in
+             os.listdir(os.path.join(TEST_DATA_DIR, 'acr', 'GE'))])
 
-    def test_object_centre(self):
-        assert self.acr_ghosting_task.centroid_com(self.dcm.pixel_array)[1] == self.centre
+        self.dcm = self.acr_ghosting_task.ACR_obj.dcm[6]
 
     def test_ghosting(self):
         assert round(self.acr_ghosting_task.get_signal_ghosting(self.dcm), 3) == self.psg
diff --git a/tests/test_acr_slice_position.py b/tests/test_acr_slice_position.py
index e3842f1..d57160d 100644
--- a/tests/test_acr_slice_position.py
+++ b/tests/test_acr_slice_position.py
@@ -4,41 +4,44 @@ import pathlib
 import pydicom
 
 from hazenlib.tasks.acr_slice_position import ACRSlicePosition
+from hazenlib.ACRObject import ACRObject
 from tests import TEST_DATA_DIR
 
 
 class TestACRSlicePositionSiemens(unittest.TestCase):
     ACR_SLICE_POSITION_DATA = pathlib.Path(TEST_DATA_DIR / 'acr')
-    centre = [128, 129]
     x_pts = [(123, 129), (123, 129)]
-    y_pts = [(40, 83), (44, 82)]
+    y_pts = [(40, 82), (44, 82)]
     dL = -0.59, -1.56
 
     def setUp(self):
         self.acr_slice_position_task = ACRSlicePosition(data_paths=[os.path.join(TEST_DATA_DIR, 'acr')])
-        self.dcm_1 = pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens', '0.dcm'))
-        self.dcm_11 = pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens', '10.dcm'))
+        self.acr_slice_position_task.ACR_obj = ACRObject(
+            [pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens', f'{i}')) for i in
+             os.listdir(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens'))])
 
-    def test_object_centre(self):
-        assert self.acr_slice_position_task.centroid_com(self.dcm_1.pixel_array)[1] == self.centre
+        self.dcm_1 = self.acr_slice_position_task.ACR_obj.dcm[0]
+        self.dcm_11 = self.acr_slice_position_task.ACR_obj.dcm[-1]
 
     def test_wedge_find(self):
         # IMAGE 1
+        img = self.dcm_1.pixel_array
         res = self.dcm_1.PixelSpacing
-        mask, _ = self.acr_slice_position_task.centroid_com(self.dcm_1.pixel_array)
-        assert (self.acr_slice_position_task.find_wedges(self.dcm_1.pixel_array, mask, res)[0] ==
+        mask = self.acr_slice_position_task.ACR_obj.mask_image(img)
+        assert (self.acr_slice_position_task.find_wedges(img, mask, res)[0] ==
                 self.x_pts[0]).all() == True
 
-        assert (self.acr_slice_position_task.find_wedges(self.dcm_1.pixel_array, mask, res)[1] ==
+        assert (self.acr_slice_position_task.find_wedges(img, mask, res)[1] ==
                 self.y_pts[0]).all() == True
 
         # IMAGE 11
+        img = self.dcm_11.pixel_array
         res = self.dcm_11.PixelSpacing
-        mask, _ = self.acr_slice_position_task.centroid_com(self.dcm_11.pixel_array)
-        assert (self.acr_slice_position_task.find_wedges(self.dcm_11.pixel_array, mask, res)[0] ==
+        mask = self.acr_slice_position_task.ACR_obj.mask_image(img)
+        assert (self.acr_slice_position_task.find_wedges(img, mask, res)[0] ==
                 self.x_pts[1]).all() == True
 
-        assert (self.acr_slice_position_task.find_wedges(self.dcm_11.pixel_array, mask, res)[1] ==
+        assert (self.acr_slice_position_task.find_wedges(img, mask, res)[1] ==
                 self.y_pts[1]).all() == True
 
     def test_slice_position(self):
@@ -55,36 +58,38 @@ class TestACRSlicePositionSiemens(unittest.TestCase):
 
 class TestACRSlicePositionGE(unittest.TestCase):
     ACR_SLICE_POSITION_DATA = pathlib.Path(TEST_DATA_DIR / 'acr')
-    centre = [253, 257]
     x_pts = [(246, 257), (246, 257)]
-    y_pts = [(82, 163), (85, 165)]
-    dL = 0.3, 0.41
+    y_pts = [(77, 164), (89, 162)]
+    dL = 0.41, 0.3
 
     def setUp(self):
         self.acr_slice_position_task = ACRSlicePosition(data_paths=[os.path.join(TEST_DATA_DIR, 'acr')])
-        self.dcm_1 = pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'GE', '0.dcm'))
-        self.dcm_11 = pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'GE', '10.dcm'))
+        self.acr_slice_position_task.ACR_obj = ACRObject(
+            [pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'GE', f'{i}')) for i in
+             os.listdir(os.path.join(TEST_DATA_DIR, 'acr', 'GE'))])
 
-    def test_object_centre(self):
-        assert self.acr_slice_position_task.centroid_com(self.dcm_1.pixel_array)[1] == self.centre
+        self.dcm_1 = self.acr_slice_position_task.ACR_obj.dcm[0]
+        self.dcm_11 = self.acr_slice_position_task.ACR_obj.dcm[-1]
 
     def test_wedge_find(self):
         # IMAGE 1
+        img = self.dcm_1.pixel_array
         res = self.dcm_1.PixelSpacing
-        mask, _ = self.acr_slice_position_task.centroid_com(self.dcm_1.pixel_array)
-        assert (self.acr_slice_position_task.find_wedges(self.dcm_1.pixel_array, mask, res)[0] ==
+        mask = self.acr_slice_position_task.ACR_obj.mask_image(img)
+        assert (self.acr_slice_position_task.find_wedges(img, mask, res)[0] ==
                 self.x_pts[0]).all() == True
 
-        assert (self.acr_slice_position_task.find_wedges(self.dcm_1.pixel_array, mask, res)[1] ==
+        assert (self.acr_slice_position_task.find_wedges(img, mask, res)[1] ==
                 self.y_pts[0]).all() == True
 
         # IMAGE 11
+        img = self.dcm_11.pixel_array
         res = self.dcm_11.PixelSpacing
-        mask, _ = self.acr_slice_position_task.centroid_com(self.dcm_11.pixel_array)
-        assert (self.acr_slice_position_task.find_wedges(self.dcm_11.pixel_array, mask, res)[0] ==
+        mask = self.acr_slice_position_task.ACR_obj.mask_image(img)
+        assert (self.acr_slice_position_task.find_wedges(img, mask, res)[0] ==
                 self.x_pts[1]).all() == True
 
-        assert (self.acr_slice_position_task.find_wedges(self.dcm_11.pixel_array, mask, res)[1] ==
+        assert (self.acr_slice_position_task.find_wedges(img, mask, res)[1] ==
                 self.y_pts[1]).all() == True
 
     def test_slice_position(self):
diff --git a/tests/test_acr_slice_thickness.py b/tests/test_acr_slice_thickness.py
index 687f234..04b7089 100644
--- a/tests/test_acr_slice_thickness.py
+++ b/tests/test_acr_slice_thickness.py
@@ -4,26 +4,27 @@ import pathlib
 import pydicom
 
 from hazenlib.tasks.acr_slice_thickness import ACRSliceThickness
+from hazenlib.ACRObject import ACRObject
 from tests import TEST_DATA_DIR
 
 
 class TestACRSliceThicknessSiemens(unittest.TestCase):
     ACR_SLICE_POSITION_DATA = pathlib.Path(TEST_DATA_DIR / 'acr')
-    centre = [128, 129]
     x_pts = [71, 181]
     y_pts = [132, 126]
     dz = 4.91
 
     def setUp(self):
         self.acr_slice_thickness_task = ACRSliceThickness(data_paths=[os.path.join(TEST_DATA_DIR, 'acr')])
-        self.dcm = pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens', '0.dcm'))
+        self.acr_slice_thickness_task.ACR_obj = ACRObject(
+            [pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens', f'{i}')) for i in
+             os.listdir(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens'))])
 
-    def test_object_centre(self):
-        assert self.acr_slice_thickness_task.centroid_com(self.dcm.pixel_array)[1] == self.centre
+        self.dcm = self.acr_slice_thickness_task.ACR_obj.dcm[0]
 
     def test_ramp_find(self):
         res = self.dcm.PixelSpacing
-        _, centre = self.acr_slice_thickness_task.centroid_com(self.dcm.pixel_array)
+        centre = self.acr_slice_thickness_task.ACR_obj.centre
         assert (self.acr_slice_thickness_task.find_ramps(self.dcm.pixel_array, centre, res)[0] ==
                 self.x_pts).all() == True
 
@@ -42,24 +43,23 @@ class TestACRSliceThicknessSiemens(unittest.TestCase):
 
 class TestACRSliceThicknessGE(unittest.TestCase):
     ACR_SLICE_POSITION_DATA = pathlib.Path(TEST_DATA_DIR / 'acr')
-    centre = [253, 255]
     x_pts = [146, 356]
     y_pts = [262, 250]
     dz = 5.02
 
     def setUp(self):
         self.acr_slice_thickness_task = ACRSliceThickness(data_paths=[os.path.join(TEST_DATA_DIR, 'acr')])
-        self.dcm = pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'GE', '10.dcm'))
+        self.acr_slice_thickness_task.ACR_obj = ACRObject(
+            [pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'GE', f'{i}')) for i in
+             os.listdir(os.path.join(TEST_DATA_DIR, 'acr', 'GE'))])
 
-    def test_object_centre(self):
-        assert self.acr_slice_thickness_task.centroid_com(self.dcm.pixel_array)[1] == self.centre
+        self.dcm = self.acr_slice_thickness_task.ACR_obj.dcm[0]
 
     def test_ramp_find(self):
         res = self.dcm.PixelSpacing
-        _, centre = self.acr_slice_thickness_task.centroid_com(self.dcm.pixel_array)
+        centre = self.acr_slice_thickness_task.ACR_obj.centre
         assert (self.acr_slice_thickness_task.find_ramps(self.dcm.pixel_array, centre, res)[0] ==
                 self.x_pts).all() == True
-
         assert (self.acr_slice_thickness_task.find_ramps(self.dcm.pixel_array, centre, res)[1] ==
                 self.y_pts).all() == True
 
diff --git a/tests/test_acr_snr.py b/tests/test_acr_snr.py
index fa90434..6d098d9 100644
--- a/tests/test_acr_snr.py
+++ b/tests/test_acr_snr.py
@@ -4,35 +4,38 @@ import pathlib
 import pydicom
 
 from hazenlib.tasks.acr_snr import ACRSNR
+from hazenlib.ACRObject import ACRObject
 from tests import TEST_DATA_DIR, TEST_REPORT_DIR
 
 
 class TestACRSNRSiemens(unittest.TestCase):
     ACR_SNR_DATA = pathlib.Path(TEST_DATA_DIR / 'acr')
-    centre = [129, 128]
     norm_factor = 9.761711312090041
-    snr = 335.27
-    sub_snr = 76.48
+    snr = 344.15
+    sub_snr = 75.94
 
     def setUp(self):
         self.acr_snr_task = ACRSNR(data_paths=[os.path.join(TEST_DATA_DIR, 'acr')],
                                    report_dir=pathlib.PurePath.joinpath(TEST_REPORT_DIR))
-        self.dcm = pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens', '6.dcm'))
-        self.dcm2 = pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens2', '7.dcm'))
+        self.acr_snr_task.ACR_obj = [ACRObject(
+            [pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens', f'{i}')) for i in
+             os.listdir(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens'))])]
+        self.acr_snr_task.ACR_obj.append(
+            ACRObject([pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens2', f'{i}')) for i in
+                       os.listdir(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens2'))]))
 
-    def test_object_centre(self):
-        assert self.acr_snr_task.centroid(self.dcm)[1] == self.centre
+        self.dcm = [i.dcm[6] for i in self.acr_snr_task.ACR_obj]
 
     def test_normalisation_factor(self):
-        SNR_factor = self.acr_snr_task.get_normalised_snr_factor(self.dcm)
+        SNR_factor = self.acr_snr_task.get_normalised_snr_factor(self.dcm[0])
         assert SNR_factor == self.norm_factor
 
     def test_snr_by_smoothing(self):
-        snr, _ = self.acr_snr_task.snr_by_smoothing(self.dcm)
+        snr, _ = self.acr_snr_task.snr_by_smoothing(self.dcm[0])
         assert round(snr, 2) == self.snr
 
     def test_snr_by_subtraction(self):
-        snr, _ = self.acr_snr_task.snr_by_subtraction(self.dcm, self.dcm2)
+        snr, _ = self.acr_snr_task.snr_by_subtraction(self.dcm[0], self.dcm[1])
         rounded_snr = round(snr, 2)
 
         print("\ntest_snr_by_subtraction.py::TestSnrBySubtraction::test_snr_by_subtraction")
@@ -44,17 +47,17 @@ class TestACRSNRSiemens(unittest.TestCase):
 
 class TestACRSNRGE(unittest.TestCase):
     ACR_SNR_DATA = pathlib.Path(TEST_DATA_DIR / 'acr')
-    centre = [253, 256]
     norm_factor = 57.12810400630368
-    snr = 39.97
+    snr = 40.19
 
     def setUp(self):
         self.acr_snr_task = ACRSNR(data_paths=[os.path.join(TEST_DATA_DIR, 'acr')],
                                    report_dir=pathlib.PurePath.joinpath(TEST_REPORT_DIR))
-        self.dcm = pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'GE', '4.dcm'))
+        self.acr_snr_task.ACR_obj = [ACRObject(
+            [pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'GE', f'{i}')) for i in
+             os.listdir(os.path.join(TEST_DATA_DIR, 'acr', 'GE'))])]
 
-    def test_object_centre(self):
-        assert self.acr_snr_task.centroid(self.dcm)[1] == self.centre
+        self.dcm = self.acr_snr_task.ACR_obj[0].dcm[6]
 
     def test_normalisation_factor(self):
         SNR_factor = self.acr_snr_task.get_normalised_snr_factor(self.dcm)
diff --git a/tests/test_acr_spatial_resolution.py b/tests/test_acr_spatial_resolution.py
index 56adce2..6b2ee64 100644
--- a/tests/test_acr_spatial_resolution.py
+++ b/tests/test_acr_spatial_resolution.py
@@ -5,6 +5,7 @@ import pydicom
 import numpy as np
 
 from hazenlib.tasks.acr_spatial_resolution import ACRSpatialResolution
+from hazenlib.ACRObject import ACRObject
 from tests import TEST_DATA_DIR, TEST_REPORT_DIR
 
 
@@ -17,23 +18,19 @@ class TestACRSpatialResolutionSiemens(unittest.TestCase):
     edge_type = 'vertical', 'downward'
     edge_loc = [5, 7]
     slope = -0.165
-    MTF50 = (1.16, 1.32)
+    MTF50 = (1.18, 1.35)
 
     def setUp(self):
         self.acr_spatial_resolution_task = ACRSpatialResolution(data_paths=[os.path.join(TEST_DATA_DIR, 'acr')],
                                                                 report_dir=pathlib.PurePath.joinpath(TEST_REPORT_DIR))
-        self.dcm = pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'SiemensMTF', '0.dcm'))
+        self.acr_spatial_resolution_task.ACR_obj = ACRObject(
+            [pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'SiemensMTF', f'{i}')) for i in
+             os.listdir(os.path.join(TEST_DATA_DIR, 'acr', 'SiemensMTF'))])
+
+        self.dcm = self.acr_spatial_resolution_task.ACR_obj.dcm[0]
         self.crop_image = self.acr_spatial_resolution_task.crop_image(self.dcm.pixel_array, self.centre[0],
                                                                       self.y_ramp_pos, self.width)
 
-    def test_object_centre(self):
-        data = self.dcm.pixel_array
-        assert self.acr_spatial_resolution_task.centroid_com(data)[1] == self.centre
-
-    def test_rotation(self):
-        data = self.dcm.pixel_array
-        assert self.acr_spatial_resolution_task.find_rotation(data) == self.rotation_angle
-
     def test_find_y_ramp(self):
         data = self.dcm.pixel_array
         res = self.dcm.PixelSpacing
@@ -68,23 +65,19 @@ class TestACRSpatialResolutionGE(unittest.TestCase):
     edge_type = 'vertical', 'upward'
     edge_loc = [5, 7]
     slope = 0.037
-    MTF50 = (0.71, 0.69)
+    MTF50 = (0.72, 0.71)
 
     def setUp(self):
         self.acr_spatial_resolution_task = ACRSpatialResolution(data_paths=[os.path.join(TEST_DATA_DIR, 'acr')],
                                                                 report_dir=pathlib.PurePath.joinpath(TEST_REPORT_DIR))
-        self.dcm = pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'GE', '10.dcm'))
+        self.acr_spatial_resolution_task.ACR_obj = ACRObject(
+            [pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'GE', f'{i}')) for i in
+             os.listdir(os.path.join(TEST_DATA_DIR, 'acr', 'GE'))])
+
+        self.dcm = self.acr_spatial_resolution_task.ACR_obj.dcm[0]
         self.crop_image = self.acr_spatial_resolution_task.crop_image(self.dcm.pixel_array, self.centre[0],
                                                                       self.y_ramp_pos, self.width)
 
-    def test_object_centre(self):
-        data = self.dcm.pixel_array
-        assert self.acr_spatial_resolution_task.centroid_com(data)[1] == self.centre
-
-    def test_rotation(self):
-        data = self.dcm.pixel_array
-        assert self.acr_spatial_resolution_task.find_rotation(data) == self.rotation_angle
-
     def test_find_y_ramp(self):
         data = self.dcm.pixel_array
         res = self.dcm.PixelSpacing
diff --git a/tests/test_acr_uniformity.py b/tests/test_acr_uniformity.py
index 2de4f1b..08d515e 100644
--- a/tests/test_acr_uniformity.py
+++ b/tests/test_acr_uniformity.py
@@ -2,34 +2,25 @@ import os
 import unittest
 import pathlib
 import pydicom
-import numpy as np
 
 from hazenlib.tasks.acr_uniformity import ACRUniformity
+from hazenlib.ACRObject import ACRObject
 from tests import TEST_DATA_DIR, TEST_REPORT_DIR
 
 
 class TestACRUniformitySiemens(unittest.TestCase):
     ACR_UNIFORMITY_DATA = pathlib.Path(TEST_DATA_DIR / 'acr')
-    centre = [129, 128]
-    piu = 68.66
-    array = np.zeros((256, 256), dtype=int)
-    array[127][126] = 1
-    array[126][127] = 1
-    array[127][127] = 1
-    array[128][127] = 1
-    array[127][128] = 1
+    piu = 67.95
 
     def setUp(self):
         self.acr_uniformity_task = ACRUniformity(data_paths=[os.path.join(TEST_DATA_DIR, 'acr')],
                                                  report_dir=pathlib.PurePath.joinpath(TEST_REPORT_DIR))
-        self.dcm = pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens', '6.dcm'))
 
-    def test_object_centre(self):
-        assert self.acr_uniformity_task.centroid_com(self.dcm.pixel_array) == self.centre
+        self.acr_uniformity_task.ACR_obj = ACRObject(
+            [pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens', f'{i}')) for i in
+             os.listdir(os.path.join(TEST_DATA_DIR, 'acr', 'Siemens'))])
 
-    def test_circular_mask(self):
-        test_circle = self.acr_uniformity_task.circular_mask([128, 128], 1, [256, 256]).astype(int)
-        assert (self.array == test_circle).all() == True
+        self.dcm = self.acr_uniformity_task.ACR_obj.dcm[6]
 
     def test_uniformity(self):
         results = self.acr_uniformity_task.get_integral_uniformity(self.dcm)
@@ -42,20 +33,20 @@ class TestACRUniformitySiemens(unittest.TestCase):
         assert rounded_results == self.piu
 
 
-# class TestACRUniformityPhilips(unittest.TestCase):
-#
+# TODO: Add unit tests for Philips datasets.
+
 class TestACRUniformityGE(unittest.TestCase):
     ACR_UNIFORMITY_DATA = pathlib.Path(TEST_DATA_DIR / 'acr')
-    centre = [253, 256]
-    piu = 84.76
+    piu = 85.17
 
     def setUp(self):
         self.acr_uniformity_task = ACRUniformity(data_paths=[os.path.join(TEST_DATA_DIR, 'acr')],
                                                  report_dir=pathlib.PurePath.joinpath(TEST_REPORT_DIR))
-        self.dcm = pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'GE', '4.dcm'))
+        self.acr_uniformity_task.ACR_obj = ACRObject(
+            [pydicom.read_file(os.path.join(TEST_DATA_DIR, 'acr', 'GE', f'{i}')) for i in
+             os.listdir(os.path.join(TEST_DATA_DIR, 'acr', 'GE'))])
 
-    def test_object_centre(self):
-        assert self.acr_uniformity_task.centroid_com(self.dcm.pixel_array) == self.centre
+        self.dcm = self.acr_uniformity_task.ACR_obj.dcm[6]
 
     def test_uniformity(self):
         results = self.acr_uniformity_task.get_integral_uniformity(self.dcm)