Intestinal tumor precise typing and prognosis prediction system based on multi-modal medical images

Abstract

This invention relates to the field of medical image processing technology, and in particular to a system for accurate subtyping and prognostic prediction of gastrointestinal tumors based on multimodal medical images. The system includes a data preprocessing module, a differential geometric supervoxel feature extraction module, a multimodal data fusion module, a subtyping model construction module, and a prognostic prediction and evaluation module. The differential geometric supervoxel feature extraction module consists of a curvature flow-driven adaptive supervoxel segmentation unit, a multi-scale geometric feature extraction unit on a Riemannian manifold, and a geodesic network tumor heterogeneity characterization unit. A global heterogeneity index is used to adaptively modulate the two-modal fusion weights, and a prognostic consistency index is used to update the weight and feature combination. The three modules are coupled in a closed loop. On a test set of 486 retrospective cohort cases, the subtyping accuracy of this invention was 89.2%, the area under the curve was 0.921, and the prognostic consistency index was 0.784.

Description

Technical Field

[0001] This invention relates to the field of medical image processing technology, and in particular to a system for accurate classification and prognosis prediction of gastrointestinal tumors based on multimodal medical images. This system uses differential geometry theory to extract features from medical images to achieve accurate molecular classification and prognosis prediction of gastrointestinal tumors. Background Technology

[0002] Gastrointestinal tumors are among the most common malignant tumors of the digestive system, and accurate identification of their molecular subtypes and prognostic assessment are crucial for developing individualized treatment plans. Traditional tumor classification mainly relies on invasive biopsies and pathological examinations, which not only cause pain for patients but also suffer from sampling errors and operator dependence. Furthermore, traditional prognostic assessments, primarily based on clinical staging and pathological indicators, struggle to accurately reflect the biological behavior of tumors and their response to individualized treatment.

[0003] In recent years, the development of radiomics has provided new insights into precise tumor subtyping and prognostic prediction. By extracting quantitative features from medical images and establishing correlation models between imaging features and molecular phenotypes and clinical prognoses, it is hoped that non-invasive tumor molecular subtyping and personalized prognostic assessment can be achieved. However, existing radiomics methods are mainly based on conventional statistical features and texture analysis, which are insufficient to fully capture the geometric structural features and heterogeneity information of tumor tissues. In addition, the effective integration and feature fusion of multimodal medical image data are also important challenges currently faced.

[0004] Furthermore, existing radiomics-based subtyping and prognostic prediction schemes generally adopt a unidirectional pipeline structure in their system architecture, consisting of four stages: feature extraction, modality fusion, model building, and prognostic assessment. Each stage only transmits feature vectors and prediction results downstream, without transmitting any information about its own reliability upstream; there is no information feedback between modules. This structure determines that the overall system performance is approximately equal to the linear sum of the individual gains of each stage: replacing the feature extraction operator with a stronger one only brings an increase in the feature level, and adding a modality only brings an increase in the amount of information; using both simultaneously does not produce additional gains. On the other hand, in existing multimodal fusion, the weights of image modalities and gene expression modalities are usually determined once during the training phase by cross-validation and remain fixed for all cases during the inference phase. This does not take into account that gene expression data is taken from tissue homogenates or local punctures, which is essentially a spatial average of the entire tumor. Its representativeness of the overall tumor decreases as the spatial heterogeneity of the tumor increases. For highly heterogeneous cases, fixed weights will systematically overestimate the reliability of gene expression modalities, and this bias cannot be eliminated by increasing the sample size or changing the classifier.

[0005] Therefore, there is an urgent need to develop a new system that can deeply mine the geometric features of multimodal medical images in order to achieve accurate classification and prognostic prediction of gastrointestinal tumors. Summary of the Invention

[0006] The purpose of this invention is to provide a precise classification and prognostic prediction system for gastrointestinal tumors based on multimodal medical images. This system extracts features from medical images by introducing differential geometry theory to achieve precise molecular classification and prognostic prediction of gastrointestinal tumors.

[0007] This invention proposes a precise classification and prognostic prediction system for gastrointestinal tumors based on multimodal medical images, comprising:

[0008] The data preprocessing module is used to acquire digestive endoscopy image datasets, CT image datasets, and gene expression datasets, and to standardize the acquired data.

[0009] A differential geometry supervoxel feature extraction module, connected to the data preprocessing module, is used to receive the standardized image data and extract radiomics features based on differential geometry theory. The differential geometry supervoxel feature extraction module includes a curvature flow-driven adaptive supervoxel segmentation unit, a multi-scale geometric feature extraction unit on a Riemannian manifold, and a geodesic network tumor heterogeneity characterization unit.

[0010] The multimodal data fusion module, connected to the differential geometry supervoxel feature extraction module, is used to fuse the radiomics features with the gene expression dataset to generate fused feature data, and adaptively modulate the fusion weights of the image modality and the gene expression modality using the global heterogeneity index output by the geodesic network tumor heterogeneity characterization unit.

[0011] The subtyping model construction module is connected to the multimodal data fusion module and is used to train a tumor molecular subtyping model based on the fused feature data;

[0012] The prognosis prediction and assessment module, connected to the subtyping model construction module, is used to predict patient prognosis and generate a prognosis assessment report based on the tumor molecular subtyping model and the fusion feature data. The prognosis prediction and assessment module uses the molecular subtype posterior probability output by the tumor molecular subtyping model as a covariate for survival analysis and feeds back the consistency index output by the prognosis prediction and assessment module to the multimodal data fusion module to update the fusion weights and feature combinations, so that the differential geometric hypervoxel feature extraction module, the multimodal data fusion module, and the prognosis prediction and assessment module form a closed loop.

[0013] Preferably, the curvature flow-driven adaptive supervoxel segmentation unit includes:

[0014] Image manifold mapping unit, used to map medical image data to a high-dimensional Riemannian manifold;

[0015] The curvature field calculation unit, connected to the image manifold mapping unit, is used to calculate the Gaussian curvature and average curvature of each point on the high-dimensional Riemannian manifold and generate a curvature field mapping.

[0016] A seed point optimization placement unit is connected to the curvature field calculation unit and is used to identify key regions and optimize the distribution of seed points based on the curvature field mapping.

[0017] The curvature flow evolution control unit, connected to the seed point optimization placement unit, is used to design curvature-dependent velocity field control of supervoxel boundary evolution and generate an initial supervoxel set.

[0018] The multi-phase supervoxel fusion unit, connected to the curvature flow evolution control unit, is used to evaluate, screen, and fuse supervoxel sets under different parameter settings, and output the optimized supervoxel partitioning results.

[0019] Preferably, the multi-scale geometric feature extraction unit on the Riemannian manifold includes:

[0020] The differential geometric invariant computation unit is used to calculate a set of differential geometric features for each hypervoxel, including curvature statistics, shape exponent, and tensor field features;

[0021] The multi-scale thermonuclear analysis unit, connected to the differential geometric invariant calculation unit, is used to construct the Laplacian matrix of the supervoxel diagram and calculate the thermonuclear signature at multiple time scales to generate multi-scale thermonuclear signature features.

[0022] The spectral decomposition and feature enhancement unit, connected to the multi-scale thermonuclear analysis unit, is used to perform spectral decomposition on the supervoxel Laplacian matrix, design frequency domain filters to enhance key features, and output the enhanced feature set.

[0023] Preferably, the geodesic network tumor heterogeneity characterization unit includes:

[0024] The geodesic network construction unit is used to identify physically adjacent hypervoxel pairs and calculate the geodesic distance between hypervoxels to construct a weighted geodesic map.

[0025] The heterogeneity multidimensional characterization unit is connected to the geodesic network construction unit and is used to calculate the feature variability and overall network topology characteristics in a local area based on the weighted geodesic map, and generate local and global heterogeneity indices.

[0026] The dynamic spatiotemporal heterogeneity monitoring unit, connected to the heterogeneity multidimensional characterization unit, is used to analyze the rate of change of heterogeneity indicators over time, establish a spatiotemporal model for predicting tumor development, and generate a spatiotemporal heterogeneity prediction report.

[0027] Preferably, the multimodal data fusion module includes:

[0028] The unified feature space mapping unit is used to normalize and standardize the features of each modality, calculate the correlation between cross-modal features, and construct a unified feature space that integrates multimodal information.

[0029] The multimodal feature optimization unit, connected to the feature space unified mapping unit, is used to evaluate the contribution of each feature to the target task, identify and remove highly redundant features, and optimize feature combinations for classification and prognosis prediction respectively.

[0030] Preferably, the fractal model construction module includes:

[0031] The training data preparation unit is used to integrate image features and molecular typing labels, and to divide the training set, validation set, and test set.

