Skip to content

Commit

Permalink
Merge pull request #185 from GSTT-CSC/fix_ghosting_script
Browse files Browse the repository at this point in the history
Updated ghosting function – fixes sampling point issue
  • Loading branch information
tomaroberts authored Dec 22, 2021
2 parents 594723e + 9f7510f commit 2c25632
Show file tree
Hide file tree
Showing 2 changed files with 97 additions and 64 deletions.
87 changes: 57 additions & 30 deletions hazenlib/ghosting.py
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ def calculate_ghost_intensity(ghost, phantom, noise) -> float:
if phantom_mean < ghost_mean or phantom_mean < noise_mean:
raise Exception(f"The mean phantom signal is lower than the ghost or the noise signal. This can't be the case ")

return 100 * (ghost_mean - noise_mean) / phantom_mean
return 100 * abs((ghost_mean - noise_mean)) / phantom_mean


def get_signal_bounding_box(array: np.ndarray):
Expand All @@ -45,14 +45,14 @@ def get_signal_bounding_box(array: np.ndarray):
if voxel > signal_limit:
signal.append(idx)


signal_column = sorted([voxel[1] for voxel in signal])
signal_row = sorted([voxel[0] for voxel in signal])

upper_row = min(signal_row)
lower_row = max(signal_row)
left_column = min(signal_column)
right_column = max(signal_column)

return left_column, right_column, upper_row, lower_row,


Expand Down Expand Up @@ -111,6 +111,7 @@ def get_background_slices(background_rois, slice_radius=5):
slices = [(np.array(range(roi[0]-slice_radius, roi[0]+slice_radius), dtype=np.intp)[:, np.newaxis], np.array(
range(roi[1]-slice_radius, roi[1]+slice_radius), dtype=np.intp))for roi in background_rois]


return slices


Expand All @@ -131,75 +132,95 @@ def get_eligible_area(signal_bounding_box, dcm, slice_radius=5):
# signal is in left half
eligible_columns = range(right_column + padding_from_box, dcm.Columns - slice_radius)
eligible_rows = range(upper_row, lower_row)
ghost_slice = np.array(
range(right_column + padding_from_box, dcm.Columns - slice_radius), dtype=np.intp)[:, np.newaxis], np.array(
range(upper_row, lower_row)
)
else:
# signal is in right half
eligible_columns = range(slice_radius, left_column - padding_from_box)
eligible_rows = range(upper_row, lower_row)
ghost_slice = np.array(
range(slice_radius, left_column - padding_from_box), dtype=np.intp)[:,
np.newaxis], np.array(
range(upper_row, lower_row))

else:
if upper_row < dcm.Rows / 2:
# signal is in top half
eligible_rows = range(lower_row + padding_from_box, dcm.Rows - slice_radius)
eligible_columns = range(left_column, right_column)
ghost_slice = np.array(
range(lower_row + padding_from_box, dcm.Rows - slice_radius), dtype = np.intp)[:,
np.newaxis], np.array(
range(left_column, right_column))
else:
# signal is in bottom half
eligible_rows = range(slice_radius, upper_row - padding_from_box)
eligible_columns = range(left_column, right_column)
ghost_slice = np.array(
range(slice_radius, upper_row - padding_from_box), dtype=np.intp)[:,
np.newaxis], np.array(
range(left_column, right_column))

return eligible_columns, eligible_rows


def get_ghost_slice(signal_bounding_box, dcm, slice_radius=5):
max_mean = 0
max_index = (0, 0)
windows = {}
arr = dcm.pixel_array

eligible_columns, eligible_rows = get_eligible_area(signal_bounding_box, dcm, slice_radius)
for idx, centre_voxel in np.ndenumerate(arr):
if idx[0] not in eligible_rows or idx[1] not in eligible_columns:
continue
else:
windows[idx] = arr[idx[0]-slice_radius:idx[0]+slice_radius, idx[1]-slice_radius:idx[1]+slice_radius]

for idx, window in windows.items():
if np.mean(window) > max_mean:
max_mean = np.mean(window)
max_index = idx


