A neural CT atlas segmentation method based on high-precision registration

Abstract

The present invention relates to the field of medical imaging technology, and in particular to a neural CT atlas segmentation method based on high-precision registration. The method comprises the following steps: obtaining a neural CT image; performing a coarse registration based on mutual information according to the CT image, and performing a multi-scale fusion based on wavelet to obtain a fused image; performing standard atlas preprocessing on a preset standard neural CT atlas library according to the fused image to obtain a preprocessed standard atlas set; adding a multi-scale feature fusion module to the pre-built basic U-Net architecture, and embedding a spatial attention module to obtain a U-Net network to be trained; performing deep registration network training on the U-Net network to be trained using the fused image and the preprocessed standard atlas set to obtain a trained U-Net network. The present invention achieves a higher precision and more robust segmentation effect through the techniques of multimodal image fusion, deep learning registration, and conditional random field optimization.

Description

Technical Field

[0001] The present invention relates to the field of medical imaging technology, and in particular to a neural CT atlas segmentation method based on high-precision registration. Background Art

[0002] NeuroCT atlas segmentation is an important research direction in the field of medical image processing. Its goal is to automatically identify and segment different tissue structures (such as gray matter, white matter, cerebrospinal fluid, etc.) in neuroCT images. Traditional CT atlas segmentation methods mainly rely on the grayscale information of the image, such as threshold segmentation, region growing and other methods. However, due to the limited grayscale contrast of CT images, the boundaries between different tissues are often blurred, resulting in low segmentation accuracy. In order to improve the segmentation accuracy, many improved methods have been proposed, among which atlas-based segmentation methods have been a hot topic in recent years. The core idea of ​​this type of method is to align the pre-constructed standard atlas to the target CT image, and then use the label information of the standard atlas to guide the segmentation of the target image.

[0003] However, the existing neural CT atlas segmentation methods have the disadvantages of low registration accuracy and lack of spatial consistency. The specific problems are as follows:

[0004] 1. Low registration accuracy: The key to the atlas-based segmentation method lies in the accuracy of registration. Traditional registration methods, such as mutual information-based registration and feature point-based registration, face many challenges when processing neuro-CT images. Due to large individual differences, there are obvious morphological differences in the brain structures of different individuals; in addition, the grayscale distribution of CT images is uneven and the noise is relatively large. These factors will affect the accuracy of registration. Insufficient registration accuracy will cause the standard atlas to be unable to accurately align to the target image, thus affecting the final segmentation result.

[0005] 2. Lack of spatial consistency: Some segmentation methods, such as those based solely on pixel classification, do not fully consider the spatial relationship between pixels, which can easily lead to isolated noise points or areas in the segmentation results and lack of spatial consistency. This makes the segmentation results not smooth and difficult to accurately reflect the real tissue structure. Summary of the invention

[0006] Based on this, it is necessary to provide a neural CT atlas segmentation method based on high-precision registration to solve at least one of the above technical problems.

[0007] To achieve the above purpose, a neural CT atlas segmentation method based on high-precision registration includes the following steps:

[0008] Step S1: Obtaining a neural CT image, an MRI-T1 image, and an MRI-T2 image; performing a coarse registration based on mutual information on the CT image, the MRI-T1 image, and the MRI-T2 image, and performing a multi-scale fusion based on wavelet to obtain a fused image;

[0009] Step S2: preprocessing the preset standard neural CT atlas library with standard atlases according to the fused image to obtain a preprocessed standard atlas set; adding a multi-scale feature fusion module to the pre-built basic U-Net architecture, and embedding the spatial attention module to obtain a U-Net network to be trained; using the fused image and the preprocessed standard atlas set to perform deep registration network training on the U-Net network to be trained to obtain a trained U-Net network; using the trained U-Net network to generate a precise registration field for the fused image to obtain a precise registration field;

[0010] Step S3: Calculate the thin plate spline deformation parameters according to the fused image and the precise registration field to obtain the thin plate spline deformation parameters; use the thin plate spline deformation parameters to perform atlas deformation migration on the preprocessed standard atlas set to obtain a deformation atlas; calculate the prior probability according to the deformation atlas to obtain a prior probability map; calculate the posterior probability according to the prior probability map to obtain a tissue probability map;

[0011] Step S4: construct a univariate potential function for the tissue probability map to obtain a univariate potential function; extract local features of the fused image to obtain a feature vector; define a binary potential function based on the fused image and the feature vector to obtain a binary potential function; construct a conditional random field energy function for the univariate potential function and the binary potential function to obtain a conditional random field energy function; generate a spectrum segmentation result based on the conditional random field energy function to obtain a spectrum segmentation result to achieve the neural CT spectrum segmentation task.

[0012] The present invention integrates the skeletal structure information of CT images and the soft tissue detail information of MRI images to generate a fused image with richer information, effectively improves the image quality and information richness, provides a more reliable data basis for subsequent registration and segmentation steps, and solves the problem of insufficient single-modal image information. High-precision nonlinear registration between fused images and standard atlases is achieved by deep learning technology. By constructing a U-Net network containing a multi-scale feature fusion module and a spatial attention module, and adopting data enhancement and joint loss function optimization strategies, a more accurate deformation field is learned, thereby more accurately describing the spatial correspondence between the fused image and the standard atlas, overcoming the limitation that traditional registration methods are difficult to handle complex deformations. Automatic deformation of the atlas and generation of tissue probability maps are achieved. The atlas is deformed using a thin plate spline interpolation method based on local structural information to better adapt to individual differences, and combined with the grayscale information of the fused image and the label information of the deformation atlas, a tissue probability map is constructed through a Gaussian mixture model and a Bayesian framework, providing a more accurate estimate of tissue category probability, laying the foundation for subsequent segmentation optimization. The conditional random field model is used to optimize the tissue probability map. By combining the unary potential function and the binary potential function based on local features, the spatial relationship and feature similarity between voxels are taken into account, which effectively improves the accuracy and robustness of the segmentation, solves the boundary blur problem caused by the similar grayscale values ​​of adjacent organs, and finally generates a more accurate and smoother segmentation result. Therefore, the present invention provides a neural CT atlas segmentation method based on high-precision registration, which comprehensively utilizes the technologies of multimodal image fusion, deep learning registration and conditional random field optimization, effectively solves the disadvantages of low registration accuracy and lack of spatial consistency in the existing neural CT atlas segmentation methods, and achieves a higher accuracy and more robust segmentation effect. Through these improvements, the present invention can better cope with the challenges of large individual differences and low grayscale contrast in neural CT images, thereby providing more accurate segmentation results and providing more reliable support for clinical applications.

[0013] Preferably, step S1 comprises the following steps:

[0014] Step S11: acquiring a neural CT image, an MRI-T1 image, and an MRI-T2 image; performing image modality preprocessing on the CT image, the MRI-T1 image, and the MRI-T2 image to obtain a preprocessed CT image, a preprocessed MRI-T1 image, and a preprocessed MRI-T2 image;

[0015] Step S12: using the pre-processed CT image as a reference image, performing coarse registration based on mutual information on the pre-processed MRI-T1 image and the pre-processed MRI-T2 image to obtain a registered MRI-T1 image and a registered MRI-T2 image;

[0016] Step S13: performing wavelet-based multi-scale fusion on the pre-processed CT image, the registered MRI-T1 image, and the registered MRI-T2 image to obtain a fused image.

[0017] The present invention effectively improves the image quality and reduces the difference between different modal images by performing preprocessing operations such as denoising, bias field correction and grayscale normalization on CT, MRI-T1 and MRI-T2 images, laying a good foundation for subsequent registration and fusion steps. The denoising operation can remove the noise in the image and improve the signal-to-noise ratio of the image; the bias field correction can eliminate the intensity unevenness caused by the magnetic field inhomogeneity in the MRI image, making the image more uniform; the grayscale normalization can unify the grayscale value range of different modal images to the same interval, which is convenient for subsequent fusion operations. The preprocessed MRI-T1 and MRI-T2 images are registered to the space of the CT image, and the preliminary alignment between the different modal images is achieved, providing spatial consistency for subsequent multi-scale fusion. The registration method based on mutual information can effectively handle the situation where the grayscale difference between different modal images is large. By maximizing the mutual information to find the optimal affine transformation parameters, the images of different modalities are aligned as much as possible in space. Through the wavelet-based multi-scale fusion method, the image information of the three modalities of CT, MRI-T1 and MRI-T2 is effectively fused into one image, making full use of the complementary information of different modal images. Wavelet transform can decompose the image into sub-bands of different scales and directions. Fusion at different scales can better preserve the image's detail information and global structure information. By using fusion rules such as weighted average and local variance maximum, the fusion coefficient can be adaptively selected according to the characteristics of different modal images, thereby obtaining a more information-rich and reliable fused image. Among them, MRI-T1 and MRI-T2 are two different weighted imaging methods of magnetic resonance imaging. The main difference between them lies in the different degrees of emphasis on signals of different tissues. In MRI-T1 weighted imaging (T1WI), T1 stands for longitudinal relaxation time. On T1WI, fat tissue usually presents high signal (bright), while tissue with high water content (such as cerebrospinal fluid) presents low signal (dark). T1WI is mainly used to display anatomical structures; in MRI-T2 weighted imaging (T2WI), T2 stands for transverse relaxation time. On T2WI, tissue with high water content (such as cerebrospinal fluid) presents high signal (bright), while fat tissue usually presents medium or low signal (relatively dark). T2WI is mainly used to display pathological changes.

[0018] Preferably, step S2 comprises the following steps:

[0019] Step S21: performing standard atlas preprocessing on a preset standard neural CT atlas library according to the fusion image to obtain a preprocessed standard atlas set;

[0020] Step S22: adding a multi-scale feature fusion module to the pre-built basic U-Net architecture and embedding a spatial attention module to obtain a U-Net network to be trained;

[0021] Step S23: using the fused image and the preprocessed standard atlas set to perform deep registration network training on the U-Net network to be trained, to obtain a trained U-Net network;

[0022] Step S24: Using the trained U-Net network to generate a precise registration field for the fused image to obtain a precise registration field.

[0023] The present invention performs preprocessing operations such as resampling and grayscale normalization on the atlas in the standard neural CT atlas library, thereby ensuring the consistency of the atlas and the fused image in voxel size, dimension and grayscale range, providing standardized input data for subsequent deep registration network training, and improving the efficiency and stability of training. By adding a multi-scale feature fusion module and a spatial attention module to the basic U-Net architecture, the feature expression ability of the network and the attention to key anatomical structures are enhanced. The multi-scale feature fusion module can combine feature information of different scales to better capture the global and local features of the image; the spatial attention module can guide the network to pay more attention to important anatomical structure areas and improve the registration accuracy. The U-Net network is trained using fused images and preprocessed standard atlas sets, so that the network learns the nonlinear mapping relationship between the fused image and the standard atlas, and can generate a high-precision deformation field. The application of data enhancement strategies increases the diversity of training data and improves the robustness and generalization ability of the model. The design of the joint loss function balances the registration accuracy and the smoothness of the deformation field, so that the generated deformation field can accurately describe the deformation and have good spatial continuity. The dynamic adjustment of the λ value enables the network to converge quickly in the early stage of training, and fine-tune in the later stage, further improving the registration accuracy. The fused image is registered using the trained U-Net network, and a precise registration field describing the fused image to the standard atlas space is generated. The precise registration field has sub-voxel accuracy and can more accurately describe the spatial correspondence between the fused image and the standard atlas, providing high-precision deformation information for subsequent atlas deformation migration.

[0024] Preferably, step S22 includes the following steps:

[0025] Step S221: extract encoder features from the pre-built basic U-Net architecture to obtain feature maps of each layer of the encoder; perform feature map channel adjustment on the feature maps of each layer of the encoder to obtain feature maps after channel adjustment;

[0026] Step S222: performing channel attention weighting on the channel-adjusted feature map to obtain an attention-weighted feature map; performing feature map splicing according to the attention-weighted feature map and a pre-built basic U-Net architecture to obtain a U-Net with a multi-scale module added;

[0027] Step S223: extract decoder features from the U-Net with the multi-scale module added to obtain feature maps of each layer of the decoder; extract local structure information from the feature maps of each layer of the decoder to obtain local structure information feature maps;

[0028] Step S224: Generate a spatial attention map for the local structure information feature map to obtain a spatial attention map; perform feature map weighting according to the U-Net with the multi-scale module added, the feature maps of each layer of the decoder, and the spatial attention map to obtain a U-Net embedded with the spatial attention module;

[0029] Step S225: Add a deformation field output layer to the U-Net embedded with the spatial attention module to obtain a U-Net network to be trained.

[0030] The present invention extracts encoder features from the pre-built basic U-Net architecture to obtain feature maps of different scales, which capture the detailed information of different levels of the image. The subsequent channel adjustment operation unifies the number of channels of the feature maps of the corresponding layers of the encoder and decoder, so that the jump connection can be carried out smoothly, thereby better integrating features of different scales and providing richer context information for the decoder. By performing global average pooling and global maximum pooling on the feature maps after channel adjustment, and combining the information of the two to generate channel attention weights, the important feature channels are effectively highlighted, the redundant feature channels are suppressed, and the feature expression efficiency of the model is improved. The feature map after attention weighting is spliced ​​with the decoder feature map, which further enriches the input information of the decoder and helps to restore the image details more accurately. The decoder part upsamples the feature map of the U-Net after adding the multi-scale module, and gradually restores the image resolution. After extracting the feature maps of each layer of the decoder, the gradient information is calculated using the Sobel operator to obtain the feature map of local structure information. These local structural information, such as edges and corners, are crucial for accurately locating the boundaries of anatomical structures and help to improve segmentation accuracy. By performing convolution and Sigmoid activation operations on the local structural information feature map, a spatial attention map is generated, which highlights the areas with rich local structural information in the image, such as organ boundaries. The spatial attention map is multiplied with the decoder feature map to achieve spatial weighting of the feature map, so that the network pays more attention to important anatomical structure areas, thereby improving the accuracy of segmentation. A 1×1×1 convolution layer and a tanh activation function are added at the end of the U-Net network as the deformation field output layer, so that the network can output a three-channel deformation field, representing the displacement in the x, y, and z directions respectively. The tanh activation function limits the output value of the deformation field to the range of [-1, 1], ensuring the numerical stability of the deformation field and helping the training convergence of the network.

[0031] Preferably, step S23 includes the following steps:

[0032] Step S231: preparing training data for the fused image and the preprocessed standard atlas set to obtain training data;

[0033] Step S232: performing image characteristic data enhancement of different modalities on the training data to obtain enhanced training data;

[0034] Step S233: inputting the enhanced training data into the U-Net network to be trained to obtain a predicted deformation field;

[0035] Step S234: deforming the preprocessed standard atlas set according to the predicted deformation field to obtain a deformed atlas; performing similarity measurement calculation on the deformed atlas and the fused image to obtain a similarity value;

[0036] Step S235: calculating the smoothness of the predicted deformation field to obtain a smoothness value;

[0037] Step S236: construct and optimize a joint loss function according to the similarity value, the enhanced training data, and the U-Net network to be trained to obtain a trained U-Net network.

[0038] The present invention forms training data by pairing fused images with preprocessed standard atlas sets, thereby ensuring the consistency and validity of input data and laying a foundation for subsequent network training. The data set is divided into a training set and a validation set, which can effectively evaluate the generalization ability of the model and prevent overfitting. By performing data enhancement operations such as random rotation, translation and scaling on the training data, the diversity of the training data is increased, various image changes that may be encountered in actual scenes are simulated, and the robustness and generalization ability of the model are effectively improved, making the model more adaptable when processing image data of different patients. At the same time, data enhancement is performed considering the characteristics of different modal images, and the structural information of the image is better retained. The enhanced training data is input into the U-Net network to be trained, and a predicted deformation field is obtained. The deformation field is the embodiment of the nonlinear mapping relationship between the fused image and the standard atlas learned by the network, and is the key to the subsequent atlas deformation and calculation of similarity measurement. The preprocessed standard atlas set is deformed using the predicted deformation field to obtain a deformed atlas, which is more spatially aligned with the fused image. By calculating the Dice coefficient (DSC) between the deformation atlas and the fused image, a similarity metric is obtained, which is used to evaluate the accuracy of the registration and is used as part of the loss function to guide the training of the network. By calculating the L2 norm of the gradient of the predicted deformation field, the deformation field smoothness value is obtained, which reflects the spatial continuity of the deformation field. A smooth deformation field is more consistent with the actual anatomical changes and avoids unreasonable distortion and deformation. By constructing a joint loss function, the similarity metric and the deformation field smoothness are combined to guide the training of the U-Net network. Minimizing the joint loss function can improve the registration accuracy and the smoothness of the deformation field at the same time. The use of the Adam optimizer can adaptively adjust the learning rate and speed up the training of the network. The strategy of dynamically adjusting the λ value can balance the importance of similarity and smoothness at different stages of training, thereby obtaining better registration results. Finally, a trained U-Net network is obtained, which can be used to generate a high-precision deformation field.

[0039] Preferably, step S3 comprises the following steps:

[0040] Step S31: extracting atlas labels from the preprocessed standard atlas set to obtain a standard atlas label set;

[0041] Step S32: Calculate the thin plate spline deformation parameters according to the fused image and the precise registration field to obtain the thin plate spline deformation parameters;

[0042] Step S33: deforming the standard atlas label set using the thin plate spline deformation parameters to obtain a deformed atlas label set;

[0043] Step S34: generating a deformation map for the standard map label set and the pre-processed standard map set to obtain a deformation map;

[0044] Step S35: Calculate the prior probability according to the deformation spectrum to obtain a prior probability map;

[0045] Step S36: constructing a grayscale distribution model of the deformation map using the grayscale information of the fused image to obtain a tissue grayscale distribution model;

[0046] Step S37: Calculate the posterior probability of the fused image according to the prior probability map and the tissue grayscale distribution model to obtain a tissue probability map.

