Radiomics is a computational imaging approach that extracts dozens or even hundreds of quantitative features from tumor images, such as texture, shape, and intensity distributions, then feeds them into machine-learning models to predict clinical outcomes. In sarcoma research, radiomics applied to MRI has already been used to estimate tumor grade, forecast metastatic relapse, and predict response to chemotherapy and radiotherapy. The promise is a non-invasive, objective window into tumor biology that outpaces what the human eye can detect.
The central problem addressed by this paper is a structural flaw in MRI: unlike CT, where Hounsfield units are physically standardized, MRI signal intensities (SIs) have no universal scale. The same patient scanned on a Siemens machine and a Philips machine will produce images with completely different raw numbers, even if the tissue is identical. This becomes critical in multicenter radiomics studies, where patient data are pooled from institutions with different scanners, coils, acquisition protocols, and field strengths. Without correcting for these differences, any features extracted across sites are not directly comparable, and any predictive model built on them may be modeling scanner differences rather than tumor biology.
Intensity harmonization techniques (IHTs) are post-processing algorithms designed to normalize or align the SI distributions across scans before radiomics feature extraction. In neuroimaging, several robust IHTs have been developed and validated, but their application to body imaging, particularly sarcomas, is far less established. Sarcomas present a special challenge because they arise in virtually any anatomical location, surrounded by different tissue types, which invalidates normalization strategies that rely on a fixed reference organ like cerebrospinal fluid.
This paper tackles, for the first time, the question of how much the choice of IHT actually changes radiomics feature values and, more critically, whether it changes the clinical conclusions drawn from those features. The authors go beyond just showing that features shift numerically. They demonstrate that the predictive models built from those features, including both unsupervised patient clustering and supervised metastasis prediction, can reach entirely different conclusions depending solely on which harmonization method was chosen.
The study population comprised 70 patients with newly diagnosed, non-metastatic, histologically confirmed high-grade soft-tissue sarcomas (STS) of the trunk wall or extremities. All patients were treated at the Bergonie Institute in Bordeaux, France, a national sarcoma reference center, between June 2006 and November 2016. The cohort was assembled by working backward from a starting pool of 163 eligible patients: 133 had received 4 to 6 cycles of anthracycline-based neoadjuvant chemotherapy followed by curative surgery, 95 had available baseline MRI, 72 had T2-weighted images without artifacts, and a final 70 had complete clinical and radiological follow-up for at least 2 years post-surgery.
Demographics and histology: Median age was 58 years (range 19-84), with 38 men (54.3%) and 32 women (45.7%). The dominant histological subtype was high-grade undifferentiated sarcoma (31/70, 44.3%), followed by synovial sarcoma and rhabdomyosarcoma (8/70 each, 11.4%), leiomyosarcoma and myxoid/round cell liposarcoma (6/70 each, 8.6%), and pleomorphic sarcoma (3/70, 4.3%). Median tumor diameter was 106 mm (range 40-273 mm), median volume was 220 cm3 (range 10.2-3,084 cm3), and 65 of 70 tumors (92.9%) were deep-seated. The lower limb was the most common location (35/70, 50%).
Treatment and outcomes: Most patients received 5 to 6 cycles of chemotherapy (52/70, 74.3%), predominantly anthracycline plus ifosfamide (64/70, 91.4%). Adjuvant radiotherapy was given in 65/70 cases (92.9%), and curative resection achieved R0 margins in 41/70 (58.5%). Histological response, defined as less than 10% viable tumor cells on post-chemotherapy surgical specimens, was achieved in only 16/70 patients (22.9%). A total of 29 of 70 patients (41.4%) experienced metastatic relapse within the first 2 years of follow-up. This was the primary outcome used to train and test supervised classification models.
MRIs were acquired on 3 different 1.5-Tesla systems: Siemens MAGNETOM Aera (41/70, 58.5%), Philips Signa (17/70, 24.3%), and General Electric Healthcare Optima (12/70, 17.1%). Repetition times ranged from 2,400 to 4,500 ms, echo times from 70 to 130 ms, and slice thicknesses from 3 to 5 mm. This scanner heterogeneity is precisely what makes harmonization necessary and what gives the study's conclusions real-world applicability.
After standardizing voxel size to 1x1x4 mm3 with b-spline interpolation and applying N4 bias field correction, a senior radiologist with 4 years of sarcoma imaging experience manually segmented the entire tumor volume, slice by slice, using LIFEx software (version 5.10). All segmentations were validated by a second radiologist with 28 years of sarcoma experience. From these segmented volumes, the authors applied 5 distinct intensity harmonization techniques to produce 5 parallel imaging datasets, plus 1 unharmonized dataset (No-IHT), for a total of 6 radiomics datasets per patient.
The five harmonization techniques: (1) IHTstd performs classical min-max normalization, rescaling each scan's voxel intensities to a 0-1 range based on the image's own minimum and maximum. (2) IHTfat normalizes by dividing all intensities by the mean signal of a manually segmented volume of at least 10 cm3 of healthy adipose tissue in that patient's scan, a natural internal reference signal. (3) IHTHM.1 applies histogram matching to a single randomly chosen patient's histogram from the dataset, aligning each scan's intensity distribution to that reference using piecewise linear mapping at 100 landmarks. (4) IHTHM.All matches each scan's histogram to the average histogram of the entire study population rather than a single patient. (5) IHTHM.All.C applies IHTHM.All first, then adds ComBat harmonization, a statistical method originally developed to remove batch effects in genomics data, using a non-parametric setting to correct for scanner-specific additive and multiplicative effects without assuming a specific probability distribution.
Radiomics feature extraction: From each of the 6 datasets, 45 three-dimensional radiomics features were extracted using LIFEx, with SIs discretized into 128 fixed bins. The feature set included 13 histogram-based features and 32 second-order texture features derived from the grey-level co-occurrence matrix (GLCM, n=7), grey-level run length matrix (GLRLM, n=11), neighborhood grey-level difference matrix (NGLDM, n=3), and grey-level zone length matrix (GLZLM, n=11).
Three layers of analysis: The authors then assessed IHT impact at three levels. First, they compared raw feature values across harmonization techniques using one-way repeated-measures ANOVA with Bonferroni-corrected Tukey post-hoc tests, and computed intraclass correlation coefficients (ICCs) to measure feature stability. Second, they performed unsupervised hierarchical clustering with Ward linkage on each of the 6 feature sets, selected 2 clusters per dataset, and tested whether those clusters predicted metastatic-relapse free survival (MFS) in multivariate Cox models adjusted for 8 clinical covariables. Third, using the same 6 feature sets, they trained a binomial elastic-net logistic regression (combining LASSO and ridge penalization) to predict 2-year metastatic relapse in 50 training patients, evaluated it with 5-times-repeated 10-fold cross-validation, and tested it on an independent 20-patient cohort. Hyperparameters alpha (mixing percentage) and lambda (regularization strength) were tuned via manual grid search.
The first and most fundamental finding is that every single one of the 45 radiomics features was significantly altered by the harmonization technique applied (p-values ranged from less than 0.0001 to 0.02 across all features). This is not a marginal numerical drift; the Bonferroni-corrected post-hoc comparisons showed that for many pairs of harmonization techniques, the majority of features were statistically different. This means that a researcher who extracts features without harmonization, or who uses a different harmonization choice than a previous study, may produce a completely different feature matrix even from the same underlying tumor images.
Pairwise differences between techniques: The largest divergence was between IHTHM.All and IHTfat, where 31 of 45 features (68.9%) differed significantly. The histogram-matching approaches diverged substantially from the fat-normalization approach in general: IHTHM.All.C versus IHTfat also showed 30/45 (66.7%) different features, and IHTHM.1 versus IHTfat showed the same 30/45. The most similar pair was IHTHM.All versus IHTHM.1, which differed on only 6/45 features (13.3%), suggesting these two histogram-matching variants produce closely related results despite using different reference histograms. IHTHM.All and IHTHM.All.C differed on 13/45 features (28.9%), showing that adding ComBat on top of histogram-matching does introduce additional feature-level changes.
Intraclass correlation coefficients: ICCs quantify how consistently a feature is reproduced across the 6 harmonization conditions. The highest ICCs (stability across methods) were observed for GLRLM_RLMNU, GLRLM_GLNU, and GLCM_Correlation, all reaching values of 0.95 or above, indicating that these particular texture features are relatively robust to harmonization choice. At the other extreme, GLZLM_ZLNU, GLZLM_LZE, HISTO_maximum, GLZLM_LZLGE, and HISTO_minimum all fell below an ICC of 0.20, meaning they are extremely sensitive to the harmonization method and essentially unreliable for cross-site comparison without careful standardization.
These results have a direct practical implication: when a published sarcoma radiomics model lists specific texture features as predictive, the reader cannot assume those features would be equally predictive if extracted after a different harmonization step. The identity of the "best" features is itself an artifact of the processing pipeline, not just of the tumor.
The second layer of analysis tested whether radiomics-based patient clustering predicted metastatic-relapse free survival (MFS). Hierarchical clustering with Ward linkage was applied to each of the 6 feature datasets, consistently producing 2 clusters per dataset, but the actual patient assignments differed substantially across harmonization conditions. The Cohen's Kappa and Baker's gamma coefficient matrices quantified how different these clusterings were from one another. The most similar pair of cluster solutions came from IHTHM.All versus IHTHM.All.C (Kappa = 0.75, Baker coefficient = 0.55), confirming that ComBat does not radically restructure the clusters when applied after IHTHM.All. The most divergent pair was No-IHT versus IHTHM.1 (Kappa = 0.18, Baker coefficient = 0.05), meaning only 18% concordance in patient assignments between the raw unharmonized data and the single-patient histogram-matched data.
Prognostic value by harmonization method: In univariate Kaplan-Meier analysis, only the IHTHM.All dataset yielded clusters with significantly different MFS (log-rank p = 0.03). The other five datasets produced clusters whose survival curves were not statistically separable at univariate level. This finding alone demonstrates that the choice of IHT determines whether a radiomics study concludes that tumor imaging captures prognostic information or not.
Multivariate Cox model results: After adjustment for 8 clinical and pathological covariables (performance status, histotype, tumor diameter, chemotherapy type and cycles, surgical margins, histological response, and adjuvant radiotherapy), 4 of the 6 datasets yielded clusters that were independently prognostic for MFS. The IHTHM.All dataset showed the strongest independent association (HR = 4.72, 95% CI 1.64-13.56, p = 0.004). IHTstd clusters were also significant (HR = 3.26, 95% CI 1.48-7.71, p = 0.007), as were No-IHT clusters (HR = 2.64, 95% CI 1.15-6.04, p = 0.02) and IHTHM.All.C clusters (HR = 2.89, 95% CI 1.19-7.05, p = 0.02). By contrast, IHTfat and IHTHM.1 clusters were not independently prognostic, meaning those harmonization approaches discarded the prognostic signal rather than preserving it.
Concordance indices for the 6 multivariate survival models ranged from 0.71 (95% CI 0.67-0.75) for IHTHM.1 to 0.75 (95% CI 0.70-0.79) for No-IHT, IHTstd, and IHTHM.All. The reference clinical-only model achieved a concordance index of 0.71 (95% CI 0.67-0.75). This indicates that only some harmonization approaches allow the radiomics features to improve upon a clinical baseline; the wrong choice adds nothing.
The supervised analysis trained elastic-net logistic regression classifiers to predict whether a patient would develop metastatic relapse within 2 years of surgery. Training used 50 patients (21/50 events, 42%) with 5-times-repeated 10-fold cross-validation; testing used an independent 20-patient cohort (8/20 events, 40%). The same random patient split was applied across all 6 feature datasets to ensure that any performance differences were attributable only to the harmonization choice, not to data partitioning.
Cross-validation performance: In the training cohort, the best cross-validated AUROCs came from IHTHM.All (0.71, 95% CI 0.66-0.76) and IHTHM.1 (0.69, 95% CI 0.64-0.74). The No-IHT dataset performed worst in cross-validation (0.57, 95% CI 0.52-0.63), barely above chance. The hyperparameter tuning revealed considerably different optimal regularization settings across datasets. For example, IHTHM.All.C selected a high lambda (0.840) suggesting aggressive feature reduction, while No-IHT used a lower lambda (0.114).
Test cohort performance: Results on the 20-patient test set showed even larger variation between methods. In descending order: IHTHM.1 achieved the best AUROC of 0.82 (95% CI 0.59-1.0), IHTfat reached 0.80 (95% CI 0.56-1.0), IHTHM.All reached 0.77 (95% CI 0.52-1.0), No-IHT reached 0.76 (95% CI 0.50-1.0), IHTHM.All.C reached 0.71 (95% CI 0.44-0.97), and IHTstd reached the lowest at 0.69 (95% CI 0.41-0.89). Notably, IHTstd was among the better performers in cross-validation but fell to last in the test cohort, while IHTHM.1 performed modestly in cross-validation but led on the test set. This reversal highlights how small test sets (n=20) make ranking unstable.
Feature selection: The number of features retained in the final models ranged from 3 (No-IHT and IHTHM.All.C) to 21 (IHTfat). The best-performing model (IHTHM.1) retained 7 features out of 48 possible (including 3 shape features). HISTO_Quartile1 and GLZLM_SZLGE were the most commonly selected features across models, appearing in 5 of 6 and 4 of 6 models respectively, suggesting some degree of cross-technique robustness for these specific features.
The discussion frames the core practical takeaway: no harmonization technique demonstrated clear superiority across all three analyses. IHTHM.All led in unsupervised prognostic clustering, while IHTHM.1 and IHTfat led in supervised test-set prediction. IHTstd was independently prognostic in the Cox model but performed worst in the test cohort. This inconsistency suggests that the "best" technique is dataset-specific and task-specific, and that treating harmonization as a fixed methodological decision rather than an empirically optimized one is scientifically inadequate.
Treating IHT as a hyperparameter: The authors propose that IHT should be treated the way voxel size and grey-level discretization are already treated in mature radiomics workflows: as a pipeline parameter that needs to be explored and selected based on the dataset at hand. Just as models are often trained with multiple feature extraction settings to identify the most predictive version, harmonization options should be explored in a training subset and the most robust choice validated prospectively. This is analogous to standard hyperparameter tuning but applied to the image processing stage rather than the model stage.
ComBat harmonization in context: IHTHM.All.C showed moderate and relatively consistent performance across both unsupervised and supervised analyses (Kappa 0.75 vs. IHTHM.All, concordance index 0.73, test AUROC 0.71), which the authors suggest may reflect a more biologically realistic normalization because it explicitly separates scanner-related technical variation from biological signal. However, they caution that the covariables supplied to the ComBat function must be carefully chosen: the clinical outcome should never be included, and for anatomically diverse sarcomas, location-specific acquisition variables may need to be added to prevent ComBat from correcting away biologically relevant differences.
The authors also note that the RAVEL algorithm, originally developed for brain MRI using cerebrospinal fluid as a control region, might be adaptable to body imaging using healthy adipose tissue as the control region in sarcoma. Alternatively, prospective harmonization at the acquisition stage through T1-mapping or T2-mapping sequences would bypass the problem entirely, though this approach cannot retroactively fix the thousands of already-collected MRIs stored in institutional archives. The review by the same group found that 17 of 31 (54.8%) MRI-based sarcoma radiomics studies did not describe their intensity harmonization method at all, a reporting gap that makes cross-study comparison essentially impossible.
The most prominent limitation is the relatively small study population of 70 patients, all from a single institution. While the authors note this is the largest study to investigate IHT in sarcoma radiomics at the time of publication (most sarcoma imaging studies do not exceed similar numbers), the 20-patient test set makes AUROC estimates extremely uncertain, as reflected in the wide confidence intervals (e.g., IHTHM.1: 0.82, 95% CI 0.59-1.0). The difference in AUROC ranking between cross-validation and the test set, where IHTstd reversed from near-top to last, is a cautionary example of how sampling variance can drive apparent method comparisons at small sample sizes.
Sequence limitation: The study analyzed only T2-weighted images, chosen because T2-WI is the most informative morphological MRI sequence for sarcomas. However, the influence of IHT on radiomics features from T1-weighted, contrast-enhanced T1-weighted, dynamic contrast-enhanced (DCE), and diffusion-weighted imaging sequences remains unknown. Each of these sequences has different tissue contrast mechanisms and different relationships to the biological parameters that drive SI variability, so findings from T2-WI cannot be assumed to transfer.
Circular reasoning in IHT validation: A fundamental methodological challenge noted by the authors is the difficulty of validating which IHT is "correct" using model performance as the criterion. This approach is only valid if MRI texture features truly have prognostic value for sarcomas (so that a technique that degrades performance is genuinely introducing noise). The authors cite prior studies confirming this relationship, which they argue justifies the approach, but acknowledge that an alternative validation strategy using phantoms or healthy-volunteer MRIs would be more rigorous in the absence of established ground truth.
Feature set scope and ComBat covariables: The 45 features were intentionally restricted to those defined by the Imaging Biomarker Standardization Initiative and available in LIFEx, excluding many features used in other software packages. Expanding the feature space would have introduced overfitting risk given the limited cohort size. Additionally, ComBat was only applied to the IHTHM.All dataset; future work should apply it to all five IHT options to more comprehensively assess its interaction with different harmonization bases. For sarcomas specifically, the tumor location may need to be included as a ComBat covariate because different body regions require different acquisition parameters, but this also risks adjusting away genuine variation in tumor phenotype by location.