A medical image segmentation method based on deep learning

Abstract

The application discloses a medical image segmentation method based on deep learning, and relates to the technical field of medical image data processing and segmentation, which comprises the following steps: obtaining a body mask based on three-dimensional gray data and constructing a radial distance map; calculating the median in the mask according to the pixel level equal-width ring band and performing logarithmic linear fitting to obtain a radial low-frequency slope sequence and an axial radial low-frequency intermodulation index; generating an intermodulation compensation field according to the slope and the distance and correcting the original image by multiplication; calculating the radial and slice direction gradients to obtain an anisotropic boundary uncertainty ratio and a local anisotropy map; and constructing a position adaptive threshold field and completing segmentation according to the anisotropic boundary uncertainty ratio and the local anisotropy map. The application can improve the stability and cross-device consistency of medical image segmentation.

Description

Technical Field

[0001] This invention relates to the field of medical image data processing and segmentation technology, and in particular to a medical image segmentation method based on deep learning. Background Technology

[0002] Medical image segmentation is a core technology in medical research. Its function is to accurately separate target tissue from the background in medical images such as CT and MRI, providing quantitative evidence for disease diagnosis, treatment planning, and prognostic assessment. With the development of imaging technology, the resolution and modal diversity of medical images have been continuously improved. However, the problem of uneven gray-level distribution is common in images. Due to the physical characteristics of imaging equipment, tissue peristalsis, or local volume effects, the same tissue exhibits significant gray-level differences in different regions. At the same time, the gradient changes of the boundary between the target and the background differ in different directions, i.e., boundary gradient anisotropy. These characteristics make accurate segmentation a long-standing challenge for the industry.

[0003] Existing segmentation techniques have significant limitations in processing 3D medical images, failing to meet the clinical demands for high-precision segmentation. Traditional thresholding methods, which separate targets from background using a single fixed threshold, are simple to operate but cannot address the prevalent uneven grayscale distribution in 3D images. Grayscale values ​​within the same target area can fluctuate significantly due to signal deviations in imaging equipment, such as inhomogeneous magnetic resonance radio frequency fields, or physiological differences in the tissue itself, leading to some areas being misclassified as background or vice versa. Region-growing methods rely on the selection of initial seed points; if the seed point is located in a grayscale transition zone, it can easily lead to overexpansion or underexpansion of the growth range, and it is difficult to accurately define boundaries with smooth grayscale transitions. While edge detection-based methods can capture local grayscale changes, they are limited by the anisotropy of gradients at 3D image boundaries. Differences in sampling resolution in different spatial directions and the directional distribution of tissue texture can cause significant differences in gradient intensity in different directions, resulting in broken and discontinuous edge detection results that fail to form a complete target contour. These technical shortcomings collectively lead to oversegmentation, undersegmentation, or boundary offsets in the segmentation results of existing methods, making it difficult to meet the stringent clinical requirements for segmentation accuracy. Summary of the Invention

[0004] The purpose of this invention is to address the problem in the prior art that it is difficult to effectively identify and correct low-frequency intensity fluctuations that vary with spatial location, resulting in inconsistent directionality and weak stability in the boundary segmentation results. Therefore, this invention proposes a medical image segmentation method based on deep learning.

[0005] To address the problems existing in the prior art, the present invention adopts the following technical solution:

[0006] A deep learning-based medical image segmentation method includes:

[0007] A body mask is generated based on 3D grayscale data, and a radial distance map is generated based on the body mask.

[0008] The median was statistically analyzed in the radial annulus within the body mask, and the logarithm of the median was linearly fitted to obtain the radial low-frequency slope sequence. The axial low-frequency intermodulation index was determined based on the radial low-frequency slope sequence.

[0009] The three-dimensional grayscale data is corrected based on the radial distance map and the radial low-frequency slope sequence to obtain the image after intermodulation compensation;

[0010] Calculate the radial gradient and slice gradient on the image after intermodulation compensation. Based on the radial gradient and slice gradient, calculate the anisotropic boundary uncertainty ratio and form a local anisotropic map.

[0011] Threshold segmentation is performed on the intermodulation-compensated image based on the radial low-frequency intermodulation index, the anisotropic boundary uncertainty ratio, the sign of the radial low-frequency slope sequence, and the local anisotropic map.

[0012] Preferably, generating a body mask based on three-dimensional grayscale data includes:

[0013] The Otsu threshold is obtained by processing slices in three-dimensional grayscale data using the Otsu threshold method.

[0014] Threshold segmentation of the slice grayscale is performed based on the Otsu method threshold to generate an initial binary region;

[0015] Select the connected region with the largest number of pixels from all connected regions of the initial binary region as the body mask of the slice.

[0016] Preferably, generating a radial distance map based on the body mask includes:

[0017] Calculate the centroid of the slice within the body mask of the slice according to the definition of geometric centroid;

[0018] Using the centroid of the slice as the center, calculate the Euclidean distance from each pixel within the slice to the centroid of the slice to obtain the radial distance map of the slice.

[0019] Preferably, the median is statistically analyzed along radial rings within the body mask, and the logarithm of the median is linearly fitted to obtain a radial low-frequency slope sequence, including:

[0020] Under the constraints of the body mask and radial distance map of each slice, a ring set is constructed according to the pixel-level equal width layer;

[0021] Calculate the median gray value within each annular band set to obtain the annular band median sequence of the slice;

[0022] The radial low-frequency slope of the slice is obtained by shifting the median sequence of the ring band in the positive domain and performing a linear least squares fitting in the logarithmic domain.

[0023] The radial low-frequency slopes of all slices are combined into a radial low-frequency slope sequence.

[0024] Preferably, determining the shaft diameter low-frequency intermodulation index based on the radial low-frequency slope sequence includes:

[0025] Calculate the discrete standard deviation of the radial low-frequency slope sequence;

[0026] Calculate the median absolute deviation of the 3D grayscale data within the body mask;

[0027] The low-frequency intermodulation index of the shaft diameter is determined by the ratio of the discrete standard deviation to the median absolute deviation.

[0028] Preferably, the three-dimensional grayscale data is corrected based on the radial distance map and the radial low-frequency slope sequence to obtain an intermodulation compensated image, including:

[0029] Based on the radial low-frequency slope sequence and radial distance map, an intermodulation compensation field of the same size as the 3D grayscale data is constructed. The formula for constructing the intermodulation compensation field is as follows:

[0030]

[0031] In the formula, For inter-compensation fields, Radial low-frequency slope, This is a pixel-level radial distance map. and These are the pixel coordinates within the slice. The slice number;

[0032] The three-dimensional grayscale data is corrected by point-by-point multiplication using the intermodulation compensation field to obtain the intermodulation compensated image.

[0033] Preferably, the radial gradient and slice gradient are calculated on the intermodulation-compensated image. Based on the radial gradient and slice gradient, the anisotropic boundary uncertainty ratio is calculated and a local anisotropic map is formed, including:

[0034] Within each slice, a radial unit vector field is defined centered on the slice's centroid. The formula for the radial unit vector field is as follows:

[0035]

[0036] In the formula, For radial unit vector fields, and These are the pixel coordinates within the slice. The slice number. and Here are the coordinates of the slice's centroid. This is a pixel-level radial distance map. To ensure numerical stability and avoid division by zero;

[0037] Spatial gradient derivative and slice directional derivative of the image after intermodulation compensation are calculated using finite difference method.

[0038] Projecting the spatial gradient derivative onto the radial unit vector field direction yields the radial gradient intensity.

[0039] The absolute value of the derivative of the slice direction is taken to obtain the gradient intensity of the slice direction;

[0040] The arithmetic mean of the gradient intensity in the slice direction and the radial direction is calculated separately within the body mask;

[0041] The anisotropic boundary uncertainty ratio is obtained by the ratio of the arithmetic mean of the gradient intensity in the slice direction to the arithmetic mean of the gradient intensity in the radial direction.

[0042] A local anisotropy map is constructed using the normalized difference between the gradient intensity in the slice direction and the gradient intensity in the radial direction.