[0047] The present invention obtains a standard atlas label set by extracting atlas label information from a preprocessed standard atlas set, providing necessary input data for subsequent label deformation and prior probability calculation. The extracted label information accurately marks the area of ​​each anatomical structure and is the basis for atlas segmentation. The thin plate spline deformation parameters are calculated based on the fused image, the precise registration field and the local gradient information of the image, so that the deformation process can better adapt to the local structural features of the image, thereby more accurately migrating the anatomical structure information of the atlas to the space of the fused image. The introduction of local structural information, especially the gradient amplitude as the control point weight, can make the deformation more fit the edge and contour of the image and avoid unreasonable distortion. The standard atlas label set is deformed using the calculated thin plate spline deformation parameters to obtain a deformed atlas label set, so that the label information of the atlas is aligned with the space of the fused image, providing important spatial position information for the subsequent probability map generation. The application of the nearest neighbor interpolation method ensures the integrity and consistency of the deformed label value. The deformed label information is combined with the preprocessed standard atlas to generate a deformed atlas, which contains both the grayscale information of the atlas and the label information aligned with the fused image, providing a more accurate reference for the subsequent prior probability calculation and grayscale distribution model construction. The deformed labels are used as a guide to interpolate and adjust the grayscale values ​​of the atlas, further improving the accuracy and consistency of the deformed atlas. The prior probability of each voxel belonging to a different tissue category is calculated based on the deformation atlas, and a prior probability map is obtained. The map provides a preliminary probability estimate of each voxel belonging to a different tissue category and provides important prior information for the subsequent posterior probability calculation. The contribution of different atlases is weighted based on the registration accuracy to improve the reliability of the prior probability. Using the grayscale information of the fused image and the label information of the deformation atlas, a Gaussian mixture model (GMM) was constructed for each tissue category. This model can more accurately describe the grayscale distribution characteristics of each tissue category, especially when the grayscale distribution is uneven or there are multiple peaks. It is more advantageous than a single Gaussian model and provides a more accurate likelihood estimate for the posterior probability calculation. Combining the prior probability map and the tissue grayscale distribution model, the posterior probability of each voxel belonging to a different tissue category is calculated to obtain a tissue probability map. This map combines the prior knowledge of the atlas and the grayscale information of the image, provides a more accurate estimate of the probability of tissue categories, and provides important input information for the subsequent conditional random field optimization.

[0048] Preferably, step S32 includes the following steps:

[0049] Step S321: performing image gradient calculation on the fused image according to the precise registration field to obtain image gradient data;

[0050] Step S322: Calculate the gradient amplitude according to the image gradient data to obtain the gradient amplitude;

[0051] Step S323: selecting control points for the precise registration field to obtain a control point set;

[0052] Step S324: Calculate the control point weights of the control point set according to the gradient amplitude to obtain weighted control points;

[0053] Step S325: Calculate thin plate spline parameters using weighted control points and precise registration fields to obtain thin plate spline deformation parameters.

[0054] The present invention obtains image gradient data by calculating the gradient of the fused image in the three directions of x, y, and z. The data reflects the degree of change of the local grayscale of the image and provides basic data for the subsequent gradient amplitude calculation and control point weight calculation. The GradientMagnitudeRecursiveGaussianImageFilter function in the SimpleITK library can efficiently calculate the image gradient, and the preprocessing of Gaussian smoothing can effectively suppress the influence of noise. The gradient amplitude is calculated according to the image gradient data, and a gradient amplitude map reflecting the intensity of the local grayscale change of the image is obtained. Areas with large gradient amplitudes usually correspond to the edges of the image or areas with rich texture details, which require more precise control during the registration process. Control points are uniformly sampled on the fused image, and these control points are used as reference points for deformation to define the thin plate spline interpolation function. Uniform sampling can ensure the uniformity of the distribution of control points and cover various areas of the image. The control point set is weighted according to the gradient amplitude to obtain weighted control points. Control points with large gradient amplitudes are given larger weights, which means that during the thin plate spline interpolation process, the areas around these control points will have larger deformations, thereby better preserving the edge and detail information of the image. The thin plate spline deformation parameters are calculated using weighted control points and precise registration fields, and the parameters describing the nonlinear deformation relationship from the fused image to the standard atlas space are obtained. The introduction of control point weights can make thin plate spline interpolation pay more attention to the local structural information of the image, especially the areas with large gradient amplitudes, thereby improving the precision and accuracy of the deformation and making the deformed atlas better aligned with the fused image.

[0055] Preferably, step S4 comprises the following steps:

[0056] Step S41: constructing a univariate potential function for the tissue probability map to obtain a univariate potential function;

[0057] Step S42: performing multi-scale Gaussian filtering on the fused image to obtain a multi-scale image; performing multi-scale image gradient calculation on the multi-scale image to obtain multi-scale gradient information;

[0058] Step S43: performing local binary pattern calculation on the multi-scale image to obtain multi-scale LBP features; constructing feature vectors for the multi-scale gradient information and the multi-scale LBP features to obtain feature vectors;

[0059] Step S44: defining a binary potential function according to the fused image and the feature vector to obtain a binary potential function;

[0060] Step S45: constructing a conditional random field energy function for the one-dimensional potential function and the two-dimensional potential function to obtain a conditional random field energy function;

[0061] Step S46: Optimizing the conditional random field energy function to obtain an optimized label graph;

[0062] Step S47: Generate a graph segmentation result for the optimized label graph to obtain a graph segmentation result.

[0063] The present invention converts the tissue probability map into a unary potential function, converts the probability information of each voxel belonging to different tissue categories into an energy value, and provides a data basis for conditional random field optimization. The higher the probability value, the lower the corresponding energy value, which conforms to the optimization principle of energy minimization. The fused image is subjected to multi-scale Gaussian filtering and the gradient information is calculated to obtain the image gradient information at different scales, capture the local structural features of the image at different scales, and provide rich multi-scale information for the subsequent construction of more distinguishable feature vectors. By calculating the local binary pattern (LBP) features of the multi-scale image and combining the multi-scale gradient information to construct a feature vector, a feature vector containing rich local texture and structural information is obtained, which can better describe the local features of the voxel and improve the accuracy of segmentation. The binary potential function is defined according to the spatial distance between the feature vector and the voxel, which effectively considers the spatial relationship and feature similarity between the voxels. Through the distance weight and the Gaussian kernel function, voxels that are close in space and have similar features are encouraged to belong to the same category, thereby improving the spatial consistency and smoothness of the segmentation result. The unary potential function and the binary potential function are combined to construct the conditional random field energy function, which comprehensively considers the individual probability information of the voxels and the spatial relationship between the voxels, and provides an objective function for the optimization of the energy function. The conditional random field energy function is optimized using the graph cut algorithm, the minimum value of the energy function is found, and the optimal label assignment for each voxel is obtained, which effectively solves the noise and isolated area problems that may exist in the segmentation results and improves the accuracy and robustness of the segmentation. The optimized label map is converted into the final atlas segmentation result, and the tissue category to which each voxel belongs is obtained, completing the segmentation task of the neural CT atlas. This step converts the numerical label into a tissue category label with more semantic information, which is convenient for subsequent clinical application and analysis.

[0064] Preferably, step S44 includes the following steps:

[0065] Step S441: extracting voxel coordinate information from the fused image to obtain voxel coordinate information; calculating the distance between voxels based on the voxel coordinate information to obtain the distance between voxels;

[0066] Step S442: Calculate the distance weight according to the distance between voxels to obtain the distance weight;

[0067] Step S443: Calculate the feature vector difference according to the feature vector and the voxel coordinate information to obtain feature vector difference data;

[0068] Step S444: performing Gaussian kernel function calculation according to the feature vector difference data to obtain Gaussian kernel similarity;

[0069] Step S445: Calculate the binary potential function value according to the Gaussian kernel similarity and the distance weight to obtain the binary potential function.

[0070] The present invention provides necessary spatial information for subsequent calculation of distance weight and feature vector difference by extracting voxel coordinate information of fused images and calculating distance between voxels. Extracting voxel coordinate information can accurately determine the position of each voxel in three-dimensional space, and calculating distance between voxels provides quantitative basis for measuring spatial relationship between voxels. Distance weight is calculated according to distance between voxels, so that the influence of spatial distance on binary potential function can be reflected. The closer the distance of voxels, the greater the distance weight, which means that in the optimization process of conditional random field, the label consistency of adjacent voxels will be more encouraged, thereby promoting the smoothness of segmentation results. By calculating the square of Euclidean distance of feature vectors between adjacent voxels, feature vector difference data is obtained, which reflects the similarity of features between adjacent voxels and provides basic data for subsequent calculation of Gaussian kernel similarity. Gaussian kernel similarity is calculated using feature vector difference data, feature difference is converted into similarity, and similarity is mapped into a suitable range by Gaussian kernel function. The introduction of Gaussian kernel function makes the influence of feature similarity on binary potential function smoother and more natural. The binary potential function value is calculated by combining Gaussian kernel similarity and distance weight, which comprehensively considers the spatial relationship and feature similarity between voxels. By multiplying the distance weight and Gaussian kernel similarity, voxels that are close in space and have similar features have lower energy, thereby encouraging them to belong to the same category in the segmentation result, improving the spatial consistency and accuracy of the segmentation result. BRIEF DESCRIPTION OF THE DRAWINGS

[0071] Figure 1 A schematic diagram of the steps of neural CT atlas segmentation based on high-precision registration;

[0072] Figure 2 Detailed implementation flow chart of step S2 in the present invention;

[0073] Figure 3 It is a schematic diagram of the detailed implementation steps of step S3 in the present invention.

[0074] The purpose, features and advantages of the present invention will be further described with reference to the accompanying drawings in conjunction with the embodiments. DETAILED DESCRIPTION

[0075] The technical method of the present invention is described clearly and completely below in conjunction with the accompanying drawings. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by technicians in this field without creative work are within the scope of protection of the present invention.

[0076] In addition, the accompanying drawings are only schematic illustrations of the present invention and are not necessarily drawn to scale. The same reference numerals in the figures represent the same or similar parts, and their repeated description will be omitted. Some of the block diagrams shown in the accompanying drawings are functional entities and do not necessarily correspond to physically or logically independent entities. The functional entities can be implemented in software form, or implemented in one or more hardware modules or integrated circuits, or implemented in different networks and / or processor methods and / or microcontroller methods.

