Method and system for three-dimensional reconstruction of consecutive levels of pathological sections based on image registration

Abstract

The present application relates to the technical field of digital pathology image processing, and discloses a pathological section continuous layer three-dimensional reconstruction method and system based on image registration, which comprises the following steps: digitizing and scanning continuous paraffin sections and performing color normalization; adopting a double strategy of feature point rigid registration and optical flow field elastic registration for interlayer registration; generating a tumor mask through convolutional neural network semantic segmentation; constructing a three-dimensional model through contour interpolation and voxelization; calculating morphological parameters such as volume and infiltration depth; and clustering and identifying heterogeneous partitions such as necrotic areas and high proliferation areas. The present application realizes high-precision three-dimensional reconstruction and heterogeneity analysis.

Description

Technical Field

[0001] This invention relates to the field of digital pathological image processing technology, specifically to a method and system for continuous-level three-dimensional reconstruction of pathological slides based on image registration. Background Technology

[0002] Tumor pathological diagnosis is the gold standard for confirming and classifying malignant tumors. Clinical pathologists determine the histological type, degree of differentiation, and extent of invasion of tumors by observing the microscopic morphological characteristics of tissue sections. However, traditional pathological diagnosis mainly relies on the observation and analysis of single-layer or a small number of two-dimensional sections, which has inherent limitations. As a three-dimensional entity, the key information of a tumor, such as its spatial growth pattern, morphology of invasion boundaries, and internal heterogeneity, is difficult to fully present with a limited number of two-dimensional sections.

[0003] Chinese invention CN118279302A discloses a three-dimensional reconstruction detection method and system for brain tumor images. This scheme uses a weighted fusion strategy to register and fuse images of different modalities for medical imaging data such as CT and MRI, uses an adaptive fuzzy inference optimization strategy to determine the tumor target region, and realizes three-dimensional morphological reconstruction of the tumor through a three-dimensional morphological constraint energy function and a ray projection algorithm.

[0004] However, the technical solution in the aforementioned document has the following shortcomings: First, the solution is designed for tomographic imaging data such as CT and MRI, and cannot be directly applied to the three-dimensional reconstruction of pathological sections. Pathological sections undergo complex non-rigid deformations during preparation, including tissue shrinkage, folding, and tearing. These deformation patterns are fundamentally different from the registration problem of CT / MRI images. Second, the weighted fusion strategy used in this solution is mainly used for information fusion of multimodal images, rather than solving the geometric registration problem between consecutive sections, and cannot effectively correct for rotation, translation, and local deformation between pathological sections. Third, the tumor region identification method based on fuzzy reasoning relies on expert knowledge to build a rule base, which is difficult to adapt to the high diversity of tumor morphology in pathological sections and lacks the ability to analyze the heterogeneous partitioning within the tumor. Fourth, the solution does not involve the calculation of three-dimensional morphological parameters specific to pathological sections, such as tumor invasion depth and spatial relationship with the basement membrane, which are important indicators for pathological staging.

[0005] Therefore, the existing technology lacks a complete technical solution designed specifically for the characteristics of pathological sections, capable of achieving high-precision interlayer registration and three-dimensional spatial morphological analysis of tumors. Summary of the Invention

[0006] To address the shortcomings of existing technologies, this invention provides a method and system for three-dimensional reconstruction of continuous layers of pathological slides based on image registration. It corrects geometric deviations between continuous slides through a dual-strategy registration scheme that combines rigid registration of feature points with elastic registration of optical flow fields. Combined with deep learning semantic segmentation and three-dimensional voxelization reconstruction technology, it achieves accurate three-dimensional morphological reconstruction and heterogeneity analysis of tumor tissues.

[0007] A first aspect of the present invention provides a method for three-dimensional reconstruction of continuous layers of pathological slides based on image registration, comprising the following steps:

[0008] S1. Digital acquisition and preprocessing steps of slides: Digitally scan the continuous paraffin slides of tumor tissue to obtain the digital image sequence of pathological slides, and perform color normalization processing on the digital image sequence of pathological slides to eliminate color shifts caused by differences in staining batches.

[0009] S2. Dual-strategy registration steps: A strategy combining rigid registration of feature points and elastic registration of optical flow fields is used to perform inter-slice registration on the color-normalized digital image sequence of pathological slides. First, feature point sets are extracted from adjacent slide images and feature point matching is performed. Based on the matched feature point pairs, a rigid transformation matrix is ​​calculated to correct the overall rotation and translation deviations between slides. Then, an optical flow field is calculated based on the rigidly registered image pairs. An elastic deformation field is generated based on the optical flow field to compensate for local non-rigid deformations, and finally, the registered slide sequence is obtained.

[0010] S3. Semantic segmentation of tumor region: A convolutional neural network is used to perform semantic segmentation of the tumor region on each layer of slice image in the registered slice sequence, outputting a tumor probability map and generating a tumor region mask through thresholding, wherein the tumor region mask identifies the spatial boundary between tumor tissue and normal tissue.

[0011] S4. Contour Interpolation and Voxelization Reconstruction Steps: Extract two-dimensional contour curves from the tumor region mask, perform shape interpolation on the two-dimensional contour curves between adjacent slice layers to generate intermediate layer contours, and convert the interpolated contour sequence into a three-dimensional voxel model.

[0012] S5. Calculation steps for three-dimensional morphological parameters: Calculate the three-dimensional morphological parameters of the tumor based on the three-dimensional voxel model. The three-dimensional morphological parameters include tumor volume, tumor surface area, maximum diameter of the tumor, depth of invasion, and boundary relationship parameters between the tumor and surrounding normal tissues.

[0013] S6. Heterogeneous Partitioning and Visualization Steps: Based on the texture features and spatial distribution features of each voxel in the three-dimensional voxel model, cluster analysis is performed to identify the heterogeneous partitioning results of the necrotic area, high-proliferation area and invasion front area inside the tumor, and a three-dimensional visualization model of the tumor and a report on its spatial distribution features are generated.

[0014] A second aspect of the present invention provides a three-dimensional reconstruction system for continuous layers of pathological slides based on image registration, for implementing the aforementioned method, comprising: a slide acquisition module, a dual-strategy registration module, a semantic segmentation module, a voxel reconstruction module, a morphological analysis module, and a heterogeneity recognition module.

[0015] Compared with the prior art, the present invention has the following beneficial effects:

[0016] First, the present invention adopts a dual-strategy registration scheme that combines rigid registration of feature points with elastic registration of optical flow field, which can effectively correct the complex geometric deformation generated during the preparation of pathological slides. Rigid registration first eliminates the overall rotation and translation deviation, while elastic registration further compensates for local non-rigid deformation. The two work together to significantly improve the registration accuracy. Compared with a single registration method, the registration error can be reduced by about 40%.

