Skip to content

Commit

Permalink
Merge branch 'main' into 378-geom-accuracy-wrong-slices
Browse files Browse the repository at this point in the history
  • Loading branch information
@mollybuckley
mollybuckley authored Jan 18, 2024
2 parents cb6edbb + c8f540c commit 4d12de9
Show file tree
Hide file tree
Showing 40 changed files with 8,460 additions and 2,714 deletions.
150 changes: 95 additions & 55 deletions hazenlib/ACRObject.py
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
import numpy as np
import cv2
import scipy
import skimage
import numpy as np


class ACRObject:
Expand All @@ -10,8 +10,6 @@ def __init__(self, dcm_list):
self.dcm_list = dcm_list
# Load files as DICOM and their pixel arrays into 'images'
self.images, self.dcms = self.sort_images()
# Store the DCM object of slice 7 as it is used often
self.slice7_dcm = self.dcms[6]
# Store the pixel spacing value from the first image (expected to be the same for all)
self.pixel_spacing = self.dcms[0].PixelSpacing
# Check whether images of the phantom are the correct orientation
Expand All @@ -20,6 +18,8 @@ def __init__(self, dcm_list):
self.rot_angle = self.determine_rotation()
# Find the centre coordinates of the phantom (circle) on slice 7 only:
self.centre, self.radius = self.find_phantom_center(self.images[6])
# Store the DCM object of slice 7 as it is used often
# self.slice7_dcm = self.dcms[6]
# Store a mask image of slice 7 for reusability
self.mask_image = self.get_mask_image(self.images[6])

Expand All @@ -35,6 +35,11 @@ def sort_images(self):
A sorted stack of dicoms
"""

# TODO: implement a check if phantom was placed in other than axial position
# This is to be able to flag to the user the caveat of measurments if deviating from ACR guidance

# x = np.array([dcm.ImagePositionPatient[0] for dcm in self.dcm_list])
# y = np.array([dcm.ImagePositionPatient[1] for dcm in self.dcm_list])
z = np.array([dcm.ImagePositionPatient[2] for dcm in self.dcm_list])
dicom_stack = [self.dcm_list[i] for i in np.argsort(z)]
img_stack = [dicom.pixel_array for dicom in dicom_stack]
Expand All @@ -54,39 +59,52 @@ def orientation_checks(self):
test_images = (self.images[0], self.images[-1])
dx = self.pixel_spacing[0]

normalised_images = [cv2.normalize(
src=image, dst=None, alpha=0, beta=255,
norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_8U
) for image in test_images]
normalised_images = [
cv2.normalize(
src=image,
dst=None,
alpha=0,
beta=255,
norm_type=cv2.NORM_MINMAX,
dtype=cv2.CV_8U,
)
for image in test_images
]

# search for circle in first slice of ACR phantom dataset with radius of ~11mm
detected_circles = [cv2.HoughCircles(
norm_image, cv2.HOUGH_GRADIENT, 1,
param1=50, param2=30,
minDist=int(180 / dx),
minRadius=int(5 / dx),
maxRadius=int(16 / dx)
) for norm_image in normalised_images]
detected_circles = [
cv2.HoughCircles(
norm_image,
cv2.HOUGH_GRADIENT,
1,
param1=50,
param2=30,
minDist=int(180 / dx),
minRadius=int(5 / dx),
maxRadius=int(16 / dx),
)
for norm_image in normalised_images
]

if detected_circles[0] is not None:
true_circle = detected_circles[0].flatten()
else:
true_circle = detected_circles[1].flatten()

if detected_circles[0] is None and detected_circles[1] is not None:
print('Performing slice order inversion to restore correct slice order.')
print("Performing slice order inversion to restore correct slice order.")
self.images.reverse()
self.dcms.reverse()
else:
print('Slice order inversion not required.')
print("Slice order inversion not required.")

if true_circle[0] > self.images[0].shape[0] // 2:
print('Performing LR orientation swap to restore correct view.')
print("Performing LR orientation swap to restore correct view.")
flipped_images = [np.fliplr(image) for image in self.images]
for index, dcm in enumerate(self.dcms):
dcm.PixelData = flipped_images[index].tobytes()
else:
print('LR orientation swap not required.')
print("LR orientation swap not required.")

def determine_rotation(self):
"""
Expand Down Expand Up @@ -122,54 +140,72 @@ def rotate_images(self):
The rotated images.
"""

return skimage.transform.rotate(self.images, self.rot_angle, resize=False, preserve_range=True)
return skimage.transform.rotate(
self.images, self.rot_angle, resize=False, preserve_range=True
)

def find_phantom_center(self, img):
"""
Find the center of the ACR phantom by filtering the uniformity slice and using the Hough circle detector.
Find the center of the ACR phantom by filtering the input slice and using the Hough circle detector.
Args:
img (np.array): pixel array of the dicom
Returns
-------
centre : tuple
Tuple of ints representing the (x, y) center of the image.
Returns:
tuple: Tuple of ints representing the (x, y) center of the image.
"""
dx, dy = self.pixel_spacing

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,
param1=50, param2=30,
minDist=int(180 / dy),
minRadius=int(180 / (2 * dy)),
maxRadius=int(200 / (2 * dx))
).flatten()

img_grad,
cv2.HOUGH_GRADIENT,
1,
param1=50,
param2=30,
minDist=int(180 / dy),
minRadius=int(180 / (2 * dy)),
maxRadius=int(200 / (2 * dx)),
).flatten()
centre = [int(i) for i in detected_circles[:2]]
radius = int(detected_circles[2])
return centre, radius

def get_mask_image(self, image, mag_threshold=0.07, open_threshold=500):
"""
"""Create a masked pixel array
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.
Args:
image (np.array): pixel array of the dicom
mag_threshold (float, optional): magnitude threshold. Defaults to 0.07.
open_threshold (int, optional): open threshold. Defaults to 500.
Returns:
np.array:
The masked image.
"""
test_mask = self.circular_mask(self.centre, (80 // self.pixel_spacing[0]), image.shape)
test_mask = self.circular_mask(
self.centre, (80 // self.pixel_spacing[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)
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)
opened_mask = skimage.morphology.area_opening(
initial_mask, area_threshold=open_threshold
)
final_mask = skimage.morphology.convex_hull_image(opened_mask)

return final_mask
Expand Down Expand Up @@ -198,7 +234,7 @@ def circular_mask(centre, radius, dims):
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
mask = (X - centre[0]) ** 2 + (Y - centre[1]) ** 2 <= radius**2

return mask

Expand Down Expand Up @@ -233,26 +269,28 @@ def measure_orthogonal_lengths(self, mask, slice_index):
horizontal_start = (horizontal, 0)
horizontal_end = (horizontal, dims[0] - 1)
horizontal_line_profile = skimage.measure.profile_line(
mask, horizontal_start, horizontal_end)
mask, horizontal_start, horizontal_end
)
horizontal_extent = np.nonzero(horizontal_line_profile)[0]
horizontal_distance = (horizontal_extent[-1] - horizontal_extent[0]) * dx

vertical_start = (0, vertical)
vertical_end = (dims[1] - 1, vertical)
vertical_line_profile = skimage.measure.profile_line(
mask, vertical_start, vertical_end)
mask, vertical_start, vertical_end
)
vertical_extent = np.nonzero(vertical_line_profile)[0]
vertical_distance = (vertical_extent[-1] - vertical_extent[0]) * dy

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
"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
Expand Down Expand Up @@ -307,10 +345,12 @@ def find_n_highest_peaks(data, n, height=1):
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_heights = peaks[1]["peak_heights"]
pk_ind = peaks[0]

peak_heights = pk_heights[(-pk_heights).argsort()[:n]] # find n highest peak amplitudes
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)
32 changes: 19 additions & 13 deletions hazenlib/HazenTask.py
View file
Original file line number Diff line number Diff line change
@@ -1,24 +1,26 @@
"""
HazenTask.py
"""