[0077] It should be understood that, although the terms "first", "second", etc. may be used herein to describe various units, these units should not be limited by these terms. These terms are used only to distinguish one unit from another unit. For example, without departing from the scope of the exemplary embodiments, the first unit may be referred to as the second unit, and similarly the second unit may be referred to as the first unit. The term "and / or" used herein includes any and all combinations of one or more of the listed associated items.

[0078] refer to Figure 1 FIG. 1 is a schematic diagram of the steps of the neural CT atlas segmentation method based on high-precision registration of the present invention. In this example, the neural CT atlas segmentation method based on high-precision registration includes the following steps:

[0079] Step S1: Obtaining a neural CT image, an MRI-T1 image, and an MRI-T2 image; performing a coarse registration based on mutual information on the CT image, the MRI-T1 image, and the MRI-T2 image, and performing a multi-scale fusion based on wavelet to obtain a fused image;

[0080] In an embodiment of the present invention, a neuro CT image, an MRI-T1 image, and an MRI-T2 image are obtained. The three modal images are preprocessed, including denoising, bias field correction, and grayscale normalization. The preprocessed CT image is used as a reference image, and the MRI-T1 and MRI-T2 images are registered to the space of the CT image using a mutual information-based registration method. Then, a wavelet-based multi-scale fusion method is used to fuse the three registered images into a fused image, which integrates the image information of the three modalities.

[0081] Step S2: preprocessing the preset standard neural CT atlas library with standard atlases according to the fused image to obtain a preprocessed standard atlas set; adding a multi-scale feature fusion module to the pre-built basic U-Net architecture, and embedding the spatial attention module to obtain a U-Net network to be trained; using the fused image and the preprocessed standard atlas set to perform deep registration network training on the U-Net network to be trained to obtain a trained U-Net network; using the trained U-Net network to generate a precise registration field for the fused image to obtain a precise registration field;

[0082] In an embodiment of the present invention, first, according to the voxel size and dimensional information of the fused image, the atlas in the preset standard neural CT atlas library is preprocessed, including resampling and grayscale normalization, to obtain a preprocessed standard atlas set. Then, a U-Net network to be trained is constructed, which includes a multi-scale feature fusion module and a spatial attention module. The fused image and the preprocessed standard atlas set are used as training data to train the U-Net network to learn the nonlinear deformation relationship from the fused image to the standard atlas. During the training process, data enhancement technology is used to improve the generalization ability of the model, and a joint loss function (including similarity measurement and deformation field smoothness) is used to optimize the network parameters. Finally, the trained U-Net network is used to align the fused image to generate a precise alignment field, which describes the precise pixel-level correspondence between the fused image and the standard atlas space.

[0083] Step S3: Calculate the thin plate spline deformation parameters according to the fused image and the precise registration field to obtain the thin plate spline deformation parameters; use the thin plate spline deformation parameters to perform atlas deformation migration on the preprocessed standard atlas set to obtain a deformation atlas; calculate the prior probability according to the deformation atlas to obtain a prior probability map; calculate the posterior probability according to the prior probability map to obtain a tissue probability map;

[0084] In an embodiment of the present invention, atlas label information is extracted from a preprocessed standard atlas set to form a standard atlas label set. The deformation parameters of the thin plate spline interpolation are calculated based on the fused image, the precise registration field and the local gradient information of the image. The standard atlas label set is deformed using the calculated thin plate spline deformation parameters to obtain a deformed atlas label set. The deformed label information is combined with the preprocessed standard atlas to generate a deformation atlas. The prior probability that each voxel belongs to a different tissue category is calculated based on the deformation atlas to obtain a prior probability map. Using the grayscale information of the fused image and the label information of the deformation atlas, a Gaussian mixture model (GMM) is constructed for each tissue category to describe its grayscale distribution. Finally, combining the prior probability and the GMM model, the posterior probability that each voxel belongs to a different tissue category is calculated to generate a tissue probability map.

[0085] Step S4: construct a univariate potential function for the tissue probability map to obtain a univariate potential function; extract local features for the fused image to obtain a feature vector; define a binary potential function based on the fused image and the feature vector to obtain a binary potential function; construct a conditional random field energy function for the univariate potential function and the binary potential function to obtain a conditional random field energy function; generate a spectrum segmentation result based on the conditional random field energy function to obtain a spectrum segmentation result, so as to realize the neural CT spectrum segmentation task;

[0086] In an embodiment of the present invention, the tissue probability map is converted into a unary potential function. The fused image is subjected to multi-scale Gaussian filtering, and the gradient information and local binary pattern (LBP) features of each scale image are calculated, and these features are combined into a feature vector. Based on the spatial distance between the feature vector and the voxel, the binary potential function is calculated, wherein the distance weight is used to adjust the influence of the spatial distance on the binary potential function, and the Gaussian kernel function is used to measure the similarity of the feature vector. The unary potential function and the binary potential function are combined to construct a conditional random field energy function. The energy function is optimized using a graph cut algorithm to obtain an optimized label map. Finally, the optimized label map is converted into the final atlas segmentation result, wherein the label of each voxel represents the tissue category to which it belongs.

[0087] Preferably, step S1 comprises the following steps:

[0088] Step S11: acquiring a neural CT image, an MRI-T1 image, and an MRI-T2 image; performing image modality preprocessing on the CT image, the MRI-T1 image, and the MRI-T2 image to obtain a preprocessed CT image, a preprocessed MRI-T1 image, and a preprocessed MRI-T2 image;

[0089] Step S12: using the pre-processed CT image as a reference image, performing coarse registration based on mutual information on the pre-processed MRI-T1 image and the pre-processed MRI-T2 image to obtain a registered MRI-T1 image and a registered MRI-T2 image;

[0090] Step S13: performing wavelet-based multi-scale fusion on the pre-processed CT image, the registered MRI-T1 image, and the registered MRI-T2 image to obtain a fused image.

[0091] In an embodiment of the present invention, a neurological CT image, an MRI-T1 image, and an MRI-T2 image are obtained. The neurological CT image, the MRI-T1 image, and the MRI-T2 image are read respectively using ITK-SNAP software. The CT image is denoised, and a Gaussian filter with a convolution kernel size of 3×3×3 and a standard deviation of 1.0 is used to perform a convolution operation on the image to smooth the image and suppress noise. The MRI-T1 image and the MRI-T2 image are respectively subjected to N4 bias field correction, and the N4BiasFieldCorrection tool in the ANTs software is used, and the shrinkfactor is set to 4, and the image is biased and corrected using B-spline interpolation to eliminate the intensity inhomogeneity in the MRI image. Grayscale normalization is performed on all images, and the grayscale value range of all images is linearly mapped to between [0, 1] to unify the grayscale range of images of different modalities. The formula is: (value-min) / (max-min), where value is the original grayscale value, and min and max are the minimum and maximum grayscale values ​​of the image, respectively. Output preprocessed CT images, MRI-T1 images, and MRI-T2 images.

[0092] The pre-processed CT image in step S11 is used as a reference image. Using Elastix software, the pre-processed MRI-T1 image and MRI-T2 image are respectively registered with the reference CT image. The registration parameters are set as follows: the transformation type is affine transformation, the similarity metric is normalized mutual information, the optimizer is a gradient descent optimizer, and the maximum number of iterations is 500. Elastix software iteratively optimizes the transformation parameters to maximize the normalized mutual information of the registered MRI-T1 image and MRI-T2 image with the reference CT image, thereby achieving spatial alignment between images of different modalities. Output the registered MRI-T1 image and MRI-T2 image.

[0093] The PyWavelets library in Python is used to perform wavelet-based multi-scale fusion on the preprocessed CT image output from step S11, the registered MRI-T1 image output from step S12, and the registered MRI-T2 image. First, the three images are decomposed into three layers of wavelets using the 'db4' wavelet basis to obtain the low-frequency subband and high-frequency subband coefficients. For the low-frequency subband coefficients, the weighted average fusion rule is adopted, and the weights are set to 0.4 (CT), 0.3 (T1), and 0.3 (T2), respectively. For the high-frequency subband coefficients, the fusion rule with the largest local variance is adopted, and the coefficient with the largest local variance is selected as the fusion coefficient. The local variance calculation window size is set to 3×3. Finally, the fused low-frequency and high-frequency coefficients are reconstructed by wavelets using the 'db4' wavelet basis to obtain the final fused image. Output the fused image.

[0094] Preferably, step S2 comprises the following steps:

[0095] Step S21: performing standard atlas preprocessing on a preset standard neural CT atlas library according to the fusion image to obtain a preprocessed standard atlas set;

[0096] Step S22: adding a multi-scale feature fusion module to the pre-built basic U-Net architecture and embedding a spatial attention module to obtain a U-Net network to be trained;

[0097] Step S23: using the fused image and the preprocessed standard atlas set to perform deep registration network training on the U-Net network to be trained, to obtain a trained U-Net network;

[0098] Step S24: Using the trained U-Net network to generate a precise registration field for the fused image to obtain a precise registration field.

[0099] As an example of the present invention, refer to Figure 2 As shown, in this example, step S2 includes:

[0100] Step S21: performing standard atlas preprocessing on a preset standard neural CT atlas library according to the fusion image to obtain a preprocessed standard atlas set;

[0101] In an embodiment of the present invention, a preset standard neuro-CT atlas library contains a plurality of standard atlases and their corresponding labels that have been registered to the same coordinate space. For the fused image, its voxel size and dimension information are obtained. A resampling operation is performed on each atlas in the standard neuro-CT atlas library, and the ResampleImageFilter function in the SimpleITK library is used to keep the voxel size and dimension of the atlas consistent with the fused image, and the interpolation method adopts B-spline interpolation. Then, the grayscale of each resampled atlas is normalized, and the grayscale value range is linearly mapped to [0, 1]. The formula is: (value-min) / (max-min), where value is the original grayscale value, and min and max are the minimum and maximum grayscale values ​​of the image, respectively. Output the preprocessed standard atlas set.

