A method and system for magnetic resonance image follow-up of amyloid-related imaging abnormalities

Abstract

This invention discloses a method and system for following up on magnetic resonance imaging (MRI) images of amyloid-related radiological abnormalities, belonging to the field of medical image processing. The method includes: acquiring multimodal MRI image data; performing bias field correction and intensity normalization based on a bias field fitting model to generate a first image set; performing global alignment and edge detection on the first image set to extract anatomical physical boundaries, generating brain region mask data; performing branch deformation field registration on the first image set within the brain parenchyma defined by the brain region mask data to generate a spatially aligned second image set; calculating the average intensity of normal brain parenchyma regions in the second image set to obtain a scaling factor, acquiring a three-dimensional residual image; performing anisotropic diffusion filtering denoising on the three-dimensional residual image, and performing differential denoising and morphological screening processing, outputting the detection results. This invention eliminates systematic artifacts caused by insufficient bias field correction.

Description

Technical Field

[0001] This invention relates to the field of medical image processing, and in particular to a method and system for following up magnetic resonance images of amyloid-related imaging abnormalities. Background Technology

[0002] In recent years, anti-amyloid monoclonal antibodies (such as lencanemab and donepemab) have been shown to help slow the progression of Alzheimer's disease, but they may cause a specific drug side effect—amyloid-related imaging abnormalities (ARIA). ARIA must be diagnosed using magnetic resonance imaging (MRI). Clinical trials and treatment monitoring protocols require longitudinal follow-up MRI scans at multiple time points during treatment. By quantitatively comparing and analyzing images of the same patient at different time points, the occurrence of ARIA can be detected in a timely manner, its severity and evolution can be assessed, and this provides crucial evidence for adjusting treatment plans (such as pausing or discontinuing medication). Therefore, automated ARIA detection technology based on longitudinal MRI is essential for ensuring patient medication safety and supporting the clinical application of anti-amyloid therapy.

[0003] However, follow-up of ARIA-E and ARIA-H based on longitudinal multimodal MRI images faces multiple technical bottlenecks, which severely limit the accuracy and clinical usability of the test results:

[0004] The problem of insufficient accuracy in bias field correction. During MRI scanning, due to the spatial inhomogeneity between the radiofrequency coil's transmitting field and the receiving sensitivity, a low-frequency, slowly changing grayscale distortion field (i.e., bias field) is superimposed on the original image, causing the same tissue type to exhibit different grayscale values ​​in different spatial regions of the image. For ARIA-E detection, it is necessary to compare voxel-by-voxel grayscale differences in images of the same patient at different time points on FLAIR sequences to identify newly developed or expanded edema / exudative areas; the grayscale inconsistency caused by the bias field will directly translate into systematic artifacts, masking the true signal changes of ARIA-E. For ARIA-H detection, the bias field will also interfere with the accurate quantification of low-signal characteristics of microbleeds on T2* or SWI sequences. While the commonly used global low-pass filtering method in the prior art is computationally simple, it inevitably damages the high-frequency edge information that carries anatomical details while filtering out low-frequency bias field components, resulting in blurred tissue boundaries and affecting the accurate delineation of lesion boundaries. In addition, the global low-pass filtering method assumes that the bias field changes uniformly and smoothly throughout the entire image space, which cannot effectively deal with the problem of local gray-level unevenness caused by the geometric differences of coils in different scanning devices.

[0005] Craniostomy presents a challenge due to its poor adaptability to atrophied brain tissue. Craniostomy is a crucial preliminary step in separating the skull, scalp, and other extracranial non-brain tissues from MRI images, preserving only the brain parenchyma. Alzheimer's patients receiving anti-amyloid therapy are predominantly elderly and typically exhibit varying degrees of brain atrophy. Severe atrophy results in extremely deep grooves and folds on the cortical surface, significantly widening the cerebrospinal fluid gaps between the brain tissue and the skull, and creating an extremely complex brain surface topology. Furthermore, edema and exudation in ARIA-E often involve the cortex and subcortical regions, while hemosiderin deposition in ARIA-H is frequently located on the brain surface. The precision of craniostomy directly impacts whether these lesions can be completely preserved. Existing brain extraction methods based on fixed thresholds or geometric templates perform poorly when faced with such complex topological changes: fixed threshold methods struggle to accurately distinguish brain tissue from cerebrospinal fluid in deep sulci; geometric template methods, due to the inability of the template shape to track the complex deformations of deep sulci, are prone to mistakenly removing normal brain tissue at the bottom of the sulci or leaving extracranial tissues such as the dura mater within the brain mask. The incorrectly removed brain tissue regions will result in the permanent loss of lesion information in subsequent analyses, while the remaining extracranial tissue will severely interfere with longitudinal subtraction operations, generating false positive signals simulating ARIA-E or ARIA-H. Summary of the Invention

[0006] One of the objectives of this invention is to provide a magnetic resonance imaging follow-up method for amyloid-related imaging abnormalities, in order to solve the problems of insufficient bias field correction accuracy, poor adaptability of craniotomy to atrophied brain tissue, interference of grayscale from new lesions in longitudinal registration, and lack of targeted artifact removal in the post-processing stage in the prior art.

[0007] This invention is achieved through the following technical solution: a magnetic resonance imaging follow-up method for amyloid-related imaging abnormalities, comprising the following steps: acquiring multimodal MRI image data; extracting white matter control points from the multimodal MRI image data and constructing a bias field fitting model; performing bias field correction and intensity normalization processing based on the bias field fitting model to generate a first image set; performing global alignment and edge detection on the first image set to extract anatomical physical boundaries; after inversely transforming the probabilistic brain mask to individual space, performing active contour evolution with the anatomical physical boundaries as the evolution stopping condition to complete skull dissection and boundary refinement processing, generating... Brain region mask data; within the brain parenchyma defined by the brain region mask data, branch deformation field registration processing is performed on the first image set to generate a spatially aligned second image set; the average intensity of normal brain parenchyma regions in the second image set is statistically analyzed to calculate the scaling factor; based on the scaling factor, linear contrast stretching and voxel-by-voxel differencing are performed on the second image set to complete residual subtraction processing and obtain a three-dimensional residual image; anisotropic diffusion filtering is performed on the three-dimensional residual image for denoising, and differential denoising and morphological screening processing are performed according to the image type to be detected by calling the corresponding differential geometric screening features, and the detection results are output.

[0008] Furthermore, the multimodal MRI imaging data included T1-weighted images, FLAIR images, and SWI images acquired by the same patient at baseline and during follow-up.

[0009] Furthermore, the formula for calculating the voxel-by-voxel difference operation is as follows:

[0010] ,in, For the follow-up images after spatial registration and linear contrast stretching, in voxels grayscale value at that location This represents the grayscale value at the same voxel location in the baseline image. The value of the three-dimensional residual image at that voxel is given.

[0011] Furthermore, the values ​​in the three-dimensional residual image whose absolute values ​​are lower than the lower limit of background tissue noise are truncated to zero. The lower limit of background tissue noise is determined based on the standard deviation of the gray values ​​of the background regions with a value of zero in the brain region mask data.

[0012] Furthermore, anisotropic diffusion filtering denoising is achieved by solving the following partial differential equation:

[0013] ;in, For the virtual time parameter of diffusion evolution, The spatial gradient vector of the three-dimensional residual image. For divergence operators, Preserve the transmission coefficient function at the edge. This is the symbol for partial differentials.

[0014] Furthermore, the formula for calculating the edge-preserving transmission coefficient function is as follows: ;in, The value of the local spatial gradient magnitude of the three-dimensional residual image at the current voxel is given. This is the gradient propagation threshold.

[0015] Further, extracting white matter control points and constructing a bias field fitting model from the multimodal MRI image data includes: constructing an axial slice grayscale histogram for each frame of the multimodal MRI image data, performing Gaussian smoothing on the axial slice grayscale histogram, and determining the average intensity of the white matter region of each axial slice using a peak finding algorithm; automatically detecting white matter control points based on the average intensity of the white matter region, and performing cubic spline fitting on the intensity coefficients of effective slices along the vertical axis; establishing a Voronoi mosaic with each white matter control point as the generation kernel, assigning the reference intensity value of the nearest white matter control point to the non-control point voxels to obtain an initial discrete estimate of the bias field; performing soap bubble smoothing iteration on the initial discrete estimate of the bias field to eliminate grayscale step jumps at the boundaries of adjacent Voronoi regions, and obtaining a converged bias field.

[0016] Furthermore, the partitioning rules for Voronoi tessellation are as follows:

[0017]

[0018] in, For three-dimensional image space domain, Let be the coordinates of any voxel position in the image space. and These are the first two white matter control points in the set. The and the first The coordinates of the white matter control points It is a Euclidean distance metric. For white matter control points The Voronoi region centered on; if Then the initial discrete estimate of the bias field is: ,in White matter control point The original image grayscale value at that location.

[0019] Furthermore, the update rule for smooth iteration of soap bubbles is as follows:

[0020]

[0021] in, This represents the current iteration number. For the white matter control point set, For voxels The set of 26 spatially adjacent voxels centered on the center. For the first Neighboring voxels at the next iteration The bias field estimate at that location, The original image grayscale value at the white matter control point;

[0022] Furthermore, for those belonging to the white matter control point set For voxels, the bias field value is fixed as the original image gray value of that point during the iteration process; for non-control point voxels, the bias field value is updated to the arithmetic mean of the surrounding 26 spatially adjacent voxels in each iteration.

[0023] Furthermore, the bias field fitting model adopts a three-dimensional tensor product cubic B-spline model, the expression of which is:

[0024]

[0025] in, To control grid nodes The spline fitting coefficients are obtained by performing least-squares fitting on the gray values ​​at each white matter control point in the white matter control point set. , , They are respectively , , One-dimensional cubic B-spline basis functions along the coordinate axes; , , This represents the number of nodes controlling the grid along each coordinate axis.

[0026] Further, bias field correction and intensity normalization processing are performed based on the bias field fitting model, including: dividing the original gray value of each frame by the converged bias field to obtain the normalized image gray value. in, For the original image in voxels grayscale value at that location The converged bias field in voxels The estimated value at that location, This is the grayscale value after intensity normalization.

[0027] Further, global alignment and edge detection are performed on the first image set to extract anatomical physical boundaries, including: downsampling the images in the first image set; performing rigid body registration and affine registration based on maximum mutual information with a standard brain template as the reference target to obtain a global transformation matrix; and applying the Laplacian operator to the images in the first image set to calculate the sum of the second-order spatial derivatives of each voxel.

[0028] ,in, The grayscale function of the three-dimensional image after intensity normalization in the first image set is given. For the Laplacian operator; detection The voxel intersection where the numerical sign is flipped is taken as the zero-crossing point, and the closed surface formed by the zero-crossing points is taken as the anatomical physical boundary.

[0029] Further, after inversely transforming the probabilistic brain mask to the individual space, active contour evolution is performed with the anatomical physical boundary as the evolution stopping condition. This includes: using the global transformation matrix and a nonlinear deformation field generated by a cross-correlation-based symmetric normalization model to inversely transform the probabilistic brain mask back to the individual space; performing binarization on the inversely transformed probabilistic brain mask to obtain an initial brain region mask; using the edge of the initial brain region mask as the initial contour, calculating the local gradient magnitude of the corresponding image in the first image set, and constructing a stopping function dependent on the local gradient magnitude; driving the initial contour to evolve along the normal direction through level set partial differential equations; when the initial contour evolves to the anatomical physical boundary, the local gradient magnitude reaches a local maximum, the value of the stopping function approaches zero, and the contour evolution is automatically frozen.

[0030] Furthermore, the stopping function, which depends on the magnitude of the local gradient, can be expressed by the following equation:

[0031] ,in, The standard deviation is The three-dimensional Gaussian smoothing kernel, The local gradient magnitude of the image after smoothing by the Gaussian smoothing kernel. This is the contrast control threshold. An exponential parameter that controls the steepness of the decay of the stopping function.

[0032] Furthermore, the level set partial differential equation can be expressed as follows:

[0033] ,in, Let the level set function be defined, and let the zero level set of the level set function define the position of the contour surface at the current moment. For the virtual time parameter of evolution, For the mean curvature term, Let be the balloon force constant. This is a convection term.

[0034] Further, the skull is dissected and its boundaries are refined to generate brain region mask data, including: performing a three-dimensional morphological closing operation on the frozen contour to fill the holes inside the contour; extracting the maximum connected component through three-dimensional connected component analysis, removing isolated tissue blocks with a volume smaller than a set threshold, and generating the brain region mask data.

[0035] Furthermore, the branch deformation field registration processing includes a first registration pathway for detecting high signal intensity in brain white matter and a second registration pathway for detecting brain microbleeds. The first registration pathway includes: using the baseline T1-weighted image as the mediating modality, establishing a synchronous cross-modal spatial mapping between the baseline T1-weighted image and the baseline FLAIR image through a symmetric normalization model; registering the follow-up FLAIR image to the follow-up T1-weighted image, and obtaining the time-dimensional displacement vector field through the registration between the baseline T1-weighted image and the follow-up T1-weighted image; performing matrix concatenation calculations on the synchronous cross-modal spatial mapping and the time-dimensional displacement vector field to synthesize a global nonlinear deformation field, and using the global nonlinear deformation field to resample the follow-up FLAIR image to align it to the coordinate space of the baseline FLAIR image. The second registration pathway includes: directly performing same-modal registration and resampling on the craniotomized follow-up SWI image and the baseline SWI image.

[0036] Furthermore, the registration in both the first and second registration paths employs a differential homeomorphism based on a symmetric normalization model, the energy functional of which is:

[0037]

[0038] ,in, and These are smooth velocity fields originating from the baseline image and the follow-up image, respectively, that meet at the midpoint of time. For differential regularization operators The defined Sobolev norm, and respectively through the velocity field and Differential homeomorphisms generated by integration along the time direction and These are reference images and moving images, respectively. For image similarity measurement functions, ' ' is a compound operator for functions.

[0039] Furthermore, the image similarity measurement function employs local cross-correlation, and the calculation formula for this local cross-correlation is as follows: ,in, For local cross-correlation similarity measurement, This is the current calculation center voxel position. For A three-dimensional local window centered on the center. and The two images to be registered are voxels in the neighborhood. grayscale value at that location and They are respectively and The average grayscale value within the three-dimensional local window.

[0040] Further, the average intensity of normal brain parenchyma regions in the second image set is statistically analyzed to calculate a scaling factor, including: within the range defined by the brain region mask data, after excluding voxels whose grayscale values ​​deviate from the normal distribution, calculating the average intensity of normal brain parenchyma regions in the baseline images and the average intensity of normal brain parenchyma regions in the follow-up images respectively; multiplying the grayscale value of each voxel in the follow-up images by a scaling factor. To complete the linear contrast stretching.

[0041] Furthermore, scaling factor It can be calculated using the following formula: ,in, This represents the average intensity of a normal brain parenchyma region. The average intensity of normal brain parenchyma regions in images taken during the follow-up period.

[0042] Further, based on the image type to be detected, the corresponding differential geometric screening features are invoked to perform differential denoising and morphological screening processing, including: when the image type to be detected is brain white matter high signal, a spherical structural element with a preset radius is constructed, and the spherical structural element is used to sequentially perform erosion and dilation operations on the brain white matter region in the three-dimensional residual image to complete the morphological opening operation, and registration edge artifacts with a width smaller than the diameter of the spherical structural element are removed; when the image type to be detected is brain microbleed, three-dimensional connected component analysis is performed on the three-dimensional residual image to extract each connected component, the volume and surface area of ​​each connected component are calculated, and the sphericity is calculated based on the volume and surface area, and connected components with sphericity lower than a preset threshold value are removed.

[0043] Furthermore, the magnetic resonance image follow-up method also includes independently performing white matter high signal segmentation on the baseline and follow-up FLAIR images, including: extracting gray-level statistical features within the white matter mask constraint range and inputting them into a linear regression prior model to calculate the segmentation threshold K; extracting seed points based on the segmentation threshold and performing iterative dilation until the hierarchical convergence rule is met; and performing connected component segmentation on the dilation result to generate a white matter high signal segmentation mask.

[0044] Furthermore, using sphericity as the differential geometric screening feature to remove artifacts in flattened or tubular vascular cross-sections includes: performing three-dimensional connected component analysis on the negative value regions corresponding to the SWI images in the three-dimensional residual images to extract each independent connected component; and calculating the volume of each connected component. and surface area The sphericity is calculated according to the following formula: ;in, Let V be the volume of the connected component. The surface area of ​​the connected component. Pi; the sphericity The range of values ​​is A preset sphericity threshold is set, and connected components with sphericity lower than the preset sphericity threshold are removed. Connected components with sphericity higher than or equal to the preset sphericity threshold are retained as suspected microbleeding lesions.

[0045] Furthermore, the preset sphericity threshold value is 0.4; wherein, when the connected component is a tubular cross-section of a small blood vessel passing through the brain, the surface area of ​​the connected component under the same volume is significantly larger than the surface area of ​​a sphere, and the sphericity is much less than 1; when the connected component is a microbleeding lesion, the shape of the connected component is close to that of a sphere, and the sphericity is close to 1.

[0046] Furthermore, the magnetic resonance image follow-up method also includes an adaptive fault-tolerant step: when the ratio of the volume of the initial brain region mask to the standard brain volume estimated based on the global transformation matrix deviates from a preset normal physiological range, it is determined that severe brain tissue atrophy exists; in response to the determination, the following adjustments are performed: recalculating local mutual information based on the tissue density distribution of the current downsampled image to update the global transformation matrix; adaptively reducing the weight coefficient of the mean curvature term in the level set partial differential equation to make the contour surface track the boundaries of the atrophic sulci with a sharper shape; adjusting the standard deviation of the Gaussian smoothing kernel in the stopping function. Or lower the contrast control threshold This expands the detection sensitivity range of the anatomical physical boundaries.

[0047] Another aspect of the present invention provides a magnetic resonance imaging follow-up system for amyloid-related imaging abnormalities, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements any of the magnetic resonance imaging follow-up methods for amyloid-related imaging abnormalities as described above.

[0048] Compared with the prior art, the present invention has the following advantages and beneficial effects:

[0049] 1. This invention employs automatic detection of white matter control points and Voronoi mosaicking to provide physically reasonable initial values ​​for the bias field in non-control point regions. This avoids the inherent one-size-fits-all smoothing problem in global low-pass filtering, making the bias field estimation locally adaptive. It can accurately capture the bias field differences at different spatial locations, and is especially suitable for non-uniform gray-level distortion caused by different devices. Furthermore, by using division operations, multiplicative gray-level distortion is completely removed from the physical level, ensuring the physical correctness of the correction. This results in the same tissue type presenting strictly consistent gray values ​​throughout the image space. Image data acquired at different time points and by different devices are transformed into comparable images with a unified gray-level standard, laying a physical foundation for gray-level consistency for longitudinal quantitative comparison across time points. This eliminates systematic false errors caused by insufficient bias field correction at the source.

[0050] 2. This invention injects the zero-crossing point of the Laplacian operator as a hard physical anatomical constraint into the active contour evolution framework of the level set. It achieves precise freezing of the contour on the real cortical surface through a nonlinear stopping function based on the local gradient magnitude. It can equally detect gray-level abrupt change locations in all directions without introducing directional bias. Its zero-crossing point corresponds precisely to the spatial location where the image gray-level gradient magnitude obtains a local maximum from the perspective of differential geometry, that is, the physical candidate location of the tissue boundary. Compared with the fixed threshold method, it has theoretical rigor and physical reliability, and overcomes the problem of insufficient adaptability of existing methods in the extraction of atrophied brain tissue.

[0051] 3. This invention provides a dual-pathway branch deformation field registration architecture for different lesion detection targets. For cross-modal registration of liquid attenuation inversion recovery sequence images, since the spatial resolution of these images is usually lower than that of T1-weighted images and the gray-white matter contrast is weak, directly performing high-precision nonlinear registration between liquid attenuation inversion recovery sequence images may lead to insufficient accuracy. By introducing high-resolution T1-weighted images as an intermediate stepping stone, synchronous cross-modal spatial mapping and time-dimensional displacement vector fields are established respectively. Then, a global nonlinear deformation field is synthesized through matrix cascading calculation, which indirectly but effectively improves the registration accuracy of liquid attenuation inversion recovery sequence images and fundamentally suppresses the generation of registration artifacts. Attached Figure Description

[0052] The accompanying drawings, which are included to provide a further understanding of embodiments of the invention and form part of this application, do not constitute a limitation thereof. In the drawings:

[0053] Figure 1 The image provided in Embodiment 1 of the present invention has undergone field correction and intensity normalization processing.

[0054] Figure 2 The image obtained after skull dissection and mask extraction is provided in Embodiment 1 of the present invention.

[0055] Figure 3 The image provided in Embodiment 1 of the present invention is an image after undergoing Laplacian isotropic edge enhancement response.

[0056] Figure 4 This is a real cerebral cortex surface image provided in Embodiment 1 of the present invention.

[0057] Figure 5 This is an image of extracted white matter provided in Embodiment 1 of the present invention. Detailed Implementation

[0058] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations.

[0059] Example 1

[0060] This embodiment discloses a method for following up on magnetic resonance images of amyloid-related imaging abnormalities. This embodiment includes the following steps:

[0061] S100: Acquire multimodal MRI images from the baseline and follow-up periods, construct a histogram parallel to the xy plane and perform Gaussian smoothing, determine the average intensity of the white matter region using a peak finding algorithm, perform bias field correction and intensity normalization, automatically detect white matter control points, fit the effective slice coefficients with cubic splines, establish a Voronoi diagram, assign non-control points to the values ​​of the nearest control points, and perform soap bubble smoothing iterative correction of voxel intensity to generate the first image set.

[0062] The baseline period refers to the time point at which MRI image data is first acquired for the same patient in a longitudinal follow-up observation study. The image data acquired at this time point serves as the reference benchmark for all subsequent longitudinal comparative analyses.

[0063] The follow-up period refers to the time points at which MRI images are recollected from the same patient at one or more pre-defined time intervals (e.g., 6 months, 12 months, or 24 months after the baseline period) to detect new occurrences, spreads, or regressions by comparing them with the baseline data.

[0064] Multimodal magnetic resonance imaging (MRI) refers to a set of image data obtained by acquiring the brain of the same patient multiple times in the same scanning session using different imaging sequence parameters. This image set includes at least: T1-weighted images (T1WI) to provide high-resolution information on brain tissue anatomy; fluid attenuation inversion recovery sequence (FLAIR) to highlight high-signal lesions in the white matter; and susceptibility-weighted images (SWI) to detect lesions containing hemosiderin deposits, such as cerebral microbleeds.

[0065] The bias field refers to a low-frequency, slowly changing gray-scale distortion field superimposed on the original image during MRI scanning due to the spatial inhomogeneity between the radio frequency coil's transmitting field and the receiving sensitivity. This bias field causes the same tissue type (e.g., white matter) to exhibit different gray-scale values ​​in different spatial regions of the image, severely interfering with subsequent longitudinal quantitative comparisons based on gray-scale differences.

[0066] Bias field correction refers to the process of accurately estimating and removing the multiplicative gray-level distortion field from a physical perspective. Its purpose is to eliminate image gray-level distribution distortion caused by the spatial non-uniformity of the sensitivity of the radio frequency coil of the scanning device.

[0067] Intensity normalization refers to the process of converting image data acquired at different times and from different devices into comparable images with a unified grayscale standard after bias field correction.

[0068] White matter control points refer to a set of voxel coordinates of white matter regions in an image that are automatically detected and confirmed to reliably represent the local sampling values ​​of the bias field at that spatial location. They are denoted as... ,in This represents the total number of control points.

[0069] The first image set refers to a collection of multimodal MRI images generated after bias field correction and intensity normalization, covering the baseline and follow-up periods. All images in this set exhibit consistent grayscale values ​​for the same tissue type across the entire image space, laying the data foundation for subsequent longitudinal quantitative comparisons across time points.