import os
import pathlib
from pydicom import dcmread

from hazenlib.logger import logger
import pathlib
import os


class HazenTask:

def __init__(self, input_data: list, report: bool = False, report_dir=None, **kwargs):
def __init__(
self, input_data: list, report: bool = False, report_dir=None, **kwargs
):
data_paths = sorted(input_data)
self.dcm_list = [dcmread(dicom)for dicom in data_paths]
self.dcm_list = [dcmread(dicom) for dicom in data_paths]
self.report: bool = report
if report_dir is not None:
self.report_path = os.path.join(str(report_dir), type(self).__name__)
else:
self.report_path = os.path.join(os.getcwd(), 'report_image',
type(self).__name__)
self.report_path = os.path.join(
os.getcwd(), "report_image", type(self).__name__
)
if report:
pathlib.Path(self.report_path).mkdir(parents=True, exist_ok=True)
else:
Expand All @@ -29,18 +31,22 @@ def init_result_dict(self) -> dict:
result_dict = {
"task": f"{type(self).__name__}",
"file": None,
"measurement": {}
"measurement": {},
}
return result_dict

def img_desc(self, dcm, properties=None) -> str:
if properties is None:
properties = ['SeriesDescription', 'SeriesNumber', 'InstanceNumber']
properties = ["SeriesDescription", "SeriesNumber", "InstanceNumber"]
try:
metadata = [str(dcm.get(field)) for field in properties]
except KeyError:
logger.warning(f"Could not find one or more of the following properties: {properties}")
metadata = [str(dcm.get(field)) for field in ['SeriesDescription', 'SeriesNumber']]
logger.warning(
f"Could not find one or more of the following properties: {properties}"
)
metadata = [
str(dcm.get(field)) for field in ["SeriesDescription", "SeriesNumber"]
]

img_desc = '_'.join(metadata).replace(' ', '_')
img_desc = "_".join(metadata).replace(" ", "_")
return img_desc
Loading

0 comments on commit 4d12de9

Please sign in to comment.