[0017] Second, this invention uses deep learning semantic segmentation to replace traditional fuzzy reasoning for tumor region identification. It can adaptively learn the morphological features of different tumor types without relying on expert knowledge to build a rule base, and has stronger robustness to the diversity of tumor morphology.

[0018] Third, this invention achieves high-resolution three-dimensional reconstruction through contour interpolation and voxelization, which can accurately calculate three-dimensional morphological parameters such as tumor volume, surface area, and invasion depth, providing quantitative basis for pathological staging and prognostic assessment.

[0019] Fourth, this invention innovatively introduces a three-dimensional heterogeneity partitioning analysis function for tumors, which can identify functional partitions such as necrotic areas, high-proliferation areas, and invasion front areas within the tumor, revealing the spatial biological characteristics of the tumor and providing three-dimensional pathological evidence for precision medical decision-making. Attached Figure Description

[0020] Figure 1 This is a flowchart of a method for three-dimensional reconstruction of continuous layers of pathological slides based on image registration, provided in an embodiment of the present invention.

[0021] Figure 2 This is an architecture diagram of a continuous-layer three-dimensional reconstruction system for pathological slides based on image registration, provided in an embodiment of the present invention. Detailed Implementation

[0022] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments.

[0023] Figure 1 A flowchart illustrating a method for continuous-slice three-dimensional reconstruction of pathological slides based on image registration, provided in an embodiment of the present invention, is shown. In one embodiment of the present invention, the method includes the following steps:

[0024] Step S1: Digital acquisition and preprocessing of pathological slides. The main purpose of this step is to obtain high-quality digital image sequences of pathological slides and eliminate interference caused by staining differences, laying the foundation for subsequent registration and segmentation.

[0025] In one specific embodiment, the tumor tissue sample is first prepared by serial paraffin sectioning. Preferably, a fully automated rotary microtome is used to perform serial sectioning of formalin-fixed paraffin-embedded tissue blocks, with a section thickness of [missing information]. A thickness of 4 μm is provided to ensure slice integrity while offering sufficient interlayer resolution. In one embodiment of the invention, the total number of layers in a continuous slice... The number of layers can be determined based on the thickness of the tumor tissue, typically between 20 and 150 layers, to ensure complete coverage of the tumor's three-dimensional structure.

[0026] After section preparation, routine H&E staining is performed. Hematoxylin staining time is controlled between 3 and 8 minutes, preferably 5 minutes, and eosin staining time is controlled between 30 and 120 seconds, preferably 60 seconds. After stained sections are dehydrated, cleared, and mounted, they are placed on glass slides for scanning.

[0027] Subsequently, the stained sections were digitally scanned using a full-slide scanner. In one embodiment of the present invention, the scanning magnification was set to 40×, corresponding to a physical resolution of approximately 0.25 μm / pixel, which is sufficient to clearly present the morphological features of cells and tissue structures. A digital image sequence of the pathological sections was obtained after the scanning was completed. ,in Indicates the first Digital images of sliced ​​layers.

[0028] Because color differences exist between different batches of staining, this invention further performs color normalization processing on the digital image sequence of pathological sections. Preferably, a staining separation method based on color deconvolution is used to achieve color normalization. This method first converts the RGB color space to an optical density space, and the conversion formula is: ,in, Indicates the optical density value. This represents the RGB channel intensity values ​​of the image. This represents the light source intensity, typically taken as 255. Based on the optical properties of H&E staining, the color deconvolution matrices of hematoxylin and eosin are established. By solving a system of linear equations, the mixed staining signal was separated into hematoxylin channels. and Yihong Channel : ,in, , , These represent the optical density values ​​for the red, green, and blue channels, respectively. Represents the residual term. Color deconvolution matrix. The elemental values ​​were determined based on the standard spectral absorption curve of H&E staining, and the specific values ​​were obtained through pre-calibration.

[0029] After separation, the intensities of the hematoxylin and eosin channels are normalized to a preset standard staining intensity range. In one embodiment of the invention, the standard intensity range is determined as follows: 100 high-quality stained slice images are randomly selected from the training set, and the mean and standard deviation of the intensity of each channel are calculated. The mean ± 2 times the standard deviation is used as the standard intensity range. The normalized image is then reconstructed into an RGB image through the inverse transform of color deconvolution, completing the color normalization process.

[0030] After the above processing, the staining differences between different slices are effectively eliminated, providing a consistent image input for subsequent registration steps. Test results of this invention show that color normalization processing can quantify the color differences between slices (in terms of the CIE Lab color space). (This indicates that) the average value decreased from 12.5 to 3.2, a decrease of approximately 74%.

[0031] Preferably, the present invention also performs quality screening on the digitized slice images. The quality screening criteria include: image sharpness assessment using the Laplacian operator to calculate the image gradient variance; images with a variance value below a preset threshold (e.g., 100) are judged as blurry images and marked as requiring rescanning; tissue integrity assessment by detecting the connectivity of tissue regions; slices with severe tissue tearing or detachment need to be re-prepared; staining quality assessment by analyzing the intensity distribution of hematoxylin and eosin channels; slices stained too darkly or too lightly can have their normalization parameters selectively adjusted or be marked as abnormal. The slice image sequence retained after quality screening ensures the quality of the basic data for 3D reconstruction.

[0032] In one embodiment of the present invention, the digital acquisition of slides also includes recording slide position information. A slide holder with a unique identifier is used to automatically record the relative position number of each slide within the tissue block during scanning, ensuring the correct spatial order of the slide sequence. Furthermore, for large slides, the scanner uses a partitioned scanning and stitching method to acquire a complete image. During stitching, feature matching of overlapping areas automatically corrects the stitching gaps, ensuring the continuity of the stitched image.

[0033] Step S2: Dual-strategy registration step. This step is the core innovation of the entire method. Through the synergistic cooperation of rigid registration of feature points and elastic registration of optical flow field, geometric deviations between consecutive sections are effectively corrected. Pathological sections undergo various deformations during preparation, including overall rotation (caused by inconsistent slide placement angles), overall translation (caused by tissue displacement on the slide), and local non-rigid deformation (caused by tissue shrinkage, folding, tearing, etc.). A single rigid registration method cannot handle local deformations, while a single elastic registration method is prone to getting trapped in local optima when large displacements exist. The dual-strategy registration scheme adopted in this invention first eliminates overall rotation and translation through rigid registration, simplifying the problem to small-displacement elastic registration, and then uses the optical flow field method to finely correct local deformations, thereby achieving high-precision registration.

[0034] In one embodiment of the present invention, the registration process is performed in an iterative manner. Taking the first... Layer slice image For a fixed image (reference image), the first Layer slice image For the floating image (the image to be registered), the registration of all adjacent slice pairs is completed sequentially, finally obtaining the registered slice sequence. .