[0070] Understandably, existing technologies often use simple global low-pass filtering to process the bias field, but this method leads to blurred anatomical edges in the image and cannot effectively address local gray-level unevenness caused by different devices. In order to accurately isolate magnetic field interference from a physical level without damaging the high-frequency edge information that carries anatomical details, this embodiment adopts a cascaded modeling strategy of control point detection—Voronoi spatial partitioning—soap bubble iterative smoothing—spline fitting to construct the mathematical model of the bias field.

[0071] Specifically, for each acquired frame of multimodal MRI image, a grayscale histogram is first constructed axially (i.e., for each two-dimensional slice parallel to the xy-plane). This histogram statistically analyzes the frequency distribution of grayscale values ​​of all voxels within the current slice. Subsequently, a Gaussian smoothing kernel with a preset standard deviation is applied to the constructed histogram for convolutional smoothing to suppress random spikes and spurious peaks caused by acquisition noise. On the smoothed histogram, a peak finding algorithm is executed to automatically locate the main peak position of the grayscale value corresponding to the white matter tissue. Since white matter usually has the highest grayscale signal in T1-weighted images, the grayscale value corresponding to this main peak is the average intensity estimate of the white matter region within the current slice, which is used as the reference intensity benchmark for that slice.

[0072] After obtaining the average white matter intensity estimates for slices along each axis, these estimates need to be smoothed and integrated along the longitudinal axis (z-axis). In this embodiment, cubic spline fitting is performed on the intensity coefficients of effective slices (i.e., slices containing sufficient white matter voxels). The cubic spline is chosen because it naturally ensures that the fitted curve is globally consistent. The smooth continuity of the curve, that is, the curve itself and its first and second derivatives are continuous, is a mathematical property that is highly consistent with the physical characteristic of the bias field changing slowly along the vertical axis.

[0073] It should be noted that the bias field is essentially a low-frequency spatial function determined by the physical properties of the magnetic field, exhibiting high smoothness and continuity. To faithfully reconstruct this physical property, three-dimensional tensor product cubic B-splines can be used as the fitting basis function. The reason for choosing cubic B-splines is that they possess local support; the spline coefficients of each control node only affect the fitting result within its finite neighborhood, making the fitting computation efficient and possessing good local adaptability. Simultaneously, cubic B-splines naturally ensure that the fitted surface has global stability. Its smooth continuity enables the stable reconstruction of the spatial distribution of magnetic field inhomogeneities across the entire globe.

[0074] For example, in this embodiment, the three-dimensional tensor product cubic B-spline fitting model of the bias field can be calculated by the following formula:

[0075]

[0076] in, To control grid nodes The spline fitting coefficients at the control point set are obtained by applying the spline fitting coefficients to the control point set. The gray values ​​at each control point are obtained by least-squares fitting. , , They are respectively , , One-dimensional cubic B-spline basis functions along the coordinate axes. Each basis function has a non-zero value in a finite local interval and a value of zero outside the interval. This reflects the local support property of B-splines. , , The number of nodes in the control grid along each coordinate axis determines the spatial resolution of the offset field reconstruction.

[0077] Understandably, by adjusting the resolution of the control grid, a balance can be struck between fitting accuracy and smoothing constraints: more nodes result in higher fitting accuracy, but also greater sensitivity to high-frequency noise; fewer nodes result in better smoothness, but may lose local details of the bias field.

[0078] While constructing the aforementioned continuous fitting model, it is also necessary to provide a reasonable initial value estimate for the bias field of the sparse control point region. To this end, this embodiment uses a Voronoi mosaic (Thieson polygon partitioning) spatial partitioning strategy to naturally divide the entire image space into several non-overlapping regions according to the nearest control point principle, and assigns the grayscale reference value of the nearest control point to each non-control point voxel.

[0079] For example, in this embodiment, the Voronoi partition weight allocation model can be expressed as follows:

[0080]

[0081] in, For the entire three-dimensional image space domain; Let be the coordinates of an arbitrary voxel in the image space; and The first one in the control point set The and the first The position coordinates of each control point; Euclidean distance metric; To control points The Voronoi region is the geometric center, and all voxels within this region are... The distance is no greater than the distance to any other control point. Based on the above partitioning, the initial discrete bias field is assigned the following value: If ,but That is, the original image grayscale value at the corresponding control point in the region is directly used as the initial value of the bias field.

[0082] Understandably, in the absence of dense sampling, the most reasonable initial approximation for estimating the bias field at any point in space is to take the gray value of the control point that is spatially closest to that point, because the spatially gradual variation of the bias field means that the closer the regions are, the more similar their magnetic field deviations will be.

[0083] It should be noted that the initial value of the bias field generated by Voronoi partitioning exhibits significant gray-level abrupt changes at the boundaries of adjacent regions, which violates the physical continuity of the bias field. To eliminate these discontinuities at artificial boundaries, spatial smoothing needs to be applied to the entire bias field while maintaining constant gray-level values ​​at control points. Inspired by the minimum surface phenomenon of soap film in physics, this embodiment models the smoothing process as solving the Laplace equation under fixed boundary conditions. The mathematical problem involves a soap film that automatically forms a smooth surface with minimum average curvature under fixed boundary constraints. Similarly, under Dirichlet boundary conditions with fixed gray values ​​at control points, by iteratively solving the discrete Laplace equation, the bias field will naturally transition between control points, forming a continuous smooth field without step jumps. In discrete three-dimensional voxel space, solving the Laplace equation is equivalent to repeatedly taking the arithmetic mean of all spatial neighbors for each voxel until convergence; this process is also known as Jacobi iteration.

[0084] For example, in this embodiment, the discretized partial differential update formula of the soap bubble smoothing iterative model can be shown as follows:

[0085]

[0086] in, This represents the current iteration number; For the white matter control point set; For voxels Centered A set of spatially adjacent voxels, which covers three types of spatial adjacency relationships of voxels: face neighbors (6), edge neighbors (12), and corner neighbors (8); For the first Neighboring voxels at the next iteration The bias field estimate at the location; This represents the original image grayscale value at the control point.

[0087] It should be noted that in the above formula, for points belonging to the control point set... The voxel's bias field value remains fixed at the original image grayscale value at that point throughout the entire iteration process. This constitutes a fixed Dirichlet boundary condition, ensuring the accuracy of the final reconstructed bias field at the known sampling locations. For non-control point voxels, their bias field values ​​are updated in each iteration to the arithmetic mean of their 26 spatial neighbors. This operation is equivalent to performing a discrete iterative solution of the Laplace equation in three-dimensional space with a maximally isotropic neighborhood. After a sufficient number of iterations and convergence, a highly smooth final bias field that accurately matches the local white matter intensity can be obtained. .

[0088] Finally, the obtained bias field is used to perform intensity normalization correction on the original image. Since the interference of the bias field on the original image is essentially a multiplicative superposition, the low-frequency gray-scale distortion caused by magnetic field inhomogeneity can be completely removed from the physical level by dividing the original image by the bias field.

[0089] For example, in this embodiment, the intensity normalization can be calculated using the following formula:

[0090]

[0091] in, For the original image in voxels The grayscale value at that location; This is the final bias field estimate after iterative convergence; The corrected and normalized image grayscale values ​​are shown below. Through this division operation, the multiplicative interference of the bias field is completely eliminated, ensuring that the same tissue type exhibits a consistent grayscale value throughout the image space. Image data acquired at different time points and using different devices are transformed into comparable images with a unified grayscale standard, laying the data foundation for subsequent longitudinal quantitative comparisons across time points. The first image set is generated by performing the above bias field correction and intensity normalization processing on all multimodal images from the baseline and follow-up periods.

[0092] In some embodiments, to more clearly illustrate the specific implementation process of establishing the Voronoi diagram and correcting voxel intensity in step S100, the steps of constructing the Voronoi diagram and correcting voxel intensity further include:

[0093] White matter control points whose intensity is automatically detected and confirmed To generate the kernel, by calculating the image space Coordinates of all non-control points coordinates of each control point Euclidean distance Each non-control point is assigned to the Thiessen polygon containing the shortest occurrence kernel. Within, and directly assign the shortest distance control point. reference strength value .

[0094] In this context, the generating kernel refers to the seed point used to drive spatial partitioning in a Voronoi mosaic, which in this embodiment is the spatial coordinates of each white matter control point. A Thiessen polygon is a convex polyhedral region centered on a given generating kernel, consisting of all spatial points closest to that kernel. Within this region, the Euclidean distance from any point to its own generating kernel is no greater than its distance to any other generating kernel. The reference intensity value refers to the grayscale value of the original image at that control point. In the context of bias field correction, this value is considered the true sampled value of the bias field at that spatial location.

[0095] Subsequently, the soap bubble smoothing iterative algorithm is initiated. This algorithm, constrained by the aforementioned discretized Laplace equation, induces curvature-driven surface tension deformation at the boundaries of each Voronoi polygon. In each iteration, the grayscale value of the non-control point voxels is updated to the arithmetic mean of its 26 spatial neighbors. This is equivalent to a Laplace iterative solution under fixed boundary conditions, smoothing the voxel intensity abrupt changes caused by discrete assignment between adjacent Voronoi cell regions. After iterative convergence, the obtained bias field transitions naturally between control points, forming a smooth spatial field consistent with the physical continuity of the magnetic field. Figure 1 The T1 and FLAIR images after field correction and intensity normalization in this embodiment are shown. Figure 1 The image preprocessing process of this application is shown. The left side is the original T1-weighted nuclear magnetic resonance image, and the right side is the T1 image after bias field correction (or gray level inhomogeneity correction). It can be seen that the gray level of the corrected image is more uniform, which lays the foundation for subsequent calculation of accurate image gradient and Laplacian second derivative.

[0096] S200: Perform cranial dissection and mask extraction on the first image set: downsample the images and perform affine transformation to obtain the global transformation matrix, apply the Laplacian operator for isotropic edge detection, and lock the anatomical physical boundary through zero-crossing points; then perform rigid body and affine registration in combination with maximum mutual information, and generate a nonlinear deformation field using a cross-correlation-based symmetric normalization model; after inversely transforming the probabilistic brain mask back to individual space binarization, use an active contour evolution mechanism based on local gradient magnitude to refine the boundary and generate brain region mask data.

[0097] Skull stripping refers to the process of removing the skull, scalp, and other extracranial non-brain tissues from MRI images, preserving only the brain parenchyma as the region of interest for subsequent analysis. This process is a crucial preliminary step for longitudinal variation detection because the presence of extracranial tissues can severely interfere with subsequent registration and subtraction analysis.

[0098] Downsampling refers to the operation of reducing the spatial resolution of the original high-resolution image by a preset scaling factor. Its purpose is to reduce the computational complexity of the subsequent global affine transformation search, so that the initial registration stage can quickly converge to the nearest neighbor region of the global optimum.

[0099] Affine transformations are a class of linear geometric transformations that preserve linearity and parallelism. They include four types of parameters: translation, rotation, scaling, and shearing, totaling 12 degrees of freedom (in three-dimensional space). This transformation is used to perform preliminary global alignment of individual images with a standard brain template in terms of overall position, orientation, and scale.

[0100] The global transformation matrix is ​​the parameterized representation of an affine transformation, which is a... The homogeneous transformation matrix encodes the global mapping from the individual image coordinate system to the standard brain template coordinate system.

[0101] The Laplace operator is an isotropic second-order differential operator, denoted as . The Laplacian operator is used to calculate the sum of the second-order partial derivatives of the image grayscale function in all directions in three-dimensional space. In brain imaging processing, the Laplacian operator can equally detect grayscale abrupt changes in all directions without introducing directional bias. Its zero-crossing point (i.e., the voxel boundary where the operator output changes from positive to negative or from negative to positive) corresponds to the physical candidate location of the tissue boundary.