[0043] Preferably, threshold segmentation is performed on the intermodulation-compensated image based on the radial low-frequency intermodulation index, the anisotropic boundary uncertainty ratio, the sign of the radial low-frequency slope sequence, and the local anisotropic map, including:

[0044] Within the body mask of each slice, the intermodulation-compensated image is processed using the Otsu method thresholding to obtain the slice baseline threshold.

[0045] The median absolute deviation of the intermodulation-compensated image is calculated within the body mask of each slice to obtain the slice intensity scale.

[0046] Based on the slice baseline threshold, slice intensity scale, axial low-frequency intermodulation index, anisotropic boundary uncertainty ratio, sign of the radial low-frequency slope sequence, and local anisotropy map, the location-adaptive threshold field is calculated. The formula for calculating the location-adaptive threshold field is as follows:

[0047]

[0048] In the formula, For location-adaptive threshold field, The slice baseline threshold, The low-frequency intermodulation index of the shaft diameter. The uncertainty ratio of the anisotropic boundary. The sign for the radial low-frequency slope. This is a local anisotropy diagram. For slice intensity scale, and These are the pixel coordinates within the slice. The slice number;

[0049] Threshold segmentation is performed on the intermodulation-compensated image based on the position-adaptive threshold field, and the segmentation result is obtained by AND operation with the body mask.

[0050] The segmentation results are evaluated by a pre-defined deep learning reviewer to obtain a score map of suspicious points.

[0051] Compared with the prior art, the beneficial effects of the present invention are:

[0052] 1. This invention uses a body mask and radial distance map as a geometric reference, statistically analyzes the median of equal-width rings within the mask, and performs logarithmic linear fitting to obtain a slice-level radial low-frequency slope sequence. The ratio of the discrete standard deviation of this sequence to the absolute deviation of the median in the body is used to form the axial low-frequency intermodulation index. Then, an intermodulation compensation field of the same size as the volume data is constructed to correct the original image. Mechanistically, it quantitatively identifies and physically consistent reverse corrections the slow intensity fluctuations that exist synchronously in the radial and slice directions, significantly reducing gain differences across devices and protocols, and unifying the contrast and background levels of similar tissues at different locations.

[0053] 2. This invention projects the spatial gradient as a radial gradient onto the image after intermodulation compensation using a radial unit vector field, and forms two orthogonal edge metrics with the slice directional derivative. The anisotropic boundary uncertainty ratio is obtained by averaging the intensity of the two directions within the body mask, and a local anisotropic map is obtained by normalizing the difference between the two directions. Finally, the position-adaptive threshold field is formed by introducing the low-frequency intermodulation index of the axis radius, the directional ratio, the slope sign, and the local distribution, using the slice baseline threshold as the anchor point. This allows the threshold to automatically shift where needed and fall back where it is stable, thereby reducing boundary shift and breakage and improving 3D coherence and robustness.

[0054] 3. This invention primarily uses deterministic statistics and geometric quantities, without introducing multiple hyperparameters that rely on experience. The computational domain is limited to within the body mask to suppress interference from the bed surface and air. Furthermore, a deep learning reviewer is used to output a question mark scoring map as a quality control prompt, enabling the system to work stably under default conditions and possessing reviewable and verifiable verification capabilities in complex scenarios. Overall, it reduces repeated manual threshold adjustments and improves consistency and usability across datasets. Attached Figure Description

[0055] The accompanying drawings, which are included to provide a further understanding of the invention and form part of this application, illustrate exemplary embodiments of the invention and, together with their description, serve to explain the invention and do not constitute an undue limitation thereof. In the drawings:

[0056] Figure 1 This is a flowchart illustrating a deep learning-based medical image segmentation method according to an embodiment of the present invention.

[0057] Figure 2 This is a schematic flowchart of a body mask generation method provided in an embodiment of the present invention;

[0058] Figure 3 This is a flowchart illustrating a method for generating radial low-frequency slope sequences according to an embodiment of the present invention.

[0059] Figure 4 This is a flowchart illustrating a method for calculating the low-frequency intermodulation index of shaft diameter according to an embodiment of the present invention.

[0060] Figure 5 This is a schematic diagram of a method for constructing a local anisotropy graph according to an embodiment of the present invention; Detailed Implementation

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

[0062] Example: This example provides a deep learning-based medical image segmentation method. See [link to example]. Figure 1 Specifically, including:

[0063] A body mask is generated based on 3D grayscale data, and a radial distance map is generated based on the body mask.

[0064] In an embodiment of the present invention, generating a body mask based on three-dimensional grayscale data and generating a radial distance map based on the body mask includes:

[0065] The Otsu threshold is obtained by processing slices in three-dimensional grayscale data using the Otsu threshold method.

[0066] Threshold segmentation of the slice grayscale is performed based on the Otsu method threshold to generate an initial binary region;

[0067] Select the connected region with the largest number of pixels from all connected regions of the initial binary region as the body mask of the slice. See [link / reference] Figure 2 ;

[0068] Specifically, 3D grayscale data is a volume of 3D image data formed by multiple consecutively acquired images and superimposed in spatial order. It uses voxels as the basic unit, and each voxel corresponds to a grayscale value that represents the strength of tissue signals. The range of grayscale values ​​is determined by the imaging device and reconstruction method. It also carries spatial scale information such as pixel spacing and layer spacing for subsequent geometric measurement and intensity calculation. A slice refers to a 2D image obtained layer by layer from 3D grayscale data along the direction perpendicular to the superposition. Each slice is composed of pixels in a regular grid, has a defined row and column size and pixel spacing, and has a defined layer spacing or thickness with adjacent slices.

[0069] Specifically, the segmentation thresholding of the 3D grayscale data is performed slice by slice using the Otsu method. This is to adaptively determine the segmentation threshold when the grayscale distribution of different slices is inconsistent, thereby stably distinguishing the foreground from the background and avoiding misclassification caused by fixed thresholds under noise, scattering, or changes in equipment imaging conditions. The initial binary region is obtained by thresholding the grayscale of the slices based on this threshold. This is to quickly form candidate human regions with minimal prior information and retain complete contour information, which is convenient for subsequent geometric and statistical calculations. The connected region with the largest number of pixels in all connected regions of the initial binary region is selected as the body mask of the slice because the human body usually occupies the largest area and has the strongest connectivity in the field of view in medical images. This selection can eliminate examination beds, markers, air areas, and scattered artifacts in one go, reducing the interference of irrelevant regions on subsequent calculations. This ensures that the subsequent centroid calculation and radial distance map generation are based on a stable and reliable anatomical range, improving the repeatability and determinism of the segmentation process and reducing computational overhead.

[0070] Specifically, the 3D grayscale data consists of multiple 2D slices arranged sequentially along the axis, each slice containing grayscale information of all pixels in the 2D plane. For a single slice, the grayscale value distribution of all pixels within the slice is first statistically analyzed to construct a grayscale histogram. The horizontal axis of the histogram represents the grayscale value range, and the vertical axis represents the number of pixels within the corresponding grayscale value range. Based on this grayscale histogram, the inter-class variance under different grayscale thresholds is calculated. The inter-class variance is calculated by dividing the pixels within the slice into a foreground pixel set and a background pixel set, using the candidate threshold as the boundary. The grayscale mean and pixel proportion of the two sets are calculated respectively, thus obtaining the inter-class variance value. All possible grayscale thresholds are iterated, and the grayscale threshold that maximizes the inter-class variance is selected as the Otsu's threshold corresponding to that slice, completing the determination of the Otsu's threshold for a single slice. All slices in the 3D grayscale data are processed sequentially according to the above process to obtain the Otsu's threshold corresponding to each slice.

[0071] Specifically, for each pixel within a slice, its grayscale value is compared with the Otsu threshold for that slice. If the pixel's grayscale value is greater than or equal to the Otsu threshold, the pixel is marked as a first value; if the pixel's grayscale value is less than the Otsu threshold, the pixel is marked as a second value. The first and second values ​​are two distinct discrete values. Through this marking operation, the slice, which originally contained continuous grayscale variations, is transformed into a binary image containing only two values. This binary image is the initial binary region.