[0035] S2-1, Feature point extraction: First, a set of feature points is extracted from adjacent slice images. This invention preferably uses the Scale Invariant Feature Transform (SIFT) algorithm for feature point extraction. This algorithm has good invariance to image rotation, scale transformation, and brightness changes, and is suitable for handling various transformations that may exist between pathological slices.

[0036] The SIFT feature point extraction process includes the following sub-steps: First, a Gaussian pyramid and a Differential Gaussian (DoG) pyramid are constructed. The location and scale of candidate feature points are determined by detecting local extrema in the DoG pyramid. In one embodiment of the invention, the Gaussian pyramid has 4 layers, and each layer has 5 scales. Then, candidate feature points are precisely located and low-contrast points are removed, retaining stable feature points. Preferably, a contrast threshold of 0.03 and an edge response threshold of 10 are set to remove unstable feature points. Finally, the principal orientation and 128-dimensional feature descriptor of the feature points are calculated.

[0037] Let the first Layer Image Extracted from There are _n_ feature points, and the set of feature points is denoted as _n_i. Each feature point Includes position coordinates ,scale , main direction and 128-dimensional descriptor vector The test results of this invention show that for a typical pathological slide image (approximately 50,000 × 30,000 pixels), approximately 50,000 to 200,000 feature points can be extracted from a single slide.

[0038] In one embodiment of the present invention, an image pyramid downsampling strategy is adopted to improve the efficiency of feature point extraction. First, rapid feature point detection is performed on a low-resolution image (downsampling factor of 8 or 16) to obtain a coarse feature point distribution. Then, fine detection is performed on the original resolution image in areas with dense feature points, ensuring detection accuracy while significantly reducing computational load. Tests show that this strategy can shorten feature point extraction time by approximately 60% while maintaining essentially the same detection quality.

[0039] Furthermore, this invention optimizes the SIFT algorithm parameters to suit the characteristics of pathological sections. Because pathological sections contain numerous repetitive texture structures (such as glandular structures and cell arrangements), they are prone to generating a large number of similar feature points, leading to matching ambiguities. Therefore, this invention increases the standard deviation of the Gaussian weighted window during descriptor generation from the default 1.5 times to 2.0 times, enabling the descriptors to capture a wider range of contextual information and improving the discriminative power of feature points. Simultaneously, for regions with significant variations in staining concentration, a feature descriptor normalization method based on gradient orientation histograms is employed to enhance robustness to brightness changes.

[0040] S2-2, Feature Point Matching and Rigid Transformation Estimation: After extracting feature points, feature point matching is performed to establish the correspondence between adjacent slices. This invention uses the Fast Nearest Neighbor Search (FLANN) algorithm to accelerate the matching process and combines it with a ratio test to screen reliable matching pairs.

[0041] Specifically, for the first Each feature point in the layer image In the Search for the nearest and second nearest neighbors of the descriptor vector in the feature point set of the layer image. Let the nearest neighbor distance be . The next nearest neighbor distance is When the ratio If the ratio is less than a preset threshold (preferably 0.7), the match is considered reliable; otherwise, the match is rejected. An initial set of matching pairs is obtained after the ratio test.

[0042] Since there may be erroneous matches in the initial matching, this invention further employs the Random Sample Consensus (RANSAC) algorithm to remove outliers. The basic idea of ​​RANSAC is to randomly sample the minimum subset from the matched point pairs to estimate the transformation model, then count the number of interior points that conform to the model, and after multiple iterations, select the model with the most interior points as the final result.

[0043] In one embodiment of the present invention, it is assumed that the obtained after screening For a matching point, let it be denoted as The rigid transformation matrix is ​​calculated based on the matching point pairs. The rigid transformation is derived from the rotation matrix. Translation vector The composition and transformation equations are as follows: ,in, Indicates the first The original coordinates of feature points in the layer image. Represents the coordinates after rigid transformation. It is a 2×2 rotation matrix. It is a 2×1 translation vector. The rotation matrix can be parameterized as: ,in, This represents the rotation angle, with a value range of [value missing]. .

[0044] The rigid transformation parameters are solved using RANSAC combined with the least squares method. In one embodiment of this invention, the number of RANSAC iterations is set to 2000, and the interior point threshold is set to 5 pixels. During the RANSAC iteration process, two pairs of matching points are randomly selected each time to solve for the transformation parameters, and then the number of points among all matching points whose transformed residuals are less than the interior point threshold is counted. After the iteration is completed, the model with the most interior points is selected, and the transformation parameters are re-estimated using the least squares method with all interior points.

[0045] Apply the obtained stiffness transformation to the first... Layered images are resampled using bilinear interpolation to obtain rigidly registered images. .

[0046] S2-3, Optical Flow Field Calculation: Even after rigid registration eliminates the overall rotation and translation between slices, local non-rigid deformations may still exist between adjacent slices. This invention uses a dense optical flow algorithm based on variational principles to calculate the optical flow field, which is used to describe pixel-level displacement relationships.

[0047] Optical flow field Represents each pixel in the image The displacement vector, where This represents the horizontal displacement component. This represents the vertical displacement component. This invention employs an improved form of the Horn-Schunck variational model, solving the optical flow field by minimizing the following energy functional: ,in, This is a data item representing the degree of compliance with optical flow constraints. The smoothing term represents the spatial smoothness of the optical flow field. This is a smoothing regularization coefficient used to balance data fidelity and smoothness.

[0048] The data term is defined as the sum of squared residuals under the assumption of constant image grayscale: ,in, For reference image, This is the image to be registered after rigid registration. The above formula represents the grayscale difference between the corresponding position of each pixel in the reference image and the image to be registered after optical flow displacement.

[0049] The smoothing term is defined as the integral of the optical flow field gradient: ,in, and These represent the spatial gradients of the optical flow components. The smoothing term constrains the spatial continuity of the optical flow field, avoiding discontinuous displacement jumps.

[0050] In one embodiment of the present invention, the smoothing regularization coefficient The adjustment is adaptive based on the tissue type. For densely structured, substantial tissue regions, a larger adjustment is used. Use values ​​(e.g., 0.5 to 1.0) to maintain a smooth optical flow field; for loosely structured interstitial regions or regions with significant tissue tearing, use smaller values. Values ​​(such as 0.01 to 0.1) are used to allow for larger local displacements. The adaptive adjustment of values ​​is achieved by analyzing the local texture complexity of the image; regions with higher texture complexity... The smaller the value.

[0051] Preferably, the present invention employs a hierarchical adaptive approach. Adjust the strategy. First, divide the image into multiple non-overlapping local regions (e.g., 128×128 pixel blocks), and calculate the texture complexity index for each region. Texture complexity is defined as the standard deviation of the gray-level gradient magnitude in a local region. Then, based on... The value determines the region value: ,in, The maximum smoothing coefficient is preferably set to 1.0. This is a texture complexity normalization constant, preferably set to the median of the global texture complexity of the image. Through this adaptive mechanism, the optical flow field in smooth regions remains continuous, while more refined displacement variations are allowed in textured regions.