[0032] The hierarchical model building unit, connected to the training data preparation unit, is used to build a two-layer classification model. The first layer distinguishes the main molecular subtypes, and the second layer refines the subtype classification.

[0033] The model optimization and verification unit, connected to the hierarchical model construction unit, is used to optimize model parameters and verify performance, and output the optimized fractal model.

[0034] Preferably, the prognosis prediction and assessment module includes:

[0035] The survival analysis model unit is used to integrate image features and survival data to build a survival prediction model.

[0036] The treatment response prediction unit, connected to the survival analysis model unit, is used to analyze the correlation between imaging features and treatment response, and to predict chemotherapy response rate and radiotherapy sensitivity.

[0037] The integrated risk assessment unit, connected to the survival analysis model unit and the treatment response prediction unit, is used to integrate survival and treatment response prediction results, classify patients into high, medium and low risk groups, and generate individualized risk assessment reports.

[0038] Preferably, the system further includes:

[0039] The system performance evaluation module, connected to the prognostic prediction and evaluation module, is used to evaluate system performance on test sets and independent external datasets, evaluate the stability of the system under different conditions, and generate a performance evaluation report.

[0040] The clinical decision support module, connected to the typing model construction module and the prognosis prediction and assessment module, is used to transform the analysis results into clinically friendly reports, integrate with the hospital PACS and electronic medical record system, collect clinical usage feedback, and continuously improve the system.

[0041] Preferably, the curvature flow-driven adaptive supervoxel segmentation unit further includes:

[0042] The supervoxel quality assurance unit, connected to the multi-phase supervoxel fusion unit, is used to calculate the consistency score between the supervoxel boundary and the image gradient, evaluate the conformity of the supervoxel with the known tissue structure, and perform adaptive re-segmentation on supervoxels that do not meet the quality standards.

[0043] Preferably, the differential geometric invariant calculation unit has the following features:

[0044] Local curvature characteristics, including mean curvature, Gaussian curvature, mean, variance, skewness, and kurtosis of principal curvatures;

[0045] Shape index features, including shape index and curvature index, discretize continuous shape features into typical categories;

[0046] Tensor field characteristics, including anisotropy, linearity, flatness, and sphericity, are used to characterize the internal structural properties of supervoxels.

[0047] Preferably, the feature space unified mapping unit is also connected to the geodesic network tumor heterogeneity characterization unit, for receiving the global heterogeneity index, and generating image modality fusion weights by monotonically increasing bounded mapping based on the difference between the global heterogeneity index and the preset heterogeneity benchmark value, and using the difference between 1 and the image modality fusion weights as gene expression modality fusion weights, and increasing the image modality fusion weights when the global heterogeneity index is higher than the heterogeneity benchmark value.

[0048] Preferably, the multimodal feature selection unit is further configured to receive the consistency index, and when the consistency index is lower than the value of the previous iteration, adjust the modulation intensity coefficient of the heterogeneity benchmark value and the bounded mapping, regenerate the image modality fusion weight, the gene expression modality fusion weight and the corresponding feature combination, and terminate the iteration and lock the fusion weight and the feature combination when the consistency index increment between two adjacent iterations is less than a preset convergence threshold or the number of iterations reaches a preset upper limit.

[0049] Preferably, the prognostic prediction and assessment module is also used to calculate the rate of change of the global heterogeneity index between two adjacent examinations and the rate of change of the risk score corresponding to the patient's prognosis within the same time interval, perform a sign consistency check on the two rates of change, and reduce the confidence of the prognostic assessment report and output a review indicator when the signs of the two rates of change are opposite.

[0050] The beneficial effects of this invention include:

[0051] 1. By introducing differential geometry theory to construct a curvature flow-driven adaptive supervoxel segmentation framework, it is possible to accurately capture tumor boundaries and internal heterogeneous regions. Using conventional gray-level co-occurrence matrix and wavelet texture feature schemes as baselines, introducing this segmentation framework alone can improve the test set classification accuracy by 6.5 percentage points and the prognostic consistency index by 0.047.

[0052] 2. Based on the multi-scale geometric feature extraction mechanism of Riemannian manifolds, a unified feature expression is established between the local differential geometric quantities of a single supervoxel and the topological quantities of the geodesic network at the whole tumor scale, which is used to characterize the geometric structural characteristics of the tumor.

[0053] 3. Using a geodesic network tumor heterogeneity characterization system, geodesic networks are used to analyze the connectivity between supravoxels within the tumor. The output global heterogeneity index depends only on the intrinsic geometry of the image manifold and not on the absolute value of gray levels, and can be used as a cross-device comparable measure of heterogeneity.

[0054] 4. Multimodal data fusion technology effectively integrates image features and gene expression information, overcoming the limitations of single-modal data. Using conventional gray-level co-occurrence matrix and wavelet texture feature schemes as baselines, introducing this fusion method alone can improve the test set typing accuracy by 4.7 percentage points and the prognostic consistency index by 0.054.

[0055] 5. The features extracted by the system are quantities with clear geometric definitions, such as mean curvature, Gaussian curvature, shape index, tensor field anisotropy, geodesic distance and network clustering coefficient. Their values ​​can be traced back to the supervoxel partitioning results and weighted geodesic maps, and the formation process of classification and prognostic conclusions can be tracked.

[0056] 6. This invention does not simply superimpose the three stages of differential geometric feature extraction, multimodal fusion, and prognostic prediction, but rather integrates them into a systematic closed loop: the global heterogeneity index output by the differential geometric supervoxel feature extraction module not only participates in modeling as a feature but also serves as a quantitative measure of gene expression modality reliability, used to adaptively modulate the fusion weights of image modality and gene expression modality; the molecular subtype posterior probability output by the genotyping model serves as a covariate for survival analysis; and the consistency index of prognostic prediction is fed back to the multimodal data fusion module to update the fusion weights and feature combinations. The mechanism lies in the fact that the global heterogeneity index given by the geodesic network is an intrinsic quantity defined on the Riemannian manifold, covering the entire tumor range and independent of the grayscale scale. It precisely measures the representativeness of gene expression data taken from local areas or homogenates to the overall tumor, thus serving as a cross-modal reliability measure. Conventional texture features are affected by window width and level, resampling, and device differences, and their values ​​lack cross-device comparability, making them unable to perform the same function. This causal relationship determines that closed-loop coupling only holds true when differential geometric feature extraction is used, and the three are mutually conditional rather than independent of each other.

[0057] 7. Due to the aforementioned closed-loop coupling, the technical effect of this invention significantly exceeds the sum of the individual effects of each module: In the control experiment, the individual introduction of differential geometric feature extraction, multimodal fusion, and closed-loop feedback resulted in classification accuracy gains of 6.5, 4.7, and 1.2 percentage points, respectively, relative to the conventional radiomics baseline. The linear superposition of the three was expected to be 12.4 percentage points, while the measured gain when all three were used simultaneously was 17.4 percentage points, with a synergistic gain of 5.0 percentage points and a nonlinear synergistic multiple of 1.40. The corresponding synergistic multiple for the prognostic consistency index was 1.36. In particular, the gain when closed-loop feedback was introduced alone was only 1.2 percentage points, close to the fluctuation range of repeated experiments, while its marginal contribution jumped to 7.1 percentage points when it coexisted with the other two. The effect of the same technical feature differed by nearly six times under different combination environments. This strong nonlinear relationship cannot be extrapolated from the effect of the feature alone.

[0058] 8. This invention improves upon existing technologies in terms of typing accuracy and prognostic prediction precision: the typing accuracy is increased from 71.8% in conventional radiomics and 80.6% in convolutional neural network methods to 89.2%; the area under the receiver operating characteristic curve (AUC) is increased from 0.742 and 0.836 to 0.921; the prognostic consistency index is increased from 0.641 and 0.673 to 0.784; on independent external datasets from different centers and devices, the typing accuracy is 86.4% and the consistency index is 0.761, with a decrease of 2.8 percentage points and 0.023 percentage points respectively relative to the internal test set, which is significantly lower than the 5.2 percentage points and 9.4 percentage points of the comparative scheme. Attached Figure Description

[0059] Figure 1 This is a schematic diagram of the overall architecture of the system of the present invention;

[0060] Figure 2 This is a schematic diagram of the differential geometry supervoxel feature extraction module of the present invention;

[0061] Figure 3 This is a flowchart illustrating the workflow of the curvature flow-driven adaptive hypervoxel segmentation unit of the present invention.

[0062] Figure 4 This is a flowchart illustrating the workflow of the multi-scale geometric feature extraction unit on the Riemannian manifold of the present invention.

[0063] Figure 5 This is a flowchart of the process of the geodesic network tumor heterogeneity characterization unit of the present invention.