[0072] Specifically, the four-neighbor connectivity criterion is used to determine connected regions within the initial binary region. This criterion considers pixels centered at a single pixel and adjacent pixels in all four directions (up, down, left, and right) as neighboring pixels. If adjacent pixels have the same value, they are considered part of the same connected region. Based on this criterion, all pixels within the initial binary region are traversed to identify all connected regions. Each connected region consists of a series of adjacent pixels with the same value. The number of pixels in each connected region is counted, and the number of pixels in all connected regions is compared. The connected region with the largest number of pixels is selected as the body mask for the corresponding slice. This process is repeated to construct the body masks for all slices in the 3D grayscale data.

[0073] Calculate the centroid of the slice within the body mask of the slice according to the definition of geometric centroid;

[0074] Using the centroid of the slice as the center, calculate the Euclidean distance from each pixel within the slice to the centroid of the slice to obtain the radial distance map of the slice;

[0075] Specifically, the centroid of the slice is calculated within the body mask based on the geometric centroid definition. This is to select a spatial reference point that is insensitive to rotation and translation to represent the position of the human body, thereby avoiding systematic bias caused by using the device coordinate origin or field boundary as a reference. The Euclidean distance from each pixel to the centroid is calculated with this centroid as the center, forming a radial distance map. This yields a scalar field that monotonically increases outward from the center, providing a uniform geometric scale for all pixels without relying on grayscale information. This ensures that subsequent annular statistics maintain consistent region division across slices of different shapes and poses, and provides a clear radial direction for the calculation of directional gradients. This, in turn, stably supports low-frequency trend estimation, directional difference measurement, and threshold field construction. At the same time, limiting the calculation range with the body mask can eliminate interference from irrelevant areas such as the bed surface and air on the centroid and distance, improving the stability and repeatability of the calculation results. Furthermore, the Euclidean distance has low computational overhead and is simple to implement, making it easy to reuse efficiently across multiple slices.

[0076] Specifically, for each slice's body mask, the two-dimensional coordinates of all pixels within the body mask are first obtained. Each pixel's two-dimensional coordinates consist of its horizontal and vertical coordinates within the slice plane. Only the coordinate information corresponding to pixels covered by the body mask is extracted, excluding pixels outside the body mask. Then, the arithmetic mean of the horizontal coordinates of all pixels within the body mask is calculated; this mean is the horizontal coordinate of the slice's centroid. Simultaneously, the arithmetic mean of the vertical coordinates of all pixels within the body mask is calculated; this mean is the vertical coordinate of the slice's centroid. The horizontal and vertical coordinates together constitute the slice's geometric centroid. The calculation of the centroid of each slice is completed sequentially according to the above process.

[0077] Specifically, using the determined centroid of the slice as a reference point, distance calculations are performed on all pixels within the slice. For a single pixel within the slice, first, the horizontal and vertical coordinates of the pixel are obtained. Then, the horizontal and vertical coordinates of the slice centroid are retrieved. The difference between the pixel's horizontal coordinate and the centroid's horizontal coordinate is calculated, as is the difference between the pixel's vertical coordinate and the centroid's vertical coordinate. These two differences are squared, and the sums are then taken as the square root to obtain the Euclidean distance from the pixel to the slice centroid. This process is repeated to calculate the Euclidean distance from all pixels within the slice to the centroid. The Euclidean distance of each pixel is used as its corresponding distance value. All pixel distance values ​​are arranged according to their coordinate positions within the slice, forming a two-dimensional distance matrix with the same size as the slice. This two-dimensional distance matrix is ​​the radial distance map of the slice.

[0078] The median was statistically analyzed in the radial annulus within the body mask, and the logarithm of the median was linearly fitted to obtain the radial low-frequency slope sequence. The axial low-frequency intermodulation index was determined based on the radial low-frequency slope sequence.

[0079] In an embodiment of the invention, the median is statistically analyzed along radial rings within the body mask, and a linear fit is performed on the logarithm of the median to obtain a radial low-frequency slope sequence. The axial low-frequency intermodulation index is then determined based on this radial low-frequency slope sequence. (See [reference]). Figure 3 and Figure 4 ,include:

[0080] Under the constraints of the body mask and radial distance map of each slice, a ring set is constructed according to the pixel-level equal width layer;

[0081] Calculate the median gray value within each annular band set to obtain the annular band median sequence of the slice;

[0082] Specifically, under the constraints of the body mask and radial distance map for each slice, a set of annular bands is constructed in pixel-level equal-width layers. This is to form a geometric hierarchy that monotonically increases from the inside out along the radial direction with the slice centroid as the center, so that pixels within the same radius band have consistent spatial meaning and facilitate direct comparison between different slices. Using equal-width layers avoids manually setting bandwidth or scale parameters, reducing subjectivity and improving the consistency of reproduction experiments. Calculating the gray-level median instead of the mean within each annular band set is to suppress the influence of a small number of high-contrast structures and metal artifacts on the statistical results, making the obtained statistics robust to noise and outliers. The annular band median sequence formed from the inside out can completely reflect the trend of gray-level slow change with radius, providing a stable data foundation for subsequent logarithmic domain linear fitting. This allows for a reliable measure of radial low-frequency change without introducing additional parameters, thus ensuring the comparability and stability of subsequent intermodulation compensation and threshold calculation. At the same time, this process is only performed within the body mask, which can eliminate interference from irrelevant areas such as air regions and bed surfaces, reduce computational overhead, and improve overall processing efficiency.

[0083] Specifically, for each slice, a ring set construction operation is performed under the constraints of its body mask and radial distance map. First, the maximum distance value of the radial distance map within the body mask of the slice is obtained, which is the maximum radial distance among the pixels covered by the body mask. Based on a preset number of layers, the maximum distance value is divided into several equally wide distance intervals. The width of each distance interval is the ratio of the maximum distance value to the number of layers, and adjacent intervals are connected end-to-end to form a continuous distance range sequence. Specifically, for each distance interval, all pixels within the body mask whose radial distance belongs to that interval are selected. These pixels together constitute a ring. Among them, the ring with a starting value of 0 in the distance interval is the innermost ring, and the rings are distributed from the inside out as the starting value of the distance interval increases. All rings obtained in the above manner are arranged in ascending order of radial distance to form the ring set of the slice. For the ring set of each slice, the median gray value corresponding to each ring is calculated sequentially. For a single ring, the gray values ​​of all pixels within the ring in the original slice are extracted, and the distribution of these gray values ​​is statistically analyzed. Arrange all grayscale values ​​in ascending order. If the number of grayscale values ​​is odd, select the middle grayscale value as the median of that ring. If the number of grayscale values ​​is even, select the arithmetic mean of the two middle grayscale values ​​as the median of that ring. Arrange the medians of each ring in the ring set from the inside out, forming an ordered numerical sequence. This sequence is the median sequence of the rings for that slice. Process all slices sequentially using the above process to obtain the median sequence of the rings for each slice.

[0084] The radial low-frequency slope of the slice is obtained by shifting the median sequence of the ring band in the positive domain and performing a linear least squares fitting in the logarithmic domain.

[0085] The radial low-frequency slopes of all slices are combined into a radial low-frequency slope sequence;

[0086] Specifically, shifting the median sequence of the ring bands into the positive domain ensures that the independent variable used for subsequent logarithmic calculations is strictly positive, thus avoiding numerical anomalies and ensuring comparability of slices at the same scale. Performing a linear least-squares fit in the logarithmic domain is beneficial because slow multiplicative fluctuations can be transformed into additive linear trends in the logarithmic domain. Linear fitting can stably characterize the slow change direction and magnitude of grayscale with radius using a single slope, and a definite solution is obtained through least squares to reduce computational uncertainty. The resulting radial low-frequency slope directly represents the systematic change of the slice from the center outwards, providing a quantitative basis for subsequent compensation mapping construction and threshold setting. Compiling the radial low-frequency slopes of all slices into a radial low-frequency slope sequence allows for unified statistical measurement and comparison at the three-dimensional level, maintaining a consistent scale across slices and providing a unique and clear scale of change for subsequent steps.