[0102] Step S22: adding a multi-scale feature fusion module to the pre-built basic U-Net architecture and embedding a spatial attention module to obtain a U-Net network to be trained;

[0103] In an embodiment of the present invention, the U-Net network to be trained is constructed using TensorFlow and Keras libraries. The basic U-Net architecture consists of an encoder and a decoder. The encoder part includes four downsampling modules, each of which consists of two 3×3×3 convolutional layers (with a step size of 1 and the same padding) and a 2×2×2 maximum pooling layer (with a step size of 2), and the convolutional layer uses the ReLU activation function. The decoder part includes four upsampling modules, each of which consists of a 2×2×2 deconvolutional layer (with a step size of 2) and two 3×3×3 convolutional layers (with a step size of 1 and the same padding), and the convolutional layer uses the ReLU activation function. A jump connection is added between the encoder and the decoder, and the output of each downsampling module in the encoder is spliced ​​with the input of the corresponding upsampling module in the decoder in the channel dimension. In the multi-scale feature fusion module, the output of each downsampling layer of the encoder is adjusted through a 1×1×1 convolutional layer, and then spliced ​​with the input of the corresponding layer of the decoder in the channel dimension. In the spatial attention module, the output of each upsampling layer of the decoder is fed into two 1×1×1 convolutional layers, and then concatenated in the channel dimension. Then, an attention map is generated through a 1×1×1 convolutional layer and a Sigmoid activation function. Finally, the attention map is multiplied by the output of the decoder layer. A 1×1×1 convolutional layer is added to the last layer of the U-Net network. The number of output channels is 3, representing the three directions of the deformation field. The tanh activation function is used to limit the output to the range of [-1, 1]. Output the U-Net network to be trained.

[0104] Step S23: using the fused image and the preprocessed standard atlas set to perform deep registration network training on the U-Net network to be trained, to obtain a trained U-Net network;

[0105] In an embodiment of the present invention, the fused image and the preprocessed standard atlas set are used to form a training data set. The training data is enhanced, including random rotation (angle range [-10°, 10°]), random translation (range [-10, 10] voxels) and random scaling (range [0.9, 1.1]). The U-Net network to be trained constructed using step S22 is input as the fused image and the standard atlas, and the output is the deformation field. The loss function is defined as the negative normalized mutual information between the deformed atlas and the fused image plus the L2 norm of the deformation field gradient, where the balance coefficient λ is set to 0.1. The Adam optimizer is used for training, the learning rate is set to 0.0001, the batch size is set to 4, and the number of training rounds is set to 100. After each round of training, the model performance is evaluated on the validation set, and the best model is saved. The trained U-Net network is output.

[0106] Step S24: using the trained U-Net network to generate a precise registration field for the fused image to obtain a precise registration field;

[0107] In an embodiment of the present invention, the fused image is input into a trained U-Net network. The network outputs a three-channel deformation field, which represents the deformation vector of the fused image to the standard atlas space. The output deformation field is scaled to match the voxel size of the fused image. The fine registration field is output.

[0108] Preferably, step S22 includes the following steps:

[0109] Step S221: extract encoder features from the pre-built basic U-Net architecture to obtain feature maps of each layer of the encoder; perform feature map channel adjustment on the feature maps of each layer of the encoder to obtain feature maps after channel adjustment;

[0110] Step S222: performing channel attention weighting on the channel-adjusted feature map to obtain an attention-weighted feature map; performing feature map splicing according to the attention-weighted feature map and a pre-built basic U-Net architecture to obtain a U-Net with a multi-scale module added;

[0111] Step S223: extract decoder features from the U-Net with the multi-scale module added to obtain feature maps of each layer of the decoder; extract local structure information from the feature maps of each layer of the decoder to obtain local structure information feature maps;

[0112] Step S224: Generate a spatial attention map for the local structure information feature map to obtain a spatial attention map; perform feature map weighting according to the U-Net with the multi-scale module added, the feature maps of each layer of the decoder, and the spatial attention map to obtain a U-Net embedded with the spatial attention module;

[0113] Step S225: Add a deformation field output layer to the U-Net embedded with the spatial attention module to obtain a U-Net network to be trained.

[0114] In the embodiment of the present invention, the encoder part of the pre-built basic U-Net architecture includes four downsampling modules, each of which consists of two 3×3×3 convolutional layers (with a step size of 1 and the same padding) and a 2×2×2 maximum pooling layer (with a step size of 2). A single-channel image of size 128×128×128 is input into the encoder. The feature maps output by each downsampling module are recorded as Feature_Map_1 (size of 64×64×64, number of channels is 64), Feature_Map_2 (size of 32×32×32, number of channels is 128), Feature_Map_3 (size of 16×16×16, number of channels is 256), and Feature_Map_4 (size of 8×8×8, number of channels is 512). In order to match the number of channels of the corresponding layer of the decoder, a 1×1×1 convolution kernel is used to adjust the channels of the feature maps of each layer of the encoder. The number of channels of Feature_Map_1 is adjusted to 512, the number of channels of Feature_Map_2 is adjusted to 256, the number of channels of Feature_Map_3 is adjusted to 128, and the number of channels of Feature_Map_4 is adjusted to 64. The feature maps after channel adjustment are output as Adjusted_Feature_Map_1 to Adjusted_Feature_Map_4.

[0115] Perform channel attention weighting on the feature maps Adjusted_Feature_Map_1 to Adjusted_Feature_Map_4 after channel adjustment. Taking Adjusted_Feature_Map_1 as an example, first perform global average pooling and global maximum pooling on it to obtain two feature vectors of size 1×1×512. Concatenate these two feature vectors to obtain a feature vector of size 1×1×1024. Then, pass the vector through two fully connected layers in sequence. The first fully connected layer has 512 neurons and uses the ReLU activation function; the second fully connected layer has 1024 neurons and does not use the activation function. Then pass the output result through the Sigmoid activation function to obtain a weight vector of size 1×1×1024. Multiply the first 512 values ​​of the weight vector with the global average pooling result, and multiply the last 512 values ​​with the global maximum pooling result, and then add the two results to obtain the final weighted feature map Weighted_Feature_Map_1. The same operation is performed on the other three feature maps to obtain Weighted_Feature_Map_2 to Weighted_Feature_Map_4. The attention-weighted feature maps Weighted_Feature_Map_1 to Weighted_Feature_Map_4 are concatenated with the feature maps of the corresponding layers in the pre-built basic U-Net architecture decoder in the channel dimension to obtain the U-Net with the multi-scale module added.

[0116] Assume that the decoder part of the U-Net with the multi-scale module added contains four upsampling modules, symmetrical to the encoder part. A single-channel image of size 128×128×128 and a single-channel atlas of the same size are input into the U-Net. Record the feature maps output by each upsampling module in the decoder part, which are recorded as Decoder_Feature_Map_1 (size 16×16×16, number of channels 512), Decoder_Feature_Map_2 (size 32×32×32, number of channels 256), Decoder_Feature_Map_3 (size 64×64×64, number of channels 128) and Decoder_Feature_Map_4 (size 128×128×128, number of channels 64). Use the Sobel operator to calculate the gradient of each decoder feature map in the x, y, and z directions respectively, and concatenate the gradient maps in the three directions in the channel dimension to obtain the local structure information feature maps Local_Structure_1 to Local_Structure_4. The number of channels of Local_Structure_1 is , and so on for the number of channels of other feature maps.

[0117] Generate spatial attention maps for local structure information feature maps Local_Structure_1 to Local_Structure_4. Taking Local_Structure_1 as an example, use a 1×1×1 convolution kernel to reduce the number of channels of Local_Structure_1 from 1536 to 512, and then pass a Sigmoid activation function to generate spatial attention map Spatial_Attention_Map_1 (size is 16×16×16, number of channels is 512). Perform the same operation on the other three local structure information feature maps to obtain Spatial_Attention_Map_2 to Spatial_Attention_Map_4. Multiply the spatial attention maps Spatial_Attention_Map_1 to Spatial_Attention_Map_4 with the corresponding decoder feature maps Decoder_Feature_Map_1 to Decoder_Feature_Map_4 by pixel, and get the weighted feature maps Attention_Weighted_Feature_Map_1 to Attention_Weighted_Feature_Map_4. Use the weighted feature maps as the new decoder feature maps, update the U-Net with the multi-scale module, and get the U-Net embedded with the spatial attention module.

[0118] Add a 1×1×1 convolutional layer after the last layer of the U-Net embedded in the spatial attention module (the output of Decoder_Feature_Map_4 after the spatial attention mechanism), with the number of output channels being 3, the convolution kernel initialized to a truncated normal distribution (mean 0, standard deviation 0.01), and the bias initialized to 0. This convolutional layer outputs a three-channel deformation field, representing the displacement in the x, y, and z directions respectively. Use the tanh activation function to limit the output value of the deformation field to the range of [-1, 1]. Output the U-Net network to be trained.

[0119] Preferably, step S23 includes the following steps:

[0120] Step S231: preparing training data for the fused image and the preprocessed standard atlas set to obtain training data;

[0121] Step S232: performing image characteristic data enhancement of different modalities on the training data to obtain enhanced training data;