[0052] A pyramid iterative scheme is employed to solve the minimization problem of the energy functional. First, a coarse optical flow is solved on a low-resolution image. Then, the image is upsampled layer by layer to a high-resolution image, refining the optical flow field. In one embodiment of this invention, the number of pyramid layers is set to 5, and the number of iterations per layer is set to 50. The Gauss-Seidel iterative method is used for the solution, and the convergence criterion is that the change in the optical flow field between two adjacent iterations is less than... pixel.

[0053] This invention further performs post-processing on the optical flow field to eliminate outliers. Median filtering is used to filter the two components of the optical flow field separately, with the filter kernel size set to 5×5 pixels, which can effectively remove isolated abnormal displacement points. For displacement vectors whose displacement amplitude exceeds a preset threshold (such as 10 times the slice spacing), they are identified as outliers and replaced with the median of the neighborhood.

[0054] S2-4, Elastic deformation and image resampling, to calculate the optical flow field. Then, it was used as an elastic deformation field to apply to the rigidly registered image. The image after deformation Obtained through the following transformation: Since the transformed coordinates usually do not fall at integer pixel positions, this invention uses bicubic interpolation for image resampling to maintain the clarity of image details.

[0055] Repeat the registration process described above, registering all adjacent slice pairs sequentially, to finally obtain the registered slice sequence. The test results of this invention show that, compared with using rigid registration or flexible registration alone, the dual-strategy registration scheme can reduce the average registration error from 8.2 pixels to 3.5 pixels, an error reduction of approximately 57%.

[0056] Step S3: Semantic segmentation of the tumor region. The purpose of this step is to accurately identify the spatial extent of the tumor region from the registered slice sequence. This invention employs a semantic segmentation method based on convolutional neural networks. Compared to methods based on fuzzy reasoning, deep learning methods can adaptively learn the complex morphological features of tumors, eliminating the need for manually constructing a rule base and exhibiting better generalization ability for different types of tumors.

[0057] In one embodiment of the present invention, the semantic segmentation network adopts an encoder-decoder architecture. The encoder is responsible for extracting multi-scale features of the image and contains 5 coding blocks, each consisting of 2 convolutional layers and 1 max pooling layer. The convolutional layers have a kernel size of 3×3, a stride of 1, and padding of 1. Each convolutional layer is followed by a batch normalization layer and a ReLU activation function. The pooling layers have a kernel size of 2×2 and a stride of 2. The number of output channels for the 5 coding blocks are 64, 128, 256, 512, and 1024, respectively.

[0058] The decoder is responsible for restoring the spatial resolution and outputting pixel-by-pixel classification results, and consists of four decoding blocks. Each decoding block first upsamples the feature map size by a factor of 2 through transposed convolution, then performs a skip connection with the feature map of the corresponding layer in the encoder, and finally fuses the features through two convolutional layers. Skip connections can fuse the fine spatial information of low-level features and the semantic information of high-level features, improving segmentation accuracy. The number of output channels for the four decoding blocks are 512, 256, 128, and 64, respectively.

[0059] The last layer of the network is a 1×1 convolutional layer with 2 output channels (corresponding to the two categories of tumor and non-tumor), followed by a Softmax function to output the probability of each pixel belonging to each category.

[0060] In one embodiment of the present invention, network training employs a combination of a binary cross-entropy loss function and a Dice loss function: ,in, For binary cross-entropy loss, For Dice's loss, and The weighting coefficient is preferably set to... , The introduction of Dice loss can alleviate the class imbalance problem and improve the segmentation performance of tumor regions.

[0061] The training dataset contains ground truth tumor regions annotated by pathologists. Data augmentation methods include random rotation, random flipping, random scaling, and color jittering. Training uses the Adam optimizer, with an initial learning rate set to... The batch size is set to 8, and the training rounds are set to 100.

[0062] During the inference phase, each slice image in the slice sequence is registered. Input a trained segmentation network, and the network outputs a tumor probability map. represents the probability that each pixel belongs to the tumor region. The probability map is converted into a binary mask through thresholding. ,in, For the first Masking of the tumor area in the layer, The segmentation threshold is preferably set to 0.5. To improve the spatial continuity of the mask, this invention further performs morphological post-processing on the mask, including opening operations to remove small noise points and closing operations to fill holes. The opening and closing operations use circular structuring elements with a radius of 5 pixels.

[0063] In one embodiment of the present invention, a sliding window strategy is employed for segmentation inference, taking into account the large size of pathological slide images. The entire slide image is divided into several overlapping image blocks (e.g., 512×512 pixels, with an overlap rate of 50%). Segmentation prediction is performed on each image block separately. Then, the prediction results of the overlapping areas are fused by weighted averaging to eliminate discontinuities at block boundaries. The weights are Gaussian distributed, with higher weights for the central region and lower weights for the edge region. The fusion formula is as follows: ,in, For the first Predicted probability map of image patches, This is the corresponding Gaussian weighted graph.

[0064] Preferably, the present invention also introduces three-dimensional consistency constraints to improve the inter-layer continuity of the segmentation. Specifically, for the segmentation results of three adjacent slices... , and Calculate the first The overlap between the layer segmentation result and the adjacent layers is considered. If the overlap is lower than a preset threshold (e.g., 0.3), then the layer segmentation result is considered lower than the threshold. The segmentation results of the layers are corrected. The correction method is as follows: pixels with a probability value higher than a low threshold (e.g., 0.3) within the union region of adjacent layer masks are included in the first layer. Layer masks are used to enhance the continuity of the three-dimensional structure. This three-dimensional consistency constraint is particularly effective in handling cases where the tumor morphology gradually changes between adjacent layers, and can avoid the breakage of the three-dimensional model caused by single-layer segmentation errors.

[0065] Step S4: Contour interpolation and voxel reconstruction. The purpose of this step is to integrate the tumor region information from the two-dimensional slice sequence into a three-dimensional voxel model. Since the spacing between adjacent slices (4 μm) is much larger than the pixel resolution within a slice (0.25 μm), directly stacking two-dimensional masks will cause obvious step-like artifacts in the Z-direction of the three-dimensional model. This invention uses contour interpolation technology to generate intermediate layer contours between adjacent slices, effectively improving the spatial smoothness of the three-dimensional model.

[0066] S4-1, Contour Extraction: First, a two-dimensional contour curve is extracted from the tumor region mask. This invention employs a boundary-tracking-based contour extraction algorithm for the binary mask. Perform boundary detection to obtain closed contour curves. Each contour curve is represented as an ordered set of points. ,in This represents the number of contour points. If there are multiple connected regions in the mask, the contours of each region are extracted separately.

