Assessing robustness of radiomic features by image perturbation

In plain words

Radiomics turns a medical scan into a long list of numbers that describe a tumour: its size, its brightness, how grainy its texture looks. A model built on these numbers only works on new scans if the numbers stay stable when the scan is taken slightly differently, for example with the patient lying a little off-centre, with a bit more image noise, or with the tumour outline drawn a little differently. A stable number is called robust.

The usual way to find stable numbers is to scan the same patient twice and see which numbers agree. That doubles the scanning, and it is often not done in routine care. The authors ask whether one scan is enough. They take that scan, make many altered copies of it by applying small, realistic changes (rotate, add noise, shift, grow or shrink the outline, and randomly redraw parts of the outline), and count how many numbers stay stable across the copies. They then check the results against the real two-scan data. The best combination of changes flags very few numbers as stable when they are not.

TL;DR

• Compares 18 chains of image perturbations with test-retest CT scans of two cohorts (31 NSCLC, 19 HNSCC) §2 §4.1.
• Robustness is the intraclass correlation coefficient ICC(1,1); a feature with ICC ≥ 0.90 counts as robust §4.5.
• The chain noise + translation + volume change + contour randomisation (NTVC) gave the fewest false-positive robust features: 3.3% (NSCLC), 10.0% (HNSCC), taking test-retest as the reference §2.3.
• The results are author-reported on small cohorts; the method is CT-only.

Problem

Radiomic models that use features which change with patient positioning, acquisition or segmentation generalise poorly §1. Test-retest imaging is the usual way to identify non-robust features, but it is rarely available for a given phenotype, and robustness does not transfer between phenotypes, modalities or computational settings (voxel size, discretisation) §1. The authors want to know whether perturbations of a single image can replace test-retest data for this purpose §1.

Method

Overview

For each CT image and its gross tumour volume (GTV) mask, the pipeline is:

1. Convert the RTSTRUCT contour to a voxel mask §4.2.
2. Apply a chain of optional perturbations: rotation, noise, translation, volume adaptation, contour randomisation. Chains are built from the five basic perturbations and every combination of their settings is generated §2 [supp. note 5].
3. Resample to isotropic spacings of 1, 2, 3 and 4 mm, after a Gaussian low-pass filter before the down-sampling §4.2 [supp. note 2].
4. Re-segment the intensity mask, exclude voxels more than three standard deviations from the ROI mean, and compute 4032 features per perturbed image §4.2 §4.4 [supp. note 3].
5. For each feature, compute ICC(1,1) of test-retest values across patients, and the perturbation ICC of each chain separately for the test and the retest image, averaged §4.5 fig. 2.
6. Call a feature robust if its ICC is at least 0.90 §4.5.
7. Compare each perturbation chain with the test-retest result feature by feature: robust under both is a true positive, non-robust under both a true negative, robust only under perturbation a false positive, robust only under test-retest a false negative §2.3 fig. 4.

# Robustness pipeline for one patient  [§2] [§4.2] [§4.5]
# Input:  test_ct, retest_ct (3-D, HU); gtv_test, gtv_retest (3-D bool)
# Output: per feature: test-retest ICC, perturbation ICC per chain

def robustness_pipeline(test_ct, retest_ct, gtv_test, gtv_retest, chains):
    # Test-retest reference, on unperturbed images  [§2.1]
    features_test = extract_features(test_ct, gtv_test)
    features_retest = extract_features(retest_ct, gtv_retest)
    icc_by_chain = {}
    # Perturbed copies of each image, one set per chain  [§4.3] [tab. 1]
    for chain in chains:
        perturbed_test = perturb_all(test_ct, gtv_test, chain)
        perturbed_retest = perturb_all(retest_ct, gtv_retest, chain)
        # Feature values of every perturbed copy  [§4.4]
        values_test = extract_each(perturbed_test)
        values_retest = extract_each(perturbed_retest)
        # Perturbation ICC: average of the two images  [§4.5]
        icc_by_chain[chain.name] = perturbation_icc(values_test, values_retest)
    return features_test, features_retest, icc_by_chain

Image processing order

The steps for one perturbed image, in the order the paper describes them, are §4.2 [supp. note 4]:

1. Rotate the image and mask in the axial plane by θ (optional). The mask is thresholded only after interpolation.
2. Add Gaussian noise (optional), with σ estimated once on the unperturbed image; round the result to integer HU.
3. Translate by a sub-voxel shift η (optional), performed together with the tri-linear interpolation onto the new grid.
4. Apply the anti-aliasing Gaussian before interpolation (see pseudocode below) and interpolate image and mask to the isotropic target spacing. Round intensities to the nearest integer HU and set the mask with a partial-volume threshold of 0.5 tab. 2.
5. Change the mask volume by a fraction τ (optional) [supp. note 4].
6. Randomise the contour with supervoxels (optional) [supp. note 4].
7. The intensity mask is copied from the morphological mask and re-segmented to the cohort range; the morphological mask is kept as it is §4.2.