[0064] Figure 6 This is a flowchart illustrating the workflow of the multimodal data fusion module of the present invention.

[0065] Figure 7 This is a flowchart illustrating the workflow of the fractal model construction module of the present invention.

[0066] Figure 8 This is a flowchart of the prognosis prediction and evaluation module of the present invention. Detailed Implementation

[0067] Please refer to Figures 1-8 The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments.

[0068] Reference Figure 1 The present invention provides a system for accurate classification and prognosis prediction of gastrointestinal tumors based on multimodal medical images, including a data preprocessing module 1, a differential geometric supervoxel feature extraction module 2, a multimodal data fusion module 3, a classification model construction module 4, and a prognosis prediction and evaluation module 5.

[0069] In one embodiment of the present invention, the data preprocessing module 1 is used to acquire a digestive endoscopy image dataset, a CT image dataset, and a gene expression dataset, and to perform standardization processing on the acquired data. Preferably, the data standardization processing includes steps such as image intensity normalization, spatial resolution unification, and noise suppression. For endoscopy images, histogram equalization can be used to enhance image contrast; for CT images, window width and window level adjustments can be used to map CT values ​​to a range suitable for displaying gastrointestinal tissues, such as setting the window width to 350 HU and the window level to 40 HU. This parameter setting can highlight soft tissue structures and is beneficial for subsequent feature extraction. In addition, for gene expression data, Z-score normalization is used to eliminate scale differences between different gene expression levels.

[0070] The differential geometry hypervoxel feature extraction module 2 is connected to the data preprocessing module 1. It receives standardized image data and extracts radiomics features based on differential geometry theory. This module is the core innovation of this invention; its design, based on differential geometry theory, enables in-depth mining of the geometric structural features of medical images. For example... Figure 2 As shown, the differential geometry supervoxel feature extraction module 2 includes a curvature flow-driven adaptive supervoxel segmentation unit 21, a multi-scale geometric feature extraction unit 22 on a Riemannian manifold, and a geodesic network tumor heterogeneity characterization unit 23.

[0071] The curvature flow-driven adaptive supervoxel segmentation unit 21 is used to segment medical images into biologically meaningful supervoxels. For example... Figure 3 As shown, the unit includes an image manifold mapping unit 211, a curvature field calculation unit 212, a seed point optimization placement unit 213, a curvature flow evolution control unit 214, and a multi-phase supervoxel fusion unit 215. The image manifold mapping unit 211 first performs multi-scale enhancement on the normalized image to strengthen tissue boundary information. Preferably, for endoscopic images, the enhancement coefficient can be set between 0.6 and 0.8; for CT images, the enhancement coefficient can be set between 0.4 and 0.6. This parameter setting is based on numerous clinical image experiments and can effectively suppress noise while preserving details.

[0072] Subsequently, the image manifold mapping unit 211 constructs a high-dimensional feature space and computes local feature descriptors, including intensity, texture, and gradient information. The feature space dimension is adaptively adjusted according to the image complexity, typically 8-12 dimensions. Finally, this unit defines a Riemannian metric tensor to establish the image manifold. The Riemannian metric tensor can be represented as:

[0073] ,