[0102] Zero crossover point refers to the three-dimensional scalar field after the Laplace operator is applied. In this context, the voxel boundary is where the numerical sign changes from positive to negative or vice versa. From the perspective of differential geometry, the zero-crossing surface corresponds precisely to the spatial location where the image gray-level gradient magnitude reaches a local maximum, which is also the anatomical physical boundary where the tissue gray-level jump is most dramatic.

[0103] Probabilistic brain mask refers to a standard spatial probability map of brain tissue constructed based on the statistical analysis of brain images of a large sample population. Each voxel has a value between 0 and 1, which represents the statistical probability that the voxel belongs to brain tissue.

[0104] Brain region mask data refers to a binary 3D mask generated after skull dissection and boundary refinement. A voxel with a value of 1 indicates a brain parenchyma region, while a voxel with a value of 0 indicates an extracranial tissue or background region. This mask is used to strictly limit the analysis scope to the brain parenchyma in subsequent steps.

[0105] Understandably, traditional threshold- or template-based brain extraction methods are prone to misremoving deep sulci or leaving extracranial tissues such as the dura mater when dealing with patients with severe brain atrophy. The fundamental reason is that atrophied sulci create extremely deep depressions on the cortical surface, and simple geometric templates or fixed thresholds cannot track these complex topological changes. To address this issue, this embodiment employs a three-layer modeling strategy: physical edge detection—dynamic contour evolution—gradient braking. First, the Laplacian operator is used to locate candidate positions of anatomical boundaries at the physical level. Then, the contour surface is driven to evolve from its initial position to the brain tissue boundary through level set equations. Finally, a nonlinear stopping function is designed to ensure that the contour is precisely frozen at the actual physical cortical boundary.

[0106] Specifically, the execution process of this step is divided into the following four stages:

[0107] Phase 1: Initial Alignment of Global Space. Individual images in the first image set are downsampled to reduce computational overhead. Then, using a standard brain template as the reference target, rigid body registration (6 degrees of freedom: 3 translations + 3 rotations) and affine registration (12 degrees of freedom) based on maximum mutual information are performed to obtain the global affine transformation matrix.

[0108] The second stage: Nonlinear deformation field generation. Based on the affine registration results, a dense nonlinear deformation field is further generated using a cross-correlation-based Symmetric Normalization (SyN) model to compensate for complex local deformations caused by individual brain morphological differences (such as sulcus depth, ventricle size, etc.). This deformation field establishes a bidirectional reversible mapping between the standard brain template space and the individual image space.

[0109] The third stage: Inverse transformation and binarization of the probabilistic brain mask. Using the inverse mapping of the aforementioned global affine transformation matrix and nonlinear deformation field, the probabilistic brain mask in standard space is inversely transformed back to the original coordinate space of the individual image. In the individual space, the inversely transformed probabilistic mask is binarized using a preset probability threshold (e.g., 0.5) to obtain the initial version of the brain region mask.

[0110] Phase Four: Boundary Refinement Based on Active Contour Evolution Mechanism. The edges of the initial binary brain mask are extracted as the initial contour, which is then driven to expand outward along the normal direction by a level set partial differential equation. When the contour evolves to the anatomical physical boundary marked by the Laplacian zero-crossing point, a stopping function constructed based on the local gradient magnitude brings the evolution rate close to zero, and the contour is precisely frozen at that point. After evolution, a three-dimensional morphological closing operation is performed to fill the holes inside the contour, and isolated tissue blocks with a volume smaller than a set threshold are removed through three-dimensional connected component analysis to generate the final brain region mask data.

[0111] In this embodiment, to more clearly illustrate the specific operational logic of isotropic edge detection using the Laplacian operator in step S200, the operation process of Laplacian operator edge detection is as follows:

[0112] Extract the image to be processed (i.e., the normalized image in the first image set generated after step S100) Three-dimensional coordinates of each voxel Calculate the sum of its second-order spatial derivatives.

[0113] For example, in this embodiment, the Laplacian isotropic edge enhancement response can be calculated using the following formula:

[0114]

[0115] in, The grayscale function of the normalized 3D image; For the Laplace operator; , , Image grayscale along , , Second-order partial derivatives in the directions of the three coordinate axes. The defined surface is the zero-crossing surface, which corresponds to the spatial location where the image gray-level gradient magnitude reaches a local maximum, i.e., the physical candidate location of the tissue boundary. This is the symbol for partial differentials.

[0116] It should be noted that there are significant gray-level abrupt changes between the brain parenchyma, skull, and cerebrospinal fluid. These abrupt changes are represented as zero-crossing points in the second derivative (Laplacian operator) of the image. As an isotropic second-order differential operator, the Laplacian operator can equally detect gray-level abrupt changes in all directions without introducing directional bias. Therefore, its zero-crossing points can serve as highly reliable candidate markers for the physical boundaries of brain tissue.

[0117] After traversal and solution In a three-dimensional scalar field, the voxel boundaries where the sign changes from positive to negative or from negative to positive are identified as zero-crossing points. The closed surface formed by these zero-crossing points is then used as the initial anatomical physical boundary of the brain tissue.

[0118] In this embodiment, to more clearly illustrate the specific implementation process of boundary refinement using the active contour evolution mechanism based on local gradient magnitude in step S200, the boundary refinement step specifically includes:

[0119] The edges of the binarized probabilistic brain mask are extracted as the initial contour, the local gradient magnitude of the corresponding original image is calculated, and a stopping function dependent on the gradient magnitude is constructed.

[0120] Here, the initial contour refers to the three-dimensional closed surface extracted from the boundary of the aforementioned binary probabilistic brain mask. This surface is implicitly represented as a higher-dimensional level set function within the level set framework. The zero level set, i.e. The position of the contour surface at the current moment is defined.

[0121] The local gradient magnitude refers to the magnitude of the first-order spatial gradient vector of the image gray-level function at each voxel location in the image space after Gaussian pre-smoothing. ,in The standard deviation is The three-dimensional Gaussian smoothing kernel, This is the image after Gaussian smoothing. The purpose of Gaussian pre-smoothing is to suppress the interference of acquisition noise on gradient calculation.

[0122] A stopping function is a nonlinear, monotonically decreasing function that depends on the magnitude of the local gradient. The design principle is as follows: In uniform regions of an image (such as the white or gray matter of the brain), the gradient magnitude is small, and the stopping function value is close to 1, allowing the contour to evolve freely; when the contour reaches the physical boundary marked by the Laplacian zero cross, the magnitude of the first-order gradient at that point reaches a local maximum, and the stopping function rapidly decays to close to 0, thus freezing the contour evolution at this point.

[0123] For example, in this embodiment, the gradient magnitude stopping function can be expressed as follows:

[0124]

[0125] in, The standard deviation is The three-dimensional Gaussian smoothing kernel; The image is Gaussian smoothed. The local gradient magnitude of the smoothed image; This is the contrast control threshold, used to set the critical value that distinguishes high gradients from low gradients. The exponential parameter controls the steepness of the decay of the function. The larger the value, the steeper the transition of the function near the threshold, and the more accurate the boundary positioning.

[0126] Understandably, when the contour is located within a uniform tissue, ,at this time The evolution equation normally drives the contour to expand outward. When the contour reaches the Laplacian zero-crossing point, it corresponds precisely to a local maximum of the first-order gradient magnitude. ,at this time The evolution rate of the contour then approaches zero, thus being precisely frozen at the actual physical boundaries of the cortex.

[0127] The initial contour is controlled to expand outward along the normal direction. This embodiment uses an active contour evolution model based on the level set method to drive the motion of the contour surface. The core idea of ​​the level set method is to implicitly represent the three-dimensional brain surface contour to be tracked as a level set function. The zero-level set, by driving the function with virtual time Evolution indirectly enables the tracking of contour surfaces. The advantage of this representation is that the contour surfaces can undergo natural topological changes (such as splitting and merging) without the need for complex topological processing in the algorithm. This is crucial for tracking the surface of the atrophic cerebral cortex, which has complex topological structures such as deep grooves and fissures.

[0128] For example, in this embodiment, the level set partial differential equation (PDE) for active contour evolution can be expressed as follows:

[0129]

[0130] in, Let be a level set function, and let its zero level set be... The current position of the contour surface is defined. Indicates extracranial regions, Indicates a region within the brain; For the virtual time parameters of evolution; The spatial gradient of the level set function. Its modulus; The MeanCurvature is the average curvature of the contour surface at the current position. Its function is to control the local smoothness of the contour and prevent the contour surface from producing excessive wrinkles and burrs due to noise interference. The Balloon Force constant provides a constant outward expansion driving force, enabling the outline to actively expand to the vicinity of the brain tissue boundary; The aforementioned velocity control function (i.e., the stopping function) based on image gradient controls the evolution speed of the contour; The advection term attracts the contour to the region where the gradient change is most significant.

[0131] It should be noted that in the above equation: the average curvature term and the balloon force constant in the first term on the right-hand side together constitute the main driving force for the contour evolution, and the stopping function Speed ​​control is applied to the driving force; the second convection term attracts the profile surface to the stopping function. The region where the gradient changes most dramatically is the actual anatomical physical boundary. When the contour evolves to the anatomical boundary where the gradient magnitude is at a local maximum (i.e., the zero-crossing point region), and The direction points towards the boundary, and the two conditions work together to make... The evolution automatically stops. Through this mechanism, the Laplace zero-crossing point is injected into the level set evolution process as a rigid physical anatomical constraint, ensuring that the final extracted brain mask boundary faithfully fits the actual anatomical cortical surface.

[0132] After the evolution is complete, a three-dimensional morphological closing operation (i.e., expansion followed by erosion) is performed to fill any small holes that may remain inside the contour due to cerebrospinal fluid filling. Then, through three-dimensional connected component analysis, the largest connected component is identified and extracted, and isolated tissue blocks with a volume smaller than a set threshold (e.g., less than 1% of the expected total brain volume) are removed to generate the final brain region mask data. Figure 2 The T1 and the image after skull dissection are shown in this embodiment. Figure 3 The image shown in this embodiment is the result of Laplacian isotropic edge enhancement response; Figure 4 This embodiment shows the final, real surface image of the cerebral cortex; wherein, Figure 2 The left side is the original T1 image; the right side is the brain image (template) after most of the non-brain tissue has been removed, providing a rough initial area for the initial step of the subsequent Laplace operation. Figure 3 The results after applying the Laplacian operator are shown. The light-dark boundary (i.e., zero-crossing point) caused by the drastic change in gray level in the figure represents the significant gray level jump position between the brain parenchyma, skull and cerebrospinal fluid. The algorithm uses these zero-crossing points as hard physical constraints to locate the highly reliable anatomical physical boundary. Figure 4 The core results of the algorithm are shown, namely the final real surface image of the cerebral cortex. The left side is the brain region mask during the evolution process; the right side is the extremely detailed brain image extracted after the contour surface evolution was driven by the level set method and precisely frozen at the Laplacian zero crossover point. It can be seen that the algorithm successfully handles complex topological structures and successfully tracks and preserves the details of cerebral cortex folds such as deep sulci and fissures.

[0133] In this embodiment, a white matter mask can be obtained based on brain region mask data and aligned to the FLAIR coordinate space. After mapping back to the individual space using a public prior probability map through inverse transformation, Gaussian smoothing and probability threshold binarization are performed. Gray-scale range correction is performed in combination with the gray-scale characteristics of T1 images. After morphological refinement, a T1 space white matter mask is generated. Subsequently, the T1 image is aligned to the FLAIR image through mutual information-driven rigid registration and symmetric normalized nonlinear registration. The white matter mask is then mapped to the FLAIR coordinate space using the obtained transform field in a nearest neighbor interpolation manner to generate FLAIR space white matter mask data.

[0134] Understandably, due to the different imaging principles of T1 and FLAIR sequences, the spatial location of the same anatomical structure varies slightly between the two sequences, while high-signal lesions in white matter are primarily more prominent in FLAIR sequences. Therefore, it is necessary to precisely align the white matter mask acquired from T1 images to FLAIR space to achieve spatial matching between the white matter region and the high signal intensity in FLAIR. Simultaneously, the acquisition of the white matter mask requires combining anatomical prior information from the prior probability map with the adaptive adjustment capability of the grayscale characteristics of T1 images to address the segmentation bias problem of a single method when individual anatomical differences are significant. Specifically, this step is executed in two stages:

[0135] Phase 1: Acquisition of the T1 spatial white matter mask. First, a publicly available prior probability map of brain white matter is selected. Using the inverse transformation matrix and inverse deformation field generated in step S200, the white matter probability map in standard space is resampled and mapped back to the patient's individual T1 image space, obtaining a preliminary white matter probability distribution map. Subsequently, the white matter probability map mapped back to the individual space is Gaussian smoothed, with the Gaussian kernel standard deviation set to 1.5 mm. The purpose of this smoothing operation is to eliminate local noise caused by registration errors between the probability map and the individual image, making the probability distribution more closely match the true spatial morphology of the individual white matter. On the smoothed probability map, a probability threshold of 0.6 is set, and regions with voxel values ​​greater than or equal to 0.6 in the probability map are marked as candidate white matter regions, generating a primary white matter mask and initially locking down the core range of the white matter.

[0136] The second stage: Aligning the white matter mask from T1 space to FLAIR space. Using the FLAIR image as the fixed image and the T1 image as the moving image, maximum mutual information is used as the similarity measure, and 6-DOF rigid body registration is performed to calculate the initial transformation matrix, eliminating the global spatial offset between the two image sequences. Based on the rigid body registration results, a symmetric normalization model is used for nonlinear deformation registration, with cross-correlation as the local similarity measure. The deformation field is iteratively optimized to accurately align the anatomical structures of the T1 image and the FLAIR image at the local level, resulting in the nonlinear deformation field. The white matter mask in T1 space is then mapped to FLAIR space using the transformation matrix and nonlinear deformation field obtained in the above steps using the nearest neighbor interpolation method. The reason for choosing nearest neighbor interpolation is that this interpolation method can avoid blurring of the mask voxel values, ensure that the binary characteristics of the mask (i.e., values ​​are only 0 or 1) are not destroyed, and accurately preserve the boundary information of the white matter region. The white matter mask mapped to FLAIR space is subjected to a three-dimensional morphological closing operation again to fill in the small holes caused by registration deviation, ensuring the spatial consistency between the mask and the white matter region in the FLAIR image, and generating the final FLAIR space white matter mask data.

[0137] S300: Branch deformation field registration based on brain region mask data: For high signal detection in brain white matter, the baseline high-resolution T1 image is used as an intermediate transition to perform cross-modal registration between T1 and FLAIR, obtain the deformation field, and resample the follow-up FLAIR image to the baseline FLAIR coordinate space; For brain microbleed detection, the same-modal registration and resampling are directly performed on the follow-up SWI and baseline SWI images after craniotomy to generate an aligned second image set.

[0138] Branch deformation field registration refers to a processing strategy that designs two independent spatial registration pathways (i.e., branch pathways) for different types of detection targets to optimize the alignment of follow-up images to the coordinate space of baseline images. The reason for using branch registration instead of a single registration pipeline is that the detection of different lesion types depends on different image modalities (white matter high signal depends on FLAIR modality, microbleeds depend on SWI modality), and there are fundamental differences in grayscale characteristics and registration difficulty between different modalities.

[0139] A deformation field is a dense displacement vector field defined in a three-dimensional image space, where a three-dimensional displacement vector is stored at each voxel location, indicating the local spatial mapping relationship from the source image coordinate space to the target image coordinate space.

[0140] Cross-modal registration refers to the registration operation that establishes a spatial correspondence between two images with different imaging contrasts (e.g., T1-weighted images and FLAIR images). Since the gray-level mapping relationship between cross-modal images is non-linear (the same tissue presents different gray-level values ​​in different modalities), similarity measures that depend on the consistency of absolute gray-level values ​​(such as mean square error) cannot be used directly. Instead, indicators that measure the statistical dependence of gray-level distribution (such as mutual information) are required.

[0141] Same-modal registration refers to the registration operation that establishes a spatial correspondence between two images with the same imaging contrast (e.g., baseline SWI and follow-up SWI). Since the gray-level mapping relationship of same-modal images is approximately linear, gray-level similarity-based metrics can be directly used.

[0142] Resampling refers to using the deformation field obtained from registration to interpolate the voxel gray values ​​in the source image into the coordinate space of the target image according to the spatial mapping relationship indicated by the deformation field, thereby generating a new image whose spatial coordinates are completely aligned with the target image.

[0143] The second image set refers to a spatially aligned image data set generated after branch deformation field registration and resampling. It includes aligned baseline and follow-up FLAIR image pairs (for white matter lesion detection) and aligned baseline and follow-up SWI image pairs (for microbleed detection). The corresponding baselines and follow-up images in this set have been accurately registered to the same coordinate space, allowing for direct voxel-by-voxel residual subtraction.

[0144] Understandably, traditional linear registration (such as rigid body transformation or affine transformation) can only compensate for global translation, rotation, and scaling, and cannot align complex nonlinear local deformations caused by brain atrophy, ventricular dilation, etc. Ordinary nonlinear registration methods face unique challenges when processing longitudinal lesion analysis: newly formed high-signal lesions (such as white matter high signal, new infarct lesions) have abnormal grayscale values, which can interfere with the similarity measurement of the registration algorithm, causing artificial local distortions in the registration field near the lesion, and subsequently generating artifacts (false positives) in the subtraction. To simultaneously address the two core issues of nonlinear alignment accuracy and lesion interference robustness, this embodiment adopts a three-layer cascaded registration strategy: rigid body / affine macroscopic alignment—differential homeomorphic nonlinear fine-tuning—local similarity measurement.

[0145] Specifically, in the macroscopic alignment stage, a global affine transformation matrix needs to be found to initially align the two images in terms of overall position, orientation, and scale. To achieve robust registration between different imaging sequences, this embodiment uses mutual information as a similarity metric. Mutual information, derived from information theory, measures the statistical dependency between the gray-level distributions of two images, rather than the absolute consistency of gray-level values. Therefore, even if the gray-level mapping relationship between two images is non-linear (such as between different imaging modalities), as long as a stable statistical dependency exists, mutual information can provide a reliable registration driving force.

[0146] For example, in this embodiment, the baseline image With follow-up images The mutual information between them can be calculated using the following formula:

[0147]

[0148] in, Let be the affine transformation matrix to be optimized, which includes translation, rotation, scaling, and shearing parameters; This indicates that the follow-up images are subjected to affine transformation. The result after mapping, ' ' is a function compound operator; and These are the grayscale values ​​of the baseline image and the transformed follow-up image, respectively; The joint probability distribution of the two images is obtained by constructing a joint gray-level histogram and estimating it using normalization. and Given their respective marginal probability distributions, the goal of registration is to search for mutual information. Transformation matrix that yields the maximum value .

[0149] Understandably, when two images are perfectly aligned, the joint probability distribution is most concentrated (i.e., the uncertainty is lowest), and the mutual information reaches its maximum. This metric is highly robust in cross-modal registration (such as between T1 and FLAIR) because it does not require the two images to have the same grayscale values, but only that there is a stable statistical dependency between them.

[0150] After completing macroscopic affine alignment, to further compensate for local nonlinear brain tissue deformation, a dense nonlinear deformation field needs to be superimposed on the affine registration result. To ensure the physical rationality of the deformation field, this embodiment uses the Symmetric Normalization (SyN) framework to construct a differential homeomorphic nonlinear registration model. The core idea of ​​Symmetric Normalization (SyN) is: instead of directly solving a unidirectional deformation field from baseline to follow-up, it simultaneously constructs two smooth velocity fields that evolve over time from two time endpoints, allowing them to meet at the midpoint of time. This bidirectional symmetric design ensures that the mapping from baseline to follow-up and the mapping from follow-up to baseline are exact inverse mappings, mathematically eliminating spatial folding and tearing phenomena and guaranteeing the topological integrity of the deformed brain tissue.

[0151] For example, in this embodiment, the symmetric normalized differential homeomorphism energy functional can be represented as follows:

[0152]

[0153] in, and From time (Corresponding to the baseline image end) and (Corresponding to follow-up image end) Starting from, at Two smooth velocity fields that meet at (midpoint of time) and change with time; For differential regularization operators The defined Sobolev norm, It is usually taken as the power form of the Laplace operator, which is used to impose a smoothness constraint on the velocity field and prevent the deformation field from exhibiting excessively severe local distortions. and respectively through the velocity field and Differential homeomorphic mappings generated by integration along the time direction, i.e. ,in It is an identity mapping; and These are the reference image and the moving image, respectively. For image similarity measurement functions, ' ' is a compound operator for functions.

[0154] It should be noted that in the above equation, the first term on the right-hand side is the regularization energy term, which ensures the global smoothness of the deformed field by penalizing the spatial irregularities of the velocity field; the second term is the data fidelity term, which drives registration by maximizing the similarity between the two deformed images. The minimization of the energy functional is achieved by minimizing the velocity field... and The joint optimization is achieved. Since the differential homeomorphism is generated by time integration of the continuous velocity field (i.e., flow mapping), its mathematical properties guarantee the invertibility and bidirectional symmetry of the mapping, thus eliminating the possibility of the Jacobian determinant in the deformable field having zero or negative values ​​(i.e., spatial folding).

[0155] Similarity measures in energy functionals To combat the interference of newly formed high-signal lesions on registration, this embodiment uses Local Cross-Correlation (CC) as a similarity metric, replacing metrics such as Mean Squared Error (MSE) which are sensitive to absolute gray values. The core idea of ​​local cross-correlation is to calculate the correlation coefficient between the gray values ​​of two images within a small local window centered on each voxel, after subtracting their respective local means. By subtracting the local means, local cross-correlation measures whether the gray-level change patterns (i.e., relative trends) of the two images are consistent within a local region, rather than whether their absolute gray values ​​are the same.

[0156] For example, in this embodiment, the local cross-correlation similarity measure can be calculated using the following formula:

[0157]

[0158] in, For local cross-correlation similarity measurement, This is the current central voxel position in the calculation; For A three-dimensional local window centered on the center (e.g., a cube neighborhood with a side length of 5 to 7 voxels). and The two images to be registered are voxels in the neighborhood. The grayscale value at that location; and They are respectively and In local window The average gray level within the range. The range of values ​​is When two images are perfectly linearly correlated within a local window, When completely unrelated, .

[0159] It should be noted that the reason why local cross-correlation measures can resist high signal absolute value bias is that the local mean is subtracted during the calculation process. and When a lesion causes an overall increase in grayscale in a localized area of ​​the follow-up image, the local mean of that area... The value of local cross-correlation also increases accordingly. The residual pattern after subtracting the mean remains consistent with the baseline image. Therefore, the value of local cross-correlation will not decrease significantly due to the presence of lesions, and the registration algorithm will not produce abnormal local deformation in order to align the lesion area. Figure 5 The extracted white matter image in this embodiment is shown; Figure 5 This displays a mask of white matter segmented from the brain with even greater precision. This high-precision segmentation of white matter structure is based on... Figure 4 This is based on the precise extraction of the physical boundaries of the cerebral cortex shown.

[0160] In this embodiment, to more clearly illustrate the specific implementation logic of the branch deformation field registration for high signal detection in white matter in step S300, the registration step specifically includes:

[0161] S310: Extract baseline T1 images and baseline FLAIR images, and establish a synchronous cross-modal spatial mapping with T1 as the main registration target through a symmetric normalization model.

[0162] Simultaneous cross-modal spatial mapping refers to the spatial correspondence established between T1 images and FLAIR images acquired at the same scanning time point. Since T1 images have the highest spatial resolution and the clearest gray-white matter boundaries, using them as the primary registration target can provide the most reliable anatomical reference frame for cross-modal registration.

[0163] Understandably, the need to establish a cross-modal mapping between T1 and FLAIR images stems from the fact that FLAIR images typically have lower spatial resolution than T1 images, and their grayscale contrast is weaker. Directly performing high-precision nonlinear registration between FLAIR images may result in insufficient registration accuracy. By introducing high-resolution T1 images as an intermediate stepping stone, the registration accuracy of FLAIR images can be indirectly improved.