[0067] S4-2, Contour Correspondence and Interpolation: To interpolate between the contours of adjacent slices, it is necessary to first establish the correspondence between contour points. This invention employs a contour correspondence method based on shape context.

[0068] Shape context descriptors characterize the shape features of a point by statistically analyzing the relative positional distribution of other points around the contour point. For a contour... Points on Its shape context descriptor Defined as a histogram of the positions of other contour points in a logarithmic polar coordinate system with the origin at that point. In one embodiment of the present invention, the polar radius direction is divided into 5 intervals, the polar angle direction is divided into 12 intervals, and the shape context descriptor has a dimension of 60.

[0069] Establish the outline of adjacent slices and When determining the point correspondence, the shape context distance of all point pairs is calculated, a cost matrix is ​​constructed, and then the optimal correspondence is solved using the Hungarian algorithm. For cases with unequal numbers of points, one-to-many or many-to-one correspondences are allowed.

[0070] After establishing the correspondence, linear interpolation or thin-plate spline interpolation is used to generate the intermediate layer profile. Assume adjacent slices... and Insertion is required between them The middle outline of the layer, then the first The coordinates of the interpolation points of the middle contour of the layer are: ,in, In one embodiment of the present invention, The value of is determined based on the desired Z-direction resolution, preferably making the Z-direction resolution close to the XY-direction resolution. For example, if the slice spacing is 4μm and the pixel resolution is 0.25μm, then Setting it to 15 makes the effective resolution in the Z direction approximately 4 / 16 = 0.25 μm.

[0071] For contour pairs with significant morphological changes, this invention employs a nonlinear interpolation method based on thin plate splines, which can produce smoother intermediate contour transitions.

[0072] S4-3, Voxelization processing, converting the interpolated contour sequence into a three-dimensional voxel model. In one embodiment of the present invention, the size and resolution of the three-dimensional voxel mesh are first determined. Let the XY plane dimensions of the slice image be... The total number of layers after contour interpolation is pixels. If the layer is 3, then the size of the three-dimensional voxel mesh is voxel resolution is ,in For pixel resolution, This represents the Z-axis resolution after interpolation.

[0073] For each contour layer, a scanline filling algorithm is used to convert the interior region of the contour into voxels. The filled binary voxel layers are stacked to form a three-dimensional voxel model. ,in This indicates that the voxel is located within the tumor region. This indicates that the voxel is located outside the tumor region.

[0074] To further improve the surface smoothness of the 3D model, this invention performs 3D Gaussian smoothing on the voxel model, and then extracts isosurfaces using the Marching Cubes algorithm to generate a 3D surface mesh model of the tumor. The kernel size of the Gaussian smoothing is set to 3×3×3 voxels, and the standard deviation is set to 1 voxel. The isosurface threshold of the Marching Cubes algorithm is set to 0.5.

[0075] In one embodiment of the present invention, for tumors with complex morphology (such as those with multiple isolated lesions), connectivity analysis is required during voxelization. A three-dimensional 26-neighborhood connectivity criterion is used to divide the voxel model into several connected components, each corresponding to an independent tumor lesion. Tiny connected components with a volume smaller than a preset threshold (e.g., 1% of the total number of primes) are identified as segmentation noise and removed, retaining the main tumor entities.

[0076] Preferably, the present invention further performs mesh optimization processing on the three-dimensional surface mesh model to improve the visualization effect. Mesh optimization includes the following sub-steps: First, the mesh surface is smoothed using the Laplacian smoothing algorithm, with the number of iterations set to 5 and the smoothing factor set to 0.5. This step can eliminate the stair-step artifacts generated by the Marching Cubes algorithm; then, a mesh simplification algorithm (such as the Quadric Error Metrics method) is used to reduce the number of mesh faces, reducing the number of faces to 30% to 50% of the original while maintaining the morphological characteristics, so as to improve the performance of subsequent rendering and interaction; finally, the normal vectors of the simplified mesh are recalculated to ensure the correctness of the surface lighting effect.

[0077] Step S5: Calculation of three-dimensional morphological parameters. The purpose of this step is to quantitatively calculate the spatial morphological characteristics of the tumor based on a three-dimensional voxel model, providing an objective basis for pathological staging and prognostic assessment. The three-dimensional morphological parameters calculated in this invention include tumor volume, tumor surface area, maximum tumor diameter, depth of invasion, and shape factor.

[0078] Tumor volume is obtained by counting the number of voxels with a value of 1 in the three-dimensional voxel model and multiplying the result by the physical volume of a single voxel. ,in, The unit is mm³ (micrometers need to be converted to millimeters). In one embodiment of the present invention, for typical pathological slide data (slide size 50000×30000 pixels, slide spacing 4μm, 50 consecutive layers), the voxel resolution is 0.25×0.25×0.25μm. 3 The volume of a single voxel is μm 3 .

[0079] The tumor surface area is obtained by calculating the total area of ​​the three-dimensional surface mesh model. Assume the surface mesh consists of... It consists of several triangular facets, and the coordinates of the vertices of each triangular facet are: The surface area is: ,in, Represents the cross product of vectors. This represents the magnitude of the vector.

[0080] The maximum diameter of a tumor is defined as the maximum distance between any two points on the tumor surface. This invention employs a rotating caliper algorithm to efficiently calculate the farthest pair of points in a 3D point set. First, the convex hull of the tumor surface mesh is calculated. Then, the rotating caliper algorithm is applied to the vertex set of the convex hull to find the farthest pair of points; the distance between these two points is the maximum diameter of the tumor. .

[0081] Invasion depth is an important indicator for tumor pathological staging, defined as the maximum depth to which tumor cells infiltrate normal tissue. The invasion depth calculation method of this invention includes the following sub-steps:

[0082] First, the location of the basement membrane is identified from the three-dimensional voxel model. The basement membrane, as the interface between epithelial and stromal tissues, can be determined by identifying the boundary line between the epithelial and stromal regions in the semantic segmentation results of each slice; this boundary line constitutes the basement membrane surface in three-dimensional space. Preferably, this invention trains a dedicated basement membrane segmentation network, or adds a basement membrane category to the tumor segmentation network.

[0083] Then, calculate the shortest distance from each point on the tumor boundary to the basement membrane surface. Let a point on the tumor boundary... The coordinates are The basement membrane surface is Then the distance from that point to the basement membrane is: This invention employs a KD-tree data structure to accelerate nearest-point lookup. Finally, the maximum distance from all tumor boundary points to the basement membrane is taken as the invasion depth. ,in, It is the surface of the tumor.

[0084] Shape factors are used to describe the overall shape characteristics of tumors. Commonly used shape factors include sphericity, compactness, and surface area to volume ratio.

[0085] Sphericity is defined as the ratio of the surface area of ​​a sphere of equal volume to the tumor to the actual surface area of ​​the tumor. The sphericity value ranges from (0, 1). The closer the value is to 1, the closer the tumor shape is to spherical. The smaller the value is, the more irregular the tumor shape is.