Perturbation definitions

• Rotation (R): in-plane rotation of image and mask about the z-axis, angle θ in [−13°, 13°] §4.3 [supp. note 5]. Image sampling is bi-linear and the image is rounded to integer HU afterwards [supp. note 4].
• Noise (N): Gaussian noise with σ equal to the noise estimated from the original image (method below) §4.3 [supp. note 4].
• Translation (T): shift of the voxel grid by a fraction η of the target grid spacing, along x, y and z; the three axes are permuted over the listed fractions §4.3 [supp. note 5].
• Volume adaptation (V): grow or shrink the mask to a target volume V_a = floor(V_0 (1 + τ)), τ in [−0.28, 0.28] §4.3 [supp. note 4].
• Contour randomisation (C): supervoxels (simple linear iterative clustering, SLIC) are selected by their overlap with the expert contour §4.3 [supp. note 4].

Chains

Eighteen chains were built: the five basic perturbations, and 13 combinations of them. Every permutation of the settings of a chain was generated. Some combinations (for example rotation + translation + volume adaptation) were not tested because the number of permutations was too large. Noise alone and contour randomisation alone were repeated 30 times tab. 1 [supp. note 5].

Per CT image the sum of the 18 counts is 565 perturbed images (sum of the column in tab. 1). The paper applies perturbations to both the test and the retest image §4.5; this total is computed from the table and is not stated by the paper.

Features

• Base set: 182 features defined in IBSI §4.4 [supp. note 3]. 118 of them need discretisation, done with a fixed bin number (8, 16, 32, 64) or a fixed bin width (6, 12, 18, 24 HU), giving 8 settings tab. S2.
• Each feature is computed at 1, 2, 3 and 4 mm, so 1008 × 4 = 4032 features per image tab. S2 §4.4.

Robustness analysis

ICC(1,1) of Shrout and Fleiss is used §4.5. The paper does not write the formula or the layout of the perturbation ICC (which factor is the subject and which the rater) [pitfall below]. The R package psych was the basis of the code §4.5.

Cohorts

• NSCLC: 31 patients from a public archive; two CT scans within 15 minutes, same scanner and protocol §4.1 [supp. note 1].
• HNSCC: 19 patients in-house; two CT scans within 4 days, using two different protocols for PET attenuation correction §4.1 [supp. note 1].
• GTV delineated by radio-oncologists in a treatment planning system §4.1.

Hyperparameters

Results

Author-reported. Setup: CT, NSCLC (31) and HNSCC (19) test-retest cohorts, 4032 IBSI-based features, ICC(1,1) with a robustness threshold of 0.90 §2 §4. Baseline is test-retest robustness in the same cohort.

Other reported results: the perturbation chains' robust fraction, and false negatives, are in fig. 3 and fig. 4; ICC distributions and CI widths are in fig. S4 and fig. S5; the test of ICC bias between test and retest is in [supp. note 6] fig. S6. The paper concludes that chained perturbations of intensity and segmentation can stand in for test-retest imaging. NTVC gave a low false-positive count in both cohorts, though RC was lowest in NSCLC §2.3 §3.

Algorithms

Anti-aliasing σ [supp. note 2, the only numbered equation, eq. 1]:

# Gaussian pre-filter width before down-sampling  [supp. note 2, eq. 1]
# Input:  beta (0 < beta <= 1), spacing_orig and spacing_target (mm)
# Output: sigma (see Pitfalls for units)

def antialias_sigma(beta, spacing_orig, spacing_target):
    # Nyquist frequency of the target grid  [supp. note 2]
    nyquist = spacing_orig / (2 * spacing_target)
    # Solve exp(-nyquist^2 * sigma^2 / 2) = beta for sigma  [supp. note 2]
    sigma = sqrt(-2 * log(beta)) / nyquist
    return sigma

Noise level and noise addition [supp. note 4, "Noise"]:

# Noise level of a CT image  [supp. note 4, "Noise"]
# Input:  image (3-D, HU), estimated once on the unperturbed image
# Output: sigma_noise (HU)

def estimate_noise(image):
    # High-pass with coiflet-1 along x and y; how the two results are
    # combined is not stated  [supp. note 4]
    diff = highpass_x_then_y(image)
    # Median absolute value scaled to a standard deviation  [supp. note 4]
    return median(abs(diff)) / 0.6754

# Add noise after rotation, then round to integer HU  [supp. note 4]
def add_noise(image, sigma_noise):
    noisy = image + normal(0, sigma_noise, image.shape)
    return round(noisy)

Volume adaptation [supp. note 4, steps 1 to 6]:

# Grow or shrink a binary ROI to a target volume  [supp. note 4]
# Input:  roi (3-D bool), tau (float, fraction)
# Output: adapted roi (3-D bool)

def adapt_volume(roi, tau):
    # tau == 0 is listed in the settings but its handling is not stated;
    # returning the input is an assumption  [supp. note 5]
    if tau == 0:
        return roi.copy()
    # Target volume, rounded down  [supp. note 4, step 2]
    volume_0 = count(roi)
    volume_target = floor(volume_0 * (1 + tau))
    # Manhattan-distance-1 structuring element  [supp. note 4, step 3]
    element = cross_element_3d()
    # Last mask that did not overshoot the target  [supp. note 4, step 4]
    previous = roi.copy()
    volume_previous = volume_0
    while True:
        # Dilate for growth, erode for shrinkage  [supp. note 4, step 5a-b]
        if tau > 0:
            current = dilate(previous, element)
        else:
            current = erode(previous, element)
        volume_current = count(current)
        # Stop when empty, or when the target is overshot  [step 5d-f]
        if volume_current == 0:
            break
        if tau > 0 and volume_current > volume_target:
            break
        if tau < 0 and volume_current < volume_target:
            break
        previous = current
        volume_previous = volume_current
    # Fix the remaining voxels from the rim  [supp. note 4, step 6]
    if volume_previous != volume_target:
        # Rim between the overshoot mask and the last good mask  [step 6b]
        rim = logical_xor(current, previous)
        # Number of voxels still to add or remove  [step 6a]
        number = abs(volume_target - volume_previous)
        # Random voxels from the rim, without replacement  [step 6c]
        chosen = choose_without_replacement(rim, number)
        # Add for growth, remove for shrinkage  [steps 6d-e]
        if tau > 0:
            previous[chosen] = True
        else:
            previous[chosen] = False
    return previous

Contour randomisation [supp. note 4, steps 1 to 7]:

# Supervoxel-based contour randomisation  [supp. note 4]
# Input:  image (3-D, HU), roi (3-D bool, expert contour),
#         spacing_mm (1, 2, 3 or 4), g_low, g_high (re-segmentation range)
# Output: randomised roi (3-D bool)

def randomise_contour(image, roi, spacing_mm, g_low, g_high):
    # Crop to 25 mm around the ROI bounding box  [supp. note 4, step 1]
    image_c, roi_c = crop_around(image, roi, margin_mm=25)
    # Widen the range by 10% on each side and clip  [supp. note 4, step 2a]
    span = g_high - g_low
    r_low = g_low - 0.1 * span
    r_high = g_high + 0.1 * span
    clipped = clip(image_c, r_low, r_high)
    # Map to the unit interval  [supp. note 4, step 2b]
    scaled = (clipped - r_low) / (r_high - r_low)
    # Estimated supervoxel count, about 0.5 cm^3 each  [supp. note 4, step 3]
    voxel_cm3 = prod(voxel_spacing_mm) / 1000
    count_est = ceil(size(image_c) * voxel_cm3 / 0.5)
    # Gaussian smoothing with sigma = voxel spacing  [supp. note 4, step 4]
    smoothed = gaussian_filter(scaled, sigma=spacing_mm)
    # SLIC: compactness 0.05, supervoxel size 0.25 to 1.5 cm^3  [step 5]
    labels = slic(smoothed, count_est, compactness=0.05)
    # Overlap of each supervoxel with the expert ROI  [supp. note 4, step 6]
    overlap = {}
    for k in unique(labels):
        # Intersection over supervoxel size; see Pitfalls  [step 6c]
        overlap[k] = count((labels == k) & roi_c) / count(labels == k)
    # Highest-overlap supervoxel is always kept  [supp. note 4, step 7a]
    best = argmax(overlap)
    selected = []
    for k in overlap:
        if k == best or overlap[k] >= 0.90:
            selected.append(k)
        elif overlap[k] < 0.20:
            continue
        # Keep with probability equal to the overlap  [step 7d]
        elif random_uniform(0, 1) <= overlap[k]:
            selected.append(k)
    # Union of selected supervoxels, then morphological closing  [step 7e]
    new_roi = binary_closing(union_of(labels, selected), cross_element_3d())
    return new_roi

Robustness of one feature §4.5:

# Robust flag and perturbation ICC  [§4.5] [fig. 2]
# Input:  feature values (patients x repeats) for test or retest images
# Output: ICC(1,1) and robust flag

def is_robust(icc):
    # Robust when ICC(1,1) reaches the threshold  [§4.5]
    return icc >= 0.90

def perturbation_icc(values_test, values_retest):
    # ICC computed separately on the perturbations of each image  [§4.5]
    icc_test = icc_1_1(values_test)
    icc_retest = icc_1_1(values_retest)
    # Average of the two  [§4.5] [fig. 2]
    return (icc_test + icc_retest) / 2