[0122] Step S233: inputting the enhanced training data into the U-Net network to be trained to obtain a predicted deformation field;

[0123] Step S234: deforming the preprocessed standard atlas set according to the predicted deformation field to obtain a deformed atlas; performing similarity measurement calculation on the deformed atlas and the fused image to obtain a similarity value;

[0124] Step S235: calculating the smoothness of the predicted deformation field to obtain a smoothness value;

[0125] Step S236: construct and optimize a joint loss function according to the similarity value, the enhanced training data, and the U-Net network to be trained to obtain a trained U-Net network.

[0126] In an embodiment of the present invention, the fused image and the preprocessed standard atlas set are paired to form training data. Assume that the preprocessed standard atlas set contains N atlases. The fused image is copied N times and paired with N atlases respectively to form N pairs of training data. Each pair of training data contains a fused image and a corresponding standard atlas. Each pair of training data is stored as a separate file, for example, in HDF5 format, for subsequent efficient reading. The training data is divided into a training set and a validation set in a ratio of 9:1. Output the paired training data, including the training set and the validation set.

[0127] Perform data enhancement on the training data. Apply the same enhancement operation to the fused image and the corresponding standard atlas to maintain the spatial correspondence between them. Use a random rotation transformation with a rotation angle range of [-10°, 10°], the rotation axis is the center of the image, and the interpolation method is trilinear interpolation. Use a random translation transformation with a translation range of [-10, 10] voxels and a mirror filling mode. Use a random scaling transformation with a scaling range of [0.9, 1.1] and a trilinear interpolation method. Considering that the fused image contains CT and MRI information, the enhancement operation of the fused image focuses more on maintaining the bone structure and soft tissue details. Perform the same enhancement operation on the standard atlas. Output the enhanced training data.

[0128] Input the enhanced training data into the U-Net network to be trained. The network input is the fused image and the standard atlas, both of which are single-channel data with a size of 128×128×128. The network output is a three-channel predicted deformation field with a size of 128×128×128, where each channel represents the displacement in the x, y, and z directions. The numerical range of the deformation field is [-1, 1], which is guaranteed by the tanh activation function of the last layer of the network. Output the predicted deformation field.

[0129] Use the predicted deformation field to deform the preprocessed standard atlas set. Use the Transformi×ImageFilter function in the SimpleITK library to deform the standard atlas with the predicted deformation field as input. The deformed atlas has the same size as the fused image, both of which are 128×128×128. Calculate the similarity between the deformed atlas and the fused image. Use the LabelOverlapMeasuresImageFilter function in the SimpleITK library to calculate the Dice coefficient (DSC). Take the deformed atlas and the fused image as input and calculate the DSC value between the two images as a similarity measure. Output the similarity value (DSC value).

[0130] Calculate the smoothness of the predicted deformation field. Use the numpy library to calculate the gradient of the deformation field in the x, y, and z directions. Then calculate the sum of the squares of the gradients and average them over all voxels to get the smoothness of the deformation field, which is the L2 norm of the gradient of the deformation field. Output the smoothness value.

[0131] A joint loss function is constructed based on the similarity value (DSC) and the smoothness value. The loss function is defined as: ;

[0132] Where λ is the weight coefficient for balancing similarity and smoothness, and its initial value is set to 0.1. Use the Adam optimizer to optimize the parameters of the U-Net network to be trained and minimize the joint loss function. The learning rate is set to 0.0001 and the batch size is set to 4. During the training process, the λ value is reduced by 0.01 every 10 rounds until the λ value is reduced to 0.01, so that smoothness is more emphasized in the early stage of training and similarity is more emphasized in the later stage. Train the model on the training set and evaluate the model performance on the validation set. Save the model with the highest DSC value on the validation set as the trained U-Net network. Output the trained U-Net network.

[0133] Preferably, step S3 comprises the following steps:

[0134] Step S31: extracting atlas labels from the preprocessed standard atlas set to obtain a standard atlas label set;

[0135] Step S32: Calculate the thin plate spline deformation parameters according to the fused image and the precise registration field to obtain the thin plate spline deformation parameters;

[0136] Step S33: deforming the standard atlas label set using the thin plate spline deformation parameters to obtain a deformed atlas label set;

[0137] Step S34: generating a deformation map for the standard map label set and the pre-processed standard map set to obtain a deformation map;

[0138] Step S35: Calculate the prior probability according to the deformation spectrum to obtain a prior probability map;

[0139] Step S36: constructing a grayscale distribution model of the deformation map using the grayscale information of the fused image to obtain a tissue grayscale distribution model;

[0140] Step S37: Calculate the posterior probability of the fused image according to the prior probability map and the tissue grayscale distribution model to obtain a tissue probability map.

[0141] As an example of the present invention, refer to Figure 3 As shown, in this example, step S3 includes:

[0142] Step S31: extracting atlas labels from the preprocessed standard atlas set to obtain a standard atlas label set;

[0143] In an embodiment of the present invention, label information corresponding to each atlas is extracted from the preprocessed standard atlas set. Each atlas in the preprocessed standard atlas set contains a corresponding label file, and the value of each voxel in the label file represents the anatomical structure label to which the voxel belongs. The label file of each atlas is read using the SimpleITK library, and the label data is stored as a numpy array. Assuming that the preprocessed standard atlas set contains N atlases, N label arrays of size 128×128×128 are obtained, constituting a standard atlas label set. The standard atlas label set is output.

[0144] Step S32: Calculate the thin plate spline deformation parameters according to the fused image and the precise registration field to obtain the thin plate spline deformation parameters;

[0145] In an embodiment of the present invention, the deformation parameters of the thin plate spline are calculated using the precise registration field and the fused image. First, the Sobel operator is used to calculate the gradient of the fused image in the x, y, and z directions, and then the gradient amplitude is calculated. Next, 1000 control points are uniformly sampled on the fused image. The average value of the gradient amplitude in the 5×5×5 neighborhood around each control point is calculated as the weight of the control point. The deformation parameters of the thin plate spline interpolation are calculated using the position of the corresponding control point in the precise registration field as the target position, as well as the position and weight of the control point in the fused image. The deformation parameters include an affine transformation matrix and non-affine transformation parameters. The thin plate spline deformation parameters are output.

[0146] Step S33: deforming the standard atlas label set using the thin plate spline deformation parameters to obtain a deformed atlas label set;

[0147] In an embodiment of the present invention, the standard atlas label set is deformed using thin plate spline deformation parameters. The Transformi×ImageFilter function in the SimpleITK library is used to deform each label map using the thin plate spline deformation parameters as input. The interpolation method uses nearest neighbor interpolation to ensure the integrity of the deformed label value. The deformed atlas label set is output, which contains N deformed label maps, each of which has a size of 128×128×128.

[0148] Step S34: generating a deformation map for the standard map label set and the pre-processed standard map set to obtain a deformation map;

[0149] In an embodiment of the present invention, a deformation atlas is generated using a deformation atlas label set and a preprocessing standard atlas set. Each preprocessing standard atlas is combined with its corresponding deformation label map. Specifically, for each voxel, if its value in the deformation label map is 0, the grayscale value of the corresponding voxel in the deformation atlas is set to 0; otherwise, the grayscale value of the corresponding voxel in the deformation atlas remains unchanged. The deformation atlas is output, which includes N deformation maps, each of which is 128×128×128 in size and includes grayscale information and deformed label information.

[0150] Step S35: Calculate the prior probability according to the deformation spectrum to obtain a prior probability map;

[0151] In an embodiment of the present invention, the prior probability that each voxel belongs to a different tissue category is calculated based on the deformation map. Assume that there are K tissue categories. For each voxel, the labels to which the voxel belongs in the N deformation maps are counted, and the frequency of occurrence of each label is calculated. The frequency of occurrence of each label is used as the prior probability that the voxel belongs to the corresponding tissue category. The prior probability map is output, with a size of 128×128×128×K, where each voxel corresponds to a vector of length K, indicating the prior probability that the voxel belongs to each tissue category.

[0152] Step S36: constructing a grayscale distribution model of the deformation map using the grayscale information of the fused image to obtain a tissue grayscale distribution model;

[0153] In an embodiment of the present invention, a tissue grayscale distribution model is constructed using a fused image and a deformation atlas. For each tissue category k, the grayscale values ​​of all voxels marked as k in the deformation atlas are extracted from the fused image. A Gaussian mixture model (GMM) is trained using the extracted grayscale values, wherein the number of Gaussian mixture components is set to 3. The GMM model is trained using the GaussianMi×ture function in the sklearn.mixture library. The tissue grayscale distribution model is output, comprising K GMM models, each corresponding to a tissue category.

[0154] Step S37: performing a posterior probability calculation on the fused image according to the prior probability map and the tissue grayscale distribution model to obtain a tissue probability map;

[0155] In the embodiment of the present invention, the posterior probability of each voxel belonging to a different tissue category is calculated based on the prior probability map, the tissue grayscale distribution model and the fused image. For each voxel i and each tissue category k, the tissue grayscale distribution model is used to calculate the probability density value of the grayscale value of voxel i in the fused image under the GMM model of tissue category k, which is recorded as P(gray_i|k). Then, P(gray_i|k) is multiplied by the prior probability P(k) to obtain the numerator of the posterior probability . Finally, add the numerators of the posterior probabilities of all tissue categories to get the denominator. Divide the numerator by the denominator to get the posterior probability P(k|gray_i) that voxel i belongs to tissue category k. Output the tissue probability map, which is 128×128×128×K, where each voxel corresponds to a vector of length K, indicating the posterior probability that the voxel belongs to each tissue category.