[0164] S320: Register the follow-up period FLAIR image to the follow-up period T1 image, and then obtain the displacement vector field in the time dimension by calculating the registration field between the baseline period T1 and the follow-up period T1 image.

[0165] The displacement vector field in the time dimension refers to a dense displacement vector field that describes the nonlinear spatial mapping relationship from the baseline T1-weighted image coordinate space to the follow-up T1-weighted image coordinate space. This vector field encodes information on local spatial deformation caused by physiological and pathological changes such as brain tissue atrophy and ventricular dilation between the two time points.

[0166] It is understandable that, since the gray-scale contrast characteristics of T1 images remain highly consistent at the two time points (same modality), and T1 images have the best spatial resolution and gray-white matter boundary sharpness, performing nonlinear registration between T1 images can obtain the most accurate temporal displacement vector field.

[0167] S330: Perform matrix concatenation calculations on cross-modal space mapping and time-dimensional displacement vector fields to synthesize a global nonlinear deformation field, and use this global nonlinear deformation field to resample FLAIR during the follow-up period.

[0168] Matrix concatenation refers to the process of mathematically combining multiple independently calculated spatial transformations (including affine transformation matrices and nonlinear deformation fields) in the correct spatial transformation combination order. Specifically, the cross-modal mapping from follow-up FLAIR to follow-up T1, the time dimension mapping from follow-up T1 to baseline T1, and the cross-modal mapping from baseline T1 to baseline FLAIR are sequentially concatenated to synthesize a global nonlinear deformation field from the follow-up FLAIR coordinate space to the baseline FLAIR coordinate space.

[0169] The global nonlinear deformation field refers to a dense displacement vector field synthesized through the aforementioned matrix cascade calculations, which directly describes the comprehensive spatial mapping relationship from the coordinate space of the follow-up FLAIR image to the coordinate space of the baseline FLAIR image. By resampling the follow-up FLAIR image using this global nonlinear deformation field, a follow-up FLAIR image that is precisely aligned with the baseline FLAIR image at the voxel level can be obtained.

[0170] It should be noted that for the registration pathway in the detection of cerebral microbleeds, since the baseline SWI and follow-up SWI are isomodal images, the grayscale mapping relationship between them is approximately linear. Therefore, affine alignment based on the aforementioned mutual information-driven method and symmetric normalized nonlinear registration based on local cross-correlation can be directly performed between the craniotomized follow-up SWI and baseline SWI images without the need for T1 images as an intermediate modality. After registration, the follow-up SWI images are resampled and, together with the output of the white matter lesion detection pathway, constitute the aligned second image set.

[0171] S400: Perform residual subtraction on the second image set: Calculate the scaling factor by statistically analyzing the average intensity of normal brain parenchyma regions, perform linear contrast stretching on baseline and follow-up images based on the scaling factor, and perform voxel-by-voxel difference operation to extract three-dimensional residual images.

[0172] The residual subtraction operation refers to the numerical calculation process of extracting the signal change between two time points by calculating the gray-level difference on a voxel-by-voxel basis between two precisely registered baseline and follow-up images in the same coordinate space. This operation is the core data extraction step for longitudinal change detection.

[0173] Normal brain parenchyma refers to the voxel region occupied by normal brain tissue, such as gray matter and white matter, after excluding known lesion areas, cerebrospinal fluid areas, and background areas in baseline and follow-up imaging.

[0174] The scaling factor is the ratio of the average grayscale value of a normal brain parenchyma region in baseline images to that in follow-up images. This factor is used to compensate for global grayscale scale differences caused by minor differences in acquisition parameters (such as gain settings, receiver bandwidth, etc.) between different scanning sessions.

[0175] Linear contrast stretching refers to the operation of linearly scaling the grayscale values ​​of an image using a scaling factor.

[0176] The scaling factor is calculated as follows: The analysis range is defined using the brain region mask data generated in step S200. Within the mask range, voxels with gray values ​​significantly deviating from the normal distribution (i.e., potential lesion areas) are excluded, and the average gray value of the remaining normal brain parenchyma voxels is calculated. The average intensity of normal brain parenchyma in the baseline images is then calculated. Mean intensity of normal brain parenchyma in follow-up images scaling factor Defined as Subsequently, the grayscale value of each voxel in the follow-up images was multiplied by the scaling factor. This completes linear contrast stretching, making the global grayscale of the two images in the normal tissue area tend to be consistent.

[0177] Understandably, the purpose of performing linear contrast stretching is to address the possibility of minor global grayscale scale differences remaining between different scanning sessions, even after the bias field correction and intensity normalization processing in step S100. Without compensation, these scale differences would be transformed into global systematic artifacts during subtraction, masking the true local lesion changes. Linear contrast stretching using a scaling factor can eliminate global grayscale shifts between different scanning sessions without altering the relative grayscale relationships within the image.

[0178] In this embodiment, to more clearly illustrate the specific mathematical process of extracting the three-dimensional residual image by performing residual subtraction in step S400, the generation of the residual image can be determined by the following formula:

[0179] After performing linear contrast stretching, the voxel intensity of the follow-up images registered in the same spatial coordinate system was analyzed. Voxel intensity of baseline images The difference is calculated for each voxel.

[0180] For example, in this embodiment, the residual image voxel The formula for generating can be calculated using the following formula:

[0181]

[0182] in, For follow-up images after spatial registration and linear contrast stretching, in voxel The grayscale value at that location; This represents the grayscale value of the baseline image at the same voxel location; This represents the value of the residual map at that voxel.

[0183] It is understandable that positive values ​​in the residual plot indicate increased signal during follow-up (such as newly developed high-signal lesions in white matter), negative values ​​indicate decreased signal (such as newly developed microbleed lesions), and ideally, the area that has not changed should have a value of zero.

[0184] It should be noted that, in order to suppress meaningless residuals generated in background areas and low signal-to-noise ratio edge areas, the following measures are taken: Values ​​with absolute values ​​below the set lower limit for background tissue noise are truncated to zero. This noise lower limit can be adaptively determined based on the standard deviation of gray values ​​in the background region outside the brain parenchyma (i.e., the region where the brain region mask value is zero) in the residual image, for example, set to 2 to 3 times the standard deviation of the background region. Through this truncation operation, only statistically significant signal change regions are retained in the residual image, thereby reducing the false positive rate of subsequent morphological screening steps.

[0185] In addition, fine segmentation of white matter hypersignal based on interpretable prior models and seed point dilation was performed across domains on FLAIR spatial white matter mask data and FLAIR images: clinically interpretable gray-scale statistical features were extracted within the constraints of the white matter mask, and the segmentation threshold was adaptively calculated by inputting them into a pre-trained linear regression prior model; based on the segmentation threshold, seed points were extracted within the white matter mask and conditional dilation was performed until the hierarchical convergence rule was met and then the dilation was stopped; connected component segmentation and quantitative statistics were performed on the dilation results, and white matter hypersignal segmentation masks were independently generated on FLAIR images at baseline and follow-up periods, respectively.

[0186] Understandably, this step can detect the signal change between two time points, but the result is a continuous grayscale difference map, which cannot directly provide the precise spatial range and boundaries of high-signal lesions in white matter. Traditional segmentation algorithms mostly rely on manually setting global or local grayscale thresholds. Due to differences in the degree of brain atrophy, demyelination level, and field strength of MRI scanners among different patients, fixed thresholds are difficult to apply universally. Existing adaptive thresholding methods often use unsupervised clustering or deep learning models to calculate segmentation thresholds, lacking clinical interpretability. To address these issues, this embodiment designs a linear regression prior model based on interpretable grayscale features to adaptively calculate the optimal segmentation threshold for different patients. Combined with seed point extraction and conditional dilation strategies, it achieves complete and fine segmentation of high-signal lesions in white matter under the spatial constraints of the white matter mask. Specifically, the execution process of this step is divided into the following three stages:

[0187] Phase 1: Adaptive calculation of segmentation threshold using an interpretable prior model. .

[0188] Within the constraints of the FLAIR spatial white matter mask data generated in step S200, six clinically significant gray-level statistical features are extracted from the current FLAIR brain parenchyma image to be processed. The six gray-level statistical features are as follows:

[0189] The first characteristic is the average gray value of the white matter region. , defined as the arithmetic mean of the gray levels of all voxels within the white matter mask, reflects the baseline gray level of the patient's normal white matter.

[0190] The second characteristic is the standard deviation of grayscale in the white matter region. , defined as the standard deviation of grayscale of all voxels within the white matter mask, reflects the degree of grayscale uniformity of normal white matter.

[0191] The third characteristic is the maximum gray value of the white matter region. , defined as the highest gray value among all voxels within the white matter mask, is used to initially determine the upper limit of gray values ​​for candidate regions that may have high white matter signals.

[0192] The fourth feature is the grayscale contrast of the white matter region. To avoid outliers caused by noise interfering with contrast calculation, an outlier removal operation is first performed on the gray values ​​of all voxels within the white matter mask: the mean gray value of all voxels within the white matter mask is calculated. and standard deviation Remove grayscale values ​​less than or greater than The voxels represent a type of voxel that is likely an extreme value caused by noise; the maximum value is calculated from the gray values ​​of the remaining valid voxels. and minimum value The final grayscale contrast of the white matter region is defined as follows: This feature reflects the true grayscale difference within the white matter region after excluding noise.

[0193] The fifth feature is the gray-scale skewness of the white matter region. , defined as the normalized value of the third central moment of the gray-level distribution of all voxels within the white matter mask, reflects the asymmetry of the gray-level distribution in the white matter region. The presence of high signals in the white matter will cause the gray-level distribution to skew towards higher values, resulting in increased skewness.

[0194] The sixth feature is the grayscale kurtosis of the white matter region. , defined as the normalized value of the fourth central moment of the gray-level distribution of all voxels within the white matter mask, reflects the steepness of the gray-level distribution in the white matter region. The concentrated presence of high signals in the white matter will cause sharp high-value peaks in the gray-level distribution, resulting in increased kurtosis.

[0195] All six features were extracted within the constraints of the white matter mask, ensuring a high correlation between the features and white matter tissue and white matter high signal intensity. Each feature corresponds to a clear clinical imaging meaning: and Describes the grayscale baseline state of normal white matter. and Describes the grayscale extreme values ​​and dynamic range of white matter regions. and These six features describe the high-order statistical morphological characteristics of grayscale distribution. They characterize the grayscale distribution of white matter regions in the current patient's FLAIR images from different dimensions, providing sufficient and interpretable input information for the regression model.

[0196] The training process of the linear regression prior model is as follows: Collect a predetermined number (e.g., more than 100 cases) of FLAIR image data with white matter hypersignal regions manually annotated by clinical radiologists. Extract the six grayscale features mentioned above for each sample. Simultaneously, based on the clinically annotated white matter hypersignal regions, the optimal threshold for each sample is calculated using the following formula. :

[0197]

[0198] in, The grayscale value of the remaining normal white matter voxels after removing clinically labeled high-signal regions of white matter from the white matter mask; This is the minimum gray value of all voxels within the high-signal region of white matter for clinical labeling. The formula takes the average of the normal white matter average gray value and the lowest gray value of the high-signal region of white matter as the optimal threshold. Its physical meaning is that the optimal segmentation threshold should be located at the midpoint between the upper edge of the normal white matter gray value distribution and the lower edge of the high-signal region of white matter, so as to achieve the maximum distinguishability between the two types of tissues.

[0199] Using the six extracted grayscale features as input variables, the optimal threshold is obtained by reverse engineering. A linear regression model was chosen as the output variable for training. Five-fold cross-validation was used during training to evaluate the model's generalization performance. The linear regression model was chosen because of its transparent model structure and clear thresholding. It can be represented as a linear weighted combination of six clinical features. The regression coefficient of each feature directly reflects the direction and magnitude of its contribution to the threshold decision, and has full interpretability.

[0200] Second stage: Extracting high-reliability white matter high-signal seed points based on the optimal threshold.