def get_ghost_slice(signal_bounding_box, dcm, slice_radius=5):
eligible_area = get_eligible_area(signal_bounding_box, dcm, slice_radius=slice_radius)
ghost_slice = np.array(
range(max_index[1] - slice_radius, max_index[1] + slice_radius), dtype=np.intp)[:, np.newaxis], np.array(
range(max_index[0] - slice_radius, max_index[0] + slice_radius)
range(min(eligible_area[1]),max(eligible_area[1])), dtype=np.intp)[:, np.newaxis], np.array(
range(min(eligible_area[0]), max(eligible_area[0]))
)
#ghost_slice returns columns, rows of indexes
return ghost_slice


def get_ghosting(dcm, plotting=False) -> dict:






def get_ghosting(dcm, report_path=False) -> dict:

bbox = get_signal_bounding_box(dcm.pixel_array)

x, y = hazenlib.get_pixel_size(dcm) # assume square pixels i.e. x=y
# ROIs need to be 10mmx10mm
slice_radius = int(10//(2*x))


signal_centre = [(bbox[0]+bbox[1])//2, (bbox[2]+bbox[3])//2]
background_rois = get_background_rois(dcm, signal_centre)
ghost_col, ghost_row = get_ghost_slice(bbox, dcm, slice_radius=slice_radius)
ghost = dcm.pixel_array[(ghost_row, ghost_col)]
ghost = dcm.pixel_array[(ghost_col, ghost_row)]
signal_col, signal_row = get_signal_slice(bbox, slice_radius=slice_radius)
phantom = dcm.pixel_array[(signal_row, signal_col)]



noise = np.concatenate([dcm.pixel_array[(row, col)] for col, row in get_background_slices(background_rois, slice_radius=slice_radius)])




eligible_area = get_eligible_area(bbox, dcm, slice_radius=slice_radius)

ghosting = calculate_ghost_intensity(ghost, phantom, noise)

if plotting:
if report_path == True:
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
x1, x2, y1, y2 = bbox

img = hazenlib.rescale_to_byte(dcm.pixel_array)
img = dcm.pixel_array
img = img.astype('float64')
#print('this is img',img)
img *= 255.0 / img.max()
#img = hazenlib.rescale_to_byte(dcm.pixel_array)
img = cv.rectangle(img.copy(), (x1, y1), (x2, y2), (255, 0, 0), 1)

for roi in background_rois:
Expand All @@ -223,25 +244,28 @@ def get_ghosting(dcm, plotting=False) -> dict:
img = cv.rectangle(img.copy(), (x1, y1), (x2, y2), (255, 0, 0), 1)

ax.imshow(img)
fig.savefig(f'{report_path}.png')
return fig, ghosting

return None, ghosting


def main(data: list, report_path=False) -> dict:
def main(data: list, report_path = False) -> dict:

results = {}
# figures = []
#figures = []
for dcm in data:
try:
key = f"{dcm.SeriesDescription.replace(' ', '_')}_{dcm.EchoTime}ms_NSA-{dcm.NumberOfAverages}"
if report_path:
report_path = key
except AttributeError as e:
print(e)
key = f"{dcm.SeriesDescription}_{dcm.SeriesNumber}"
try:
fig, results[key] = get_ghosting(dcm, report_path)
if report_path:
fig.savefig(key + '.png')
fig, results[key] = get_ghosting(dcm, report_path=True)



except Exception as e:
print(f"Could not calculate the ghosting for {key} because of : {e}")
Expand All @@ -251,5 +275,8 @@ def main(data: list, report_path=False) -> dict:
return results



if __name__ == "__main__":
main([os.path.join(sys.argv[1], i) for i in os.listdir(sys.argv[1])])


74 changes: 40 additions & 34 deletions tests/test_ghosting.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,17 +15,17 @@ class TestGhosting(unittest.TestCase):
PADDING_FROM_BOX = 30
SLICE_RADIUS = 5
PE = 'ROW'
ELIGIBLE_GHOST_AREA = range(5, SIGNAL_BOUNDING_BOX[0] - PADDING_FROM_BOX), range(
SIGNAL_BOUNDING_BOX[2], SIGNAL_BOUNDING_BOX[3])
ELIGIBLE_GHOST_AREA = range(5, SIGNAL_BOUNDING_BOX[0] - PADDING_FROM_BOX), range(SIGNAL_BOUNDING_BOX[2], SIGNAL_BOUNDING_BOX[3])

SIGNAL_SLICE = np.array(range(SIGNAL_CENTRE[0] - SLICE_RADIUS,
SIGNAL_CENTRE[0] + SLICE_RADIUS), dtype=np.intp)[:, np.newaxis], np.array(
range(SIGNAL_CENTRE[1] - SLICE_RADIUS, SIGNAL_CENTRE[1] + SLICE_RADIUS), dtype=np.intp)
SIGNAL_SLICE = np.array(range(SIGNAL_CENTRE[0] - SLICE_RADIUS, SIGNAL_CENTRE[0] + SLICE_RADIUS), dtype=np.intp)[:, np.newaxis], np.array(
range(SIGNAL_CENTRE[1] - SLICE_RADIUS, SIGNAL_CENTRE[1] + SLICE_RADIUS), dtype=np.intp)

GHOST_SLICE = np.array(
range(min(ELIGIBLE_GHOST_AREA[1]), max(ELIGIBLE_GHOST_AREA[1])), dtype=np.intp)[:, np.newaxis], np.array(
range(min(ELIGIBLE_GHOST_AREA[0]), max(ELIGIBLE_GHOST_AREA[0])))

GHOSTING = (None, 0.11803264099090763)

GHOST_SLICE = np.array(range(13 - SLICE_RADIUS, 13 + SLICE_RADIUS), dtype=np.intp)[:, np.newaxis], np.array(
range(283 - SLICE_RADIUS, 283 + SLICE_RADIUS)
)
GHOSTING = (None, 0.6115742783493184)

def setUp(self):
self.file = str(TEST_DATA_DIR / 'ghosting' / 'GHOSTING' / 'IM_0001.dcm')
Expand All @@ -47,14 +47,10 @@ def test_calculate_ghost_intensity(self):
phantom=np.asarray([100]),
noise=np.asarray([5]))

assert -5.0 == hazen_ghosting.calculate_ghost_intensity(ghost=np.asarray([5]),
phantom=np.asarray([100]),
noise=np.asarray([10]))

def test_get_signal_bounding_box(self):
(left_column, right_column, upper_row, lower_row,) = hazen_ghosting.get_signal_bounding_box(
self.dcm.pixel_array)

def test_get_signal_bounding_box(self):
(left_column, right_column, upper_row, lower_row,) = hazen_ghosting.get_signal_bounding_box(self.dcm.pixel_array)
assert (left_column, right_column, upper_row, lower_row) == self.SIGNAL_BOUNDING_BOX

def test_get_signal_slice(self):
Expand Down Expand Up @@ -85,38 +81,48 @@ class TestCOLPEGhosting(TestGhosting):
PADDING_FROM_BOX = 30
SLICE_RADIUS = 5
ELIGIBLE_GHOST_AREA = range(SIGNAL_BOUNDING_BOX[0], SIGNAL_BOUNDING_BOX[1]), range(
SLICE_RADIUS, SIGNAL_BOUNDING_BOX[2] - PADDING_FROM_BOX)
SLICE_RADIUS, SIGNAL_BOUNDING_BOX[2] - PADDING_FROM_BOX)