[0086] Step S6: Heterogeneous Partitioning and Visualization. The innovative purpose of this step is to reveal the spatial heterogeneity within the tumor, identify different functional partitions, and provide three-dimensional pathological information for precision medicine. As a heterogeneous population, tumors often contain different subregions, including rapidly proliferating tumor cell-rich areas, necrotic areas due to hypoxia, and infiltration front areas that invade surrounding tissues. These functional partitions have different biological behaviors and prognostic significance.

[0087] S6-1, Voxel Feature Extraction. First, a multi-dimensional feature vector is extracted for each voxel in the 3D voxel model. The features extracted in this invention include the following categories:

[0088] Gray-scale features include the gray-scale value at the corresponding voxel position, the mean gray-scale value and standard deviation within the local window.

[0089] Texture features: Texture parameters calculated based on the gray-level co-occurrence matrix, including energy, local texture entropy, contrast, and homogeneity. In one embodiment of the present invention, the gray-level co-occurrence matrix is ​​calculated within an 11×11×11 neighborhood window around the voxel, and the texture parameters in 13 directions are statistically analyzed and averaged.

[0090] Cell density features: Cell density estimates for local regions are obtained through cell detection algorithms. This invention trains a cell detection model based on a fully convolutional network. The model takes a sliced ​​image region as input and outputs the cell center point detection results, from which the number of cells per unit area is calculated as the estimated cell density.

[0091] Proliferative activity characteristics: The proliferative activity of local areas is estimated using a Ki-67 positivity rate prediction model. If the slide contains Ki-67 immunohistochemical staining results, the proportion of positive cells can be calculated directly; if only H&E stained slides are available, this invention trains a regression network to predict the Ki-67 positivity rate from H&E images.

[0092] Based on the above characteristics, each voxel is represented as a single... 3D feature vector ,in As a feature dimension, in one embodiment of the present invention .

[0093] S6-2, Cluster Analysis and Heterogeneous Partitioning: Cluster analysis based on voxel feature vectors is used to divide the tumor interior into different functional partitions. This invention preferably employs density-based spatial clustering (DBSCAN) or Gaussian mixture model (GMM) clustering algorithms.

[0094] Taking the DBSCAN algorithm as an example, this algorithm can identify clusters of arbitrary shapes without requiring a pre-specified number of clusters. The core parameters of the algorithm include the neighborhood radius. and minimum points In one embodiment of the present invention, Set to the median of the Euclidean distance in the feature space. Set as ,in For feature dimensions.

[0095] The clustering results divide tumor voxels into several clusters, each cluster corresponding to a heterogeneous partition. This invention performs semantic annotation on the partitions based on the characteristic statistics of each cluster:

[0096] Necrotic areas are characterized by low cell density, high texture entropy, and low grayscale values, corresponding to regions within the tumor that have undergone necrosis due to hypoxia. Necrotic areas are typically located in the central part of the tumor and exhibit morphological features such as structural disorder, fragmented or dissolved nuclei. This invention identifies necrotic areas by detecting regions with cell density below a threshold (e.g., 50 cells / mm²) and texture entropy above a threshold (e.g., 4.0).

[0097] High-proliferative regions are characterized by high cell density, high Ki-67 predictive values, and high texture energy, corresponding to areas of rapid tumor cell division and proliferation. Cells in high-proliferative regions are typically densely packed, with large, deeply stained nuclei and numerous mitotic figures. This invention identifies high-proliferative regions by detecting areas with a Ki-67 predictive positivity rate higher than a threshold (e.g., 30%) and a cell density higher than a threshold (e.g., 500 cells / mm²).

[0098] Invasion front zone: Characterized by moderate cell density, high texture contrast, and location near the tumor boundary, corresponding to the area where the tumor invades surrounding normal tissue. Tumor cells in the invasion front zone exhibit a pattern of single-cell or small cell nest infiltration, interspersed with surrounding normal tissue, demonstrating high spatial heterogeneity. This invention identifies the invasion front zone by detecting areas located within a certain distance (e.g., 200 μm) from the tumor boundary with texture contrast exceeding a threshold.

[0099] Preferably, the present invention also supports user-defined heterogeneity partitioning criteria. Pathologists can adjust clustering parameters and partitioning thresholds based on the biological characteristics of specific tumor types, or introduce additional biomarker prediction results (such as HER2 expression prediction, PD-L1 expression prediction, etc.) to refine heterogeneity partitioning. This flexible configuration mechanism enables the system to adapt to the needs of different clinical application scenarios.

[0100] S6-3, 3D Visualization and Report Generation: Finally, a 3D visualization model of the tumor and a report on its spatial distribution characteristics are generated. The 3D visualization model uses different colors to identify different heterogeneous partitions. Users can interactively rotate, zoom, and cut the model to observe the spatial structure of the tumor from any perspective. In one embodiment of this invention, the VTK visualization toolkit is used to achieve 3D rendering.

[0101] The spatial distribution characteristics report includes the following: basic tumor information (volume, surface area, maximum diameter, depth of invasion), heterogeneity partition statistics (volume percentage of each partition, spatial distribution map), a summary table of three-dimensional morphological parameters, and representative cross-sectional images of the tumor. The report is output in PDF format for pathologists' reference.

[0102] This invention also provides a registration quality feedback optimization mechanism, which guides the adjustment of the pre-registration parameters through post-evaluation indicators to form a closed-loop optimization.

[0103] Registration quality score Defined as the structural similarity index (SSIM) between adjacent slices after registration: SSIM is the structural similarity index, with a value range of [0, 1]. The closer the value is to 1, the better the registration quality.

[0104] when When the value falls below a preset threshold (preferably 0.85), a parameter adjustment mechanism is triggered. Parameter adjustment strategies include: increasing the upper limit of the number of feature point detections, decreasing the ratio test threshold to obtain more matching points, and reducing the optical flow field smoothing coefficient to allow for larger local displacements. After adjusting the parameters, registration step S2 is re-executed, iteratively optimizing until... The requirements are met or the maximum number of iterations is reached.

[0105] Figure 2 This diagram illustrates the architecture of a continuous-slice 3D reconstruction system for pathological slides based on image registration, provided in an embodiment of the present invention. This system is used to implement the methods described in the foregoing method embodiments and includes the following modules:

[0106] The slide acquisition module is configured to digitally scan sequential paraffin sections of tumor tissue to obtain a sequence of digital images of pathological slides, and to perform color normalization processing on the sequence of digital images of pathological slides. In one embodiment of the present invention, the slide acquisition module includes a full-slide scanner interface submodule and a color normalization submodule. The full-slide scanner interface submodule supports interface with mainstream full-slide scanners (such as Aperio, Hamamatsu, Leica, etc.) and can read digital slide files in formats such as SVS, NDPI, and SCN. The color normalization submodule implements the color deconvolution normalization algorithm described in the method embodiment.

