jrkerns / pylinac

An image analysis library for medical physics
https://pylinac.readthedocs.io/en/latest/
MIT License
157 stars 99 forks source link

ACR geometric distortions incorrect in v3.27 #521

Open mitch-1211 opened 1 month ago

mitch-1211 commented 1 month ago

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 commented 1 month ago

MR_1.3.12.2.1107.5.2.50.177688.zip

mitch-1211 commented 1 month ago

image

jrkerns commented 3 weeks ago

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