SIGNAL_SLICE = np.array(range(SIGNAL_CENTRE[0] - SLICE_RADIUS, SIGNAL_CENTRE[0] + SLICE_RADIUS), dtype=np.intp)[:,np.newaxis],\
np.array(range(SIGNAL_CENTRE[1] - SLICE_RADIUS, SIGNAL_CENTRE[1] + SLICE_RADIUS),dtype=np.intp)

GHOST_SLICE = np.array(
range(min(ELIGIBLE_GHOST_AREA[1]), max(ELIGIBLE_GHOST_AREA[1])), dtype=np.intp)[:, np.newaxis], np.array(
range(min(ELIGIBLE_GHOST_AREA[0]), max(ELIGIBLE_GHOST_AREA[0])))

SIGNAL_SLICE = np.array(range(SIGNAL_CENTRE[0] - SLICE_RADIUS, SIGNAL_CENTRE[0] + SLICE_RADIUS), dtype=np.intp)[:,
np.newaxis], np.array(range(SIGNAL_CENTRE[1] - SLICE_RADIUS, SIGNAL_CENTRE[1] + SLICE_RADIUS),
dtype=np.intp)
GHOST_SLICE = np.array(range(197 - SLICE_RADIUS, 197 + SLICE_RADIUS), dtype=np.intp)[:, np.newaxis], np.array(
range(135 - SLICE_RADIUS, 135 + SLICE_RADIUS))
PE = "COL"
GHOSTING = (None, 0.317544586211758)
GHOSTING = (None, 0.015138960417776908)