[0156] Preferably, step S32 includes the following steps:

[0157] Step S321: performing image gradient calculation on the fused image according to the precise registration field to obtain image gradient data;

[0158] Step S322: Calculate the gradient amplitude according to the image gradient data to obtain the gradient amplitude;

[0159] Step S323: selecting control points for the precise registration field to obtain a control point set;

[0160] Step S324: Calculate the control point weights of the control point set according to the gradient amplitude to obtain weighted control points;

[0161] Step S325: Calculate thin plate spline parameters using weighted control points and precise registration fields to obtain thin plate spline deformation parameters.

[0162] In an embodiment of the present invention, a fused image and the SimpleITK library are used to calculate image gradient data. First, the GradientMagnitudeRecursiveGaussianImageFilter function in the SimpleITK library is used to perform Gaussian smoothing with a sigma value of 1.0, and then the gradient of each voxel in the x, y, and z directions is calculated. The gradient values ​​in the three directions are stored as three independent arrays, Gradient_x, Gradient_y, and Gradient_z, respectively. The size of each array is the same as that of the fused image, which is 128×128×128. Output image gradient data Gradient_x, Gradient_y, and Gradient_z.

[0163] Calculate the gradient magnitude using the image gradient data Gradient_x, Gradient_y, and Gradient_z. For each voxel, calculate its gradient magnitude using the formula: `Magnitude=sqrt(Gradient_x^2+Gradient_y^2+Gradient_z^2)`. Store the calculated gradient magnitude into an array Gradient_Magnitude of size 128×128×128. Output the gradient magnitude Gradient_Magnitude.

[0164] Sample control points evenly on the fused image. Set the sampling interval to 4 voxels, and sample a control point every 4 voxels in the x, y, and z directions. Store the 3D coordinates of the sampled control points in an array Control_Points, where each row of the array represents the 3D coordinates of a control point. Since the size of the fused image is 128×128×128 and the sampling interval is 4, the number of control points is Output control point set Control_Points.

[0165] Calculate the weighted control points using the gradient magnitude Gradient_Magnitude and the control point set Control_Points. For each control point, calculate the average of the gradient magnitudes in the 5×5×5 neighborhood around it. Specifically, for each control point, find its position in the fused image, and then extract the gradient magnitude in the 5×5×5 neighborhood centered on that position. Calculate the average of these gradient magnitudes as the weight of the control point. Store the weight of each control point in an array Control_Point_Weights, the size of the array is 32768, the same as the number of control points. Output the weighted control points Control_Point_Weights.

[0166] The thin plate spline deformation parameters are calculated using the control point set Control_Points, the weighted control points Control_Point_Weights and the fine registration field. The thin plate spline interpolation function is constructed using the weighted control points and the positions of the corresponding control points in the fine registration field. The coordinates of the control points in the fused image are used as the source control points, and the displacements of the corresponding control points in the fine registration field are used as the displacements of the target control points. The thin plate spline interpolation function is weighted using the control point weights so that the areas with larger gradient amplitudes have a greater impact on the deformation. The thin plate spline interpolation function is solved to obtain the deformation parameters, including the affine transformation matrix (4×4 matrix) and the non-affine transformation parameters (each control point corresponds to a three-dimensional vector). Output the thin plate spline deformation parameters.

[0167] Preferably, step S4 comprises the following steps:

[0168] Step S41: constructing a univariate potential function for the tissue probability map to obtain a univariate potential function;

[0169] Step S42: performing multi-scale Gaussian filtering on the fused image to obtain a multi-scale image; performing multi-scale image gradient calculation on the multi-scale image to obtain multi-scale gradient information;

[0170] Step S43: performing local binary pattern calculation on the multi-scale image to obtain multi-scale LBP features; constructing feature vectors for the multi-scale gradient information and the multi-scale LBP features to obtain feature vectors;

[0171] Step S44: defining a binary potential function according to the fused image and the feature vector to obtain a binary potential function;

[0172] Step S45: constructing a conditional random field energy function for the one-dimensional potential function and the two-dimensional potential function to obtain a conditional random field energy function;

[0173] Step S46: Optimizing the conditional random field energy function to obtain an optimized label graph;

[0174] Step S47: Generate a graph segmentation result for the optimized label graph to obtain a graph segmentation result.

[0175] In an embodiment of the present invention, a unary potential function is constructed using a tissue probability map. The size of the tissue probability map is 128×128×128×K, where K is the number of tissue categories. For each voxel, the value on the corresponding tissue probability map is a vector of length K, which represents the probability that the voxel belongs to each category. The probability value is converted to an energy value using the formula: `Unary_Potential=-ln(Probability)`. Wherein, ln is the natural logarithm function. The probability vector of each voxel is converted into an energy vector to form a unary potential function. The size of the unary potential function is the same as that of the tissue probability map, which is also 128×128×128×K. Output the unary potential function Unary_Potential.

[0176] Perform multi-scale Gaussian filtering on the fused image. Use Gaussian filters with standard deviations of 1, 2, and 3 to convolve the fused image, and get three multi-scale images Multiscale_Image_1, Multiscale_Image_2, and Multiscale_Image_3, all of which are 128×128×128 in size. Use the GradientMagnitudeRecursiveGaussianImageFilter function in the SimpleITK library to calculate the gradients of the three multi-scale images in the x, y, and z directions, and get nine gradient maps Gradient_x_1, Gradient_y_1, Gradient_z_1, Gradient_x_2, Gradient_y_2, Gradient_z_2, Gradient_x_3, Gradient_y_3, and Gradient_z_3, all of which are 128×128×128 in size. Output multi-scale gradient information.

[0177] Calculate the local binary pattern (LBP) for the multiscale images Multiscale_Image_1, Multiscale_Image_2, and Multiscale_Image_3. Use the local_binary_pattern function in the skimage.feature library, set the radius to 1, the number of sampling points to 8, and the method to 'uniform' to calculate the LBP feature maps LBP_1, LBP_2, and LBP_3 of the three multiscale images, respectively, with a size of 128×128×128. Combine the nine gradient maps and the three LBP feature maps into a feature vector. For each voxel, combine the 12 eigenvalues ​​into a 12-dimensional feature vector. Output the feature vector, with a size of 128×128×128×12.

[0178] The binary potential function is defined using the eigenvectors. For each pair of adjacent voxels i and j, the square of the Euclidean distance between their eigenvectors is calculated: ,in and Represent the 12-dimensional feature vectors of voxels i and j respectively. Calculate the binary potential function value: , where σ is the standard deviation of the Gaussian kernel function, set to 5, and w is the weight coefficient, set to 1. Output the binary potential function Pairwise_Potential.

[0179] The conditional random field energy function is constructed using the unary potential function Unary_Potential and the binary potential function Pairwise_Potential obtained in step S44. The energy function is defined as: ,in is a univariate potential function, is the binary potential function, represents the label of voxel i.

[0180] Use the graph cut algorithm to optimize the conditional random field energy function. Use the PyMa×flow library to implement the graph cut algorithm. Convert the energy function to a graph cut problem, where each voxel corresponds to a node and an edge connects adjacent voxels. The weight of the edge is determined by the binary potential function. The weight of the node is determined by the unary potential function. Use the graph cut algorithm to find the minimum value of the energy function and obtain the optimal label assignment for each voxel, that is, the optimized label map Optimized_Label_Map, with a size of 128×128×128. Output the optimized label map.

[0181] The optimized label map Optimized_Label_Map obtained in step S46 is used to generate the final atlas segmentation result. The label value of each voxel in the label map is converted into the corresponding tissue category label. For example, if the label value is 1, it means that the voxel belongs to the gray matter of the brain; if the label value is 2, it means that the voxel belongs to the white matter of the brain, and so on. The final atlas segmentation result is output with a size of 128×1128×128, where the value of each voxel represents the tissue category to which it belongs.

[0182] Preferably, step S44 includes the following steps:

[0183] Step S441: extracting voxel coordinate information from the fused image to obtain voxel coordinate information; calculating the distance between voxels based on the voxel coordinate information to obtain the distance between voxels;

[0184] Step S442: Calculate the distance weight according to the distance between voxels to obtain the distance weight;

[0185] Step S443: Calculate the feature vector difference according to the feature vector and the voxel coordinate information to obtain feature vector difference data;

[0186] Step S444: performing Gaussian kernel function calculation according to the feature vector difference data to obtain Gaussian kernel similarity;

[0187] Step S445: Calculate the binary potential function value according to the Gaussian kernel similarity and the distance weight to obtain the binary potential function.

[0188] In an embodiment of the present invention, for a fused image of size 128×128×128, the coordinate information of each voxel is extracted. Three arrays X, Y and Z are used to store the coordinate values ​​of each voxel in the x, y and z directions respectively. The sizes of arrays X, Y and Z are all 128×128×128. The values ​​in array X increase from 0 to 127, indicating the coordinates of the voxels in the x direction. Arrays Y and Z are similar, representing the coordinates of the voxels in the y and z directions respectively. Calculate the distance between each pair of adjacent voxels. Consider 6-neighborhood connectivity, that is, each voxel is connected to its directly adjacent voxels in the x, y and z directions. For each voxel i, calculate its distance to the 6 adjacent voxels. Use the Euclidean distance formula to calculate the distance: ,in and are the coordinates of voxels i and j respectively. Output the distance between voxels Distance.

[0189] The distance weight is calculated using the distance between voxels Distance. The distance weight is calculated using the inverse function: To avoid the distance being 0, the weight of the distance being 0 is set to 1. Output the distance weight Distance_Weight.