Paper vs code

No official code.

Searches for this card: the paper contains no code link (its only URLs are for PyWavelets and Python) §4.2 [supp. note 4]; no project page exists. A GitHub repository search for "perturbation radiomic" returned only third-party repositories, which are not used. Searches of the authors' GitHub user account and of the OncoRay organisation found no repository for this method. The OncoRay organisation's oncoray/mirp package (EUPL-1.2) is the lab's IBSI-compliant radiomics package; its README does not mention perturbations or robustness, and the paper does not link it, so it is not used here.

Pitfalls

• Overlap formula is printed inverted. The paper writes the supervoxel overlap as voxels of the supervoxel divided by voxels of the intersection [supp. note 4, step 6c], which is ≥ 1, but states 0 ≤ η ≤ 1. The pseudocode uses intersection divided by supervoxel size. Confirm this.
• Anti-aliasing σ is printed in a form that is not real for β < 1. The paper's last line, σ = −2 (d2/d1) √(2 ln β) [supp. note 2], is negative under the root. Solving its own equation gives σ = √(−2 ln β) / ω_N with ω_N = d1 / (2 d2) [supp. note 2, eq. 1]. The units are not stated; the derivation uses the original voxel grid, so σ is taken in original-voxel units (inferred from the derivation).
• β is used twice. The anti-aliasing filter parameter (0.93 in tab. 2) and the SLIC compactness (0.05) [supp. note 4] are different quantities.
• τ = 0 is not handled. Several chains list τ = 0.0 [supp. note 5]. The described loop neither dilates nor erodes at τ = 0 and has no break condition for it. The pseudocode returns the input mask; the paper does not say this.
• Volume adaptation edge cases. The paper does not say what happens when the rim holds fewer voxels than the number still to change, or whether the final step uses the overshoot mask or the last valid mask (the pseudocode follows the text in using both).
• Noise estimate. The paper gives the coiflet-1 high-pass and the 0.6754 constant, but not how the x and y responses are combined (sum, product, or maximum). The noise σ is taken from the original image before any step [supp. note 4].
• Translation units. η is a fraction of the target grid spacing, applied per axis, so each axis is permuted independently [supp. note 5] §4.3. The paper uses "isotropic voxel spacing" in the same sentence; the target grid is the reading the paper's interpolation description supports.
• Mask thresholding order. Rotation of the mask is not thresholded until after interpolation [supp. note 4]. A naive implementation that thresholds earlier will give different masks.
• Random seeds. Noise, supervoxel selection, rim selection and the Moran's I / Geary's C approximation all use random draws; no seed is stated. Exact reproduction needs the seeds.
• Moran's I / Geary's C are approximated from 100 random voxels, repeated until the standard error falls below 0.002 [supp. note 3]. The stopping rule is stochastic and not otherwise specified.
• ICC design is not stated. The paper cites Shrout and Fleiss and says the psych package code was adapted §4.5, but does not say which factor is the target and which is the rater for the perturbation ICC, or whether the single-measure form is used. Read the code if you need this; the code is not public (see Paper vs code).
• Test-retest is itself confounded. HNSCC retests used two different PET protocols, with exposure about four times higher on average and two reconstruction kernels; NSCLC used one protocol §3 [supp. note 1]. The paper suggests this may explain the higher false-negative rate in NSCLC.
• Small cohorts. 31 and 19 patients. The CI widths are large (test-retest HNSCC average 0.35 §3). The headline percentages can move a lot with small changes in the reference.
• Robust fraction depends on the spacing and β choices. The main results use β = 0.93 tab. 2. The β sweep [supp. note 2] shows that HNSCC robustness changes strongly with β; choose and report β explicitly.
• Chain counts. Table 1 counts agree with the setting lists in [supp. note 5] for all 18 chains (for example, RV: 6 rotation angles × 5 volume fractions = 30). The text's "NTCV" in §3 is a typo for NTVC.
• Function names. The paper names skimage.segmentation.slic_superpixels for SLIC, with scikit-image 0.13.1 [supp. note 4]. Check the function name and call signature (supervoxel count, size bounds) in the scikit-image version you use; the paper's name is not confirmed here.
• Feature averaging. For GLCM and GLRLM, features are computed per direction and averaged, not computed on an averaged matrix [supp. note 3].
• Perturbation counts are not uniform. Chain sizes range from 27 to 40 images tab. 1; the paper keeps this to limit sample-size effects on ICC. Compare chains only within the same cohort.

Lineage

• Builds on: IBSI, the Image Biomarker Standardisation Initiative, for feature definitions and image processing §4.2 §4.4 (not a card here).
• Builds on: Chang's noise estimate and SLIC supervoxels [supp. note 4] (not cards here).
• Compared against: test-retest ICC from the same two cohorts §2.
• Superseded by: not stated in this paper.