Enhancing PET/CT Radiomics Robustness Through Graph Signal Processing

1. Introduction

Prostate cancer (PCa) is one of the leading causes of cancer-related death among men. In autopsy studies, bone metastases are present in approximately 90.1% of men who died from metastatic prostate cancer [1,2,3]. Tumor cells detach from the primary tumor through epithelial-to-mesenchymal transition (EMT) and prepare optimal proliferative conditions through exosome release. After prior intravasation, tumor cells reach bone tissue and begin forming metastases through either the destruction (osteolytic metastases) or disorganized formation (osteoblastic metastases) of bone tissue, causing complications that may lead to physical, metabolic, and neurological effects [4]. For adequate prevention, it is essential to have an accurate diagnosis to identify an appropriate interventional or therapeutic pathway. However, it is not uncommon for a patient’s clinical profile, due to factors related to complex etiology or possible intra-subject heterogeneity, to fall within a clinical context where confounders are numerous and the discriminative elements between physiological and pathological conditions are limited or not easily identifiable through visual inspection of images derived from conventional diagnostic examinations [5]. In this complex scenario, radiomics has emerged as a computational imaging approach that has taken on a leading role in clinical research, owing to its effectiveness in providing concrete support tools from both diagnostic and therapeutic perspectives [6,7,8,9,10,11]. Radiomics is based on a set of procedures that allow the extraction of quantitative information or “features” from medical images [12,13], which are used for predictive modelling [10,14,15],with the aim of discriminating the clinical condition of a patient. A standard radiomics pipeline includes preprocessing, segmentation of the volume of interest [16], feature extraction [17], and model development with machine learning algorithms [18]. However, the workflow of a standard radiomics analysis is characterized by intrinsic limitations [19,20]. Scanner characteristics, including hardware and acquisition parameters, define the voxel grid on which images are reconstructed, while segmentation determines which voxels are included in the analysis. Therefore, differences in scanner hardware and acquisition settings, together with the inter/intra operator variability introduced by manual segmentation, directly impact radiomics features, as they determine the voxel grid that defines the data and the intensity scales associated with each voxel. Although research efforts have aimed to automate segmentation and statistical methods such as ComBat [21,22] have been proposed to mitigate scanner-induced variability, none of these approaches fully resolve the intrinsic dependence of radiomics features on acquisition, segmentation, and voxel discretization. It is also common that radiomics features fail to capture the complex morpho-functional patterns underlying specific clinical conditions, limiting their ability to fully represent the lesions’ phenotype [23,24]. The construction of robust models is therefore often precluded by the intrinsic nature of the input data. Although radiomics originated as a stand-alone technique characterized by standard procedures in the literature, the data and objects processed during the analysis are sufficiently versatile to enable the integration of approaches and results suitable for completely different scientific and operational fields. This study introduces the concept of point clouds into radiomics analysis. A point cloud is an unstructured set of points that may derive from or represent volumes and specific spatial structures [25]. The connection with the underlying structure is latent but can be made explicit by applying appropriate processing techniques. This scenario allows, with appropriate precautions, the mapping of data into a completely different analytical domain compared to that in which the target volume is defined, reducing the dependence of the representation on the original grid. To make the new representation robust and informative, we applied a strategy to compensate for the controlled loss of information caused by extracting the point cloud from the volumes. The solution we adopted to recover the underlying geometric information consisted of using the extracted point clouds as input to construct meshes [26] that accurately approximated the volumes, allowing us to integrate graph signal processing (GSP) [27] tools for feature extraction. In this study, 50 patients, for a total of 92 lesions, were selected to evaluate the pipeline. For each patient, standard radiomics features were extracted in accordance with the Image Biomarker Standardisation Initiative (IBSI) [28], along with the features derived from the proposed analysis. The study was structured around three primary objectives: (i) to validate the additional information content of the proposed features; (ii) to verify the robustness of GSP predictors to the variability introduced by segmentation and to batch effects; and (iii) to assess their ability to build robust predictive models, either as stand-alone inputs or in combination with classical features. For the first objective, a patient-level Spearman correlation analysis was performed [29]. For the second objective, classical and GSP features were recomputed after introducing perturbations to lesion segmentation masks, and percentage differences were analyzed using descriptive statistics to evaluate sensitivity to segmentation variability. Batch effect was further investigated at the feature level by stratifying patients by scanner type and testing group differences using the Kruskal–Wallis test [30]. For the third objective, the classification performances were estimated using Receiver Operating Characteristic (ROC) analysis and accuracy. The Results section reports the outcomes of the analyses and classification performances, and the Discussion addresses the strengths and limitations of the proposed framework.

2. Materials and Methods

In this study, each lesion was segmented using masks derived from computed tomography (CT) scans. From the resulting segmented volumes, IBSI-compliant radiomics features [28] were extracted from the corresponding positron emission tomography (PET) images. Subsequently, a point cloud was extracted from the volume, the mesh was constructed, and GSP analysis was applied to compute the proposed features. From each lesion, 206 classical features (comprising 49 intensity, 21 shape, and 136 textures features) and 77 GSP features were extracted. At the end of the analysis, the two sets of features and the combined one were used to perform a patient-level correlation analysis, to verify the robustness of GSP predictors to the variability introduced by segmentation and to batch effects, and to train an LDA classifier and an SVM classifier; for each classifier, three different models were trained evaluating their predictive performance.

