Median vs. mean roi values

[Devin Baillie]

Why are ROI values implemented as medians instead of means?  Is there a reason to prefer one over the other for this application?

[James Kerns]

The calculation is over a numpy array with NaNs in it. At the time I wrote the module, nanmean wasn't available, just nanmedian. Given that ROI arrays are usually normally distributed, it pretty close.

[Devin Baillie]

Makes sense, thanks for the response.  The issue I'm having is that some of my CT data comes in as integers, and taking the median forces the result to be an integer.  Not a major issue, but is there any chance of getting this changed in a future version?

As well, I'm not sure the assumption of normality is accurate, particularly for air on some CT scanners due to the range being bounded at -1000.

Thanks,
Devin

[James Kerns]

For backwards stability it will likely stay the same as the default. I will put in a ticket for a new setting.

[Devin Baillie]

Thanks.

[James Han]

Sorry I'm new to being a developer and little to none experience with Github - I'll try to follow the developer's guide when I get some time.

But I rarely do changes to source code, so not sure I will get to it.

I implemented a "mean" version of my ROIs instead of median for CT phantom analysis.  

Honestly, I like median but my group had some dissenting voices and the quoted TG mentions "mean" and doesn't mention median  :(
So I thought I'd share what I did.

Code changes:

FYI: I use virtual environment and analyzing gammex_phantom clinically.

Gammex code (Custom code for Gammex.py) was described in the forums (Thanks for that btw!!!)
and it inherits CheeseModule & CheesePhantomBase (cheese.py)
which inherits from CatPhanBase, CatPhanModule (ct.py)
which inherits from DiskROI (roi.py)

Lib\site-packages\pylinac\core\ROI.py (Maybe this is all that is needed for the actual mean to surface, but I made other changes below):

    def pixel_value(self) -> float:
        """The mean pixel value of the ROI.""" 
        masked_img = self.circle_mask()
        return float(np.mean(masked_img))

Lib\site-packages\pylinac\ct.py (I think this is for CatPhan504 cause I started using this for testing but then later I noticed my CheesePhantom/Gammex didn't see the change to mean so I went to ROI.py):

def rois_to_results(dict_mapping: dict[str, DiskROI]) -> dict[str, ROIResult]:
    """Converts a dict of HUDiskROIs to a dict of ROIResults. This is for dumping to simple data formats for results_data and RadMachine"""
    flat_dict = {}
    for name, roi in dict_mapping.items():
        flat_dict[name] = ROIResult(
            name=name,
            value=np.mean(roi.pixel_values), # Correctly calculates the mean from the pixel array
            stdev=roi.std,
            difference=roi.value_diff if hasattr(roi, "value_diff") else None,
            nominal_value=roi.nominal_val if hasattr(roi, "nominal_val") else None,
            passed=roi.passed if hasattr(roi, "passed") else None,
        )
    return flat_dict

and:

class ROIResult(BaseModel):
    name: str = Field(description="The region the ROI was sampled from.")
    value: float = Field(description="The measured HU value. NOTE: MODIFIED to MEAN not MEDIAN")
    stdev: float = Field(description="The pixel value standard deviation of the ROI.")
    difference: float | None = Field(description="The difference between the measured and nominal values.")
    nominal_value: float | None = Field(description="The nominal HU value.")
    passed: bool | None = Field(description=" Whether the ROI passed.")

Lib\site-packages\pylinac\cheese.py
    def results(self, as_list: bool = False) -> str | list[str]:
        """Return the results of the analysis as a string. Use with print().

        Parameters
        ----------
        as_list : bool
            Whether to return as a list of strings vs single string. Pretty much for internal usage.
        """
        results = [
            f" - {self.model} Phantom Analysis - ",
            " - HU Module - ",
        ]
        results += [
            f"ROI {name} Mean: {roi.pixel_value:.2f}, stdev: {roi.std:.1f}"
            for name, roi in self.module.rois.items()
        ]
        if as_list:
            return results
        else:
            return "\n".join(results)

and

    def _quaac_datapoints(self) -> dict[str, QuaacDatum]:
        results_data = self.results_data(as_dict=True)
        data = {}
        data["Phantom roll"] = QuaacDatum(
            value=results_data["phantom_roll"],
            unit="degrees",
        )
        for roi_num, roi_data in results_data["rois"].items():
            data[f"ROI {roi_num}"] = QuaacDatum(
                value=roi_data["mean"], 
                unit="HU",
            )
        return data

BTW - the code's roll uncertainty is pretty impressive.  I intentionally did an angular sensitivity via software to see how the code would find ROIs and was impressed when it already has some robust rotational analysis within 10 degrees I think.

[James Han]

Some results of code change Mean vs Median HU on gammex phantom:

Rod Type Substitute Median HU Mean HU
CT Solid Water 9 9.27
LN-300 Lung -682 -683.25
CT Solid Water 9 8.59
IB Inner Bone 212 212.26
CT solid Water 10.5 10.71
LV1 Liver 82 81.41
B200 Bone Material 229 228.24
LN-450 Lung -529 -530.57
BR-12 Breast -37 -36.48
SB3 Cortical Bone 1176 1176.22
BRN-SR2 Brain 33 32.81
AP6 adipose -79 -78.42
CB2-50% CaCO3 780 780.11
CT Solid Water 11 11.56
CB2-30% CaCO3 436 436.48
Water Insert 3 0.74

Difference (respectively): mean - median
0.27
1.25
0.41
0.26
0.21
0.59
0.76
1.57
0.52
0.22
0.19
0.58
0.11
0.56
0.48
2.26