[0107] Preferably, the slice acquisition module further includes an image quality assessment submodule, used to automatically detect the sharpness, staining quality, and tissue integrity of the slice images, and generate warning prompts for slice images that do not meet the quality standards. The image quality assessment submodule assesses sharpness by calculating the Brenner gradient function value of the image, assesses staining quality by analyzing the intensity distribution of H&E staining channels, and assesses tissue integrity by detecting the connectivity of tissue regions.

[0108] The dual-strategy registration module is configured to perform inter-slice registration on the color-normalized digital image sequence of pathological slides using a combination of rigid registration of feature points and elastic registration of optical flow fields, thereby obtaining a registered slide sequence. In one embodiment of the present invention, the dual-strategy registration module includes a feature point extraction submodule, a feature point matching submodule, a rigid transformation estimation submodule, an optical flow field calculation submodule, and an image resampling submodule. The functions of each submodule are consistent with the descriptions of the corresponding steps in the method embodiment.

[0109] The dual-strategy registration module also includes a registration quality assessment submodule, which calculates the registration quality score and triggers parameter adjustments. This submodule comprehensively evaluates the registration effect by calculating multiple indicators, such as the structural similarity index of adjacent slices after registration, normalized mutual information, and consistency of overlapping regions. When the registration quality score falls below a threshold, this submodule automatically adjusts the registration parameters and triggers re-registration, achieving closed-loop optimization of registration quality.

[0110] The semantic segmentation module is configured to use a convolutional neural network to perform semantic segmentation of the tumor region on the registered slice sequence, generating a tumor region mask. In one embodiment of the present invention, the semantic segmentation module includes a deep learning inference engine, a segmentation network model, and a post-processing submodule. The deep learning inference engine supports GPU-accelerated inference and is compatible with model formats such as ONNX and TensorRT. The segmentation network model is a pre-trained encoder-decoder network, which can load corresponding model weights according to different tumor types. The post-processing submodule implements thresholding and morphological filtering.

[0111] Preferably, the semantic segmentation module supports switching between segmentation models for multiple tumor types. The system has built-in pre-trained segmentation models for common tumor types (such as colorectal cancer, breast cancer, lung cancer, gastric cancer, etc.), and users can select the corresponding model based on actual cases. In addition, the system provides a model fine-tuning function, allowing users to fine-tune the pre-trained model using the institution's labeled data to adapt it to the institution's staining style and pathological diagnostic standards.

[0112] The voxel reconstruction module is configured to extract two-dimensional contour curves from the tumor region mask and perform inter-layer contour interpolation, converting the interpolated contour sequence into a three-dimensional voxel model. In one embodiment of the present invention, the voxel reconstruction module includes a contour extraction submodule, a contour correspondence submodule, a contour interpolation submodule, and a voxelization submodule. The contour correspondence submodule implements a contour point matching algorithm based on shape context. The voxelization submodule implements scan line filling and Marching Cubes surface reconstruction.

[0113] The voxel reconstruction module also includes a 3D model optimization submodule, used for post-processing optimization of the generated 3D voxel model and surface mesh model. Optimization includes removing isolated noise points, filling micro-holes, smoothing surface meshes, and simplifying the number of mesh patches to improve the visual effect of the 3D model and the accuracy of subsequent analysis.

[0114] The morphological analysis module is configured to calculate the three-dimensional morphological parameters of tumors based on a three-dimensional voxel model. In one embodiment of the present invention, the morphological analysis module includes a volume calculation submodule, a surface area calculation submodule, a maximum diameter calculation submodule, an invasion depth calculation submodule, and a shape factor calculation submodule. The implementation of each submodule is consistent with the description of the corresponding steps in the method embodiment. The morphological analysis module outputs structured morphological parameter data for report generation.

[0115] Preferably, the morphological analysis module also supports the calculation of custom morphological parameters. Users can define new morphological indicators and their calculation formulas through configuration files, and the system automatically calculates and outputs the custom parameters based on the configuration. This scalable design enables the system to adapt to the needs of different clinical studies and supports the exploration of new morphological biomarkers.

[0116] The heterogeneity identification module is configured to perform cluster analysis based on a three-dimensional voxel model to identify heterogeneous partitioning within the tumor and generate a three-dimensional visualization model of the tumor and a report on its spatial distribution characteristics. In one embodiment of the invention, the heterogeneity identification module includes a feature extraction submodule, a cluster analysis submodule, a partitioning annotation submodule, a three-dimensional rendering submodule, and a report generation submodule. The three-dimensional rendering submodule uses VTK to achieve interactive three-dimensional visualization. The report generation submodule integrates the analysis results into a structured PDF report.

[0117] The 3D rendering submodule offers rich interactive features, including: model rotation, scaling, and translation; arbitrary angle cross-section display; independent display and hiding of heterogeneous partitions; transparency adjustment; distance and angle measurement tools, etc. Users can observe the 3D morphology and internal structure of tumors from any perspective using mouse and keyboard, obtaining intuitive spatial pathological information.

[0118] The report generation submodule automatically generates an analysis report based on a preset report template. The report includes: basic case information, slide processing parameters, 3D reconstruction quality indicators, a summary table of tumor morphology parameters, heterogeneity partition statistics and distribution maps, representative 3D view screenshots, and an area for pathologist-editable diagnostic opinions. The report is output in PDF format for easy archiving and sharing.

[0119] This invention's system can be deployed on workstations or cloud servers in hospital pathology departments, providing services to users through a web interface or client application. The system supports batch processing of multiple cases and provides database storage and query functions, facilitating the management and retrospective analysis of historical data.

[0120] In one embodiment of the present invention, the system employs a distributed computing architecture to support the processing of large-scale pathological slide data. The system includes a task scheduling server, multiple computing nodes, and a shared storage server. The task scheduling server is responsible for receiving analysis tasks submitted by users and allocating them to idle computing nodes. The computing nodes execute specific image processing and analysis algorithms, and the shared storage server provides high-speed data access services. By horizontally scaling the number of computing nodes, the system can linearly increase processing throughput, meeting the batch processing needs of large medical institutions.

[0121] The efficacy of this invention was validated on a test set containing 120 tumor samples. The test set included 30 colorectal cancer, 30 breast cancer, 30 lung cancer, and 30 gastric cancer samples, with each sample containing 20 to 80 serial sections.

[0122] Regarding registration accuracy, the dual-strategy registration scheme of this invention reduces the average registration error from 7.8 pixels to 3.2 pixels compared to the affine registration method; compared to using the elastic registration method alone, the average registration error is reduced from 5.6 pixels to 3.2 pixels, and the registration accuracy is improved by approximately 43%.