def setUp(self):
self.file = str(TEST_DATA_DIR / 'ghosting' / 'PE_COL_PHANTOM_BOTTOM_RIGHT' / 'PE_COL_PHANTOM_BOTTOM_RIGHT.IMA')
self.dcm = pydicom.read_file(self.file)
self.file = str(TEST_DATA_DIR / 'ghosting' / 'PE_COL_PHANTOM_BOTTOM_RIGHT' / 'PE_COL_PHANTOM_BOTTOM_RIGHT.IMA')
self.dcm = pydicom.read_file(self.file)


class TestAxialPhilipsBroomfields(TestGhosting):
SIGNAL_BOUNDING_BOX = (217, 299, 11, 93)
SIGNAL_CENTRE = [(SIGNAL_BOUNDING_BOX[0] + SIGNAL_BOUNDING_BOX[1]) // 2,
(SIGNAL_BOUNDING_BOX[2] + SIGNAL_BOUNDING_BOX[3]) // 2]
(SIGNAL_BOUNDING_BOX[2] + SIGNAL_BOUNDING_BOX[3]) // 2]
BACKGROUND_ROIS = [(258, 264), (194, 264), (130, 264), (66, 264)]
PADDING_FROM_BOX = 30
SLICE_RADIUS = 5
ELIGIBLE_GHOST_AREA = range(SLICE_RADIUS, SIGNAL_BOUNDING_BOX[0] - PADDING_FROM_BOX), range(
SIGNAL_BOUNDING_BOX[2], SIGNAL_BOUNDING_BOX[3])
SIGNAL_SLICE = np.array(range(SIGNAL_CENTRE[0] - SLICE_RADIUS, SIGNAL_CENTRE[0] + SLICE_RADIUS), dtype=np.intp)[:,
np.newaxis], np.array(range(SIGNAL_CENTRE[1] - SLICE_RADIUS, SIGNAL_CENTRE[1] + SLICE_RADIUS),
dtype=np.intp)
GHOST_SLICE = np.array(range(11 - SLICE_RADIUS, 11 + SLICE_RADIUS), dtype=np.intp)[:, np.newaxis], np.array(
range(16 - SLICE_RADIUS, 16 + SLICE_RADIUS))
SIGNAL_BOUNDING_BOX[2], SIGNAL_BOUNDING_BOX[3])
SIGNAL_SLICE = np.array(range(SIGNAL_CENTRE[0] - SLICE_RADIUS, SIGNAL_CENTRE[0] + SLICE_RADIUS), dtype=np.intp)[:,np.newaxis],\
np.array(range(SIGNAL_CENTRE[1] - SLICE_RADIUS, SIGNAL_CENTRE[1] + SLICE_RADIUS),dtype=np.intp)
GHOST_SLICE = np.array(
range(min(ELIGIBLE_GHOST_AREA[1]), max(ELIGIBLE_GHOST_AREA[1])), dtype=np.intp)[:, np.newaxis], np.array(
range(min(ELIGIBLE_GHOST_AREA[0]), max(ELIGIBLE_GHOST_AREA[0])))

GHOSTING = (None, 0.6244366401836245)
GHOSTING = (None, 0.007246960909896829)

def setUp(self):
self.file = str(TEST_DATA_DIR / 'ghosting' / 'GHOSTING' / 'axial_philips_broomfields.dcm')
self.dcm = pydicom.read_file(self.file)
self.file = str(TEST_DATA_DIR / 'ghosting' / 'GHOSTING' / 'axial_philips_broomfields.dcm')
self.dcm = pydicom.read_file(self.file)








0 comments on commit 2c25632

Please sign in to comment.