[0201] Based on the optimal threshold K obtained in the first stage, voxels with gray values ​​strictly greater than K within the white matter mask are marked as candidate seed points, and noisy seed points are filtered out through 3D morphological opening operation and neighborhood gray value verification.

[0202] Among them, 3D morphological opening operation filtering refers to selecting 2×2×2 spherical structural elements to perform morphological erosion and dilation operations in sequence. First, the erosion operation removes isolated seed points with a volume ≤ 3 voxels, and then the dilation operation restores the structure of the real seed points.

[0203] The neighborhood grayscale verification operation refers to calculating the average grayscale value within the 3×3×3 neighborhood of each candidate seed point. If the average grayscale value of a candidate seed point is ≤K, it means that the seed point only has a high grayscale value of its own and there is no continuous high signal area around it. It is likely to be local noise or edge artifacts, so it is removed. Only seed points with an average grayscale value of neighborhood greater than K are retained to ensure that there is a continuous high signal environment around the seed point.

[0204] The third stage: Three-dimensional conditional dilation operation based on hierarchical convergence rules.

[0205] This stage uses the highly reliable seed points extracted in the second stage as a foundation to perform a three-dimensional conditional dilation operation. Relying on the grayscale continuity of the high-signal region in white matter, the seed points are gradually expanded outward to cover the entire high-signal region of white matter. At the same time, to avoid introducing false positives in non-white matter regions due to over-dilation, or missing high-signal voxels at the edges due to premature convergence, this stage designs a layered convergence rule to dynamically control the dilation iteration process, ultimately generating a high-signal segmentation mask that fits the actual anatomical morphology.

[0206] The input to this stage of the algorithm is a FLAIR brain parenchyma image, an aligned white matter mask, and a set of WMH seed points. The output is a high-signal segmentation mask that has undergone dilation and smoothing. The specific execution logic is as follows:

[0207] First, the structuring element for the dilation operation is determined, as its selection directly affects the efficiency of dilation and the smoothness of the segmentation edges. This embodiment uses a 3×3×3 spherical structuring element, with all voxels within it assigned a value of 1. The advantage of the spherical structuring element is that it allows the dilation operation to expand uniformly in all directions of three-dimensional space, avoiding deformation of the segmented region due to excessive dilation in a single direction. Simultaneously, the 3×3×3 size balances dilation speed and detail preservation, ensuring both efficient execution of the dilation operation and preventing the omission of edge voxels in high-signal areas, thus preserving a complete edge morphology foundation for subsequent smoothing processing.

[0208] After selecting the structural elements, conditional dilation iteration is initiated. During the iteration, the segmentation region is not expanded indiscriminately. Instead, a dual-conditional screening is performed on the 3×3×3 neighborhood voxels of the currently dilated high-signal region. Only voxels that simultaneously meet the following two conditions are included in the segmentation region: First, the voxel must be located within the white matter mask, i.e., its voxel value is 1. This constraint strictly limits the dilation operation to within the white matter tissue, avoiding invasion into non-white matter regions such as cerebrospinal fluid and gray matter, thus eliminating false positives at the spatial level. Second, the voxel's FLAIR gray value is ≥ K-5. Here, the gray value threshold is slightly lower than the optimal threshold K. The core reason is that the gray value of edge voxels in white matter high-signal regions is usually slightly lower than that of the core region, but still falls within the range of true high signal. This setting effectively avoids missing edge high-signal voxels and improves the integrity of the segmentation.

[0209] After each dilation iteration, the volume change ΔV of the high-signal segmentation mask is calculated. ΔV is calculated as the difference between the segmentation volume after the current iteration and the volume of the previous iteration. The volume change is the core indicator for determining whether dilation has converged. Its numerical trend directly reflects whether the dilation operation has covered all real high-signal regions. When ΔV tends to stabilize, it indicates that the dilation has nearly completely covered the high-signal region, and iteration can be considered for termination.

[0210] To accommodate high-signal regions of white matter of varying sizes and avoid premature convergence or overexpansion caused by fixed convergence rules, this embodiment employs a convergence rule of "hierarchical threshold + priority," with the specific determination process as follows:

[0211] First, calculate the total volume V of the high-signal segmentation mask after the current iteration. tota； If V total For regions with fewer than 500 voxels, corresponding to small-volume high-signal regions, the volume change is inherently small. Using a relative threshold could lead to premature convergence of expansion, missing some edge voxels. Therefore, only "ΔV ≤ 10 voxels for three consecutive iterations" is used as the convergence condition to ensure complete coverage of these small-volume high-signal regions. If V totalFor voxels greater than or equal to 500, corresponding to medium / large volume high-signal regions, a relative threshold is preferentially used for determination, i.e., when ΔV ≤ V total When the value reaches 0.5%, it indicates that the volume change has decreased to an extremely low level, and the dilation operation basically covers the entire high-signal region. At this point, the iteration is terminated directly. If this relative threshold condition is not met, the iteration continues until the fixed threshold condition of "ΔV ≤ 10 voxels for 3 consecutive times" is met, balancing the integrity of the segmentation and computational efficiency. When any of the above convergence conditions is met, the dilation iteration is stopped immediately. At this point, the obtained segmentation mask basically covers the entire high-signal concentration region, avoiding the problem of missed detection in small volume regions and effectively avoiding false positives due to excessive dilation in medium / large volume regions.

[0212] After convergence iteration, the initial high-signal segmentation mask is smoothed to eliminate jagged artifacts at the edges, making the segmentation boundaries more closely resemble the true anatomical morphology of high-signal white matter. This stage employs 3D Gaussian smoothing with a standard deviation of 0.8 mm. The core of Gaussian smoothing is to perform convolution operations between a Gaussian kernel and the mask image, applying a weighted average to the edge pixels. A standard deviation of 0.8 mm achieves a balance between smoothing effect and boundary preservation: effectively eliminating jagged artifacts at the edges without blurring the boundaries of high-signal areas due to excessive smoothing, thus preserving the overall spatial structure and morphological characteristics of the high-signal areas.

[0213] S500: Performs differential denoising and morphological screening on 3D residual images to output detection results: Denoises residual images of brain microbleeds using 3D anisotropic diffusion filtering: Performs 3D connected domain analysis and calculates sphericity assessment index to remove flat artifacts and vascular cross sections.

[0214] Three-dimensional anisotropic diffusion filtering is an adaptive denoising method based on partial differential equations. Its core idea comes from the physical analogy of the heat conduction equation: the image grayscale is regarded as a temperature field, and the grayscale values ​​are thermally diffused in space to tend to be smooth. However, unlike standard isotropic thermal diffusion (i.e., Gaussian smoothing), the anisotropic diffusion model introduces a conduction coefficient function that depends on the local gradient, making the diffusion process directionally selective—allowing sufficient diffusion in uniform noise regions with small gradients (equivalent to strong smoothing), while suppressing diffusion at real edges with large gradients (equivalent to maintaining sharpness), thus achieving the goal of denoising without blurring edges.

[0215] Morphological opening is a composite operation that sequentially performs morphological erosion and dilation. First, the erosion operation is performed to eliminate slender structures whose width is smaller than the diameter of the structuring element. Then, the dilation operation is performed to restore the original size of larger structures that were reduced but not completely eliminated by the erosion step.

[0216] Three-dimensional connected component analysis refers to the operation of aggregating all spatially connected foreground voxels into independent connected components in a three-dimensional binary image based on a preset voxel adjacency relationship (such as 6-connected, 18-connected, or 26-connected), and assigning a unique identifier to each connected component.

[0217] Sphericity is a dimensionless three-dimensional geometric topological index that measures how close the shape of a three-dimensional connected region is to an ideal sphere. Its definition is based on a corollary of the isoperimetric inequality.

[0218] The detection result refers to the set of three-dimensional connected regions in the residual image that is finally retained after all the above differential denoising and morphological screening processes. These regions are identified as candidate regions for suspected new lesions (white matter high signal lesions or microbleed lesions) and can be further reviewed by clinical radiologists.

[0219] Specifically, firstly, the three-dimensional residual image Anisotropic diffusion filtering is performed. This filtering process is achieved by solving a diffusion partial differential equation that depends on the local gradient.

[0220] For example, in this embodiment, the partial differential equation for three-dimensional anisotropic diffusion filtering can be expressed as follows:

[0221]

[0222] in, The virtual time parameter for diffusion evolution controls the intensity of the filter. The larger the value, the stronger the smoothing effect; This is the spatial gradient vector of the residual image at the current location; It is a divergence operator; The edge propagation coefficient function is preserved for factors that depend on the magnitude of the local gradient.

[0223] For example, in this embodiment, the edge-preserving conduction coefficient function can be expressed as follows:

[0224]

[0225] in, The magnitude of the local spatial gradient of the residual image at the current voxel; The gradient propagation threshold is a core parameter that distinguishes noise from edges and needs to be set according to the noise level of the residual map.

[0226] It should be noted that in the above formula: when the voxel is located in a uniform noise region, The exponent term approaches zero, at which point... The diffusion equation degenerates into a standard thermal equation, performing sufficient isotropic smoothing and effectively eliminating random noise jumps. When the voxel is located at the boundary of the actual lesion, When the exponent term takes a maximum negative value, at this time... In this direction, diffusion conduction is almost completely blocked, thus sharply preserving the boundary features of small lesions. Through this adaptive conduction coefficient design, the anisotropic diffusion filter achieves an automatic balance between powerful denoising in uniform regions and precise preservation at edges. The selection of the value is crucial: if it is too large, the edge protection will be insufficient and the edge of the lesion will be blurred; if it is too small, the noise suppression will be insufficient and the false positive rate will increase.

[0227] After anisotropic diffusion denoising, differential morphological screening was initiated on the residuals of cerebral microbleeds. The specific screening steps are as follows:

[0228] In the three-dimensional connected component analysis results of the SWI residual image, each independent connected component is extracted. .

[0229] In this context, a connected component refers to an independent region composed of all spatially directly or indirectly connected foreground voxels, based on a predefined three-dimensional voxel adjacency relationship (e.g., 26-connected), after binarization of the SWI residual image. Each connected component corresponds to a potential signal variation region in the residual image, which may correspond to a real microbleed lesion or a cross-sectional artifact of a vascular structure.

[0230] For a single connected component Calculate its volume and the actual surface area surrounding the connected component. .

[0231] Among them, volume This refers to the sum of the physical volumes of all voxels within the connected region, obtained by summing the physical volumes of all voxels within the region: Surface area It refers to the surface boundary of the connected region. The actual physical area is determined by the boundary area infinitesimal element. The integral yields: .

[0232] For example, in this embodiment, the sphericity of the three-dimensional connected domain can be calculated by the following formula:

[0233]

[0234] in, For connected components Volume; For connected components Surface area; Pi; It is an index of sphericity.

[0235] It should be noted that sphericity The range of values ​​is It reaches its maximum value only when the connected region is a perfect sphere. Its physical meaning is based on a corollary of the isoperimetric inequality: among all three-dimensional closed regions with the same volume, a sphere has the smallest surface area. Therefore, sphericity quantifies the degree to which a connected region deviates from a sphere by comparing the relationship between the actual volume and surface area of ​​the region.

[0236] It is understandable that real microbleeding lesions are formed in three-dimensional space by the uniform deposition of blood components in all directions, and their shape tends to be spherical. The value is relatively high, usually close to 1. Small blood vessels running through the brain appear as elongated tubular structures in three-dimensional space, with a cross-sectional area much smaller than their longitudinal length. This results in a significantly larger surface area than a sphere for the same volume, leading to a higher sphericity. .

[0237] Set a preset sphericity critical value (e.g.) ), directly remove Components with sphericity values ​​below a threshold are identified as residual images of blood vessel cross-sections that are elongated or flattened. When the sphericity of a connected region is below the threshold, it is identified as a cross-section of a blood vessel structure and removed; when the sphericity is above the threshold, it is retained as a suspected microbleeding lesion.