[0123] In terms of segmentation accuracy, on a test set annotated by three pathology experts, the semantic segmentation method of this invention achieved an average Dice coefficient of 0.923 and an average IoU of 0.861, which is better than the traditional threshold-based segmentation method (Dice 0.812, IoU 0.724).

[0124] Regarding the three-dimensional reconstruction effect, according to the evaluation of pathology experts, the three-dimensional tumor model reconstructed by this invention can accurately reflect the spatial morphological characteristics of the tumor and has good consistency with the actual tumor morphology seen during surgery.

[0125] Regarding heterogeneous partitioning, the overlap rate between the necrotic areas identified by this invention and the necrotic areas labeled by immunohistochemical Ki-67 / HIF-1α joint staining reached 0.87, and the overlap rate between the high-proliferation areas and the Ki-67 high-expression areas reached 0.82, demonstrating the biological validity of the heterogeneous partitioning results.

[0126] The above description is merely a preferred embodiment of the present invention and does not limit the patent scope of the present invention. Any equivalent structural transformations made based on the inventive concept of the present invention and the contents of the specification and drawings of the present invention, or direct or indirect applications in other related technical fields, are similarly included within the patent protection scope of the present invention.

Claims

1. A method for continuous-layer three-dimensional reconstruction of pathological sections based on image registration, characterized in that, Includes the following steps: S1. Digital acquisition and preprocessing steps of slides: Digitally scan the continuous paraffin slides of tumor tissue to obtain the digital image sequence of pathological slides, and perform color normalization processing on the digital image sequence of pathological slides to eliminate color shifts caused by differences in staining batches. S2. Dual-strategy registration steps: A strategy combining rigid registration of feature points and elastic registration of optical flow fields is used to perform inter-slice registration on the color-normalized digital image sequence of pathological slides. First, feature point sets are extracted from adjacent slide images and feature point matching is performed. Based on the matched feature point pairs, a rigid transformation matrix is ​​calculated to correct the overall rotation and translation deviations between slides. Then, an optical flow field is calculated based on the rigidly registered image pairs. An elastic deformation field is generated based on the optical flow field to compensate for local non-rigid deformations, and finally, the registered slide sequence is obtained. S3. Semantic segmentation of tumor region: A convolutional neural network is used to perform semantic segmentation of the tumor region on each layer of slice image in the registered slice sequence, outputting a tumor probability map and generating a tumor region mask through thresholding, wherein the tumor region mask identifies the spatial boundary between tumor tissue and normal tissue. S4. Contour Interpolation and Voxelization Reconstruction Steps: Extract two-dimensional contour curves from the tumor region mask, perform shape interpolation on the two-dimensional contour curves between adjacent slice layers to generate intermediate layer contours, and convert the interpolated contour sequence into a three-dimensional voxel model. S5. Calculation steps for three-dimensional morphological parameters: Calculate the three-dimensional morphological parameters of the tumor based on the three-dimensional voxel model. The three-dimensional morphological parameters include tumor volume, tumor surface area, maximum diameter of the tumor, depth of invasion, and boundary relationship parameters between the tumor and surrounding normal tissues. S6. Heterogeneous Partitioning and Visualization Steps: Based on the texture features and spatial distribution features of each voxel in the three-dimensional voxel model, cluster analysis is performed to identify the heterogeneous partitioning results of the necrotic area, high-proliferation area and invasion front area inside the tumor, and a three-dimensional visualization model of the tumor and a report on its spatial distribution features are generated.

2. The method according to claim 1, characterized in that, In step S1, the color normalization process uses a staining separation method based on color deconvolution to decompose the digital image of the pathological slide into hematoxylin channels and eosin channels, and normalizes the staining intensity of each channel to a preset standard staining intensity range.

3. The method according to claim 1, characterized in that, In step S2, the feature point set is extracted using the scale-invariant feature transformation algorithm, the feature point matching is screened using fast nearest neighbor search combined with ratio test, and the rigid transformation matrix is ​​solved using the least squares estimation after removing erroneous matching point pairs using the random sampling consensus algorithm.

4. The method according to claim 1, characterized in that, In step S2, the optical flow field is calculated using a dense optical flow algorithm based on the variational principle. The optical flow field is solved by minimizing the energy functional containing data terms and smoothing terms, wherein the regularization coefficient of the smoothing term is adaptively adjusted according to the tissue type.

5. The method according to claim 1, characterized in that, In step S3, the convolutional neural network adopts an encoder-decoder architecture. The encoder contains multiple convolutional layers and pooling layers to extract multi-scale features, and the decoder restores spatial resolution through upsampling and skip connections. The network output is a pixel-wise classification result with the same size as the input image.

6. The method according to claim 1, characterized in that, In step S4, the contour interpolation uses a contour correspondence method based on shape context to establish the correspondence between adjacent slice contour points, and generates the intermediate layer contour through linear interpolation or nonlinear interpolation based on thin plate splines.

7. The method according to claim 1, characterized in that, In step S5, the method for calculating the invasion depth includes: first, identifying the tumor boundary and basement membrane positions from the three-dimensional voxel model; then, calculating the shortest distance from each point on the tumor boundary to the basement membrane; and taking the maximum value of the shortest distance as the invasion depth.

8. The method according to claim 1, characterized in that, In step S6, the cluster analysis of heterogeneous partitions adopts a density-based spatial clustering algorithm. The clustering features include voxel gray values, local texture entropy, estimated cell density, and predicted Ki-67 positivity rate.

9. The method according to claim 1, characterized in that, It also includes a registration quality feedback optimization step: calculate the registration quality score of the registration slice sequence, and when the registration quality score is lower than a preset threshold, adjust the feature point detection parameters or the optical flow field smoothing coefficient, and re-execute step S2.

10. A three-dimensional reconstruction system for continuous slices of pathological tissue based on image registration, used to implement the method of any one of claims 1 to 9, characterized in that, include: The slide acquisition module is configured to digitally scan continuous paraffin sections of tumor tissue to obtain digital image sequences of pathological slides, and to perform color normalization processing on the digital image sequences of pathological slides. The dual-strategy registration module is configured to perform inter-layer registration on the color-normalized digital image sequence of pathological slides by combining a strategy of rigid registration of feature points and elastic registration of optical flow field, so as to obtain a registered slide sequence. The semantic segmentation module is configured to use a convolutional neural network to perform semantic segmentation of the tumor region on the registered slice sequence and generate a tumor region mask. The voxel reconstruction module is configured to extract two-dimensional contour curves from the tumor region mask and perform interlayer contour interpolation, converting the interpolated contour sequence into a three-dimensional voxel model. The morphological analysis module is configured to calculate the three-dimensional morphological parameters of tumors based on a three-dimensional voxel model. The heterogeneity identification module is configured to perform cluster analysis based on a three-dimensional voxel model to identify the heterogeneous partitioning results within the tumor and generate a three-dimensional visualization model of the tumor and a report on its spatial distribution characteristics.

Images