[0087] Specifically, the minimum value in the median sequence of the annular zones is determined by iterating through all values ​​in the sequence. A translation constant is then set based on this minimum value. The translation constant must satisfy the condition that all values ​​in the translated sequence are greater than zero, and the translation constant must be a fixed, small positive value to avoid subsequent logarithmic operations failing due to zero or negative values. Each value in the median sequence of the annular zones is added to the translation constant to obtain the translated median sequence of the annular zones, completing the positive domain translation. A logarithmic domain transformation is then performed on the translated median sequence of the annular zones. For each translated value, the natural logarithm is performed to obtain the corresponding logarithmic value. All logarithmic values ​​are arranged in the original order of the median sequence of the annular zones to form the logarithmic median sequence of the annular zones. Simultaneously, the radial distance mean for each annular zone is obtained. This mean is the arithmetic mean of the upper and lower boundary values ​​of the distance interval to which each annular zone belongs, forming a radial distance mean sequence in order of the annular zones from the inside out. Linear least squares fitting is then performed based on the logarithmic median sequence of the annular zones and the radial distance mean sequence. The linear fitting model is defined as follows: the logarithmic value equals the intercept term plus the radial low-frequency slope multiplied by the radial distance mean, where the radial low-frequency slope is the core parameter to be solved. By calculating the covariance of the two sequences, the variance of the radial distance mean sequence, and the mean of the two sequences, and substituting these statistics into the linear least squares fitting formula, the intercept term and radial low-frequency slope that minimize the sum of squares of the fitting error are obtained. The sum of squares of the fitting error is the sum of the squares of the differences between each value in the logarithmic ring median sequence and the predicted value of the fitting model; the resulting slope parameter is the radial low-frequency slope of that slice. Following the axial arrangement order of the slices in the 3D grayscale data, the radial low-frequency slopes calculated for each slice are sequentially arranged. Starting from the radial low-frequency slope corresponding to the first slice at the axial starting position, the radial low-frequency slopes of all subsequent slices are sequentially connected to form a complete ordered numerical sequence, which is the radial low-frequency slope sequence.

[0088] Calculate the discrete standard deviation of the radial low-frequency slope sequence;

[0089] Calculate the median absolute deviation of the 3D grayscale data within the body mask;

[0090] The low-frequency intermodulation index of the shaft diameter is determined based on the ratio of the discrete standard deviation to the median absolute deviation.

[0091] Specifically, the discrete standard deviation is first calculated for the sequence composed of the radial low-frequency slopes of each slice. This is to quantify the overall fluctuation of the slope in the slice direction. The larger the value, the more obvious the difference in radial trend between different slices, which is more likely to produce inconsistent boundary performance in the slice direction and radial direction. Using only the discrete standard deviation is affected by the overall grayscale scale and noise level, making it difficult to directly compare between different devices and different individuals. Therefore, the median absolute deviation of the 3D grayscale data is calculated in the body region defined by the body mask. This median absolute deviation is used as a robust intensity scale benchmark to reflect the typical grayscale fluctuation intensity of the body data itself. The ratio of the discrete standard deviation to the median absolute deviation is used to obtain a dimensionless full-volume scale. This scale can automatically cancel out gain differences and unit differences, reduce the interference of abnormally high-intensity structures and metal artifacts on the results, and adaptively scale with the statistical characteristics of the data itself. In the subsequent threshold field calculation, this scale is used as a weighting factor, so that the threshold offset only plays a role when the fluctuation is significant, and naturally falls back to the slice baseline threshold when the fluctuation is slight. This improves the stability and cross-scene consistency of the segmentation results without adding additional parameters.

[0092] Specifically, a constructed radial low-frequency slope sequence is obtained, which contains the radial low-frequency slopes corresponding to each slice in the 3D grayscale data. All slope values ​​in the sequence are traversed, and the arithmetic mean of the sequence is calculated. The arithmetic mean is the ratio of the sum of all slope values ​​to the number of values ​​in the sequence. For each slope value, the difference between it and the arithmetic mean is calculated, and each difference is squared to obtain a series of squared differences. The arithmetic mean of all squared differences is calculated, and then the square root operation is performed on the arithmetic mean. The result is the discrete standard deviation of the radial low-frequency slope sequence. The grayscale values ​​of all pixels covered by the body mask in the 3D grayscale data are collected to form a grayscale value set. This set only contains the original grayscale information corresponding to pixels within the body mask of each slice, excluding the grayscale values ​​of pixels outside the body mask. All values ​​in the grayscale value set are arranged in ascending order. If the number of values ​​in the set is odd, the value in the middle position is selected as the median of the grayscale value set; if the number of values ​​is even, the arithmetic mean of the two middle values ​​is selected as the median. For each value in the grayscale value set, the difference between it and the median is calculated, and the absolute value of each difference is taken to obtain a series of absolute deviation values. All absolute deviation values ​​are arranged in ascending order, and the median selection rule is repeated to obtain the median of the absolute deviation value set. This median is the median absolute deviation of the 3D grayscale data within the body mask. The discrete standard deviation of the previously calculated radial low-frequency slope sequence is used as the numerator, and the median absolute deviation of the 3D grayscale data within the body mask is used as the denominator. A division operation is performed, and the quotient obtained is the radial low-frequency intermodulation index. This index quantifies the relative relationship between the dispersion of the radial low-frequency slope and the background grayscale fluctuation of the image, reflecting the intensity level of radial low-frequency intermodulation.

[0093] The three-dimensional grayscale data is corrected based on the radial distance map and the radial low-frequency slope sequence to obtain the image after intermodulation compensation;

[0094] In an embodiment of the present invention, three-dimensional grayscale data is corrected based on a radial distance map and a radial low-frequency slope sequence to obtain an image after intermodulation compensation, including:

[0095] Based on the radial low-frequency slope sequence and radial distance map, an intermodulation compensation field of the same size as the three-dimensional grayscale data is constructed.

[0096] Specifically, based on the radial low-frequency slope sequence and radial distance map, an intermodulation compensation field of the same size as the 3D grayscale data is constructed. The formula for constructing the intermodulation compensation field is as follows:

[0097]

[0098] In the formula, For inter-compensation fields, Radial low-frequency slope, This is a pixel-level radial distance map. and These are the pixel coordinates within the slice. The slice number. It is an exponential function;

[0099] Specifically, medical imaging generates a multiplicative low-frequency term during illumination and reconstruction that varies slowly with the radius and slice size. The observed grayscale can be considered as the product of this low-frequency term and the intrinsic grayscale of the tissue. Taking the logarithm of the multiplicative term transforms it into an additive term, whose first-order approximation along the radius within the slice is a constant term plus the product of the radial low-frequency slope and the radial distance map. Therefore, when performing inverse correction in the intensity domain, the negative exponential value of this additive term should be multiplied. That is, the intermodulation compensation field is equal to the negative radial low-frequency slope of the exponential function multiplied by the radial distance map. Thus, as the radius increases, the intensity of the periphery is automatically suppressed or enhanced based on the sign of the radial low-frequency slope. The radial distance map at the center position is equal to zero, making the compensation factor one, ensuring that the central grayscale is not changed. The exponential mapping naturally ensures that the compensation factor is positive, which conforms to the physical constraint of non-negative grayscale and mainly acts on the low-frequency trend without destroying the high-frequency edge. This calculation form is consistent with the general expression of Bell-Lambert's law, which describes that the signal changes exponentially with the path length. Therefore, exponential inversion can obtain intensity correction consistent with physical laws. Based on the coordinates of all pixels in the 3D grayscale data, the intermodulation compensation field value of each pixel is calculated sequentially, ultimately forming a 3D compensation field matrix with the same size as the 3D grayscale data. This matrix is ​​the intermodulation compensation field, and the value of each element is determined by the above calculation method.

[0100] The three-dimensional grayscale data is corrected by point-by-point multiplication using the intermodulation compensation field to obtain the intermodulation compensated image.