2.1. The Dataset

The dataset consisted of PET/CT images acquired from 50 patients, for a total of 92 lesions, with two different scanners and stored in DICOM format. For each patient, an average of 290 slices was acquired. Each lesion was associated with an identifying label, for a total of 75 benign bone lesions and 17 malignant lesions. The dataset used in this study was derived from a previously conducted study and consists of a subset data that have been fully and irreversibly anonymized from Cannizzaro Hospital of Catania (Italy). In particular, a destructive anonymization procedure was applied, ensuring that no subject can be re-identified. The original study obtained informed consent from all participants, including consent for the use of anonymized data for secondary research purposes. Due to its origin and privacy constraints, the dataset is not publicly available, but it can be shared by the authors upon reasonable request and subject to appropriate data protection agreements. The original study from which the dataset was derived should be cited as follows [31].

2.2. Preprocessing

The entire pipeline is designed to operate on PET images, which by construction are inherently noisy and characterized by very low spatial resolution, due to the technological limitations of the scanner components and low-count phenomena. To correctly extract the diagnostic content from the images, we developed a specific filtering pipeline, aimed at mitigating the intrinsic issues of PET data.

An adaptive Exponentially Weighted Moving Average (EWMA) filter has been implemented to improve the local Signal to Noise Ratio (SNR). The filter output is then convolved with a Sharpening kernel. Since the targets of interest generally stand out from the surrounding tissue, the filtering pipeline is designed to emphasize these transitions while simultaneously reducing residual noise.

2.2.1. The Adaptive EWMA Filter

The EWMA filter is a particular variant of the classic moving average widely used in numerical signal analysis. The main issue with the standard moving average, which makes it unsuitable in this context, concerns the inability to robustly control the amount of smoothing. This could seriously compromise the result, with a high possibility of introducing excessive blurring, thereby eliminating diagnostic information. The EWMA variant partially solves this problem thanks to the possibility of actively adjusting the amount of smoothing applied during filtering. Let us consider the one-dimensional sequence $x_n$and assume it is passed through an EWMA filter. The output $y_n$will have the following form:

$$y_n=\alpha x_n+\left(1-\mathsf{\alpha }\right)y_{n-1},\hspace{1em}0<\alpha <1$$

(1)

The smoothing action is weighted through the parameter $\alpha$. A low weight factor emphasizes the contribution of the previous output, reducing the influence of the current input, while a high weight factor, conversely, emphasizes the current input and reduces the contribution of the previous output. Low weight factors therefore correspond to very aggressive filtering, whereas high values favor fidelity to the input sequence, performing much lighter smoothing. Based on the one-dimensional formulation, we developed a two-dimensional extension of the filter in the image domain.. Let $x[n_1,n_2]$be the input image;

the output $y[n_1,n_2]$ will take the following form:

$$\{\begin{array}{c}{y}^{\left(h\right)}\left[n_1,n_2\right]=\left(1-\mathsf{\alpha }\right)\hspace{0.17em}{y}^{\left(h\right)}\left[n_1,n_2-1\right]+\alpha \hspace{0.17em}x\left[n_1,n_2\right]\\ y\left[n_1,n_2\right]=\left(1-\mathsf{\alpha }\right)\hspace{0.17em}y\left[n_1-1,n_2\right]+\alpha \hspace{0.17em}{y}^{\left(h\right)}\left[n_1,n_2\right]\end{array},0<\alpha <1$$

(2)

Filtering is performed in two steps: first, a row-wise filtering is performed producing the temporary output $y^{\left(h\right)}[n_1,n_2]$; subsequently, this is used as input for column-wise filtering, producing the final output $y[n_1,n_2]$. Although the EWMA filter in this form allows global filtering control, it does not allow control or modulation of local smoothing. This could be critical, especially in clinical contexts where local gradients are usually high due to anatomical complexity. Therefore, we developed an adaptive EWMA, capable of actively modulating the smoothing parameter based on the local image gradient. The idea is to apply strong smoothing in flat, low-gradient areas (usually less informative) and almost no smoothing in high-gradient areas rich in details, such as edges and corners. This approach improves local SNR while preserving high spatial frequency structures and avoiding blurring the target boundaries that will be segmented. To estimate the image gradient, Sobel filters were used: two FIR kernels that allow gradient estimation. Their matrix definitions are:

$$S_{\mathrm{row}}\left[n_1,n_2\right]=\left[\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right]$$

(3)

$$S_{\mathrm{column}}\left[n_1,n_2\right]=\left[\begin{array}{ccc}-1& -2& -1\\ 0& 0& 0\\ 1& 2& 1\end{array}\right]$$

(4)

Convolving the original image with $S_{\mathrm{r}\mathrm{o}\mathrm{w}}$and $S_{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{m}\mathrm{n}}$estimates the gradient components along rows and columns, respectively. Let $G_x[n_1,n_2]$be the row component and $G_y[n_1,n_2]$the column component. The gradient is computed as:

$$G\left[n_1,n_2\right]=\sqrt{G_x^2\left[n_1,n_2\right]+G_y^2\left[n_1,n_2\right]}$$

(5)

The gradient is then normalized to [0,1] to make smoothing control scale-invariant. The modulation of the weight parameter $\alpha$is performed using a Gaussian kernel applied to the normalized gradient:

$$\mathsf{\alpha }\left[n_1,n_2\right]={\mathsf{\alpha }}_{min}+\left({\mathsf{\alpha }}_{max}-{\mathsf{\alpha }}_{min}\right)\left(1-\exp \left(-{\displaystyle \frac{G_{\mathrm{norm}}^2\left[n_1,n_2\right]}{2{\mathsf{\sigma }}^2}}\right)\right)$$

(6)

This relationship describes the behavior discussed above. The parameter $\sigma$determines the weight parameter’s sensitivity to the gradient, while ${\alpha }_{\mathrm{m}\mathrm{i}\mathrm{n}}$and ${\alpha }_{\mathrm{m}\mathrm{a}\mathrm{x}}$indicate the minimum and maximum weight values, corresponding to the maximum and minimum smoothing.

2.2.2. The Sharpening Filter

The second filtering block is a Sharpening kernel. It is a FIR filter that extracts the high spatial frequency components of the image (e.g., edges, corners) and adds them to the input image, in this case the output of the adaptive EWMA. The result coincides with the original image enriched with details. We adopted a 3x3 Sharpening kernel to avoid excessively amplifying residual noise, which would negate the EWMA effect. A basic 3x3 kernel has the following form:

$$\mathrm{Sharp}\left[n_1,n_2\right]=\left[\begin{array}{ccc}0& -1& 0\\ -1& 5& -1\\ 0& -1& 0\end{array}\right]$$

(7)

Let us consider the matrices:

$$\delta \left[n_1,n_2\right]=\left[\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right]$$

(8)

$$\mathrm{Lap}\left[n_1,n_2\right]=\left[\begin{array}{ccc}0& -1& 0\\ -1& 4& -1\\ 0& -1& 0\end{array}\right]$$

(9)

where $\delta [n_1,n_2]$represents the spatial impulse and $\mathrm{L}\mathrm{a}\mathrm{p}[n_1,n_2]$represents a Laplacian kernel. The latter is a high-pass filter that extracts high spatial frequency components by subtracting from the current element the average of neighboring elements. If the area is flat, the filter output is zero; if the current element stands out from its context, the filter detects it. Comparing the 7 with the 8 and 9, the Sharpening kernel can be rewritten as the sum of the spatial impulse and the Laplacian kernel. Let $x[n_1,n_2]$be the input image, then the output $y[n_1,n_2]$can be written as:

$$y\left[n_1,n_2\right]=x\left[n_1,n_2\right]\mathrm{*}\left(\delta \left[n_1,n_2\right]+\mathrm{Lap}\left[n_1,n_2\right]\right)=x\left[n_1,n_2\right]+x_H\left[n_1,n_2\right]$$

(10)

where $x_H[n_1,n_2]$is the high spatial frequency component of the input image.

2.3. Standard Radiomics Features

Classical radiomic features are quantitative descriptors used to characterize specific properties of segmented volumes. According to IBSI guidelines [28] three main classes of standardized descriptors are defined:

• Intensity features: This category includes all first-order statistical metrics that can be extracted from the gray level distribution related to the segmented volume.
• Shape features: This class of metrics describes the geometric characteristics of the target.
• Texture features: This last category of features includes metrics derivable from the spatial distribution of specific patterns associated with voxel intensity.

2.4. Segmentation and Classical Feature Extraction

The pipeline was applied to PET images, as these provided information for the functional patterns of interest, whereas CT was used for anatomical reference [32,33] to create the masks delineating the lesions. For feature extraction, the Medical Imaging Toolbox of MATLAB [34] was employed; specifically, upon loading and filtering the images, the metadata were extracted and combined with the volumetric data to create a radiomics object, from which the standard features were extracted using the toolbox’s built-in functions with default parameters.

2.5. Point Cloud Extraction and Mesh Construction

From the segmented volume, a pseudorandom subset of voxels was extracted, while strictly controlling the size of the point cloud. Through preliminary experiments, a maximum size of 600 points was determined as the optimal trade-off between the goodness of the geometric representation and computational times. Automatic subsampling was applied to point clouds exceeding this limit, reducing the number of points to 600. Once the point cloud was extracted, we applied the Delaunay algorithm [35] to construct the mesh. It is a widely used algorithm in Computer Vision due to its ability to build robust and reliable meshes. Delaunay triangulation takes the points as input and connects them with tetrahedra such that no point of the point cloud lies inside the circumscribed sphere of any tetrahedron. This property allows the construction of regular and well-conditioned tetrahedra, avoiding distortion of local geometric curvature while reducing dependence on the regular grid on which the target volume is defined. As a toy case, Figure 1 shows an example of a spherical point cloud (a) and its corresponding Delaunay triangulation (b).

2.5.1. Mesh-Based Shape Features

Let us consider only the external triangles of a Delaunay triangulation obtained by considering only the extreme points from the original point cloud. Let $\mathcal{T}\left(v_1,v_2,v_3\right)$denote a generic triangle of the mesh defined by vertices $v_1,v_2,v_3$. Let vectors $\mathrm{a},\mathrm{b}\in {\mathbb{R}}^3$be defined as:

$$a=v_2-v_1,\hspace{1em}b=v_3-v_1$$

(11)

We define now the unit normal:

$$\widehat{n}={\displaystyle \frac{a\times b}{\Vert a\times b\Vert }}$$

(12)

Let $\theta$(in degrees) be the angle between $\widehat{n}$and the mean normal $\widehat{n_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}}$:

$$\mathsf{\theta }=\arccos \left(\widehat{n}\cdot \widehat{n_{\mathrm{mean}}}\right)\cdot {\displaystyle \frac{180}{\mathsf{\pi }}}$$

(13)

From the Delaunay triangulation, the boundary triangles were identified, and for each triangle the normal vector was computed together with its angular dispersion relative to the mean normal vector. From the multivariate distribution built from the coordinates of all triangle normals, univariate statistics such as mean, median, percentile-based descriptors, maximum, minimum, variance, skewness, and kurtosis were computed for each coordinate distribution. The same analysis was applied to the angular distribution. These statistical descriptors, particularly dispersion and shape metrics, provided complementary measures of global surface orientation and geometric irregularity. To identify the presence of preferential directions, the statistical entropy of angular distribution was computed according to the following formulation:

$$\mathcal{H}=-{\sum }_{i}p\left(i\right)\hspace{0.17em}\log p\left(i\right)$$

(14)

where $p\left(i\right)$ represents the probability that an observation falls into the $i$-th bin of the normalized histogram. A high entropy indicates the absence of preferential orientations, suggesting a more isotropic surface configuration, whereas lower entropy values reflect the presence of dominant directional patterns.

2.5.2. Connectivity Features

Let us consider the complete mesh, including internal triangle connections. We can now define $G=\left(V,E\right)$ as the graph associated with the point cloud, where $\mathrm{V}$ denotes the set of points and E the set of connections (edges) between them. The adjacency matrix $\mathrm{A}$is defined as the binary matrix that describes the connections between nodes [36]. Each element of the matrix takes the value 1 if the corresponding pair of nodes is connected, 0 otherwise:

$${\mathit{A}}_{\mathit{i}\mathit{j}}=\{\begin{array}{c}:1\hspace{1em}\mathrm{i}\mathrm{f}\left(v_i,v_j\right)\in E\\ 0\hspace{1em}otherwise\end{array}$$

(15)

Considering the generic node $v_i$, we can define its degree as the number of nodes to which it is directly connected [28]:

$$d_i=\sum _{j=1}^N{A}_{ij}$$

(16)

where $d_i$denotes the degree and $N$ the number of nodes. From this last definition, it is possible to introduce the Laplacian operator [36] as the difference between the diagonal matrix which contains the degrees of all nodes composing the graph, and the adjacency matrix. Let us consider the endomorphism reported below:

$$Lv=\mathsf{\lambda }v$$

(17)

$\mathsf{\lambda }$and $v$denote respectively the generic eigenvalue and eigenvector of the Laplacian. From the eigenvalue and eigenvector decomposition it is possible to infer some interesting properties of the graph such as connectivity and its topological complexity [28]. The eigenvectors correspond to the basis for signals associated with the nodes of the graph; they indicate the modes of variation supported by the topology. The eigenvalues instead encode the spatial frequencies associated with the eigenvectors. Eigenvectors with strong nodal variability will be associated with high eigenvalues and vice versa. The first eigenvalue is always equal to 0 and corresponds to the constant eigenvector and is always supported by all nodes [36]. From the Delaunay mesh, we computed the graph Laplacian and its eigenvalue–eigenvector decomposition. We used the extreme position statistical metrics on the distribution of eigenvalues (excluding the null eigenvalue) to estimate the topological complexity of the graph, while the dispersion and orientation metrics provided estimates of the average density of connections and possible inhomogeneities.

2.5.3. Features Derived from Transforms

From the mesh representation, it is possible to define a nodal signal associated with the graph. From a radiomics perspective, since the point cloud underlying the graph is a subset of the voxels composing the segmented volume, we assigned each point the gray level (radiotracer activity) it had in the original image. Within the graph analysis, we then applied transformations that allowed the extraction of information regarding the radiotracer activity while accounting for the topological structure on which the signal was defined.

Features Derived from The Graph Fourier Transform

The Graph Fourier Transform (GFT) [37] is the graph equivalent of the classical Fourier transform. The bases, in this case, are not sinusoids but the eigenvectors of the Laplacian. Since “transforming” a signal means projecting it onto a new basis where each new element is calculated by performing a similarity measure (correlation) between the signal and the generic basis element, the GFT assumes the following synthesis formulation:

$$G={\mathsf{\Phi }}^{\mathrm{\top }}f$$

(18)

where $G$ denotes the transform, $f$is the node signal and $\mathsf{\Phi }$is the matrix of Laplacian eigenvectors. The GFT is characterized by both frequency and spatial localization. The coefficients establish the similarity with the modes of variation supported by the graph. Now consider the energy spectral density [37] obtained from $G.$ The i-th coefficient can be computed according to the following formulation:

$$\mathrm{ES}{\mathrm{D}}_{\mathrm{i}}={\left|G_i\right|}^2,i=1,\dots ,N$$

(19)

where $G_i$represents the i-th spectral coefficient. Figure 2 shows an example of signal on graph extracted from a lesion in the dataset (a) and the corresponding energy spectral density (b). For greater graphical clarity, the connections between the nodes have not been shown.

From this representation, we extracted features that quantified the total energy and its distribution. The energy was normalized by the total number of nodes to compare the total metabolic activity across graphs of different sizes:

$$E={\displaystyle \frac{1}{N}}\sum _{i=1}^N\mathrm{ES}{\mathrm{D}}_{\mathrm{i}}$$

(20)

For the computation of the metrics, it was preferable to filter out the constant component (first eigenvalue). We separated low frequencies from high frequencies [38] by defining a threshold value ${\lambda }_{\mathrm{t}\mathrm{h}}$:

$$\mathcal{L}=\{\hspace{0.17em}i\mid {\mathsf{\lambda }}_i\le {\mathsf{\lambda }}_{\mathrm{th}}\hspace{0.17em}\}$$

(21)

$$\mathcal{H}=\{\hspace{0.17em}i\mid {\lambda }_i>{\lambda }_{\mathrm{th}}\hspace{0.17em}\}$$

(22)

We evaluated the energies associated with the two frequency bands by summing the spectral coefficients corresponding to the bands of interest:

$$E_{\mathrm{LF}}={\displaystyle \frac{1}{N}}\sum ^{i\in \mathcal{L}}\mathrm{ES}{\mathrm{D}}_{\mathrm{i}}$$

(23)

$$E_{\mathrm{HF}}={\displaystyle \frac{1}{N}}\sum ^{i\in \mathcal{H}}\mathrm{ES}{\mathrm{D}}_{\mathrm{i}}$$

(24)

By normalizing these two quantities by the total energy, we obtained a relative measure of energy partition between the two bands. To separate the frequency range of the energy density spectrum derived from the GFT, we used a threshold equal to 30% of the maximum eigenvalue. In order to show the rationale behind this operational approach, three lesions from different patients in the dataset were considered. For each lesion, the high-frequency energy component was calculated, selecting the spectral threshold at 30% and 50% of the maximum eigenvalue, respectively. The results of the procedure are reported in Figure 3.

From the figure, in particular from lesion 2, it can be seen that the threshold greatly influences the results. In this scenario, the choice of a threshold equal to the 50% of the maximum eigenvalue tends to concentrate most of the energy in the low-frequency band, reducing the contribution of high-frequency components and limiting the sensitivity to localized variations, which are often indicators of potential metabolic peaks in the signal. A threshold equal to 30% allows obtaining a balanced spectral decomposition, where both global trends and local fluctuations of the signal are meaningfully represented.

Features Derived from The Graph Wavelet Transform

The second type of transform we considered was the Graph Wavelet Transform (GWT) [39]. This approach allows investigating how the nodal signal varies with respect to its context across different spatial scales. For the analysis, we defined weighting factors for the graph before applying the decomposition. We built weighting factors that reflected the Euclidean distances between connected nodes, modulating them through a Gaussian kernel [39]:

$$W_{ij}=\{\begin{array}{c}\exp \left(-{\displaystyle \frac{|x_i-x_j{|}_2^2}{2{\mathsf{\sigma }}^2}}\right)\\ 0otherwise\end{array}$$

(25)

where $W_{ij}$denotes the weight assigned to the edge between node $i$and node $j$. If the edge does not exist, no weight is assigned. Then we introduced the matrix $D$, a diagonal matrix where the generic diagonal element is the sum of all weights associated with the corresponding node and we used it to build the normalized Laplacian [36,39]:

$$L=I-D^{-1/2}W{D}^{-1/2}$$

(26)

where $I$is the identity matrix and $W$is the weight matrix. We used the normalized Laplacian for the Graph Wavelet Transform to enable a more interpretable multi-scale analysis, as normalization mitigates the influence of varying node degrees and weighted edges, producing spectral components that are stable and comparable across nodes and scales. As in the classical Wavelet transform, it was necessary to identify a Wavelet kernel to use for the decomposition. Following Hammond et al., we adopted one of the proposed spectral kernels, defined as:

$$g\left(s\mathsf{\lambda }\right)=s\mathsf{\lambda }\hspace{0.17em}{e}^{-s\mathsf{\lambda }}$$

(27)

where $\lambda$is an eigenvalue of the normalized Laplacian and $s$represents a specific scale. It is a band-pass filter that emphasizes specific frequency bands. Small scales correspond to broader bands and a theoretical peak shifted toward higher frequencies, allowing the capture of local variations that manifest over a relatively extended frequency support. Large scales correspond to narrower bands and a theoretical peak shifted toward lower frequencies, enabling the filter to selectively highlight patterns that manifest at a macroscale. To apply the wavelet, we selected three spatial scales. Figure 4 shows the shape of the kernel for the chosen scales.

To apply the Wavelet transform, the node signal must be projected onto the eigenvectors of the normalized Laplacian [39]. The GFT is subsequently weighted by the kernel. At last, the result is reprojected onto the nodal space through the eigenvector matrix:

$${\mathsf{\psi }}_s=U\hspace{0.17em}g\left(s\mathsf{\lambda }\right)\hspace{0.17em}{U}^{\mathrm{\top }}f$$

(28)

where ${\psi }_s$denotes the Wavelet transform, $U$the eigenvector matrix, and $f$the nodal signal. Positive nodal coefficients reflect an increase of the nodal signal with respect to the local context determined by the scale, whereas negative nodal coefficients reflect a decrease. For each scale, the decomposition was performed, resulting in three distinct distributions representing the radiotracer activity on graph at three different spatial resolutions. From the Wavelet distributions, we used the statistical metrics listed in section 2.5.1 to extract information on trend and dispersion patterns at different scales. Figure 5 shows an example of Wavelet decomposition applied to the signal in Figure 2a.

2.6. Validation Analysis

The procedures described above were repeated for all lesions, and at the end of the process, the resulting data were used to assess the reliability of the proposed methodology through informative, robustness and predictive validation.

2.7. Informative Validation

To evaluate the added informational value of the new features, a correlation analysis was conducted using Spearman’s correlation. From the correlation matrix, the distribution of the average correlation with the classical features was subsequently derived for each proposed feature, and a left-tailed Wilcoxon signed-rank test [40] was performed to test the median of the distribution against different reference values.

2.8. Robustness Validation

The dependence on segmentation was evaluated by recomputing the classical and GSP metrics extracted from a subset of three lesions after perturbing the original mask used for target extraction. The percentage differences between the classical and GSP feature vectors were analyzed using a descriptive statistical approach to assess robustness to the variability introduced by mask perturbation. The dependence of classical and GSP features on the batch effect was evaluated at the patient level. For each feature, the distribution of its observations was considered. This distribution was split into two groups based on scanner membership, and the Kruskal-Wallis test was performed to assess the presence of batch effect.

2.9. Predictive Validation

To assess the discriminative power of the proposed features, we implemented two different Machine Learning methods: Linear Discriminant Analysis (LDA) [41] and Support Vector Machine (SVM) [42]. For each algorithm, three models were constructed and trained on the proposed features, the classical features, and a combination of both, respectively, and their performances were compared in terms of accuracy and Receiver Operating Characteristic (ROC) curves.

2.9.1. Performance Metrics

The predictive performance of a model can be quantified using several metrics. Considering a positive-negative dichotomy, the most relevant metrics are:

• True Positive (TP): Total positive observations correctly classified. Normalized by the total number of positives, this gives the True Positive Rate (TPR), also called sensitivity.
• True Negative (TN): Total negative observations correctly classified. Normalized by the total number of negatives, this gives the True Negative Rate (TNR), also called specificity.
• False Positive (FP): Total negative observations incorrectly classified as positive. Normalized by the total number of negatives, this gives the False Positive Rate (FPR), also called fall-out.
• False Negative (FN): Total positive observations incorrectly classified as negative. Normalized by the total number of positives, this gives the False Negative Rate (FNR).

These metrics can be combined to estimate the accuracy, which represents the overall discriminative ability of the model. It is defined as the ratio of correctly classified observations to the total number of observations:

$$\mathrm{Accuracy}={\displaystyle \frac{TP+TN}{TP+TN+FP+FN}}$$

(29)

Class assignment is determined by comparing the score to a threshold:

$$\{\begin{array}{c}{z}_i\ge \theta \Rightarrow \widehat{c_i}=0\\ z_i<\theta \Rightarrow \widehat{c_i}=1\end{array}$$

(30)

where $z_i$denotes the score corresponding to the $i$-th observation, and $\widehat{c_i}$denotes the predicted class. For an established threshold, the classifier will be characterized by a specific TPR and a specific FPR; this pair represents a point in the FPR-TPR plane that describes the classifier. By varying the threshold, a new model will be obtained, characterized by different TPR and FPR corresponding to a new point in the plane. The set of points obtained by varying the threshold is called the ROC curve and represents a fundamental element for evaluating the classifier’s performance. The area under curve (AUC) indicates the capability of the classifier to separate the two classes.

2.9.2. Feature Selection

To reduce dimensionality and avoid overfitting, we implemented the LASSO algorithm for feature selection [43]. Starting from an observation $y_i$characterized by a set of features and a class to be predicted, linear regression is based on the hypothesis that the observation can be obtained through a linear combination of the features plus an error term ${\epsilon }_i$.

Let $x_i^{\mathrm{\top }}=\left(1,x_{i1},x_{i2},\dots ,x_{ip}\right)$ be the row vector representing the features associated with the observation, and ${\mathsf{\beta }}^{\mathrm{\top }}=\left({\mathsf{\beta }}_0,{\mathsf{\beta }}_1,{\mathsf{\beta }}_2,\dots ,{\mathsf{\beta }}_p\right)$ the vector of coefficients weighting the features. Linear regression is defined as:

$$y_i=\sum _{j=0}^p{\mathsf{\beta }}_j{x}_{ij}+{\mathsf{\epsilon }}_i$$

(31)

where ${\epsilon }_i$represents the unknown error term. The prediction obtainable using a generic vector ${\beta }^{\mathrm{\top }}$is:

$$\widehat{y_i}=\sum _{j=0}^p{\beta }_j{x}_{ij}$$

(32)

The prediction error is defined as the difference between the observed and predicted value. Linear regression aims to find the optimal vector $\widehat{\beta }$that minimizes the error over the $n$observations of the dataset:

$$\widehat{\mathsf{\beta }}=\arg \underset{\mathsf{\beta }}{\min }\sum _{i=1}^n{\left(y_i-\sum _{j=0}^p{\mathsf{\beta }}_j{x}_{ij}\right)}^2$$

(33)

To perform feature selection, criteria, often statistical or based on the elements of $\beta$, are used to eliminate less informative features based on their contribution to the model in predictive terms. LASSO allows feature selection by introducing a penalization term into the minimization problem:

$$\widehat{\mathsf{\beta }}=\arg \underset{\mathsf{\beta }}{\min }\sum _{i=1}^n{\left(y_i-\sum _{j=0}^p{\mathsf{\beta }}_j{x}_{ij}\right)}^2+\mathsf{\lambda }\sum _{j=0}^p\left|{\mathsf{\beta }}_j\right|$$

(34)

where $\lambda$is the penalty parameter. If many coefficients have large values, the penalty term has a stronger impact. This procedure can shrink several coefficients to zero, facilitating the identification of features to remove. For binary classification tasks, logistic regression is used instead. In the linear approach, considering data generation, $y_i$can be treated as a realization of a random variable $Y_i$, given a fixed error value. It can be shown that:

$$E\left[Y_i\right|x_i]={\mathsf{\beta }}_0+{\mathsf{\beta }}_1{x}_{i1}+{\mathsf{\beta }}_2{x}_{i2}+\dots +{\mathsf{\beta }}_p{x}_{ip}$$

(35)

This expression identifies the deterministic part of the model. For dichotomous outcomes [37], this expected value corresponds to the conditional probability that $Y_i$ belongs to class 1:

$$E\left[Y_i\right|x_i]=P\left(Y_i=1\right|x_i)={\displaystyle \frac{1}{1+e^{-\left({\mathsf{\beta }}_0+{\mathsf{\beta }}_1{x}_{i1}+{\mathsf{\beta }}_2{x}_{i2}+\dots +{\mathsf{\beta }}_p{x}_{ip}\right)}}}$$

(36)

Let $P(Y_i\mid x_i)$denote the probabilities that the model explains the true value of the observation; the final model is constructed from the product of probabilities for each observation:

$$L\left(\mathsf{\beta }\right)=\prod _{i=1}^{n}P\left(Y_i\right|x_i)=\prod _{i=1}^n{p}_i^{Y_i}{\left(1-p_i\right)}^{1-Y_i},\hspace{1em}{p}_i=P\left(Y_i=1\right|x_i)$$

(37)

Taking the logarithm yields:

$$\mathrm{l}\left(\mathsf{\beta }\right)=\log L\left(\mathsf{\beta }\right)=\sum _{i=1}^n\left(Y_i\log p_i+\left(1-Y_i\right)\log \left(1-p_i\right)\right)$$

(38)

LASSO introduces a penalty term to this function as well:

$$\widehat{\mathsf{\beta }}=\arg \underset{\mathsf{\beta }}{\min }-\sum _{i=1}^n\left(Y_i\log p_i+\left(1-Y_i\right)\log \left(1-p_i\right)\right)+\mathsf{\lambda }\sum _{j=0}^p\left|{\mathsf{\beta }}_j\right|$$

(39)

2.9.3. LDA Model Development

Let us consider the vector $x_i=\left[\hspace{0.17em}{x}_{i1},\hspace{0.25em}{x}_{i2},\hspace{0.25em}\dots ,\hspace{0.25em}{x}_{id}\hspace{0.17em}\right]$, identifying the features associated with the i-th observation, and the vector $w=\left[\hspace{0.17em}{w}_1,\hspace{0.25em}{w}_2,\hspace{0.25em}\dots ,\hspace{0.25em}{w}_d\hspace{0.17em}\right]$ representing a direction in the multidimensional space $R^d$. Linear Discriminant Analysis consists of projecting the vector $x_i$ onto the direction identified by $w$, obtaining the score $z_i$ representative of the observation. To ensure that the classification is reliable, it is appropriate to identify the optimal vector $w$ that guarantees maximum separability of the classes. The optimal vector corresponds to the direction that guarantees maximum separability from the centers of the distributions and minimum intra-class variance.

2.9.4. SVM Model Development

The Support Vector Machine algorithm identifies an optimal separating hyperplane that maximizes the minimum distance between the decision boundary and the training samples of each class. The separating hyperplane is defined as:

$$w^{T}x+b=0$$

(40)

where $w$denotes the normal vector to the hyperplane and $b$the bias term. The optimal hyperplane is obtained by solving the following convex optimization problem:

$$\underset{w,b}{\min }{\displaystyle \frac{1}{2}}{\Vert w\Vert }^2$$

(41)

subject to:

$$y_i\left(w^T{x}_i+b\right)\ge 1\forall i$$

(42)

where $y_i\in \{-1,+1\}$are the class labels. Once the optimal hyperplane is estimated, a new sample $x$is classified according to:

$$f\left(x\right)=\mathrm{sign}\left(w^{T}x+b\right)$$

(43)

3. Results

The results are presented according to the three main objectives: informative validation of the proposed features, robustness validation and evaluation of predictive modeling.

3.1. Correlation Results

The correlation Analysis was performed setting a threshold of 0.6 to define a high correlation. Figure 6 reports the correlation matrix (a), where each element above the threshold is highlighted with a rectangular red marker, and the average correlation of each GSP feature with all classic features (b).

According to a Spearman correlation threshold of 0.6, 28.6% of GSP features were not correlated with any classical feature, indicating that these features provide additional information beyond the classical radiomics metrics. The median of the distribution in Figure 6b was compared with the reference values reported in Table 1 and the corresponding p-values were calculated.

3.2. Robustness Results

3.2.1. Dependence on Segmentation

For each of the three lesions selected for the analysis, the mean, median, and standard deviation of the percentage difference vectors of the features were computed. Table 2 reports the results; the lesions are labeled (a), (b), and (c) for the first, second, and third lesion, respectively.

3.2.2. Dependence on Batch Effect

Given the high number of comparisons to assess the presence of a batch effect and in order to control for false positives, the p-value of each comparison was corrected (q-value) using the Benjamini–Hochberg False Discovery Rate (FDR) correction [44]. In Table 3, the percentages of GSP and classical features for which the comparison yielded a q-value lower than the reference value $\alpha =0.05$are reported, respectively.

3.3. Predictive Modelling Results

For each model, feature selection was performed using LASSO and, in order to avoid data leakage, it was performed internally within the k-fold cross-validation, where observations from the same patient were assigned to the same train/test fold. For each model development, the classifier was implemented and trained using stratified 5-fold cross-validation. Since the results, in terms of accuracy and ROC curve, depended on how the data were split into train/test sets, the training procedure was repeated 10 times for each model.

3.3.1. LDA Results

This section reports the results obtained by training the LDA classifier. Figure 7 shows the mean ROC curves (a), obtained by averaging the curves across cross-validation runs, and the distributions of accuracies for the three models (b).

Table 4 reports the mean accuracies and mean AUC values for the three models. The metrics were accompanied by the 95% confidence interval.

3.3.2. SVM Results

This section reports the results obtained by training the SVM classifier. Figure 8 shows the mean ROC curves (a), obtained by averaging the curves across cross-validation runs, and the distributions of accuracies for the three models (b).

Table 5 reports the mean accuracies and mean AUC values for the three models. The metrics were accompanied by the 95% confidence interval.

4. Discussion

This study introduces an exploratory approach aimed at validating an innovative methodology in the context of radiomics analysis. The results of the correlation study suggest that the proposed graph-based analysis captures discriminative patterns not entirely encoded by standard voxel-based descriptors, while the model performance highlights an apparent predictive contribution of the proposed features. In particular, the models trained exclusively on GSP features exhibited higher accuracy and AUC compared to the other two models, suggesting that the information captured through the new representation may enable the construction of models with improved class-separation capacity compared to those obtained using classical approaches, showing no evident advantage from combining it with the standard informational domain. We have verified that the GSP metrics exhibit greater robustness to variability introduced by segmentation and to batch effects; therefore, the results obtained may reflect these improvements. However, we fully acknowledge that the study is based on a relatively small cohort (which represents an inherent limitation, particularly for machine learning–based approaches). Furthermore, the evaluation does not include direct comparisons with more recent or advanced approaches, such as deep learning–based methods or modern radiomics frameworks. This choice was primarily driven by the limited size of the available dataset, which we believe would not support a fair, stable, and meaningful training and evaluation of more complex models without a high risk of overfitting. For these reasons, we emphasize that the results should be interpreted as preliminary and hypothesis-generating rather than definitive. Future studies will investigate the discriminative capabilities of the models on larger datasets, incorporating complementary techniques to classical GSP, such as graph neural networks (GNN) [45,46] and topological signal processing (TSP) [26,47]. These approaches will be introduced to provide advanced additional information beyond that captured by the proposed features, with the goal of maximizing their discriminative power within this translational framework.

5. Conclusions

The study introduces a translational approach that enables the implementation of innovative data processing and analysis techniques in radiomics studies, beyond the standards proposed in the literature. The introduced metrics aim to describe the geometry, structure, and characteristics of the metabolic activity of the lesion; these aspects may be particularly relevant for distinguishing benign from malignant lesions, as they reflect differences in structural organization and metabolic heterogeneity. The interpretability of the metrics is therefore closely linked to the surface, structural, and metabolic aspects of the lesion. The type of information described is more detailed and sector-specific from a qualitative perspective compared to the information associated with classical predictors, which describe much more general characteristics of the examined targets. The results indicate the ability to capture morpho-functional characteristics of the targets that are not directly accessible through classical approaches, providing a benefit to the predictive capabilities of the models. Future studies will focus on acquiring larger experimental cohorts and implementing highly specific processing strategies, with the aim of validating the clinical reliability of the proposed approach.