pylinac icon indicating copy to clipboard operation
pylinac copied to clipboard

ACR geometric distortions incorrect in v3.27

Open mitch-1211 opened this issue 1 year ago • 3 comments

Describe the bug v3.27 fixed the following issue: "Phantoms that were not fully filled and/or had air gaps at the top would sometimes cause the geometric distortion analysis to fail. This was caused by the air gap not completing the full circle of the phantom. The algorithm no longer relies on this assumption and is robust to these air gaps for the geometric distortion analysis."

  1. This issue can still occur for particular datasets (attached example dataset)
  2. The geometric distortion in the diagonal directions is reported incorrectly for all datasets (nominal value is 190mm) Geometric Distortions: {'horizontal': '191.41mm', 'vertical': '24.41mm', 'negative diagonal': '134.81mm', 'positive diagonal': '135.74mm'} To Reproduce Steps to reproduce the behavior:
  1. Use the attached DICOM dataset
  2. Perform an ACRMRILarge analysis with default settings

Expected behavior

  1. Geometric distortions are returned for all directions
  2. Geometric distortion values are all close to 190 mm

Notes For problem 2) If we consider the function _setup_rois in the class GeometricDistortionModule:

  def _setup_rois(self) -> None:
       """This is mostly for plotting purposes. This is why we use FWXMProfile
       instead of FWXMProfilePhysical. The lines to plot should be in pixel coordinates, not physical.
       We convert to physical just for the field width calculation."""
       px_to_cut_off = int(round(5 / self.mm_per_pixel))
       self.profiles = {}
       bin_image = self.image.as_binary(threshold=np.percentile(self.image, 60))
       bin_image = ndimage.binary_fill_holes(bin_image).astype(float)
       # calculate horizontal
       data = bin_image[int(self.phan_center.y), :]
       # cutoff 3mm from the search area
       f_data = fill_middle_zeros(data, cutoff_px=px_to_cut_off)
       prof = FWXMProfile(values=f_data)
       line = Line(
           Point(prof.field_edge_idx(side="left"), self.phan_center.y),
           Point(prof.field_edge_idx(side="right"), self.phan_center.y),
       )

       self.profiles["horizontal"] = {
           "width (mm)": prof.field_width_px * self.mm_per_pixel,
           "line": line,
       }
       # calculate vertical
       data = bin_image[:, int(self.phan_center.x)]
       f_data = fill_middle_zeros(data, cutoff_px=px_to_cut_off)
       prof = FWXMProfile(values=f_data)
       line = Line(
           Point(self.phan_center.x, prof.field_edge_idx(side="left")),
           Point(self.phan_center.x, prof.field_edge_idx(side="right")),
       )
       self.profiles["vertical"] = {
           "width (mm)": prof.field_width_px * self.mm_per_pixel,
           "line": line,
       }
       # calculate negative diagonal
       # calculate slope equation intercept
       # b = y - (+1)x
       b = self.phan_center.y - self.phan_center.x
       xs = np.arange(0, self.image.shape[1])
       ys = xs + b
       coords = ndimage.map_coordinates(bin_image, [ys, xs], order=1, mode="mirror")
       f_data = fill_middle_zeros(coords, cutoff_px=px_to_cut_off)
       # pixels are now diagonal and thus spacing between pixels is now the hypotenuse
       prof = FWXMProfile(values=f_data)
       line = Line(
           Point(
               xs[int(round(prof.field_edge_idx(side="left")))],
               ys[int(round(prof.field_edge_idx(side="left")))],
           ),
           Point(
               xs[int(round(prof.field_edge_idx(side="right")))],
               ys[int(round(prof.field_edge_idx(side="right")))],
           ),
       )
       self.profiles["negative diagonal"] = {
           "width (mm)": prof.field_width_px * self.mm_per_pixel,
           "line": line,
       }
       # calculate positive diagonal
       # calculate slope equation intercept
       # b = y - (-1)x
       b = self.phan_center.y + self.phan_center.x
       ys = -xs + b
       coords = ndimage.map_coordinates(bin_image, [ys, xs], order=1, mode="mirror")
       f_data = fill_middle_zeros(coords, cutoff_px=px_to_cut_off)
       prof = FWXMProfile(values=f_data)
       line = Line(
           Point(
               xs[int(round(prof.field_edge_idx(side="left")))],
               ys[int(round(prof.field_edge_idx(side="left")))],
           ),
           Point(
               xs[int(round(prof.field_edge_idx(side="right")))],
               ys[int(round(prof.field_edge_idx(side="right")))],
           ),
       )
       self.profiles["positive diagonal"] = {
           "width (mm)": prof.field_width_px * self.mm_per_pixel,
           "line": line,
       }

We see the comment

pixels are now diagonal and thus spacing between pixels is now the hypotenuse

It looks as though the hypotenuse is never actually calculated i.e. the pixel width is used, not the diagonal pixel distance. In this example we are returned Geometric Distortions: {'horizontal': '191.41mm', 'vertical': '24.41mm', 'negative diagonal': '134.81mm', 'positive diagonal': '135.74mm'} which if we calculate sqrt(135.74^2 + 135.74^2) we get 191.96 mm which much closer to our expected value of 190 mm.

Potentially need to update self.mm_per_pixel for the diagonals to be the hypotenuse distance

mitch-1211 avatar Oct 15 '24 06:10 mitch-1211

image

mitch-1211 avatar Oct 15 '24 06:10 mitch-1211

Thanks for the write-up. I've been out of town but it's added to my list!

jrkerns avatar Oct 28 '24 14:10 jrkerns