[0190] The eigenvector difference is calculated using the eigenvector and voxel coordinate information. For each pair of adjacent voxels i and j, the square of the Euclidean distance between their eigenvectors is calculated: ,in and Represents the 12-dimensional feature vectors of voxels i and j respectively. Output feature vector difference data Feature_Difference.

[0191] The Gaussian kernel similarity is calculated using the feature vector difference data Feature_Difference. The Gaussian kernel function is defined as: , where σ is the standard deviation of the Gaussian kernel function, set to 5. Output Gaussian kernel similarity Gaussian_Similarity.

[0192] Use Gaussian kernel similarity Gaussian_Similarity and distance weight Distance_Weight to calculate the binary potential function value. The binary potential function is defined as: . Output the binary potential function Pairwise_Potential.

[0193] Therefore, the embodiments should be regarded as illustrative and non-restrictive from all points, and the scope of the present invention is limited by the appended claims rather than the above description, and it is therefore intended that all changes falling within the meaning and range of equivalent elements of the application documents are included in the present invention.

[0194] The above description is only a specific embodiment of the present invention, so that those skilled in the art can understand or implement the present invention. Various modifications to these embodiments will be apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to the embodiments shown herein, but should conform to the widest scope consistent with the principles and novel features invented herein.

Claims

1. A neural CT atlas segmentation method based on high-precision registration, characterized in that: The following steps are involved: Step S1: Obtaining a neural CT image, an MRI-T1 image, and an MRI-T2 image; performing a coarse registration based on mutual information on the CT image, the MRI-T1 image, and the MRI-T2 image, and performing a multi-scale fusion based on wavelet to obtain a fused image; Step S2: preprocessing the preset standard neural CT atlas library with standard atlases according to the fused image to obtain a preprocessed standard atlas set; adding a multi-scale feature fusion module to the pre-built basic U-Net architecture, and embedding the spatial attention module to obtain a U-Net network to be trained; using the fused image and the preprocessed standard atlas set to perform deep registration network training on the U-Net network to be trained to obtain a trained U-Net network; using the trained U-Net network to generate a precise registration field for the fused image to obtain a precise registration field; Step S3: Calculate the thin plate spline deformation parameters according to the fused image and the precise registration field to obtain the thin plate spline deformation parameters; Using thin plate spline deformation parameters to perform atlas deformation migration on the preprocessed standard atlas set, a deformation atlas is obtained; Calculate the prior probability according to the deformation spectrum to obtain a prior probability map; The posterior probability is calculated based on the prior probability map to obtain the tissue probability map; Step S4: constructing a univariate potential function for the tissue probability map to obtain a univariate potential function; extracting local features for the fused image to obtain a feature vector; defining a binary potential function based on the fused image and the feature vector to obtain a binary potential function; constructing a conditional random field energy function for the univariate potential function and the binary potential function to obtain a conditional random field energy function; The atlas segmentation result is generated according to the conditional random field energy function to obtain the atlas segmentation result to achieve the neural CT atlas segmentation task.

2. The neural CT atlas segmentation method based on high-precision registration according to claim 1, characterized in that: Step S1 includes the following steps: Step S11: acquiring a neural CT image, an MRI-T1 image, and an MRI-T2 image; performing image modality preprocessing on the CT image, the MRI-T1 image, and the MRI-T2 image to obtain a preprocessed CT image, a preprocessed MRI-T1 image, and a preprocessed MRI-T2 image; Step S12: using the pre-processed CT image as a reference image, performing coarse registration based on mutual information on the pre-processed MRI-T1 image and the pre-processed MRI-T2 image to obtain a registered MRI-T1 image and a registered MRI-T2 image; Step S13: performing wavelet-based multi-scale fusion on the pre-processed CT image, the registered MRI-T1 image, and the registered MRI-T2 image to obtain a fused image.

3. The neural CT atlas segmentation method based on high-precision registration according to claim 1, characterized in that: Step S2 includes the following steps: Step S21: performing standard atlas preprocessing on a preset standard neural CT atlas library according to the fusion image to obtain a preprocessed standard atlas set; Step S22: adding a multi-scale feature fusion module to the pre-built basic U-Net architecture and embedding a spatial attention module to obtain a U-Net network to be trained; Step S23: using the fused image and the preprocessed standard atlas set to perform deep registration network training on the U-Net network to be trained, to obtain a trained U-Net network; Step S24: Using the trained U-Net network to generate a precise registration field for the fused image to obtain a precise registration field.

4. The neural CT atlas segmentation method based on high-precision registration according to claim 3 is characterized in that: Step S22 includes the following steps: Step S221: extract encoder features from the pre-built basic U-Net architecture to obtain feature maps of each layer of the encoder; perform feature map channel adjustment on the feature maps of each layer of the encoder to obtain feature maps after channel adjustment; Step S222: performing channel attention weighting on the channel-adjusted feature map to obtain an attention-weighted feature map; performing feature map splicing according to the attention-weighted feature map and a pre-built basic U-Net architecture to obtain a U-Net with a multi-scale module added; Step S223: extract decoder features from the U-Net with the multi-scale module added to obtain feature maps of each layer of the decoder; extract local structure information from the feature maps of each layer of the decoder to obtain local structure information feature maps; Step S224: Generate a spatial attention map for the local structure information feature map to obtain a spatial attention map; perform feature map weighting according to the U-Net with the multi-scale module added, the feature maps of each layer of the decoder, and the spatial attention map to obtain a U-Net embedded with the spatial attention module; Step S225: Add a deformation field output layer to the U-Net embedded with the spatial attention module to obtain a U-Net network to be trained.

5. The neural CT atlas segmentation method based on high-precision registration according to claim 3 is characterized in that: Step S23 includes the following steps: Step S231: preparing training data for the fused image and the preprocessed standard atlas set to obtain training data; Step S232: performing image characteristic data enhancement of different modalities on the training data to obtain enhanced training data; Step S233: inputting the enhanced training data into the U-Net network to be trained to obtain a predicted deformation field; Step S234: deforming the preprocessed standard atlas set according to the predicted deformation field to obtain a deformed atlas; performing similarity measurement calculation on the deformed atlas and the fused image to obtain a similarity value; Step S235: calculating the smoothness of the predicted deformation field to obtain a smoothness value; Step S236: construct and optimize a joint loss function according to the similarity value, the enhanced training data, and the U-Net network to be trained to obtain a trained U-Net network.

6. The neural CT atlas segmentation method based on high-precision registration according to claim 1, characterized in that: Step S3 includes the following steps: Step S31: extracting atlas labels from the preprocessed standard atlas set to obtain a standard atlas label set; Step S32: Calculate the thin plate spline deformation parameters according to the fused image and the precise registration field to obtain the thin plate spline deformation parameters; Step S33: deforming the standard atlas label set using the thin plate spline deformation parameters to obtain a deformed atlas label set; Step S34: generating a deformation map for the standard map label set and the pre-processed standard map set to obtain a deformation map; Step S35: Calculate the prior probability according to the deformation spectrum to obtain a prior probability map; Step S36: constructing a grayscale distribution model of the deformation map using the grayscale information of the fused image to obtain a tissue grayscale distribution model; Step S37: Calculate the posterior probability of the fused image according to the prior probability map and the tissue grayscale distribution model to obtain a tissue probability map.

7. The neural CT atlas segmentation method based on high-precision registration according to claim 6, characterized in that: Step S32 includes the following steps: Step S321: performing image gradient calculation on the fused image according to the precise registration field to obtain image gradient data; Step S322: Calculate the gradient amplitude according to the image gradient data to obtain the gradient amplitude; Step S323: selecting control points for the precise registration field to obtain a control point set; Step S324: Calculate the control point weights of the control point set according to the gradient amplitude to obtain weighted control points; Step S325: Calculate thin plate spline parameters using weighted control points and precise registration fields to obtain thin plate spline deformation parameters.

8. The neural CT atlas segmentation method based on high-precision registration according to claim 1, characterized in that: Step S4 includes the following steps: Step S41: constructing a univariate potential function for the tissue probability map to obtain a univariate potential function; Step S42: performing multi-scale Gaussian filtering on the fused image to obtain a multi-scale image; performing multi-scale image gradient calculation on the multi-scale image to obtain multi-scale gradient information; Step S43: performing local binary pattern calculation on the multi-scale image to obtain multi-scale LBP features; constructing feature vectors for the multi-scale gradient information and the multi-scale LBP features to obtain feature vectors; Step S44: defining a binary potential function according to the fused image and the feature vector to obtain a binary potential function; Step S45: constructing a conditional random field energy function for the one-dimensional potential function and the two-dimensional potential function to obtain a conditional random field energy function; Step S46: Optimizing the conditional random field energy function to obtain an optimized label graph; Step S47: Generate a graph segmentation result for the optimized label graph to obtain a graph segmentation result.

9. The neural CT atlas segmentation method based on high-precision registration according to claim 8, characterized in that: Step S44 includes the following steps: Step S441: extracting voxel coordinate information from the fused image to obtain voxel coordinate information; calculating the distance between voxels based on the voxel coordinate information to obtain the distance between voxels; Step S442: Calculate the distance weight according to the distance between voxels to obtain the distance weight; Step S443: Calculate the feature vector difference according to the feature vector and the voxel coordinate information to obtain feature vector difference data; Step S444: performing Gaussian kernel function calculation according to the feature vector difference data to obtain Gaussian kernel similarity; Step S445: Calculate the binary potential function value according to the Gaussian kernel similarity and the distance weight to obtain the binary potential function.

Images