[0101] Specifically, the intermodulation compensation field provides an inverse vector of the linear trend in the logarithmic domain in exponential form, which means a slowly changing gain field. Multiplying this gain field point by point with the 3D grayscale data at voxel locations is equivalent to subtracting the low-frequency additive drift from the observed grayscale in the logarithmic domain, thus achieving a direct inversion of multiplicative inhomogeneities. This point-by-point multiplication method has three advantages: First, the compensation factor is always positive, which can maintain the physical non-negativity and dimensional consistency of grayscale and avoid the numerical instability that is prone to occur in division. Second, the compensation field changes slowly with space and only affects low-frequency components, which can correct the overall brightness gradient while preserving high-frequency edges and texture details as much as possible, providing a more stable contrast background for subsequent gradient calculation and thresholding. Third, the point-by-point operation with the same size as the 3D grayscale data avoids intensity mixing across regions, ensuring that each voxel is only affected by the trend correction of its current position, reducing artifact diffusion near the boundary and improving the consistency and repeatability of the whole-volume processing.

[0102] Specifically, the process involves acquiring the constructed intermodulation compensation field and the original 3D grayscale data. The intermodulation compensation field is a 3D compensation field matrix of the same size as the 3D grayscale data, with each element corresponding to the intermodulation compensation value of a pixel. The original 3D grayscale data contains the original grayscale information of all pixels within the slice. For each pixel in the 3D grayscale data, its coordinates in 3D space are determined. Based on these coordinates, the corresponding values ​​in the intermodulation compensation field and the corresponding grayscale values ​​in the original 3D grayscale data are retrieved. The intermodulation compensation field values ​​and the original grayscale values ​​are then multiplied point-by-point to obtain the corrected grayscale value of the pixel. Following the coordinate order of all pixels in the 3D grayscale data, the point-by-point multiplication correction operation is performed on each pixel sequentially. All corrected grayscale values ​​are then recombine according to their coordinate positions to form a 3D image matrix of the same size as the original 3D grayscale data. This matrix is ​​the intermodulation-compensated image.

[0103] Calculate the radial gradient and slice gradient on the image after intermodulation compensation. Based on the radial gradient and slice gradient, calculate the anisotropic boundary uncertainty ratio and form a local anisotropic map.

[0104] In an embodiment of the present invention, the radial gradient and slice gradient are calculated on the intermodulation-compensated image. Based on the radial gradient and slice gradient, the anisotropic boundary uncertainty ratio is calculated and a local anisotropic map is formed. See [link to relevant documentation]. Figure 5 ,include:

[0105] Define a radial unit vector field centered on the centroid of each slice;

[0106] Specifically, a radial unit vector field is defined within each slice, centered on the slice's centroid. The formula for the radial unit vector field is as follows:

[0107]

[0108] In the formula, For radial unit vector fields, and These are the pixel coordinates within the slice. The slice number. and Here are the coordinates of the slice's centroid. This is a pixel-level radial distance map. To ensure numerical stability and avoid division by zero, To select the largest operator, choose and The larger of the two;

[0109] Specifically, the radial unit vector field is obtained by normalizing the vector pointing from the centroid of the slice to the pixel position according to the Euclidean norm. It belongs to the unit direction vector required in the theory of directional derivatives and can be used to project the spatial gradient onto the physical direction from the center outward and calculate the directional derivative along the radius. This construction sets the third component to zero in the slice plane, geometrically excluding the slice direction component and ensuring that only radial changes in the plane are measured. The denominator takes the larger value of the radial distance and the numerically stable constant to avoid division by zero at the centroid and maintain computability in the whole domain, while ensuring that the vector magnitude does not exceed one to meet the normalization requirement. Using the slice centroid as a reference point can obtain a directional reference that is insensitive to translation and rotation, so that the radial definition is consistent under different slices and different body positions. When the distance is measured in pixels or millimeters, the numerator and denominator have the same dimensions, and the normalized result is a dimensionless quantity, which is convenient for dot product with the gradient and obtaining comparable radial gradient strength. This definition not only conforms to the common formula of Euclidean geometry and directional derivatives, but is also numerically stable and can be directly connected with subsequent gradient projection and anisotropy measurement. Calculate the radial unit vector of all pixels within the slice sequentially, ultimately forming a two-dimensional vector field with the same size as the slice. This vector field is the radial unit vector field of that slice. Repeat the above operation for all slices to complete the construction of the radial unit vector field for each slice.

[0110] Spatial gradient derivative and slice directional derivative of the image after intermodulation compensation are calculated using finite difference method.

[0111] Projecting the spatial gradient derivative onto the radial unit vector field direction yields the radial gradient intensity.

[0112] The absolute value of the derivative of the slice direction is taken to obtain the gradient intensity of the slice direction;

[0113] Specifically, the spatial gradient derivative and slice directional derivative of the intermodulation-compensated image are calculated using finite difference to obtain the rate of change of local grayscale with spatial position using a first-order approximation. Finite difference has the advantages of simple implementation, numerical stability, and uniform application on regular grids, and can suppress high-order noise interference while preserving edge positions. Projecting the spatial gradient derivative in the radial unit vector field direction is to convert the change in any direction into the change in a single geometric direction from the centroid of the slice outward, thereby obtaining a boundary strength measure related to the radius, which is convenient for direct comparison between different slices and under different body positions. Taking the absolute value of the slice directional derivative is to eliminate the difference in sign caused by the different grayscale order inside and outside the tissue, so that the boundary strength in this direction is determined only by the change amplitude and is not affected by polarity. The above processing together establishes two comparable intensity quantities in orthogonal directions, one reflecting the edge sharpness along the radial direction in the plane, and the other reflecting the edge sharpness in the slice stacking direction, providing a reliable, unified, and verifiable physical quantity basis for the subsequent construction of anisotropic boundary uncertainty ratio and local anisotropy map based on the intensity relationship of the two directions.

[0114] Specifically, the intermodulation-compensated image is a three-dimensional matrix containing multiple slices arranged along the axis, with pixels within each slice distributed according to two-dimensional coordinates. For each pixel in the image, the finite difference method is used to calculate its horizontal and vertical gradient components in the slice plane, forming the spatial gradient derivative. The horizontal gradient component is obtained by dividing the gray-level difference between the pixel and its right-hand neighbor by the horizontal pixel spacing, and the vertical gradient component is obtained by dividing the gray-level difference between the pixel and its lower neighbor by the vertical pixel spacing. If a pixel is located at the edge of a slice without adjacent pixels, one-sided difference is used, i.e., the gray-level difference between the pixel and its nearest neighbor on the same side is used for calculation. Simultaneously, the slice direction derivative of the pixel in the axial direction is calculated by dividing the gray-level difference between the pixel and its neighboring upper slice pixel at the same coordinate by the axial slice spacing. If it is the first or last slice without adjacent slices, one-sided difference is used for calculation.

[0115] Specifically, the spatial gradient derivative is projected onto the radial unit vector field to obtain the radial gradient intensity. For each pixel, the horizontal and vertical components of its corresponding spatial gradient derivative are retrieved to form a two-dimensional gradient vector; simultaneously, the radial unit vector of that pixel in the radial unit vector field, which contains both horizontal and vertical components, is retrieved. The dot product of the two-dimensional gradient vector and the radial unit vector is calculated, i.e., the horizontal components are multiplied, the vertical components are multiplied, and then the sum is obtained to get the projection result. The absolute value of this projection result is taken to obtain the radial gradient intensity of that pixel. This process is repeated for all pixels to form a radial gradient intensity matrix of the same size as the intermodulation-compensated image.

[0116] Specifically, the slice orientation gradient intensity is obtained by taking the absolute value of the slice orientation derivative. For each pixel, its slice orientation derivative, calculated using finite difference, is retrieved; this derivative represents the axial grayscale change rate. The absolute value of this slice orientation derivative is then performed to eliminate directional influence, retaining only the intensity of the change. The calculation results for all pixels are arranged according to their coordinate positions in the 3D image, forming a slice orientation gradient intensity matrix of the same size as the intermodulation compensated image. The value of each element in the matrix represents the slice orientation gradient intensity of the corresponding pixel.

[0117] The arithmetic mean of the gradient intensity in the slice direction and the radial direction is calculated separately within the body mask;

