Catphan 604 Contrast

[crcrewso]

I'm looking into the CTP515 low contrast visibility section. We're trying to get signal to noise directly to compare with previous tools. Does the algorithm take the HU values without modification or is there a normalization of some kind to account for the scale starting at -1000? Is there a reference for how this decision was made?

[jrkerns]

HU values are not modified beyond the typical loading process. It uses the mean of the two surrounding background ROIs when comparing to the ROI. Any chance the HU values are negative like #539 ? (I've not looked at that one closely yet)

[crcrewso]

I'm looking at the algorithms Doselab uses and it seems they're all rescaled up by 1000 to every HU, and then reported as percentages. The consistency of noise and agreement with Doselab seems to be improved by rescaling. The area we used to analyze for this purpose was non-traditional and the signal could swing from -50 to +10.

What are your thoughts?

Has there been any thought in reducing the size of the ROI's for contrast and for HU by some small amount to remove the edges and use median in place of mean?

[jrkerns]

Ah. I can see why that's in place although 1000 seems arbitrary. If there's a paper that mentions doing that I'd consider doing it. I've had to bend over backwards for Doselab before and it's not my favorite. That being said, let me look at #539 first and see if I can do anything since that appears related. Worst case, I'll see if I can get a script workaround together.

"Has there been any thought in reducing the size of the ROI's for contrast and for HU by some small amount to remove the edges and use median in place of mean?" In general, when it comes to changing ROI size I've always been bitten because it's just on the right side of the line for some datasets and just on the wrong side for others (but of course this is customizable). As for median/mean, that choice was arbitrary. I'm fine using either.

[crcrewso]

I'll spend some time looking for papers to see if/how other sources scale. I believe the +1000 factor comes from HU's being defined as -1000 as air and 0 as water, whereas most photography image quality measurements are using black = 0 and everything else is positive. Are there any IP concerns with trying to agree with Image Owl's algorithm (There's a reference in the Catphan manual about it being a possible reference algorithm).

Automated computer analysis
To assist our customers, The Phantom Laboratory has worked with Image Owl, Inc. to develop an automated QA cloud-based service.

Image Owl

Our traditional contrast measurement for tracking drifting image quality was completely arbitrary. We compared the Poly ROI to the background of the phantom at a little further out from the central pixel from where you have the 604 background region 2 defined. We've found that the background at the same radius as the density inserts can often have the most significant artifacts of the module. This was chosen because DoseLab's ability to analyze the low contrast module on CBCT's taken from our iX's were worse than useless. It does so many things oddly that I'm surprised when I find it using a reasonable solution.

On a side-note the existence of values below -1000 HU in a CT have also been a head scratcher but I've never looked too deep into the calibration and reconstruction algorithms. I've included my code as an FYI.

import os
import pylinac
import statistics
import numpy as np
from pylinac import CatPhan604
from pylinac import core
from pylinac.ct import CTP404CP604, CTP486, CTP515, CTP528CP604


class SCACatphanCTP404(CTP404CP604):

# set smaller radius to remove potential edge artifacts
    roi_radius_mm = 4 
      # nominal diameter of inserts are 1 cm
    roi_settings = {
      # Copy/pasted directly from ct.py's CTP404CP604
    }
    background_roi_settings = {
        "1": {
            "angle": -30,
# background moved further away from the center, this should probably be static
            "distance": roi_dist_mm + 2 * roi_radius_mm, 
            "radius": roi_radius_mm,
        },
        "2": {
            "angle": -210,
            "distance": roi_dist_mm + 2 * roi_radius_mm,
            "radius": roi_radius_mm,
        },
    }


class SCACatphan604(CatPhan604):
    modules = {
        SCACatphanCTP404: {"offset": 0},
        CTP486: {"offset": -80},
        CTP528CP604: {"offset": 40},
        CTP515: {"offset": -40},
    }

def CustomAnalysis(Bin_File):

# Setup
    ctSet = SCACatphan604.from_zip(BIN_FILE)
    ctSet.analyze()
    tempDict = ctSet.results_data(as_dict=True)
  
# Get values
    signal_array = ctSet.ctp404.rois["Poly"].pixel_values + 1000
    background_array = ctSet.ctp404.background_rois["2"].pixel_values + 1000
    signal_mean = round(statistics.mean(signal_array))
    signal_std = statistics.stdev(ctSet.ctp404.rois["Poly"].pixel_values)
    background_val = round(statistics.mean(background_array))
    background_std = statistics.stdev(background_array)
  
# Calculations - Manually added for clarity  
    contrast = abs(signal_mean - background_val) / (abs(signal_mean + background_val))
    noise = np.sqrt(signal_std**2 + background_std**2) / np.sqrt(signal_mean**2 + background_val**2)

# Making results available with the rest of the results_data dictionary object
    tempDict["ctp404"]["Contrast"] = {}
    tempDict["ctp404"]["Contrast"]["signal"] = signal_mean
    tempDict["ctp404"]["Contrast"]["background"] = background_val
    tempDict["ctp404"]["Contrast"]["contrast"] = contrast * 100 # Doselab reports percent
    tempDict["ctp404"]["Contrast"]["noise"] = noise * 100 # Doselab reports percent
    tempDict["ctp404"]["Contrast"]["cnr"] = contrast / noise
    tempDict["ctp404"]["Contrast"]["snr"] = abs(signal_mean) / signal_std
    tempDict["status"] = "success"
    
    return tempDict

[jrkerns]

Thanks for the context! If 1000 is the norm then I'm willing to add that as an option.

[crcrewso]

Do you have an idea for what an implementation of this option might look like?