[0238] It should be noted that, through this geometric morphology judgment criterion based on rigorous mathematical definition, the algorithm can forcibly decouple and isolate the vascular structure from the actual bleeding point from the geometric dimension of physical three-dimensional space, without relying on unreliable grayscale threshold judgment, thus fundamentally improving the specificity of microbleed detection.

[0239] In addition, to enhance the robustness and adaptability of the method when dealing with images of patients with severe brain atrophy, this embodiment also includes an adaptive fault-tolerant mechanism adapted to technical scenarios where severe brain atrophy leads to abnormal boundary overlap:

[0240] When the ratio of the initial brain region mask volume extracted in step S200 to the standard brain volume estimated by downsampling deviates from the normal physiological range, it is determined that there is excessive local boundary deformation caused by severe brain tissue atrophy.

[0241] The initial brain region mask volume refers to the sum of the physical volumes of all voxels with a value of 1 in the initial brain region mask obtained after inversely transforming the probabilistic brain mask back to the individual space and binarizing it in the third stage of step S200. The standard brain volume refers to the expected normal brain tissue volume of the individual estimated based on the scaling parameters of the standard brain template and the affine transformation matrix during the downsampling and affine alignment process in the first stage of step S200. The normal physiological range refers to the normal range (e.g., between 0.75 and 1.10) of the ratio of brain volume to expected volume established based on brain volume statistics from a large sample population; values ​​exceeding this range are considered to indicate significant abnormalities in brain tissue volume.

[0242] Understandably, patients with severe brain atrophy experience a significant reduction in brain volume, with extremely deep sulci and gyri forming on the cortical surface, corresponding to highly curvature local structures. In this situation, the standard active contour evolution mechanism may face the following problem: if the curvature smoothing term has an excessive weight, the contour surface will be over-smoothed, making it impossible to deeply trace these atrophied sulci and gyri, resulting in the erroneous removal of brain tissue at the bottom of the sulci.

[0243] At this point, before generating the nonlinear deformation field, a punitive constraint logic is initiated: local mutual information is recalculated based on the tissue density distribution of the current downsampled image to update the estimation accuracy of the affine transformation matrix for the shrunken region; the smoothing weight term controlling the contour curvature in the active contour evolution mechanism (i.e., the average curvature term in the aforementioned level set evolution equation) is adaptively reduced. The weighting coefficients allow the contour surface to more sharply track the boundaries of atrophic grooves and gyri, rather than being overly smoothed out; and the influence range of the anatomical physical boundaries calculated by Laplacian edge detection is dynamically reduced (i.e., in the stopping function). Medium-enlarged Gaussian pre-smoothing kernel Standard deviation Or lower the contrast control threshold This allows the evolutionary profile to be forcibly traced deep into the recessed sulci regions of the brain.

[0244] It should be noted that through this adaptive fault-tolerant mechanism, when a severely shrunken area is detected, the curvature term weight is adaptively reduced and the sensitivity range of edge detection is adjusted. This allows the contour surface to track the shrunken boundary more deeply with a sharper shape, avoiding the tissue at the shrunken area being erroneously cut off, thus achieving a dynamic balance between smoothness and boundary tracking accuracy.

[0245] Example 2

[0246] This embodiment discloses a magnetic resonance imaging (MRI) image follow-up system for amyloid-related imaging abnormalities. This system can be integrated into an electronic device, such as a terminal or server. The terminal can be a medical imaging workstation, a smart diagnostic terminal, a tablet computer, a laptop computer, or a personal computer. The server can be a single server, a server cluster consisting of multiple servers, a dedicated high-performance computing server deployed within a hospital radiology department's local area network, or a cloud-based medical image processing platform.

[0247] In some embodiments, the magnetic resonance imaging follow-up system for amyloid-related imaging abnormalities can also be integrated into multiple electronic devices. For example, the lesion follow-up detection system based on spatial alignment and residual analysis can be integrated into multiple servers, with multiple servers implementing the method in Embodiment 1 of this application.

[0248] In some embodiments, the server may also be implemented as a terminal.

[0249] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method of magnetic resonance image follow-up of amyloid-related imaging abnormalities, characterized in that, The magnetic resonance image follow-up method includes: Acquire multimodal nuclear magnetic resonance imaging data, extract white matter control points from the multimodal nuclear magnetic resonance imaging data and construct a bias field fitting model, perform bias field correction and intensity normalization processing based on the bias field fitting model, and generate a first image set; Global alignment and edge detection are performed on the first image set to extract anatomical physical boundaries. After inversely transforming the probabilistic brain mask to the individual space, active contour evolution is performed with the anatomical physical boundary as the evolution stopping condition to complete the skull dissection and boundary refinement, and generate brain region mask data. Within the brain parenchyma defined by the brain region mask data, branch deformation field registration processing is performed on the first image set to generate a spatially aligned second image set. The average intensity of normal brain parenchyma regions in the second image set is statistically analyzed to calculate the scaling factor. Based on the scaling factor, linear contrast stretching and voxel-by-voxel difference operations are performed on the second image set to complete the residual subtraction processing and obtain a three-dimensional residual image. Anisotropic diffusion filtering is applied to the three-dimensional residual image for denoising. Based on the type of image to be detected, the corresponding differential geometric screening features are called to perform differential denoising and morphological screening processing, and the detection results are output.

2. The method for magnetic resonance imaging follow-up of amyloid-related imaging abnormalities according to claim 1, characterized in that, Extracting white matter control points from the multimodal nuclear magnetic resonance imaging data and constructing a bias field fitting model includes: For each frame of the multimodal MRI image data, an axial slice grayscale histogram is constructed by axially slicing the image. After Gaussian smoothing of the axial slice grayscale histogram, the average intensity of the white matter region of each axial slice is determined by a peak finding algorithm. Based on the average intensity of the white matter region, white matter control points are automatically detected, and cubic spline fitting is performed on the intensity coefficients of effective slices along the vertical axis. A Voronoi mosaic is established using each white matter control point as the generation kernel. Non-control point voxels are assigned the reference intensity value of the nearest white matter control point to obtain the initial discrete estimate of the bias field. A soap bubble smoothing iteration is performed on the initial discrete estimate of the bias field to eliminate gray-level step jumps at the boundaries of adjacent Voronoi regions, thereby obtaining the converged bias field.

3. The method for magnetic resonance imaging follow-up of amyloid-related imaging abnormalities according to claim 1, characterized in that, The step of performing global alignment and edge detection on the first image set to extract anatomical physical boundaries includes: The images in the first image set are downsampled, and rigid body registration and affine registration based on maximum mutual information are performed with a standard brain template as the reference target to obtain the global transformation matrix; The sum of the second spatial derivatives of each voxel is calculated by applying the Laplacian operator to the images in the first image set: ； in, The grayscale function of the three-dimensional image after intensity normalization in the first image set is given. For the Laplace operator; Detection The voxel intersection where the numerical sign is flipped is taken as the zero-crossing point, and the closed surface formed by the zero-crossing points is taken as the anatomical physical boundary.

4. The method for magnetic resonance imaging follow-up of amyloid-related imaging abnormalities according to claim 3, characterized in that, The process of inversely transforming the probabilistic brain mask to the individual space, and then performing active contour evolution with the anatomical physical boundary as the evolutionary stopping condition, includes: Using the global transformation matrix and the nonlinear deformation field generated by the cross-correlation-based symmetric normalization model, the probabilistic brain mask is inversely transformed back to the individual space, and the inversely transformed probabilistic brain mask is binarized to obtain the initial brain region mask. Using the edges of the initial brain region mask as the initial contour, calculate the local gradient magnitude of the corresponding image in the first image set, and construct a stopping function that depends on the local gradient magnitude. The initial profile evolves along the normal direction through a level set partial differential equation. When the initial contour evolves to the anatomical physical boundary, the local gradient magnitude reaches a local maximum, the value of the stopping function approaches zero, and the contour evolution is automatically frozen.

5. The method for magnetic resonance imaging follow-up of amyloid-related imaging abnormalities according to claim 1 or 4, characterized in that, The process of completing skull dissection and boundary refinement to generate brain region mask data includes: Perform a 3D morphological closing operation on the frozen contour to fill the holes inside the contour; The brain region mask data is generated by extracting the largest connected component through three-dimensional connected component analysis and removing isolated tissue blocks with a volume smaller than a set threshold.

6. The method for magnetic resonance imaging follow-up of amyloid-related imaging abnormalities according to claim 1, characterized in that, The branched deformation field registration process includes a first registration pathway for detecting high signal intensity in brain white matter and a second registration pathway for detecting brain microbleeds. The first registration path includes: Using the baseline T1-weighted image as the intermediate mode, a synchronous cross-modal spatial mapping between the baseline T1-weighted image and the baseline FLAIR image is established through a symmetric normalization model. The FLAIR images during the follow-up period were registered to the T1-weighted images during the follow-up period, and the time-dimensional displacement vector field was obtained by registering the T1-weighted images during the baseline period with the T1-weighted images during the follow-up period. The synchronous cross-modal space mapping and the time-dimensional displacement vector field are matrix concatenated to synthesize a global nonlinear deformation field. The global nonlinear deformation field is then used to resample the follow-up FLAIR image to align it to the coordinate space of the baseline FLAIR image. The second registration path includes: Same-modal registration and resampling were performed directly on the follow-up SWI images after craniotomy and the baseline SWI images.

7. The method for magnetic resonance imaging follow-up of amyloid-related imaging abnormalities according to claim 1, characterized in that, The step of calculating the scaling factor by statistically analyzing the average intensity of normal brain parenchyma regions in the second image set includes: Within the range defined by the brain region mask data, after excluding voxels whose gray values ​​deviate from the normal distribution, the average intensity of the normal brain parenchyma region in the baseline image and the average intensity of the normal brain parenchyma region in the follow-up image are calculated respectively. Multiply the gray value of each voxel in the follow-up images by a scaling factor. To complete the linear contrast stretching.

8. The method for magnetic resonance imaging follow-up of amyloid-related imaging abnormalities according to claim 1, characterized in that, The step of calling the corresponding differential geometric screening features based on the image type to be detected to perform differential denoising and morphological screening processing includes: When the image type to be detected is high signal in brain white matter, after performing linear contrast stretching on the FLAIR images in the second image set based on the scaling factor, within the white matter mask constraint range aligned to FLAIR space, clinically interpretable gray-level statistical features are extracted and input into a pre-trained linear regression prior model to adaptively calculate the segmentation threshold. Based on the segmentation threshold, seed points are extracted within the white matter mask, and conditional dilation is performed until the hierarchical convergence rule is met, at which point the dilation stops. Three-dimensional connected component segmentation and morphological opening are performed on the dilation results to remove isolated noise points. When the image type to be detected is brain microbleed, three-dimensional connected component analysis is performed on the three-dimensional residual image to extract each connected component, the volume and surface area of ​​each connected component are calculated, and the sphericity is calculated based on the volume and surface area. Connected components with sphericity lower than a preset threshold are removed.

9. The method for magnetic resonance imaging follow-up of amyloid-related imaging abnormalities according to claim 1 or 4, characterized in that, The magnetic resonance image follow-up method also includes an adaptive fault-tolerant step: When the ratio of the volume of the initial brain region mask to the standard brain volume estimated based on the global transformation matrix deviates from the preset normal physiological range, it is determined that there is severe brain tissue atrophy. In response to the determination, the following adjustments are performed: The local mutual information is recalculated based on the tissue density distribution of the current downsampled image to update the global transformation matrix; The weighting coefficient of the mean curvature term in the level set partial differential equation is adaptively reduced so that the contour surface tracks the boundary of the shrinkage groove with a sharper shape. Adjusting the standard deviation of the Gaussian smoothing kernel in the stopping function Or lower the contrast control threshold This expands the detection sensitivity range of the anatomical physical boundaries.

10. A magnetic resonance imaging follow-up system for amyloid-related imaging abnormalities, characterized in that, The magnetic resonance imaging follow-up system includes: processor; A memory storing a computer program that, when executed by a processor, implements a magnetic resonance imaging follow-up method for amyloid-related imaging abnormalities as described in any one of claims 1 to 9.

Images