[0118] The anisotropic boundary uncertainty ratio is obtained by the ratio of the arithmetic mean of the gradient intensity in the slice direction to the arithmetic mean of the gradient intensity in the radial direction.

[0119] Specifically, the arithmetic mean of the gradient intensity in the slice direction and the radial direction is calculated within the body mask. This is to globally aggregate the boundary intensities in the two orthogonal directions within the effective anatomical region. By using large-sample averaging, the influence of pixel-level noise and local anomalies is suppressed, so that the obtained average value can represent the overall sharpness level in that direction and has cross-slice comparability. Furthermore, the ratio of the arithmetic mean of the gradient intensity in the slice direction to the arithmetic mean of the gradient intensity in the radial direction is determined. This is because when the average values ​​in the two directions are comparable, a ratio close to 1 indicates that the overall boundary uncertainty is approximately balanced in the two directions, while a ratio deviating from 1 quantitatively reflects the degree and direction of directional ambiguity. At the same time, using a ratio instead of a difference can eliminate the influence of overall gain changes and contrast scaling on the measurement, resulting in dimensionless results. This ensures stable comparability under different devices, different subjects, and different acquisition conditions, and provides a reliable global adjustment factor for the subsequent calculation of the position-adaptive threshold field.

[0120] Specifically, the arithmetic mean of the slice orientation gradient intensities is calculated within the body mask. All pixels in the 3D image are traversed, and for each pixel, it is determined whether it lies within the corresponding slice body mask, retaining only pixels covered by the body mask. The slice orientation gradient intensity values ​​of these pixels are collected, forming a slice orientation gradient intensity set. The sum of all values ​​in this set is calculated, and then the sum is divided by the number of values ​​in the set, i.e., the total number of pixels within the body mask. The result is the arithmetic mean of the slice orientation gradient intensities within the body mask. The arithmetic mean of the radial orientation gradient intensities is also calculated within the body mask. Using the same pixel selection rules as for the slice orientation gradient intensities, only pixels within the body mask are selected, and the radial orientation gradient intensity values ​​of these pixels are collected, forming a radial orientation gradient intensity set. The sum of all values ​​in this set is calculated, and then the sum is divided by the number of values ​​in the set (consistent with the total number of pixels within the body mask). The result is the arithmetic mean of the radial orientation gradient intensities within the body mask. The anisotropic boundary uncertainty ratio is calculated based on the above two arithmetic means. The arithmetic mean of the gradient intensities along the slice direction within the body mask is used as the numerator, and the arithmetic mean of the gradient intensities along the radial direction within the body mask is used as the denominator. A division operation is performed, and the resulting quotient is the anisotropic boundary uncertainty ratio. This ratio reflects the relative difference in boundary gradient intensities along the slice direction and the radial direction, quantifying the anisotropic characteristics of boundary uncertainty.

[0121] A local anisotropy map is constructed using the normalized difference between the gradient intensity in the slice direction and the gradient intensity in the radial direction;

[0122] Specifically, a local anisotropy map is constructed using the normalized difference between the gradient intensity in the slice direction and the gradient intensity in the radial direction. This is to provide a bounded and interpretable measure of directionality using edge information from two orthogonal directions at the same pixel. The normalized difference is expressed as the ratio of the difference to the sum of the two gradients. Its value is naturally limited to between -1 and +1. The sign of the value directly indicates which direction is more blurred or sharper, and the absolute value reflects the degree of inconsistency between the two directions. This ratio form can simultaneously offset the effects of overall gain variation and contrast scaling, making the results obtained from different slices and under different acquisition conditions comparable. At the same time, this construction only relies on the two gradient intensities calculated in the previous step, without introducing new parameters or external priors. This facilitates the formation of a smooth and stable spatial distribution at the pixel level, which can be used as a local adjustment factor for the subsequent threshold field, thereby refining the processing of easily confused regions while maintaining global consistency.

[0123] Specifically, the calculated slice orientation gradient intensity matrix and radial orientation gradient intensity matrix are obtained. Both matrices are consistent with the image size after intermodulation compensation. Each element corresponds to the slice orientation gradient intensity and radial orientation gradient intensity of the same coordinate pixel in the 3D image, and only pixels within the body mask have valid gradient intensity values. For each pixel in the 3D image, the difference between the slice orientation gradient intensity and the radial orientation gradient intensity is calculated. The value of the pixel in the slice orientation gradient intensity matrix is ​​retrieved, and its value in the radial orientation gradient intensity matrix is ​​subtracted to obtain the gradient intensity difference of the pixel. The difference can be positive or negative, reflecting the relative magnitude of the slice orientation and radial orientation gradient intensities. The gradient intensity sum value corresponding to each pixel is calculated. The slice orientation gradient intensity value and the radial orientation gradient intensity value of the pixel are added to obtain the gradient intensity sum value. To avoid the sum being zero, which would cause abnormal subsequent division operations, a preset numerical stability constant is introduced. This constant is a small positive value, which is added to the gradient intensity sum value to obtain the corrected sum value as the denominator for the normalization operation. Dividing the gradient intensity difference by the corrected sum yields the local anisotropy value for that pixel. This value is constrained to a range of -1 to +1, where a positive value indicates dominance of the gradient intensity in the slice direction and a negative value indicates dominance in the radial direction. The local anisotropy values ​​of all pixels are arranged sequentially according to their coordinate positions in the 3D image, forming a 3D matrix with dimensions identical to the intermodulation-compensated image. This matrix is ​​the local anisotropy map, where the value of each element quantifies the gradient intensity anisotropy characteristics at the corresponding location.

[0124] Threshold segmentation is performed on the intermodulation-compensated image based on the radial low-frequency intermodulation index, the anisotropic boundary uncertainty ratio, the sign of the radial low-frequency slope sequence, and the local anisotropic map.

[0125] In embodiments of the present invention, threshold segmentation is performed on the intermodulation-compensated image based on the radial low-frequency intermodulation index, the anisotropic boundary uncertainty ratio, the sign of the radial low-frequency slope sequence, and the local anisotropic map, including:

[0126] Within the body mask of each slice, the intermodulation-compensated image is processed using the Otsu method thresholding to obtain the slice baseline threshold.

[0127] The median absolute deviation of the intermodulation-compensated image is calculated within the body mask of each slice to obtain the slice intensity scale.

[0128] Specifically, the Otsu method thresholding is applied to the intermodulation-compensated image within the body mask of each slice. This is to adaptively determine the grayscale segmentation baseline of the slice in a statistical domain containing only the effective anatomical region, thereby avoiding histogram interference from the background and air. The resulting slice baseline threshold can reflect the true contrast conditions of the slice and serve as the anchor point for the subsequent position-adaptive threshold field. Furthermore, the median absolute deviation of the intermodulation-compensated image is calculated within the same body mask to describe the grayscale fluctuation intensity of the slice with a robust scale, suppressing the influence of isolated highlights or metallic artifacts on scale evaluation. This ensures that the resulting slice intensity scale is dimensionless with grayscale and comparable to different devices and acquisition protocols. Consequently, in the subsequent threshold offset calculation, the dimensionless adjustment term is converted into an absolute offset consistent with grayscale, ensuring the stability, interpretability, and cross-slice consistency of the threshold setting.

[0129] Specifically, for each slice in the intermodulation-compensated image, the slice baseline threshold is calculated within its corresponding body mask. First, the body mask of the slice is extracted, and the pixel range covered by the mask is determined. Only the grayscale values ​​of these pixels in the intermodulation-compensated image are selected to form the grayscale set within the mask for that slice, excluding the grayscale information of pixels outside the body mask. A grayscale histogram is constructed based on the grayscale set within the mask, with the horizontal axis representing the grayscale value range and the vertical axis representing the number of pixels in each range. All possible grayscale thresholds are iterated, and the grayscale set within the mask is divided into a foreground pixel subset and a background pixel subset using each candidate threshold as a boundary. The grayscale mean and pixel percentage of each subset are calculated. The inter-class variance is calculated based on the grayscale mean and pixel percentage of the two subsets, reflecting the degree of separation between the two subsets. The candidate threshold that maximizes the inter-class variance is selected and determined as the slice baseline threshold for that slice. All slices in the intermodulation-compensated image are processed sequentially according to the above procedure to obtain the slice baseline threshold corresponding to each slice.