[0074] in, For point Riemannian metric tensor elements at the location. For the Kronecker function (when The value is 1 when the time is right, and 0 otherwise. For image At point Along the first Gradient components in each direction, This is a weighting coefficient used to adjust the degree of influence of the gradient on the metric, typically set between 0.1 and 0.5. This metric definition elongates the distance of high-gradient regions (such as tissue boundaries) on the manifold, resulting in more accurate boundary localization in subsequent segmentation.

[0075] The curvature field calculation unit 212 calculates the Gaussian curvature and mean curvature of each point based on the image manifold data, generating a curvature field mapping. and mean curvature The calculation formula is:

[0076] ,

[0077] ,

[0078] in, For Gaussian curvature, For the second basic form of determinant, The determinant is the first basic form. For the mean curvature, For matrix traces, and These are the first and second fundamental forms of the manifold, respectively. First fundamental form It is a quadratic form describing a metric on a manifold, and the second fundamental form This describes the curvature of the manifold in the embedded space. To suppress the influence of noise, the curvature field is smoothed at multiple scales, which can be achieved using a Gaussian kernel function.

[0079] ,

[0080] in, The curvature value after smoothing. The image domain (i.e., the area covered by the entire image). For point The original curvature value at that point, The standard deviation is The Gaussian function is expressed as: , For vectors The Euclidean norm, The value is usually set to 1-3 pixels, a range that can effectively smooth out noise while preserving important structural information.

[0081] The seed point optimization placement unit 213 identifies key regions and optimizes the distribution of seed points based on curvature field mapping. Preferably, this unit uses the rate of curvature change as a key indicator, increasing the seed point density by 50% in regions with a high rate of curvature change (such as tumor boundaries). Increasing the seed point density can improve the segmentation accuracy of these key regions, especially for tumor boundaries and invasive areas with complex shapes.

[0082] The curvature flow evolution control unit 214 is designed to control the evolution of the supervoxel boundary using a curvature-dependent velocity field, generating an initial set of supervoxels. The evolution process can be described as follows:

[0083] ,

[0084] in, Supervoxel boundary curve Over time rate of change, For Gaussian curvature and mean curvature velocity function, Let be the boundary normal vector (a vector of unit length pointing outwards from the hypervoxel). Preferably, the velocity function can be defined as:

[0085] ,

[0086] in, This is a velocity function, representing the basic velocity parameter, typically set to 0.5-1.0. and The curvature influence coefficient controls the degree of influence of Gaussian curvature and mean curvature on the evolution rate, and is usually set in the range of 0.3-0.7. and These are the absolute values ​​of the Gaussian curvature and the mean curvature, respectively. This velocity field design slows down the evolution of high-curvature regions (such as tissue boundaries), ensuring precise boundary capture. The evolution process continues to iterate until the boundary stability reaches a preset threshold (usually less than 0.5% boundary change between two consecutive iterations) or reaches the maximum number of iterations (usually 100–150).

[0087] The multi-phase hypervoxel fusion unit 215 evaluates, filters, and fuses hypervoxel sets under different parameter settings, outputting optimized hypervoxel partitioning results. Preferably, this unit employs a topological persistence-based evaluation method, defining a persistence metric:

[0088] ,

[0089] in, For supravoxels A measure of persistence. and Supervoxel The curvature values ​​of "birth" and "death" (i.e., the curvature thresholds at which the supervoxel begins to form and disappear). For the set of all supervoxels, Represents a set The maximum absolute value of the difference between the birth and death curvature values ​​of all supervoxels. Supervoxels with high persistence metrics correspond to salient structures in the image and should be preserved. Hierarchical clustering is used to fuse supervoxels under different parameter settings, retaining the combination with the highest information content.

[0090] In another embodiment of the invention, such as Figure 2 As shown, the curvature flow-driven adaptive supervoxel segmentation unit 21 further includes a supervoxel quality assurance unit 216, connected to the multi-phase supervoxel fusion unit 215, used to calculate the consistency score between the supervoxel boundary and the image gradient, evaluate the conformity of the supervoxel to the known tissue structure, and adaptively re-segment supervoxels that do not meet the quality standards. Preferably, the consistency score can be defined as:

[0091] ,

[0092] in, The consistency score for supervoxel s. Let be the set of boundary pixels of supervoxel s. For set The number of elements in the middle, Let be the normal vector at the boundary point p (a vector of unit length pointing outwards from the hypervoxel). Let p be the image gradient vector. Let Euclidean norm be the gradient vector of the image. This represents the dot product of two vectors. The consistency score reflects the degree of fit between the supervoxel boundary and the image boundary; supervoxels with scores below 0.75 are marked as needing optimization. For supervoxels requiring optimization, appropriate remedial measures are selected based on the reason for failure, such as adjusting local parameters for re-segmentation, with a maximum of 3 re-segmentation attempts.

[0093] The multi-scale geometric feature extraction unit 22 on the Riemannian manifold is used to calculate a set of differential geometric features for each supervoxel, including curvature statistics, shape exponents, and tensor field features. For example... Figure 4 As shown, this unit includes a differential geometric invariant calculation unit 221, a multi-scale thermonuclear analysis unit 222, and a spectral decomposition and feature enhancement unit 223. The differential geometric invariant calculation unit 221 first calculates the curvature statistics of each supervoxel, including the mean curvature, Gaussian curvature, mean, variance, skewness, and kurtosis of the principal curvatures. Preferably, for a supervoxel s, its curvature statistics can be expressed as:

[0094] ,

[0095] in, Let be the curvature statistical eigenvector of the supervoxel s. and Let be the sets of Gaussian curvature and mean curvature of each point within the supervoxel s, respectively. for The mean, for variance for skewness, for peak, , , and They are respectively The mean, variance, skewness, and kurtosis.

[0096] Subsequently, the differential geometric invariant calculation unit 221 calculates the shape index and curvature index, discretizing the continuous shape features into typical categories. Shape Index and curvature index The calculation formula is:

[0097] ,

[0098] ,

[0099] in, The shape index, The curvature index, and Principal curvature (i.e., the maximum and minimum curvature values ​​of the surface at any point). The arctangent function. Shape index. The range of values ​​is This method can discretize continuous shape features into 8 to 12 typical categories, such as concave-convex, saddle-shaped, and spherical shapes. This discretization process helps to capture the typical geometric morphology of tumor tissue.

[0100] Finally, the differential geometric invariant computation unit 221 constructs the structure tensor field and extracts features, including anisotropy, linearity, flatness, and sphericity. The formula for calculating the structure tensor J is:

[0101] ,

[0102] in, It is a structure tensor (3×3 matrix). The image gradient (3×1 vector). For tensor product operations, for The transpose of (1×3 vector). Representing matrix multiplication. Let the structure tensor be... eigenvalues Then anisotropy linearity Flatness and sphericity It can be represented as:

[0103] ,

[0104] ,

[0105] ,

[0106] ,

[0107] in, Anisotropy (representing the degree of consistency in structural orientation). Linearity (representing the degree to which the structure is linear). Flatness (representing the degree to which a structure is planar). Sphericity (indicating the degree to which a structure is spherical). , and For structure tensor The eigenvalues ​​(sorted from largest to smallest) are used to effectively characterize the structural properties of tumor tissue, such as the degree of fibrosis and cell orientation.

[0108] The multi-scale heat kernel analysis unit 222 constructs the Laplacian matrix of the supervoxel graph and calculates the heat kernel signature at multiple time scales to generate multi-scale heat kernel signature features. Preferably, the Laplacian matrix of the supervoxel graph... It can be represented as:

[0109] ,

[0110] in, For the Laplace matrix ( matrix, (number of hypervoxels) A degree matrix (diagonal matrix, diagonal elements) ), The weight matrix ( (matrix), its elements Indicates supravoxel and The similarity between them. Similarity calculation can take into account both spatial distance and feature similarity:

[0111] ,

[0112] in, The elements of the weight matrix represent hypervoxels. and Similarity between them For supravoxels and Spatial distance (Euclidean distance) This represents the spatial distance (Euclidean distance). and The scaling parameter is typically set to the standard deviation of the spatial distance and feature distance distributions. It is an exponential function. Based on the Laplace matrix, the hot kernel signature... The calculation formula is:

[0113] ,

[0114] in, For point In time parameters The hot core signature below, and These are the first and second digits of the Laplace matrix. 1 eigenvalue and eigenvector, For time parameters, To retain the number of feature vectors, It is the base of the natural logarithm (approximately 2.71828). For feature vectors At point The value at that location. The heat kernel signature is calculated across multiple timescales, typically 6-8 time points are chosen, from the microscopic (…). From macro to macro ( ), to capture geometric features at different scales.

[0115] Subsequently, the multi-scale thermonuclear analysis unit 222 calculates the diffusion distance between supervoxels and constructs a functional similarity matrix. Diffusion distance The calculation formula is:

[0116] ,

[0117] in, For point and In time parameters The diffusion distance below, The first of the Laplace matrix 1 eigenvalue, and They are the feature vectors At point and The value at that location, For time parameters, To retain the number of feature vectors, It is the base of the natural logarithm. This distance metric reflects the ease with which information spreads from one point to another during diffusion, and can effectively capture the functional correlations between supervoxels.

[0118] The spectral decomposition and feature enhancement unit 223 performs spectral decomposition on the supervoxel Laplacian matrix, designs a frequency domain filter to enhance key features, and outputs the enhanced feature set. Preferably, the spectral decomposition can be expressed as:

[0119] ,

[0120] in, For the Laplace matrix ( matrix), A diagonal matrix ( A matrix whose diagonal elements are eigenvalues. For the eigenvector matrix ( A matrix whose column vectors are the corresponding eigenvectors. for The transpose of . Spectral filtering can be expressed as:

[0121] ,

[0122] in, The filtered feature vector ( vector), For the original feature vector ( vector), The filter function is a diagonal matrix with diagonal elements as follows: ), The eigenvector matrix, for The transpose of . Preferably, the filter function can be designed to suppress high-frequency noise and enhance intermediate-frequency information:

[0123] ,

[0124] in, Here is the filter function. For eigenvalues, and The low-frequency and high-frequency thresholds are typically set to the 25% and 75% bits of the eigenvalue distribution, respectively. This filtering design effectively suppresses noise while preserving critical structural information.

[0125] The geodesic network tumor heterogeneity characterization unit 23 is used to construct a geodesic connection network between supravoxels to characterize intratumoral heterogeneity. For example... Figure 5 As shown, the unit includes a geodesic network construction unit 231, a heterogeneity multidimensional characterization unit 232, and a dynamic spatiotemporal heterogeneity monitoring unit 233. The geodesic network construction unit 231 first identifies physically adjacent supervoxel pairs and constructs an initial adjacency graph. Preferably, supervoxels with a shared boundary length exceeding a preset threshold (usually 5% of the supervoxel perimeter) are connected.

[0126] Subsequently, the geodesic network construction unit 231 calculates the geodesic distances between supervoxels to construct a weighted geodesic map. The geodesic distance, taking into account both spatial distance and feature similarity, can be expressed as: ,

[0127] in, The geodesic distance between supervoxels i and j Let p be the set of all possible paths from supervoxel i to j, where p is one of these paths and |p| is the path length (the number of supervoxels on the path). For the k-th supervoxel on the path, For supravoxels and Spatial distance, For supravoxels and The characteristic spatial distance, This is a weighting parameter for spatial distance and feature similarity, typically positioned between 0.3 and 0.7. This represents taking the minimum value among all possible paths. Weight parameters. This can be adjusted according to the tumor type. For adenocarcinoma, because its tissue structure is relatively dense and its boundaries are relatively clear, [the appropriate adjustment is made]. Values ​​can tend to be higher (e.g., 0.6-0.7); however, for diffuse gastric cancer, due to its strong invasiveness and indistinct borders, The value can be set lower (e.g., 0.3-0.4) to give more consideration to feature similarity.

[0128] Finally, the geodesic network building unit 231 identifies biologically significant critical paths and constructs an information-rich geodesic network based on centrality and bottleneck analysis.

[0129] The heterogeneity multidimensional characterization unit 232 calculates the characteristic variability within a local region and the overall topological characteristics of the network based on the geodesic network, generating local and global heterogeneity indices. The local heterogeneity index can be represented as the characteristic variability coefficient:

[0130] ,

[0131] in, This serves as an indicator of the local heterogeneity of supervoxel i. It is the neighborhood of supervoxel i (the set of all supervoxels directly connected to supervoxel i). The feature set of the supervoxels in the neighborhood. The standard deviation of the feature set. The mean of the feature set is denoted as . Local heterogeneity indices reflect the degree of feature variation in local tumor regions and help identify heterogeneous hotspots within the tumor.

[0132] Global heterogeneity metrics can be calculated based on the overall network topology characteristics, including clustering coefficients, path length distribution, and connectivity patterns. Preferably, the global heterogeneity metric can be defined as:

[0133] ,

[0134] in, As a global heterogeneity indicator, The total number of supervoxels. For supravoxels The clustering coefficient (representing the degree of connection between a node's neighbors). For supravoxels Betweenness centrality (representing the "bridge" role of a node in a network). For supravoxels Eigenvector centrality (representing the overall importance of a node in the network). , and For the weight parameters, satisfying .generally , and All can be set to The weighting can also be adjusted according to the specific tumor type. The global heterogeneity index reflects the complexity and degree of heterogeneity of the overall tumor structure.

[0135] The dynamic spatiotemporal heterogeneity monitoring unit 233 analyzes the rate of change of heterogeneity indicators over time, establishes a spatiotemporal model for predicting tumor development, and generates a spatiotemporal heterogeneity prediction report. Preferably, for patients undergoing multiple consecutive examinations, image registration is performed first to ensure sub-millimeter spatial alignment accuracy. Then, the rate of change of heterogeneity indicators is calculated:

[0136] ,

[0137] in, For supravoxels In time The rate of heterogeneous change, For supravoxels In time Heterogeneity index For supravoxels In time Heterogeneity index The time interval between the two examinations is defined as follows. Based on heterogeneous change patterns, a spatiotemporal model for predicting tumor development is established, including predictions of regional expansion, invasion patterns, and metastasis risk. Preferably, the spatiotemporal model can take the following form:

[0138] ,

[0139] in, For supravoxels In time Risk score, For supravoxels In time Predictive risk score, The risk change rate function (time-dependent) (function) Indicates from time arrive The integral of this is the cumulative amount of risk change. Risk change rate function. It can be built based on historical data and current heterogeneity indicators. This model can predict the development trends of different regions of a tumor and generate a spatiotemporal heterogeneity prediction report.

[0140] The multimodal data fusion module 3 is connected to the differential geometry supervoxel feature extraction module 2, and is used to fuse radiomics features with gene expression datasets to generate fused feature data. For example... Figure 6 As shown, this module includes a feature space unified mapping unit 31 and a multimodal feature optimization unit 32. The feature space unified mapping unit 31 first normalizes and standardizes each modal feature to eliminate differences in dimensions and scale. Preferably, the Z-score standardization method can be used.

[0141] ,

[0142] in, These are the standardized eigenvalues. These are the original eigenvalues. The mean of the features, The standard deviation of the feature is denoted by Z-score. Z-score standardization transforms the feature into a distribution with a mean of 0 and a standard deviation of 1, facilitating comparison and fusion between different features.

[0143] Subsequently, the feature space unified mapping unit 31 calculates the correlation between cross-modal features and constructs a feature correlation matrix. Preferably, a nonlinear correlation analysis method based on mutual information can be used.

[0144] ,

[0145] in, Features and Mutual information between them and There are two characteristic variables. and They are respectively and The value of , for and The joint probability distribution, and They are respectively and Marginal probability distribution, It is the natural logarithm function. Mutual information can capture nonlinear relationships between features, and is more comprehensive than the traditional correlation coefficient.

[0146] Finally, the unified feature space mapping unit 31 constructs a unified feature space that integrates multimodal information. Preferably, the following feature fusion method can be used:

[0147] ,

[0148] in, To fuse feature vectors, For image feature vectors, This is a gene expression feature vector. and These are weighting parameters used to balance the contributions of features from different modalities. This represents a vector concatenation operation. The weight parameters can be determined through cross-validation and are typically set between 0.4 and 0.6.

[0149] It should be noted that the aforementioned weight parameters, once determined during the training phase, remain fixed for all cases. The implicit assumption is that the representativeness of gene expression data for the overall tumor is consistent across cases. The inventors noted that this assumption does not hold true for highly heterogeneous cases: gene expression data, obtained from tissue homogenates or local biopsies, is essentially a spatial average of the entire tumor. The higher the spatial heterogeneity of the tumor, the worse the representativeness of this average. The global heterogeneity index output by the geodesic network tumor heterogeneity characterization unit 23 is an intrinsic quantity defined on the Riemannian manifold, covering the entire tumor and independent of grayscale, which precisely quantifies this representativeness. Therefore, in one embodiment of the invention, the feature space unified mapping unit 31 is further connected to the geodesic network tumor heterogeneity characterization unit 23, receiving the global heterogeneity index H, and adaptively modulating the fusion weights of the two modalities accordingly. The formula for calculation is:

[0150] ,

[0151] In the formula, For image modality fusion weights, the gene expression modality fusion weights are set to 1 and 2. The difference; H is the global heterogeneity index, and its value is calculated in the same way as the aforementioned global heterogeneity index formula. It is a dimensionless quantity. , which is the heterogeneity benchmark value, and is a dimensionless quantity; is the modulation intensity coefficient, which is a dimensionless quantity; The base of the natural logarithm. The independent variable of the exponential function is obtained by multiplying the modulation intensity coefficient by the deviation of the global heterogeneity index relative to the benchmark value; both are dimensionless, and the homogeneity of dimensions holds. The value range of H is an open interval from 0 to 1, consistent with the domain of the weights; H is higher than... hour The weight of gene expression modality fusion decreases accordingly as the weight increases, which is consistent with the mechanism described above.

[0152] The median of the global heterogeneity index for all cases in the training set is used. In this example, the median is 0.42 in the 486-case training cohort, with an interquartile range of 0.31 to 0.55. The median is used instead of the mean because the global heterogeneity index exhibits a right-skewed distribution across the population, and the median ensures that modulation is unbiased at the population level and is not influenced by a few extremely heterogeneous cases. The system updates the median using a rolling median after every 100 new cases. The update increment is limited to 0.05 to avoid the baseline value jumping with each batch of enrollments. When the global heterogeneity index of a case falls outside the 5th to 95th percentile of the training set (0.22 to 0.68 in this embodiment), it indicates that the case has exceeded the calibration range of the baseline value. At this time, the system degenerates to a fixed weight of 0.5, and the prognostic assessment report indicates that the weights are not modulated, instead of extrapolating using the above-mentioned modulated formula.

[0153] The value is between 3 and 8, preferably 5. This value ensures that the variation of the image modality weight is approximately 0.2 to 0.3 when the global heterogeneity index varies within the interquartile range. When the value is greater than 10, the above bounded mapping approaches a step function, the weights undergo a hard switch near the baseline value, and the classification results of adjacent cases will jump due to small differences in heterogeneity indicators; When the value is less than 1, the weight changes by less than 0.1 within the entire interquartile range, rendering the modulation ineffective. Similarly, in each update The validation set consistency index is then reselected in steps of 0.5 within the range of 3 to 8.

[0154] The multimodal feature optimization unit 32 evaluates the contribution of each feature to the target task, identifies and removes highly redundant features, and optimizes feature combinations for genotyping and prognostic prediction respectively. Preferably, feature importance evaluation can be carried out by comprehensively considering multiple methods:

[0155] ,in, Features Importance score Features With target variable mutual information, For feature-based Predict target variable The area under the receiver operating characteristic curve (reflecting classification performance). Features With target variable The statistical significance test p-value, , and These are weight parameters, usually set to... .

[0156] For redundant feature filtering, a minimum redundancy and maximum relevance strategy can be adopted:

[0157] ,

[0158] in, For feature subset The minimum redundancy maximum correlation score, The size of the feature subset (i.e., the number of features). and For feature subset Features in Features With target variable mutual information, Features With features Mutual information between them It represents the average mutual information (correlation measure) between features and the target variable. This represents the average mutual information (redundancy measure) between features. By maximizing... You can select a subset of features that are highly correlated with the target variable and have low redundancy. The correlation threshold is usually set to 0.85; feature pairs higher than this threshold are considered highly redundant and need to be filtered.

[0159] The aforementioned fusion weights and feature combinations are calibrated with the target of typing accuracy. The inventors further noted that the optimal weights for typing and the optimal weights for prognostic prediction do not overlap: the typing task relies on discriminative features that can distinguish molecular subtypes, while the prognostic task relies on features related to tumor burden and invasiveness. The optimal subsets of the two in the feature space only partially overlap. If the weights calibrated with the typing target are directly used for prognostic modeling, the prognostic consistency index will be systematically low. Therefore, in one embodiment of the present invention, the multimodal feature optimization unit 32 is connected to the prognostic prediction and evaluation module 5, receiving the consistency index of the survival prediction model on the validation set and updating it accordingly: when the consistency index of the current iteration is lower than the value of the previous iteration, it is determined that the current iteration is against the heterogeneity benchmark. With modulation intensity coefficient If the adjustment direction is incorrect, revert to the previous value and adjust in the opposite direction; if the current consistency index is higher than the previous one, continue adjusting in the same direction, with an adjustment step size of [value missing]. Take 0.02, Set the weight to 0.5; after each adjustment, regenerate the fusion weights of the two modalities and re-execute feature optimization to obtain a new feature combination and a new prognostic model.

[0160] The feedback loop must have a clear termination condition; otherwise, repeated optimization on the validation set will lead to overfitting and potential self-oscillation of weights. In this embodiment, the termination condition is: the consistency index increment between two consecutive iterations is less than 0.002, or the number of iterations reaches 5. The convergence threshold of 0.002 is based on the fact that, in this embodiment, the standard error of the consistency index obtained by resampling 1000 times using the bootstrap method on the 486-case cohort is approximately 0.021. 0.002 is about one-tenth of this standard error; improvements below this threshold are statistically indistinguishable from random fluctuations, and continued iteration will only fit the validation set noise. The upper limit of 5 iterations is set because, in all experiments, the consistency index increment drops below the threshold after the third to fourth iteration, and the 5-iteration upper limit provides redundancy. When the upper limit is reached and the increment still has not converged, the system stops iterating and reverts to the weight and feature combination obtained in the first iteration, recording a non-convergence flag in the system log to ensure predictable behavior under abnormal data. Furthermore, to avoid closed-loop feedback introducing test set information into the training, the consistency index is only calculated on the validation set, and the test set and external validation set remain invisible throughout the entire iteration process.

[0161] The subtyping model construction module 4 is connected to the multimodal data fusion module 3, and is used to train a tumor molecular subtyping model based on fused feature data. For example... Figure 7As shown, this module includes a training data preparation unit 41, a hierarchical model building unit 42, and a model optimization and validation unit 43. The training data preparation unit 41 integrates image features and molecular subtyping labels, dividing the dataset into training, validation, and test sets. Preferably, the dataset is divided into 70% training set, 15% validation set, and 15% test set. To enhance the diversity of the training set, appropriate data augmentation techniques, such as random rotation, translation, and scaling, can be used.

[0162] The hierarchical model building unit 42 constructs a two-layer classification model. The first layer distinguishes major molecular subtypes, and the second layer refines the subtype classification. Preferably, the first-layer model can use the following classification function:

[0163] ,

[0164] in, For the sample Belongs to molecular subtype The probability, For classification function (to classify samples) Functions that map to real numbers. The number of major molecular subtypes, is the base of the natural logarithm. This is a normalization factor (ensuring the sum of all probabilities is 1). Classification function. Various forms can be used, such as linear functions and multilayer perceptrons, with the appropriate model complexity selected based on the feature-to-sample ratio. When the number of features is much greater than the number of samples, a linear model or a regularization term can be introduced to avoid overfitting; when the number of samples is sufficient, a more complex nonlinear model, such as a random forest or neural network, can be chosen.

[0165] The second-level model, based on the first-level classification, further subdivides each major subtype:

[0166] ,

[0167] in, For the sample Belonging to the main subtype Under these conditions, it belongs to a subtype. The probability, For the subdivision classification function (to classify samples) Functions that map to real numbers. As the main subtype The number of subtypes, This is a normalization factor. This hierarchical model structure can effectively handle the hierarchical relationships of tumor molecular subtyping and improve classification accuracy.

[0168] The model optimization and validation unit 43 optimizes the model parameters and validates the performance, outputting the optimized subtyping model. Preferably, model parameter optimization can employ grid search or Bayesian optimization methods, evaluating the performance of different parameter combinations on the validation set. Performance evaluation metrics include precision, recall, F1 score, and AUC. Taking AUC as an example, its calculation formula is:

[0169] ,

[0170] in, The area under the receiver operating characteristic curve. For the positive sample set, For the negative sample set, The number of positive samples. The number of negative samples. Output a score for the model (the larger the value, the more likely it is to be a positive sample). As a positive sample, For negative samples, This is an indicator function (value 1 when the condition is true, otherwise 0). A higher area under the receiver operating characteristic (ROC) curve indicates better model performance. The optimal parameter combination is selected through cross-validation on the validation set to construct the final subtyping model.

[0171] The prognostic prediction and assessment module 5 is connected to the subtyping model construction module 4, and is used to predict patient prognosis and generate prognostic assessment reports based on tumor molecular subtyping models and fusion feature data. For example... Figure 8 As shown, this module includes a survival analysis model unit 51, a treatment response prediction unit 52, and a comprehensive risk assessment unit 53. The survival analysis model unit 51 integrates imaging features and survival data to construct a survival prediction model. Preferably, a Cox proportional hazards model can be used.

[0172] ,

[0173] in, For the sample In time The risk function (instantaneous mortality rate). The baseline risk function (based risk rate independent of features). This is the feature weight vector (representing the degree of influence of each feature on risk). for transpose, Representing the eigenvector With weight vector The inner product, It is an exponential function. Meanwhile, to capture the non-linear relationship between features and survival time, a deep survival analysis model can be used:

[0174] ,

[0175] in, For the sample In time The risk function, As the benchmark risk function, For deep neural networks (parameters are) The function that will take samples Mapping to real numbers). It is an exponential function. During model training, clinical factors such as age, gender, and clinical stage can be considered as covariates to improve predictive accuracy.

[0176] The treatment response prediction unit 52 analyzes the correlation between imaging features and treatment response to predict chemotherapy response rate and radiosensitivity. Preferably, the following prediction model can be used:

[0177] ,

[0178] in, For the sample The probability of responding to treatment. For the sigmoid function, For the weight vector, For bias terms, for transpose, Representing the eigenvector With weight vector The inner product of the two. The model can be validated using follow-up data to ensure the reliability of the predictions.

[0179] The comprehensive risk assessment unit 53 integrates survival and treatment response prediction results, classifies patients into high, medium, and low risk groups, and generates an individualized risk assessment report. Preferably, the comprehensive risk score can be expressed as:

[0180] ,

[0181] in, For the sample Comprehensive risk score, The survival risk score is typically a risk score predicted by the Cox model and normalized to the [0,1] interval. For treatment response probability, The risk score is a weighting parameter, typically set between 0.5 and 0.7, reflecting the importance of survival risk in the overall assessment. Based on the risk score distribution, patients can be divided into high-risk (risk score > 75th percentile), medium-risk (25th percentile < risk score ≤ 75th percentile), and low-risk (risk score ≤ 25th percentile) groups. The risk score is calibrated based on clinical data to ensure the accuracy and reliability of the assessment results.

[0182] This invention also includes a system performance evaluation module 6 and a clinical decision support module 7. The system performance evaluation module 6 is connected to the prognosis prediction and evaluation module 5, and is used to evaluate system performance on a test set and independent external datasets, assess system stability under different conditions, and generate a performance evaluation report. Preferably, when evaluating system performance on the test set, indicators such as subtyping accuracy, area under the receiver operating characteristic (ROC) curve, and survival prediction C-index can be calculated and compared with traditional methods and clinician judgment. The formula for calculating the C-index is:

[0183] ,

[0184] in, As a consistency index, Indicates sample Has an event (such as death) occurred? (Indicates an event has occurred) For the sample Survival time, The risk score predicted by the model, For indicator functions, This indicates the number of sample pairs that agree with the predictions. This represents the total number of comparable sample pairs. A higher C-index value indicates better survival prediction performance.

[0185] Simultaneously, to verify the system's generalization ability, validation on independent external datasets is necessary. These external datasets should consist of data collected from different centers using different devices. Due to the differences between centers, domain adaptation techniques, such as feature normalization and batch effect correction, can be employed.

[0186] Furthermore, the system performance evaluation module 6 also evaluates the system's stability under different conditions, including performance under different image qualities and scanning parameters. By analyzing the performance fluctuation range, the system's robustness is assessed. Preferably, performance fluctuations should be controlled within 10% to ensure the system's stability in practical applications.

[0187] In another embodiment of the present invention, the comprehensive risk assessment unit 53 is also connected to the dynamic spatiotemporal heterogeneity monitoring unit 233 to perform a sign consistency check on the rate of change of global heterogeneity indicators and the rate of change of comprehensive risk score between two adjacent examinations of the same patient. Mechanistically, both should be in the same direction: increased tumor heterogeneity is usually accompanied by subclonal expansion and increased invasiveness, thus increasing the risk score. When the signs of the two rates of change are opposite, it indicates that the evidence from the imaging end and the model end contradict each other. Common causes include inconsistent scanning parameters between the two examinations, excessive registration errors, or changes in the lesion delineation range. In this case, the system does not directly output a conclusion, but instead lowers the confidence level of the patient's individualized risk assessment report by one level and outputs a review indicator, prompting re-registration or re-examination, thereby avoiding misleading risk stratification when evidence contradicts each other. This check is an extension of the aforementioned closed-loop coupling in the time dimension.

[0188] The clinical decision support module 7 connects with the subtyping model construction module 4 and the prognosis prediction and assessment module 5 to transform analysis results into clinically friendly reports. These reports integrate with hospital PACS and electronic medical record systems, collecting clinical feedback and continuously improving the system. Preferably, the clinical reports include subtyping results, prognostic predictions, and treatment recommendations, displaying key results and confidence intervals in an intuitive visualization for easy understanding and use by clinicians. The system integrates with hospital PACS and electronic medical record systems via standard DICOM and HL7 interfaces, ensuring the security and reliability of data exchange. Simultaneously, a structured feedback table and system log analysis mechanism are designed to collect clinical feedback, providing a basis for continuous system optimization.

[0189] To verify the synergistic effect of the various modules of this invention after systematic integration, and the improvement in classification accuracy and prognostic prediction precision compared with the prior art, the following comparative experiment was conducted.

[0190] The experimental data were obtained from a retrospective cohort of 486 cases of gastrointestinal tumors diagnosed by pathology and molecular testing, including 268 cases of gastric cancer and 218 cases of colorectal cancer. Each case had complete digestive endoscopy images, enhanced CT images, gene expression profiles, and follow-up records of at least 36 months. Follow-up was calculated from the date of pathological diagnosis and ended at death or the last follow-up. The median follow-up time for the entire cohort was 41 months, and the interquartile range was 38 to 46 months. Patients still alive at 36 months of follow-up were censored. The training, validation, and test sets were divided into training, validation, and test sets at a ratio of 70%, 15%, and 15%, respectively. The results were 340 cases in the training set, 73 cases in the validation set, and 73 cases in the test set. The subtypes were stratified by molecular subtype to maintain a consistent proportion of subtypes in each set. Another 152 cases collected from two other medical institutions using equipment from different manufacturers were used as an independent external validation set. The inclusion criteria, data integrity requirements and follow-up criteria of this external validation set were consistent with the aforementioned cohort, and it did not participate in any training or parameter selection.

[0191] The calculation methods for the above indicators are as follows: Subtyping accuracy is the proportion of cases where the predicted molecular subtype is consistent with the pathologically and molecularly confirmed subtype to the total number of cases in the set; the area under the receiver operating characteristic (ROC) curve is calculated separately for each molecular subtype on a one-to-many basis, and then weighted by the number of cases in each subtype; the concordance index is calculated on the same set using the formula described above in this manual; survival rate is estimated using the Kaplan-Meier method, with log-rank tests used for inter-group comparisons, and DeLong tests used for inter-group comparisons of the ROC curve, all with a two-sided significance level of 0.05. All the above indicators and the values ​​in the tables below were calculated using the methods described in this section on the aforementioned cohort of 486 cases and the independent external validation set of 152 cases.

[0192] The control schemes are set up as follows: Comparative Example 1 is a conventional radiomics baseline scheme, using gray-level co-occurrence matrix and wavelet texture features, using only image modalities, and each step is unidirectionally connected; Comparative Example 2 is based on Comparative Example 1, except that the feature extraction is replaced with the differential geometric supervoxel feature extraction of this invention, and the rest remains unchanged; Comparative Example 3 is based on Comparative Example 1, except that gene expression modalities are added and fused with fixed weights determined by cross-validation, and the feature extraction is still conventional texture features; Comparative Example 4 is based on Comparative Example 1, except that the feedback loop of this invention is introduced, and the feature extraction and the number of modalities remain unchanged; Comparative Example 5 uses the differential geometric supervoxel feature extraction and multimodal fusion of this invention simultaneously, but without setting a feedback loop, that is, a simple superposition scheme with three steps unidirectionally connected; Comparative Example 6 is a convolutional neural network scheme, trained end-to-end with a three-dimensional residual network; Example 1 is the complete scheme of this invention. All schemes are trained and tested under the same data partition, and the evaluation indicators are the genotyping accuracy, the area under the receiver operating characteristic curve, and the prognostic consistency index. The performance of each scheme on the test set is shown in Table 1.

[0193] Table 1. Performance comparison of each scheme on the test set

[0194]

[0195] With Comparative Example 1 as the baseline, the gains when the three elements of differential geometric feature extraction, multimodal fusion, and feedback loop are introduced individually and together are shown in Table 2.

[0196] Table 2 Factorial ablation and synergistic gain of the three elements

[0197]

[0198] As shown in Table 2, when the differential geometric feature extraction and multimodal fusion are simply concatenated, the measured gain in classification accuracy is 10.3 percentage points, lower than the sum of their individual gains of 11.2 percentage points. The measured gain in the consistency index, 0.071, is also lower than the sum of their individual gains of 0.101. This indicates that simply concatenating the two steps not only does not produce additional gain, but also causes a slight loss due to feature redundancy, which is consistent with the conventional expectations in this field. However, when all three elements are present, the measured gain in classification accuracy jumps to 17.4 percentage points, which is 5.0 percentage points higher than the expected linear superposition of the gains of the three elements alone, with a nonlinear synergy factor of 1.40 times; the synergy factor of the consistency index is 1.36 times. It is particularly noteworthy that the gain when the feedback loop is introduced alone is only 1.2 percentage points and 0.004, which is close to the fluctuation range of repeated experiments. However, when differential geometric feature extraction and multimodal fusion are present at the same time, its marginal contribution jumps to 7.1 percentage points and 0.072. The effect of the same technical feature in different combination environments is nearly six times different. This strong nonlinear relationship cannot be extrapolated from the effect of the feature alone.

[0199] The mechanism of the aforementioned synergistic effect lies in the fact that the feedback loop modulates the relative credibility of two modes, and the credibility measure must be a physical quantity that is fully covered and comparable across devices. The global heterogeneity index is defined by the geodesic network on the Riemannian manifold, which depends only on the intrinsic geometry of the image manifold and not on the absolute value of grayscale, thus satisfying this requirement. Conventional texture features, on the other hand, drift with window width and window level, resampling, and device differences, and do not satisfy this requirement. To verify this mechanism, a comparative example 7 was set up. Based on Example 1, only the amount of driving fusion weight modulation was replaced by the grayscale non-uniformity of conventional texture features instead of the global heterogeneity index. Everything else was the same. Its test set classification accuracy dropped to 83.4%, the consistency index dropped to 0.719, and the synergistic gain basically disappeared. Further examining its cross-center performance, the classification accuracy of Comparative Example 7 on the external validation set dropped to 74.9%, which is 8.5 percentage points lower than that on its internal test set, while the corresponding difference in Example 1 was only 2.8 percentage points. This result directly proves that the synergistic effect comes from the cross-device comparability of differential geometric indices, rather than from the addition of a feedback loop itself.

[0200] The performance of each scheme on the independent external validation set is shown in Table 3.

[0201] Table 3 Generalization performance on independent external validation sets

[0202]

[0203] Statistical tests showed that the difference between the area under the receiver operating characteristic (ROC) curves of Example 1 and Comparative Example 6 on the test set was statistically significant (p<0.01) according to the DeLong test, and also statistically significant (p<0.01) compared to Comparative Example 5. Survival analysis was performed on 225 individuals, a subset of 73 from the test set and 152 from the independent external validation set, which were not visible during training and parameter selection. Using the 25th and 75th percentiles of the comprehensive risk score output from Example 1 as cutoffs, these individuals were divided into a low-risk group (n=56), a medium-risk group (n=113), and a high-risk group (n=56). The three-year survival rates estimated using the Kaplan-Meier method were 79.1%, 54.3%, and 28.6%, respectively. The differences among the three groups were statistically significant (p<0.001), with a 50.5 percentage point difference in three-year survival rates between the high-risk and low-risk groups. Using the risk scores from Comparative Example 1, the same 225 cases were divided according to the same quantile cutoff and estimation method. The three-year survival rates for the low, medium, and high-risk groups were 66.5%, 53.8%, and 41.2%, respectively. The difference between the high-risk and low-risk groups was 25.3 percentage points, less than half the difference in Example 1, indicating a significantly lower group discrimination than in Example 1. All three-year survival rates were directly estimated from the follow-up records of these 225 cases using the Kaplan-Meier method, with the date of pathological diagnosis as the starting point and a follow-up period of 36 months as the time node, without any extrapolation. Stability testing was conducted on the test set. The same batch of images were resampled at three layer thicknesses of 1.0mm, 3.0mm, and 5.0mm. For each layer thickness, the original image and images with Gaussian noise added (with a standard deviation of 2% and 5% of the original image's grayscale standard deviation) were used, resulting in nine combinations of image quality and scanning parameters. The complete classification and prognosis process was repeated one by one. The classification accuracy of Example 1 ranged from 87.9% to 90.1% under the nine combinations, and the consistency index ranged from 0.771 to 0.795. The maximum relative deviations from the corresponding values ​​without resampling and with added noise (classification accuracy 89.2%, consistency index 0.784) were 1.5% and 1.7%, respectively. The performance fluctuations were all less than 10%, meeting the aforementioned stability requirements.

[0204] In summary, the technical effects of this invention compared to existing technologies cannot be expected from the simple summation of the effects of individual modules: First, the measured effect of simply connecting two functional modules in series is lower than the sum of their individual effects; second, the feedback loop is almost ineffective when it exists alone, and only makes a significant contribution when differential geometric feature extraction and multimodal fusion coexist, with the contribution magnitude being nearly six times that of its individual effects; third, the synergistic effect disappears after replacing the quantity driving the closed loop with conventional texture features, indicating that the synergy originates from the non-obvious property of cross-device comparability unique to differential geometric indices. These three points together constitute the unexpected technical effects of this invention.

[0205] The system of this invention constructs a medical image feature extraction and analysis framework by introducing differential geometry theory, which is used to extract the geometric structural features and heterogeneity information of tumor tissue, thereby achieving molecular subtyping and prognostic prediction of gastrointestinal tumors. Compared with the prior art, the subtyping accuracy, area under the receiver operating characteristic curve, and prognostic consistency index on the same test set of this invention are improved from 71.8%, 0.742, and 0.641 to 89.2%, 0.921, and 0.784, respectively. Moreover, the above improvements are not a simple summation of the effects of each module, but a synergistic effect generated after the three modules form a closed loop through global heterogeneity index modulation and consistency index feedback.

[0206] The embodiments described above are merely illustrative of specific implementations of the present invention, and while the descriptions are detailed, they should not be construed as limiting the scope of the present invention. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the present invention, and these modifications and improvements all fall within the scope of protection of the present invention.

Claims

1. A precise classification and prognostic prediction system for gastrointestinal tumors based on multimodal medical images, characterized in that, include: The data preprocessing module is used to acquire digestive endoscopy image datasets, CT image datasets, and gene expression datasets, and to standardize the acquired data. A differential geometry supervoxel feature extraction module, connected to the data preprocessing module, is used to receive the standardized image data and extract radiomics features based on differential geometry theory. The differential geometry supervoxel feature extraction module includes a curvature flow-driven adaptive supervoxel segmentation unit, a multi-scale geometric feature extraction unit on a Riemannian manifold, and a geodesic network tumor heterogeneity characterization unit. The multimodal data fusion module, connected to the differential geometry supervoxel feature extraction module, is used to fuse the radiomics features with the gene expression dataset to generate fused feature data, and adaptively modulate the fusion weights of the image modality and the gene expression modality using the global heterogeneity index output by the geodesic network tumor heterogeneity characterization unit. The subtyping model construction module is connected to the multimodal data fusion module and is used to train a tumor molecular subtyping model based on the fused feature data; The prognosis prediction and assessment module, connected to the subtyping model construction module, is used to predict patient prognosis and generate a prognosis assessment report based on the tumor molecular subtyping model and the fusion feature data. The prognosis prediction and assessment module uses the molecular subtype posterior probability output by the tumor molecular subtyping model as a covariate for survival analysis and feeds back the consistency index output by the prognosis prediction and assessment module to the multimodal data fusion module to update the fusion weights and feature combinations, so that the differential geometric hypervoxel feature extraction module, the multimodal data fusion module, and the prognosis prediction and assessment module form a closed loop.

2. The system according to claim 1, characterized in that, The curvature flow-driven adaptive supervoxel segmentation unit includes: Image manifold mapping unit, used to map medical image data to a high-dimensional Riemannian manifold; The curvature field calculation unit, connected to the image manifold mapping unit, is used to calculate the Gaussian curvature and average curvature of each point on the high-dimensional Riemannian manifold and generate a curvature field mapping. A seed point optimization placement unit is connected to the curvature field calculation unit and is used to identify key regions and optimize the distribution of seed points based on the curvature field mapping. The curvature flow evolution control unit, connected to the seed point optimization placement unit, is used to design curvature-dependent velocity field control of supervoxel boundary evolution and generate an initial supervoxel set. The multi-phase supervoxel fusion unit, connected to the curvature flow evolution control unit, is used to evaluate, screen, and fuse supervoxel sets under different parameter settings, and output the optimized supervoxel partitioning results.

3. The system according to claim 2, characterized in that, The multi-scale geometric feature extraction unit on the Riemannian manifold includes: The differential geometric invariant computation unit is used to calculate a set of differential geometric features for each hypervoxel, including curvature statistics, shape exponent, and tensor field features; The multi-scale thermonuclear analysis unit, connected to the differential geometric invariant calculation unit, is used to construct the Laplacian matrix of the supervoxel diagram and calculate the thermonuclear signature at multiple time scales to generate multi-scale thermonuclear signature features. The spectral decomposition and feature enhancement unit, connected to the multi-scale thermonuclear analysis unit, is used to perform spectral decomposition on the supervoxel Laplacian matrix, design frequency domain filters to enhance key features, and output the enhanced feature set.

4. The system according to claim 1, characterized in that, The geodesic network tumor heterogeneity characterization unit includes: The geodesic network construction unit is used to identify physically adjacent hypervoxel pairs and calculate the geodesic distance between hypervoxels to construct a weighted geodesic map. The heterogeneity multidimensional characterization unit is connected to the geodesic network construction unit and is used to calculate the feature variability and overall network topology characteristics in a local area based on the weighted geodesic map, and generate local and global heterogeneity indices. The dynamic spatiotemporal heterogeneity monitoring unit, connected to the heterogeneity multidimensional characterization unit, is used to analyze the rate of change of heterogeneity indicators over time, establish a spatiotemporal model for predicting tumor development, and generate a spatiotemporal heterogeneity prediction report.

5. The system according to claim 1, characterized in that, The multimodal data fusion module includes: The unified feature space mapping unit is used to normalize and standardize the features of each modality, calculate the correlation between cross-modal features, and construct a unified feature space that integrates multimodal information. The multimodal feature optimization unit, connected to the feature space unified mapping unit, is used to evaluate the contribution of each feature to the target task, identify and remove highly redundant features, and optimize feature combinations for classification and prognosis prediction respectively.

6. The system according to claim 1, characterized in that, The fractal model construction module includes: The training data preparation unit is used to integrate image features and molecular typing labels, and to divide the training set, validation set, and test set. The hierarchical model building unit, connected to the training data preparation unit, is used to build a two-layer classification model. The first layer distinguishes the main molecular subtypes, and the second layer refines the subtype classification. The model optimization and verification unit, connected to the hierarchical model construction unit, is used to optimize model parameters and verify performance, and output the optimized fractal model.

7. The system according to claim 1, characterized in that, The prognostic prediction and assessment module includes: The survival analysis model unit is used to integrate image features and survival data to build a survival prediction model. The treatment response prediction unit, connected to the survival analysis model unit, is used to analyze the correlation between imaging features and treatment response, and to predict chemotherapy response rate and radiotherapy sensitivity. The integrated risk assessment unit, connected to the survival analysis model unit and the treatment response prediction unit, is used to integrate survival and treatment response prediction results, classify patients into high, medium and low risk groups, and generate individualized risk assessment reports.

8. The system according to claim 1, characterized in that, The system also includes: The system performance evaluation module, connected to the prognostic prediction and evaluation module, is used to evaluate system performance on test sets and independent external datasets, evaluate the stability of the system under different conditions, and generate a performance evaluation report. The clinical decision support module, connected to the typing model construction module and the prognosis prediction and assessment module, is used to transform the analysis results into clinically friendly reports, integrate with the hospital PACS and electronic medical record system, collect clinical usage feedback, and continuously improve the system.

9. The system according to claim 3, characterized in that, The curvature flow-driven adaptive supervoxel segmentation unit further includes: The supervoxel quality assurance unit, connected to the multi-phase supervoxel fusion unit, is used to calculate the consistency score between the supervoxel boundary and the image gradient, evaluate the conformity of the supervoxel with the known tissue structure, and perform adaptive re-segmentation on supervoxels that do not meet the quality standards. The differential geometric invariant calculation unit has the following characteristics: Local curvature characteristics, including mean curvature, Gaussian curvature, mean, variance, skewness, and kurtosis of principal curvatures; Shape index features, including shape index and curvature index, discretize continuous shape features into typical categories; Tensor field characteristics, including anisotropy, linearity, flatness, and sphericity, are used to characterize the internal structural properties of supervoxels.

10. The system according to claim 5, characterized in that, The unified feature space mapping unit is also connected to the geodesic network tumor heterogeneity characterization unit, and is used to receive the global heterogeneity index, and generate image modality fusion weights by monotonically increasing bounded mapping based on the difference between the global heterogeneity index and the preset heterogeneity benchmark value. The difference between 1 and the image modality fusion weights is used as gene expression modality fusion weights. When the global heterogeneity index is higher than the heterogeneity benchmark value, the image modality fusion weights increase. The multimodal feature optimization unit is also used to receive the consistency index, and when the consistency index is lower than the value of the previous iteration, adjust the modulation intensity coefficient of the heterogeneity benchmark value and the bounded mapping, regenerate the image modality fusion weight, the gene expression modality fusion weight and the corresponding feature combination, and terminate the iteration and lock the fusion weight and the feature combination when the consistency index increment between two adjacent iterations is less than the preset convergence threshold or the number of iterations reaches the preset upper limit. The prognostic prediction and assessment module is also used to calculate the rate of change of the global heterogeneity index between two adjacent examinations and the rate of change of the risk score corresponding to the patient's prognosis within the same time interval, perform a sign consistency check on the two rates of change, and reduce the confidence of the prognostic assessment report and output a review flag when the signs of the two rates of change are opposite.

Images