[0130] Specifically, the median absolute deviation of the intermodulation-compensated image is calculated within the body mask of each slice to determine the slice intensity scale. For a single slice, the grayscale values ​​of all pixels covered by the intermodulation-compensated image are extracted based on its body mask, forming a mask-in-mask compensation grayscale set. The number of elements in this set is the same as the number of pixels within the body mask. All values ​​in the mask-in-mask compensation grayscale set are arranged in ascending order. If the number of values ​​is odd, the value in the middle position is selected as the median of the set; if the number of values ​​is even, the arithmetic mean of the two middle values ​​is selected as the median. The difference between each value in the mask-in-mask compensation grayscale set and the median is calculated. The absolute value of each difference is taken to obtain a series of absolute deviation values, forming an absolute deviation set. The values ​​in the absolute deviation set are arranged in ascending order, and the median selection rule is repeated to obtain the median of the absolute deviation set. This median is the slice intensity scale of the slice. All slices are processed in the above manner to obtain the slice intensity scale corresponding to each slice.

[0131] The location-adaptive threshold field is calculated based on the slice baseline threshold, slice intensity scale, axial low-frequency intermodulation index, anisotropic boundary uncertainty ratio, sign of the radial low-frequency slope sequence, and local anisotropic map.

[0132] Specifically, the formula for calculating the location-adaptive threshold field is as follows:

[0133]

[0134] In the formula, For location-adaptive threshold field, The slice baseline threshold, The low-frequency intermodulation index of the shaft diameter. The uncertainty ratio of the anisotropic boundary. The sign for the radial low-frequency slope. This is a local anisotropy diagram. For slice intensity scale, and These are the pixel coordinates within the slice. The slice number;

[0135] Specifically, the position-adaptive threshold field uses the slice baseline threshold as an anchor point and introduces four types of quantities—global scale and directionality information, spatial distribution, and intensity scale—to complete a clear offset correction: the axial low-frequency intermodulation index is used as a weight for the fluctuation intensity across slices; the greater the fluctuation, the more significant the threshold offset. The magnitude and direction of the deviation from isotropy are obtained by subtracting one from the anisotropic boundary uncertainty ratio; this quantity is positive when the slice direction is relatively ambiguous and negative otherwise, thus determining the polarity along which the offset should be increased or suppressed. The sign term of the radial low-frequency slope is used to determine whether the systematic brightness trend from the center to the periphery should be canceled upwards or downwards, achieving reverse correction of low-frequency multiplicative drift. The local anisotropy map provides the location where the threshold needs to be strengthened or relaxed at the pixel level, so that the offset only occurs in areas with prominent inconsistencies in direction. The slice intensity scale maps the aforementioned dimensionless product into an absolute offset with the same dimension as grayscale, ensuring that it can be directly added to the slice baseline threshold without destroying dimensional consistency. The above product and summation structure conforms to the principle of dimension balancing and the physical inversion idea of ​​log-linear fluctuations. When the low-frequency intermodulation index of the axis diameter is close to zero, or the uncertainty ratio of the anisotropy boundary is close to one, or the sign term is zero, or the local anisotropy map is close to zero, the offset naturally decays to zero, and the threshold returns to the slice baseline threshold. Thus, it maintains a steady state when the fluctuations and anisotropy are not significant, and adaptively amplifies and corrects when the problem is significant, which ensures interpretability and improves the consistency and robustness of the segmentation.

[0136] Specifically, for each pixel in the 3D image, its slice number is determined, and the slice baseline threshold corresponding to that slice is retrieved. This threshold is calculated using the Otsu method for intermodulation compensation within the slice's body mask. The radial low-frequency intermodulation index is retrieved; this index is the ratio of the discrete standard deviation of the radial low-frequency slope sequence to the median absolute deviation of the 3D grayscale data within the body mask. Simultaneously, the anisotropic boundary uncertainty ratio is retrieved; this ratio is the ratio of the arithmetic mean of the slice directional gradient intensity to the arithmetic mean of the radial gradient intensity within the body mask. The radial low-frequency slope corresponding to the slice to which the pixel belongs is determined, and the sign of the slope is determined: if the radial low-frequency slope is greater than 0, the sign is 1; if the radial low-frequency slope is less than 0, the sign is -1; if the radial low-frequency slope is equal to 0, the sign is 0. The value of the pixel in the local anisotropy map and the slice intensity scale of the slice to which the pixel belongs are retrieved; this scale is the median absolute deviation calculated within the slice's body mask for intermodulation compensation. Calculate the anisotropic boundary uncertainty ratio minus an adjustment term. Then, multiply the radial low-frequency intermodulation index, the adjustment term, the radial low-frequency slope sign, the local anisotropy map value, and the slice intensity scale sequentially to obtain the adaptive adjustment component. Add the slice baseline threshold to the adaptive adjustment component to obtain the position-adaptive threshold for that pixel. This process is repeated for all pixels in the 3D image, ultimately forming a 3D threshold matrix with the same size as the 3D image. This matrix is ​​the position-adaptive threshold field, used to provide the final threshold for threshold segmentation for each pixel.

[0137] Threshold segmentation is performed on the intermodulation-compensated image based on the position-adaptive threshold field, and the segmentation result is obtained by AND operation with the body mask.

[0138] Specifically, threshold segmentation is performed on the intermodulation-compensated image based on the location-adaptive threshold field. This involves using a decision boundary at each pixel that matches the statistical state and physical background of that location. This allows the threshold to be automatically adjusted up or down in areas with significant directional inconsistencies to suppress false positives and false negatives, and to naturally fall back to the slice baseline threshold in areas with insignificant changes to avoid over-segmentation. Subsequently, an AND operation is performed with the body mask to strictly limit the decision result to the effective anatomical range, eliminate irrelevant components such as the bed surface, air, and areas outside the field of view, eliminate the influence of external intensity fluctuations on the result, and maintain the consistency between and within slices. This achieves a balance between global stability and local finesse without adding new parameters, improving the stability, noise resistance, and cross-slice consistency of the segmentation result.

[0139] Specifically, the position-adaptive thresholding field is a three-dimensional threshold matrix of the same size as the intermodulation-compensated image, with each pixel corresponding to a position-adaptive threshold. The intermodulation-compensated image contains corrected three-dimensional grayscale information. For each pixel in the three-dimensional image, its grayscale value in the intermodulation-compensated image is retrieved, and simultaneously, the corresponding threshold in the position-adaptive thresholding field is retrieved. The pixel's grayscale value is compared with the corresponding position-adaptive threshold. If the grayscale value is greater than or equal to the threshold, the pixel is determined to be a target pixel and marked as the first value; if the grayscale value is less than the threshold, the pixel is determined to be a non-target pixel and marked as the second value. The first and second values ​​are distinct discrete values. Through the above comparison operation, the intermodulation-compensated image is transformed into a three-dimensional binary matrix containing only two values. This matrix is ​​the preliminary result of threshold segmentation. The body mask is a three-dimensional binary matrix, where pixels inside the body mask are marked as the first value, and pixels outside the body mask are marked as the second value, consistent with the numerical definitions in the preliminary result of threshold segmentation. A bitwise AND operation is performed on corresponding pixels in the two matrices. If the corresponding pixel value is the first value in both matrices, the result is the first value; if the corresponding pixel value is the second value in either matrix, the result is the second value. The results of the bitwise AND operations on all pixels are combined to form a 3D binary matrix with the same size as the original 3D grayscale data. This matrix is ​​the final segmentation result. In the segmentation result, pixels corresponding to the first value represent the segmented target region, and pixels corresponding to the second value represent the background region. All target region pixels are located within the body mask area.

[0140] The segmentation results are evaluated by a pre-set deep learning reviewer to obtain a question mark scoring map.

[0141] Specifically, the pre-defined deep learning reviewer serves as the result quality control module, evaluating the segmentation results to generate a questionable scoring map. This scoring map provides a visual score at the pixel or small region level, reflecting potential problems such as abnormal morphology and structure, weak boundary gradients, insufficient cross-slice coherence, and isolated small connected regions. Thus, without changing the threshold settings and segmentation results, it provides clear and traceable prompts for manual review, improving quality control efficiency and clinical usability.

[0142] Specifically, the preset deep learning reviewer is a trained neural network model. This model takes the segmentation result and the intermodulation-compensated image as input and outputs a suspicious point score for the corresponding pixel. The model training process is based on labeled sample data, which includes pixel samples known to be correctly and incorrectly segmented. Training is completed by minimizing the error between the predicted suspicious point score and the true label of the sample, enabling the model to identify potential error regions in the segmentation result. The input format is adapted for the segmentation result and the intermodulation-compensated image. The 3D binary matrix of the segmentation result and the 3D grayscale image after intermodulation compensation are aligned by pixel position and merged into multi-channel input data. One channel contains the binary information of the segmentation result, and the other channel contains the grayscale information of the corresponding position, ensuring that the input data retains the original 3D spatial coordinate relationship. The adapted input data is sent to the preset deep learning reviewer to perform inference calculations. The reviewer extracts local and global features of the input data through a multi-layer neural network structure. Local features include the grayscale changes and segmentation status of pixels and their neighborhoods, while global features include the morphology and spatial distribution patterns of the target region. By employing feature fusion and nonlinear transformation, a suspicion score is calculated for each pixel. This score quantifies the likelihood of an error in the segmentation result for that pixel; a higher score indicates a greater probability of deviation in the segmentation result. The suspicion scores of all pixels output by the deep learning reviewer are collected and arranged according to their coordinate positions in 3D space, forming a 3D scoring matrix with the same size as the original 3D grayscale data. This matrix is ​​the suspicion score map. The value of each element in the suspicion score map corresponds to the degree of suspicion of the segmentation result for that pixel, providing a quantitative basis for optimizing the segmentation results or for manual review.

[0143] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and inventive concept of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. A medical image segmentation method based on deep learning, characterized in that, include: A body mask is generated based on 3D grayscale data, and a radial distance map is generated based on the body mask. The median was statistically analyzed in the radial annulus within the body mask, and the logarithm of the median was linearly fitted to obtain the radial low-frequency slope sequence. The axial low-frequency intermodulation index was determined based on the radial low-frequency slope sequence. The three-dimensional grayscale data is corrected based on the radial distance map and the radial low-frequency slope sequence to obtain the image after intermodulation compensation; Calculate the radial gradient and slice gradient on the image after intermodulation compensation. Based on the radial gradient and slice gradient, calculate the anisotropic boundary uncertainty ratio and form a local anisotropic map. Threshold segmentation is performed on the intermodulation-compensated image based on the radial low-frequency intermodulation index, the anisotropic boundary uncertainty ratio, the sign of the radial low-frequency slope sequence, and the local anisotropic map.

2. The medical image segmentation method based on deep learning according to claim 1, characterized in that, Generate a body mask based on 3D grayscale data, including: The Otsu threshold is obtained by processing slices in three-dimensional grayscale data using the Otsu threshold method. Threshold segmentation of the slice grayscale is performed based on the Otsu method threshold to generate an initial binary region; Select the connected region with the largest number of pixels from all connected regions of the initial binary region as the body mask of the slice.

3. The medical image segmentation method based on deep learning according to claim 2, characterized in that, Generate a radial distance map based on the body mask, including: Calculate the centroid of the slice within the body mask of the slice according to the definition of geometric centroid; Using the centroid of the slice as the center, calculate the Euclidean distance from each pixel within the slice to the centroid of the slice to obtain the radial distance map of the slice.

4. The medical image segmentation method based on deep learning according to claim 1, characterized in that, Within the body mask, the median is statistically analyzed along radial annular bands, and a linear fit is performed on the logarithm of the median to obtain the radial low-frequency slope sequence, including: Under the constraints of the body mask and radial distance map of each slice, a ring set is constructed according to the pixel-level equal width layer; Calculate the median gray value within each annular band set to obtain the annular band median sequence of the slice; The radial low-frequency slope of the slice is obtained by shifting the median sequence of the ring band in the positive domain and performing a linear least squares fitting in the logarithmic domain. The radial low-frequency slopes of all slices are combined into a radial low-frequency slope sequence.

5. The medical image segmentation method based on deep learning according to claim 1, characterized in that, The low-frequency intermodulation index of the shaft diameter is determined based on the radial low-frequency slope sequence, including: Calculate the discrete standard deviation of the radial low-frequency slope sequence; Calculate the median absolute deviation of the 3D grayscale data within the body mask; The low-frequency intermodulation index of the shaft diameter is determined by the ratio of the discrete standard deviation to the median absolute deviation.

6. The medical image segmentation method based on deep learning according to claim 1, characterized in that, The three-dimensional grayscale data is corrected based on the radial distance map and the radial low-frequency slope sequence to obtain the intermodulation compensated image, including: Based on the radial low-frequency slope sequence and radial distance map, an intermodulation compensation field of the same size as the 3D grayscale data is constructed. The formula for constructing the intermodulation compensation field is as follows: In the formula, For inter-compensation fields, Radial low-frequency slope, This is a pixel-level radial distance map. and These are the pixel coordinates within the slice. The slice number; The three-dimensional grayscale data is corrected by point-by-point multiplication using the intermodulation compensation field to obtain the intermodulation compensated image.

7. The medical image segmentation method based on deep learning according to claim 1, characterized in that, Calculate the radial and slice gradients on the intermodulation-compensated image. Based on the radial and slice gradients, calculate the anisotropic boundary uncertainty ratio and generate a local anisotropic map, including: Within each slice, a radial unit vector field is defined centered on the slice's centroid. The formula for the radial unit vector field is as follows: In the formula, For radial unit vector fields, and These are the pixel coordinates within the slice. The slice number. and The coordinates of the slice's centroid. This is a pixel-level radial distance map. To ensure numerical stability and avoid division by zero; The spatial gradient derivative and slice directional derivative of the image after intermodulation compensation are calculated using finite difference calculus. Projecting the spatial gradient derivative onto the radial unit vector field direction yields the radial gradient intensity. The absolute value of the derivative of the slice direction is taken to obtain the gradient intensity of the slice direction; The arithmetic mean of the gradient intensity in the slice direction and the radial direction is calculated separately within the body mask; The anisotropic boundary uncertainty ratio is obtained by the ratio of the arithmetic mean of the gradient intensity in the slice direction to the arithmetic mean of the gradient intensity in the radial direction. A local anisotropy map is constructed using the normalized difference between the gradient intensity in the slice direction and the gradient intensity in the radial direction.

8. The medical image segmentation method based on deep learning according to claim 1, characterized in that, Threshold segmentation of the intermodulation-compensated image is performed based on the radial low-frequency intermodulation index, the anisotropic boundary uncertainty ratio, the sign of the radial low-frequency slope sequence, and the local anisotropic map, including: Within the body mask of each slice, the intermodulation-compensated image is processed using the Otsu method thresholding to obtain the slice baseline threshold. The median absolute deviation of the intermodulation-compensated image is calculated within the body mask of each slice to obtain the slice intensity scale. Based on the slice baseline threshold, slice intensity scale, axial low-frequency intermodulation index, anisotropic boundary uncertainty ratio, sign of the radial low-frequency slope sequence, and local anisotropy map, the location-adaptive threshold field is calculated. The formula for calculating the location-adaptive threshold field is as follows: In the formula, For location-adaptive threshold field, The slice baseline threshold, The low-frequency intermodulation index of the shaft diameter. The uncertainty ratio of the anisotropic boundary. The sign for the radial low-frequency slope. This is a local anisotropy diagram. For slice intensity scale, and These are the pixel coordinates within the slice. The slice number; Threshold segmentation is performed on the intermodulation-compensated image based on the position-adaptive threshold field, and the segmentation result is obtained by AND operation with the body mask. The segmentation results are evaluated by a pre-defined deep learning reviewer to obtain a score map of suspicious points.

Images