Multi-time-point neural image comparative analysis method and system

Abstract

The invention relates to the technical field of medical image processing, in particular to a multi-time-point neural image comparative analysis method and system.The method comprises the following steps that a first three-dimensional neural image matrix and a second three-dimensional neural image matrix are obtained, rigid transformation matrix registration is calculated, and a third three-dimensional neural image matrix is generated; constructing a three-dimensional elastic grid model and adjusting node coordinates based on a mutual information criterion; generating a three-dimensional displacement vector field; constructing a Jacobian matrix and calculating determinant values; executing natural logarithmic transformation to generate a microscopic volume change rate map; and calculating a logarithm Jacobian mean value and determining a volume change state. According to the method, pose deviation is eliminated by calculating a rigid transformation matrix, a vector field reflecting fine displacement is constructed by utilizing a mutual information criterion, voxel-level quantitative analysis is realized in combination with a Jacobian matrix and logarithmic transformation, and minimal lesion recognition precision and diagnosis efficiency are improved in cooperation with a standard brain anatomy template.

Description

Technical Field

[0001] This invention relates to the field of medical image processing technology, and in particular to a method and system for comparative analysis of neural images at multiple time points. Background Technology

[0002] The field of medical image processing technology refers to the technical scope of acquiring, displaying, storing, transmitting, and post-processing digital image signals acquired by medical imaging equipment using computer hardware and software systems. Among these, the traditional multi-time-point neural image comparison analysis method involves medical staff acquiring brain scan data from patients at different time points using a magnetic resonance imaging (MRI) scanner. The DICOM format image files are loaded onto the monitor of a reading workstation. The operator manually traces the outlines of lesions or specific neural structures layer by layer on the tomographic scan image using a mouse. The volume index is obtained by calculating the total number of voxels within the traced area, and the volume values ​​from the two time points are entered into a paper or electronic medical record form for comparison.

[0003] Current technology relies on manual drawing of lesion outlines layer by layer on tomographic images. This is limited by the operator's subjective experience and hand stability, resulting in differences in the delineated boundaries and low reproducibility. Manual image reading cannot eliminate spatial misalignment caused by differences in patient scanning position, leading to inconsistent comparative analysis benchmarks. The calculation method based on the total number of voxels in the delineated area can only reflect overall macroscopic differences and cannot capture the details of microscopic deformation or local atrophy within nerve tissue. Furthermore, the purely manual processing of massive tomographic data is time-consuming and lengthy, severely restricting the efficiency of clinical diagnosis and the implementation of large-scale disease progression analysis. Summary of the Invention

[0004] The purpose of this invention is to overcome the shortcomings of existing technologies and to propose a multi-time point neural image comparison and analysis method and system.

[0005] To achieve the above objectives, the present invention employs the following technical solution: a multi-time-point neural image comparison and analysis method, comprising the following steps: S1: Use a magnetic resonance imaging (MRI) instrument to acquire a first three-dimensional neural image matrix and a second three-dimensional neural image matrix, calculate the rigid transformation matrix of the second three-dimensional neural image matrix relative to the first three-dimensional neural image matrix, register the second three-dimensional neural image matrix based on the rigid transformation matrix, and generate a third three-dimensional neural image matrix. S2: Construct a three-dimensional elastic mesh model using the first three-dimensional neural image matrix and the third three-dimensional neural image matrix, adjust the node coordinates of the three-dimensional elastic mesh model based on the mutual information criterion, calculate the displacement vector of the third three-dimensional neural image matrix, and generate a three-dimensional displacement vector field. S3: Calculate the axial partial derivative of the three-dimensional displacement vector field, construct the Jacobian matrix based on the axial partial derivative, calculate the determinant of the Jacobian matrix, perform a natural logarithmic transformation on the determinant, and generate a microscopic volume change rate map. S4: Map the standard brain anatomical structure template onto the first three-dimensional neural image matrix to delineate the region of interest, extract the logarithmic Jacobian determinant value of the region of interest in the microscopic volume change rate map, calculate the arithmetic mean of the logarithmic Jacobian determinant value, compare the arithmetic mean with the lesion threshold, and determine the volume change state of the region of interest.

[0006] As a further aspect of the present invention, step S1 specifically comprises: S11: At the baseline time point, the original image data of the subject's brain is acquired using a magnetic resonance imaging device, and the original image data is processed by skull dissection and grayscale non-uniformity correction to generate a first three-dimensional neural image matrix. At the follow-up time point, brain image data is acquired using the same imaging parameters and the same preprocessing is performed to generate a second three-dimensional neural image matrix. S12: Extract the geometric center coordinates of the first three-dimensional neural image matrix and the second three-dimensional neural image matrix, calculate the translation vector and rotation matrix that can make the two geometric centers coincide and the principal axis directions aligned based on the least squares method, and combine the translation vector and rotation matrix to construct a rigid transformation matrix. S13: Based on the spatial mapping relationship defined by the rigid transformation matrix, perform inverse coordinate transformation on each voxel coordinate in the second three-dimensional neural image matrix, and use the trilinear interpolation algorithm to calculate the gray value of the transformed coordinate position to generate the third three-dimensional neural image matrix.

[0007] As a further aspect of the present invention, step S2 specifically comprises: S21: Define a set of uniformly distributed grid control points in the image space of the first three-dimensional neural image matrix, construct a three-dimensional elastic grid model using B-spline basis functions, and map the third three-dimensional neural image matrix as a floating image to the coordinate system defined by the model; S22: Calculate the joint gray-level histogram between the first three-dimensional neural image matrix and the third three-dimensional neural image matrix affected by mesh deformation, calculate the similarity measure between the two images according to the mutual information criterion, and iteratively adjust the spatial coordinates of multiple control points in the three-dimensional elastic mesh model using the gradient descent optimization algorithm until the similarity measure reaches the convergence condition. S23: Obtain the offset of the control point coordinates after optimization and convergence relative to the initial grid position, and calculate the non-rigid deformation component of each voxel point in the image space using the B-spline tensor product formula to generate a three-dimensional displacement vector field.

[0008] As a further aspect of the present invention, step S3 specifically comprises: S31: Traverse each voxel node in the three-dimensional displacement vector field, and use the central difference method to calculate the rate of change of the displacement vector at the node in the X-axis, Y-axis and Z-axis directions of the spatial rectangular coordinate system, and obtain nine axial partial derivatives; S32: Construct a Jacobian matrix using the nine axial partial derivatives, calculate the determinant value of the Jacobian matrix to quantify the local volume expansion or contraction factor of the voxel point between two time points, and obtain the determinant value. S33: Perform a logarithmic operation with the natural constant as the base on the determinant value, map the volume shrinkage change to a negative value and the volume expansion change to a positive value, and construct a microscopic volume change rate map based on the operation result.

[0009] As a further aspect of the present invention, step S4 specifically comprises: S41: Obtain a standard brain anatomical structure template including preset anatomical partition labels, calculate the affine transformation parameters and nonlinear deformation parameters of the standard brain anatomical structure template relative to the first three-dimensional neural image matrix, and use the parameters to inversely deform the standard brain anatomical structure template and resample it to the coordinate space where the first three-dimensional neural image matrix is ​​located. S42: Generate a binary mask matrix based on the anatomical partition labels in the mapped standard brain anatomy template, perform a Boolean AND operation between the binary mask matrix and the microscopic volume change rate map, delineate the region of interest, and extract the logarithmic Jacobian determinant values ​​of all effective voxels within the region; S43: Calculate the logarithmic Jacobian determinant value of all voxels in the region of interest and its arithmetic mean. Compare the arithmetic mean with a preset lesion threshold. If the mean exceeds the threshold range, generate an abnormal state identifier to determine the volume change status of the region of interest.

[0010] As a further aspect of the present invention, the process of constructing the rigid transformation matrix includes: Obtain the voxel grayscale sets of the first three-dimensional neural image matrix and the second three-dimensional neural image matrix, construct a rigid registration objective function based on normalized mutual information, and use a multi-resolution pyramid strategy to perform hierarchical downsampling processing on the image; At each resolution level, the Powell optimization algorithm is used to search for the optimal six-degree-of-freedom rigid parameters, which include translation along the three coordinate axes and rotation angles around the three coordinate axes. The convergent parameter combination at the highest resolution level is converted into a four-by-four homogeneous transformation matrix to establish the rigid transformation matrix.

[0011] As a further aspect of the present invention, the node coordinate adjustment process of the three-dimensional elastic mesh model includes: Obtain the probability distribution of the reference image in the current iteration step. Floating image probability distribution and the joint probability distribution of the two under the current mesh deformation. According to the following formula: ; Calculate the mutual information objective function as a registration similarity measure; in, Representative reference image With floating images Mutual information values ​​between them Represents the gray level of the reference image. Represents the grayscale level of the floating image. Represents the numerical value of the joint probability distribution. Represents the edge probability distribution of the reference image. Represents the edge probability distribution of the floating image; The negative gradient direction of the mutual information objective function relative to the coordinates of the grid control points is calculated. The coordinates of the control points are updated along the negative gradient direction with an adaptive step size. The histogram and gradient are repeatedly calculated until the rate of change of the objective function is lower than the preset convergence threshold, thereby generating a non-rigid deformation grid.

[0012] As a further aspect of the present invention, the numerical calculation process of the determinant of the Jacobian matrix includes: Obtain the three-dimensional displacement vector field in coordinates displacement components at the location , and Construct the local deformation gradient tensor matrix as shown below: ; in, The value representing the determinant of the Jacobian matrix is... Displacement vector Quantity in Axial partial derivative of the direction, Displacement vector Quantity in Axial partial derivative of the direction, Displacement vector Quantity in The axial partial derivatives of the direction, and the remaining elements in the matrix represent the cross partial derivatives; The determinant value of the matrix is ​​calculated to obtain the local volume change factor at the coordinate point. When the value is greater than 1, it represents local volume expansion. When the value is less than 1 but greater than 0, it represents local volume contraction. The determinant value is generated.

[0013] As a further aspect of the present invention, the process of determining the volume change state of the region of interest includes: Obtain the target atrophy detection threshold for the hippocampal region and the target expansion detection threshold for the ventricle region, and calculate the average of the logarithmic Jacobian determinant values ​​of all voxels in the hippocampal region and the average of the logarithmic Jacobian determinant values ​​of all voxels in the ventricle region. Determine if the average number of hippocampal regions is less than the atrophy detection threshold; if so, generate diagnostic prompts. Determine whether the average number of ventricular regions is greater than the expansion detection threshold. If so, generate diagnostic prompts and generate a report on the volume change status of the region of interest based on the above judgment results.

[0014] A multi-time-point neural image comparison and analysis system is provided, the system being used to implement the aforementioned multi-time-point neural image comparison and analysis method, the system comprising: The image acquisition and rigid registration module is used to acquire a first three-dimensional neural image matrix and a second three-dimensional neural image matrix using a magnetic resonance imaging instrument, calculate the rigid transformation matrix of the second three-dimensional neural image matrix relative to the first three-dimensional neural image matrix through geometric center alignment and iterative optimization algorithms, and perform spatial resampling and registration operations on the second three-dimensional neural image matrix based on the rigid transformation matrix to generate a third three-dimensional neural image matrix. The elastic mesh and vector field generation module is used to construct a three-dimensional elastic mesh model based on B-splines using the first three-dimensional neural image matrix and the third three-dimensional neural image matrix as input data, iteratively adjust the node coordinates of the three-dimensional elastic mesh model to match local features of the image based on the maximum mutual information criterion, and calculate the displacement vector of the third three-dimensional neural image matrix relative to the reference image by interpolation to generate a three-dimensional displacement vector field. The micro-volume change analysis module is used to calculate the axial partial derivative of the three-dimensional displacement vector field using the central difference algorithm, construct a Jacobian matrix describing the local deformation gradient based on the axial partial derivative, calculate the determinant of the Jacobian matrix and perform a natural logarithmic transformation on it, thereby generating a micro-volume change rate map. The region status assessment and determination module is used to map the standard brain anatomical structure template to the first three-dimensional neural image matrix space through nonlinear transformation to delineate the region of interest. It uses masking technology to extract the logarithmic Jacobian determinant value of the corresponding region in the microscopic volume change rate map, calculates its arithmetic mean and compares it with the preset lesion threshold, and outputs the volume change status of the region of interest.

[0015] Compared with the prior art, the advantages and positive effects of the present invention are as follows: In this invention, spatial registration of images at different time points is performed by calculating a rigid transformation matrix to eliminate pose deviations caused by differences in scanning position. The mutual information criterion is used to drive the deformation of the elastic mesh model to construct a three-dimensional vector field reflecting subtle displacements within the tissue. Combined with the microscopic volume change rate map generated by the Jacobian matrix determinant and natural logarithmic transformation, voxel-level quantitative analysis of the local expansion or contraction state of neural structures is achieved. With the help of a standard brain anatomy template, the region of interest is automatically delineated, which improves the accuracy of identifying minute lesions and the efficiency of automated diagnostic analysis while ensuring the consistency of anatomical positioning. Attached Figure Description

[0016] Figure 1 This is a flowchart of the multi-time-point neural image comparison and analysis method of the present invention; Figure 2 This is a flowchart of the rigid transformation matrix construction and image registration process of the present invention; Figure 3 This is a flowchart of the three-dimensional elastic mesh model adjustment and displacement field generation process of the present invention; Figure 4 This is a flowchart of the process for generating the microscopic volume change rate spectrum of the present invention; Figure 5 This is a flowchart for determining the volume change state of the region of interest in this invention. Detailed Implementation

[0017] To make the objectives, technical solutions, and advantages of this invention clearer, the software-based technical solution is described in detail below with reference to system architecture diagrams and embodiments. It should be understood that the specific embodiments described herein are only for explaining the technical solutions of this invention and do not constitute a limitation on the scope of protection.

[0018] In the description of this invention, the system architecture relationships or data processing flows indicated by terms such as "layer," "module," "interface," "data flow," "client," and "server" are all defined based on the architecture diagram or flowchart corresponding to the embodiments. This way of describing is only used to clearly illustrate the logical relationships between the elements in the technical solution, and not to limit the physical deployment form. The term "multiple" includes two or more technical units, including but not limited to multiple data nodes, processing threads, service instances, or functional components and other scalable elements. The specific number is determined according to the actual business scenario and needs to be specifically specified.

[0019] Please see Figure 1 and Figure 2 This invention provides a technical solution: a method for comparative analysis of neural images at multiple time points, comprising the following steps: S1: Use a magnetic resonance imaging (MRI) scanner to acquire a first three-dimensional neural image matrix and a second three-dimensional neural image matrix. Calculate the rigid transformation matrix of the second three-dimensional neural image matrix relative to the first three-dimensional neural image matrix. Based on the rigid transformation matrix, register the second three-dimensional neural image matrix to generate a third three-dimensional neural image matrix.

[0020] S11: At the baseline time point, the original image data of the subject's brain is acquired using a magnetic resonance imaging (MRI) instrument, and the original image data is processed by skull dissection and grayscale non-uniformity correction to generate a first three-dimensional neural image matrix. At the follow-up time point, brain image data is acquired using the same imaging parameters and the same preprocessing is performed to generate a second three-dimensional neural image matrix.

[0021] S12: Extract the geometric center coordinates of the first three-dimensional neural image matrix and the second three-dimensional neural image matrix. Calculate the translation vector and rotation matrix that can make the two geometric centers coincide and the principal axis directions aligned based on the least squares method. Combine the translation vector and rotation matrix to construct a rigid transformation matrix.

[0022] S13: Based on the spatial mapping relationship defined by the rigid transformation matrix, perform inverse coordinate transformation on the coordinates of each voxel in the second three-dimensional neural image matrix, and use the trilinear interpolation algorithm to calculate the gray value of the transformed coordinate position to generate the third three-dimensional neural image matrix.

[0023] The process of constructing a rigid transformation matrix includes: Obtain the voxel grayscale sets of the first and second three-dimensional neural image matrices, construct a rigid registration objective function based on normalized mutual information, and perform hierarchical downsampling processing on the images using a multi-resolution pyramid strategy.

[0024] At each resolution level, the Powell optimization algorithm is used to search for the optimal six-degree-of-freedom rigid parameters. The six-degree-of-freedom rigid parameters include the translation along the three coordinate axes and the rotation angles around the three coordinate axes. The parameter combination that converges at the highest resolution level is converted into a four-by-four homogeneous transformation matrix to establish the rigid transformation matrix.

[0025] In specific clinical implementation scenarios, a superconducting magnetic resonance scanner with a magnetic field strength of 3.0T was first selected, and a 64-channel head and neck combined coil was configured. At the baseline time point... The subject underwent an initial brain scan. The scan sequence employed a three-dimensional magnetized fast gradient echo sequence, with specific imaging parameters strictly set as follows: repetition time 2300 ms, echo time 2.98 ms, reversal time 900 ms, flip angle 9 degrees, and field of view [missing information]. The voxel isotropic resolution is After the scan is complete, the acquired DICOM format data is imported into an image processing workstation and parsed to generate the first three-dimensional neural image matrix. At the follow-up time point ,in and The time interval was 12 months, and the same subject was scanned twice using the same imaging parameters to generate a second three-dimensional neuroimaging matrix. .

[0026] After acquiring the raw data, skull dissection was performed using a brain extraction tool. The algorithm was set with a high score of 0.5 and a vertical gradient of 0 to remove signals from non-brain tissues such as the skull, scalp, and meninges, retaining only the brain parenchyma. Next, the N4ITK algorithm was used for grayscale non-uniformity correction, with a B-spline fitting distance of 200 mm, a convergence threshold of 0.001, and an iteration count of [50, 50, 50, 50] to eliminate bias field effects caused by radio frequency field non-uniformity. After preprocessing, the program first traverses... and For all non-zero voxels, calculate the geometric center coordinates of the two images respectively. and Using the principle of least squares, a preliminary translation vector is constructed. This makes the centroids of the two images coincide in space, and this state is used as the initial solution for rigid registration.

[0027] The aforementioned N4ITK algorithm refers to an improved algorithm based on nonparametric non-uniform intensity normalization (N3), which iteratively optimizes the B-spline bias field to correct low-frequency gray-level non-uniformity in MRI images. ITK indicates that it is implemented based on the InsightSegmentationandRegistrationToolkit.

[0028] Subsequently, the system enters the stage of fine-tuning the rigid transformation matrix. A multi-resolution pyramid strategy is employed, with four pyramid layers. At each layer, the image undergoes Gaussian smoothing and downsampling. The standard deviations for the four pyramid layers are 4.0, 2.0, 1.0, and 0.5, respectively, to prevent the optimization process from getting trapped in local minima. A registration objective function based on normalized mutual information is constructed. The optimization process uses the Powell algorithm, which does not require calculating the derivative of the objective function and searches for the optimal solution in the six-degree-of-freedom parameter space using the conjugate direction method. The six parameters specifically include the translation amounts along the X, Y, and Z axes. and rotation angles around the X, Y, and Z axes After convergence at the top level, i.e., the original resolution level, the final rigid transformation parameters are obtained. In this embodiment, the calculated translation parameters are: The rotation parameters are Convert these six parameters to homogeneous transformation matrix Using the spatial mapping relationship defined by this matrix, for Each voxel coordinate in Perform the inverse coordinate transformation and calculate its coordinates in the coordinate system. Corresponding position in space Since the calculated coordinates are usually non-integer, a trilinear interpolation algorithm is used. This algorithm calculates the gray value at a given location by weighted averaging of the voxel gray values ​​within a 26-neighborhood of the target coordinates. This generates a third 3D neural image matrix that is perfectly aligned with the first 3D neural image matrix space. Experimental data show that, after applying the above-mentioned multi-resolution and N4ITK correction strategies, the average target overlap rate of rigid registration reaches 98.5% within the whole-brain mask, providing an accurate anatomical benchmark for subsequent microscopic deformation analysis.

[0029] The Powell algorithm mentioned above is a direct search method for multidimensional unconstrained optimization problems. It finds the extreme points of a function by performing a one-dimensional search in a set of conjugate directions. It is particularly suitable for scenarios where the derivative of the objective function is difficult to calculate or is discontinuous.

[0030] Please see Figure 1 and Figure 3 S2: Construct a three-dimensional elastic mesh model using the first and third three-dimensional neural image matrices, adjust the node coordinates of the three-dimensional elastic mesh model based on the mutual information criterion, calculate the displacement vector of the third three-dimensional neural image matrix, and generate a three-dimensional displacement vector field.

[0031] S21: Define a set of uniformly distributed grid control points in the image space of the first three-dimensional neural image matrix, construct a three-dimensional elastic grid model using B-spline basis functions, and map the third three-dimensional neural image matrix as a floating image to the coordinate system defined by the model.

[0032] S22: Calculate the joint gray-level histogram between the first three-dimensional neural image matrix and the third three-dimensional neural image matrix affected by mesh deformation. Calculate the similarity measure between the two images based on the mutual information criterion. Iteratively adjust the spatial coordinates of multiple control points in the three-dimensional elastic mesh model using the gradient descent optimization algorithm until the similarity measure reaches the convergence condition.

[0033] S23: Obtain the offset of the control point coordinates after optimization and convergence relative to the initial grid position, and calculate the non-rigid deformation component of each voxel point in the image space using the B-spline tensor product formula to generate a three-dimensional displacement vector field.

[0034] The node coordinate adjustment process for a 3D elastic mesh model includes: Obtain the probability distribution of the reference image in the current iteration step. Floating image probability distribution and the joint probability distribution of the two under the current mesh deformation. According to the following formula: ; Calculate the mutual information objective function as a registration similarity measure. .

[0035] in, Representative reference image With floating images Mutual information values ​​between them Represents the gray level of the reference image. Represents the grayscale level of the floating image. Represents the numerical value of the joint probability distribution. Represents the edge probability distribution of the reference image. This represents the edge probability distribution of the floating image.

[0036] Calculate the negative gradient direction of the mutual information objective function relative to the coordinates of the grid control points, update the coordinates of the control points along the negative gradient direction with an adaptive step size, and repeat the calculation of the histogram and gradient until the rate of change of the objective function is lower than the preset convergence threshold, thereby generating a non-rigid deformation grid.

[0037] Building upon rigid alignment, we further capture the microscopic, non-rigid deformations of brain tissue over time. Firstly, in the first three-dimensional neural imaging matrix... Defined image spatial domain Within this model, a three-dimensional elastic mesh is constructed. The control points of the mesh are evenly distributed, and the spacing between the control points is set to [value missing]. The choice of this spacing is based on prior knowledge of the brain atrophy gradient; a 5mm spacing effectively balances computational complexity and local deformation capture capability. A continuous displacement field model is constructed using cubic B-spline functions as tensor product basis functions. The third three-dimensional neural imaging matrix is ​​then used... As a floating image, it is mapped to the coordinate system defined by the elastic grid.

[0038] The cubic B-spline function mentioned above refers to a basis function composed of piecewise polynomials. It has second-order continuous differentiability, which can guarantee the smoothness of the deformation field. It is often used in the construction of free-form deformation models in medical image registration.

[0039] During the optimization phase, the core task is to adjust the three-dimensional coordinates of the mesh control points so that the deformed... and To maximize mutual information between the two images, the calculation process first involves calculating the joint gray-level histogram of the two images within the current overlapping region. The number of bins for each gray level is set to 64 to ensure sufficient statistical significance in the histogram statistics. The histogram is then smoothed using the Parzen window method combined with a Gaussian kernel function, and the probability distribution of the reference image is estimated. Floating image probability distribution and joint probability distribution The mutual information objective function, used as a measure of registration similarity, is calculated according to the following formula. : ; in, Representative reference image With floating images Mutual information values ​​between them Represents the gray level of the reference image. Represents the grayscale level of the floating image. This represents the joint probability distribution of grayscale values ​​occurring simultaneously in the reference image and the floating image at corresponding spatial locations. The edge probability distribution representing the grayscale values ​​of the reference image. The edge probability distribution representing the grayscale values ​​of a floating image.

[0040] To maximize this measure, a stochastic gradient descent optimization algorithm is employed. The algorithm calculates the negative gradient direction of the objective function with respect to the displacement parameters of each grid control point. An initial learning rate is set. An adaptive step-rate decay strategy is employed, meaning that if the objective function value does not improve after five consecutive iterations, the learning rate is halved until it falls below the stopping threshold. .

[0041] In the specific implementation calculation, for coordinates as In a local region, the initial mutual information value was 0.85. After 50 iterations of optimization, the control point positions were fine-tuned, resulting in a high degree of overlap in local tissue texture, and the mutual information value increased to 1.12. The optimization process converged when the rate of change of the objective function was less than [a certain value]. Alternatively, when the maximum number of iterations (500) is reached, the final coordinates of all control points are locked. Then, based on the optimized control point displacements, the B-spline tensor product formula is used to calculate the coordinates of each voxel in the image space. Non-rigid deformation components The resulting three-dimensional displacement vector field covers the entire brain region, with a resolution consistent with the original image. Each vector in this field precisely describes the brain tissue from the reference time point. At the follow-up time point The physical displacement trajectory in spatial coordinates is represented. By introducing a bending energy regularization term and setting its weight to 0.01 to constrain mesh deformation, the generated displacement field is ensured to have topology-preserving properties, avoiding mesh folding. Experimental comparison data show that using a 5mm mesh spacing with a B-spline model, compared to the traditional fluid model, the registration accuracy at fine structures such as the hippocampus is reduced from 1.2mm to 0.6mm, improving the reliability of subsequent volume change analysis.

[0042] Please see Figure 1 and Figure 4 S3: Calculate the axial partial derivatives of the three-dimensional displacement vector field, construct the Jacobian matrix based on the axial partial derivatives, calculate the determinant of the Jacobian matrix, perform a natural logarithmic transformation on the determinant, and generate a microscopic volume change rate map.

[0043] S31: Traverse each voxel node in the three-dimensional displacement vector field, and use the central difference method to calculate the rate of change of the displacement vector at the node in the X-axis, Y-axis and Z-axis directions of the spatial rectangular coordinate system, and obtain nine axial partial derivatives.

[0044] S32: Construct the Jacobian matrix using nine axial partial derivatives, calculate the determinant of the Jacobian matrix to quantify the local volume expansion or contraction of a voxel between two time points, and obtain the determinant value.

[0045] S33: Perform a logarithmic operation with the natural constant as the base on the determinant value, map the volume shrinkage change to a negative value and the volume expansion change to a positive value, and construct a microscopic volume change rate map based on the operation results.

[0046] The numerical calculation process of the determinant of the Jacobian matrix includes: Obtain the three-dimensional displacement vector field in coordinates displacement components at the location , and Construct the local deformation gradient tensor matrix as shown below: ; in, The value representing the determinant of the Jacobian matrix. Displacement vector Quantity in Axial partial derivative of the direction, Displacement vector Quantity in Axial partial derivative of the direction, Displacement vector Quantity in The axial partial derivatives of the direction, and the remaining elements in the matrix represent the cross partial derivatives.

[0047] The determinant of the matrix is ​​calculated to obtain the local volume change factor at the coordinate point. When the value is greater than 1, it represents local volume expansion. When the value is less than 1 but greater than 0, it represents local volume contraction. The determinant value is generated.

[0048] The macroscopic displacement vector field is transformed into a quantitative map reflecting changes in the microscopic volume of tissues. The program first loads the whole-brain three-dimensional displacement vector field generated in step S2, and then traverses each voxel node in the field. The central difference method is used to calculate the rate of change of the displacement vector along the three principal axes in space to eliminate first-order truncation errors and improve calculation accuracy. Specifically, the partial derivatives in the X direction are calculated... Using the formula Perform calculations, where This represents the step size of one voxel unit. Similarly, the step size can be calculated. There are a total of nine axial partial derivative components.

[0049] Based on these nine partial derivatives, construct The Jacobian matrix. The determinant of this matrix is ​​numerical. It has a clear physical meaning: it quantifies the expansion or contraction factor of the tiny tissue volume represented by the voxel within a time interval. The local deformation gradient tensor matrix is ​​constructed according to the following formula, and the determinant value is calculated: ; in, The value representing the determinant of the Jacobian matrix. Represents determinant operations. Displacement vector Quantity in Axial partial derivative of the direction, Displacement vector Quantity in Axial partial derivative of the direction, Displacement vector Quantity in The axial partial derivatives of the direction, and the remaining elements in the matrix represent the cross partial derivatives.

[0050] like This indicates local volume expansion; if This indicates local volume contraction; if This indicates that the volume remains unchanged. To ensure the statistical distribution of the data conforms to a normal distribution, and to assign equal weights to contraction and expansion, for all valid... Value, that is In the case of natural constants, the execution is performed. Logarithmic operations with base 0, calculate At this point, negative values ​​intuitively represent contraction, positive values ​​represent expansion, and 0 represents no change.

[0051] The central difference method mentioned above is a numerical differentiation method that approximates the derivative of a target point by using the difference in function values ​​between two neighboring points before and after the target point. Its truncation error is on the order of the square of the step size.

[0052] Taking the actual calculation of a voxel point in the hippocampus region as an example: Assume that the partial derivative matrix obtained by the displacement component of this point and its neighborhood difference calculation through step S2 is, i.e., the matrix part in the above formula is: Calculate its determinant value to obtain According to calculations, This indicates that the voxel shrank to 94.1% of its original volume during the follow-up period. Further logarithmic transformation yielded... This value of -0.0608 is then written into the corresponding coordinate position of the microscopic volume change rate map. The final generated microscopic volume change rate map is a three-dimensional scalar field with the same dimensions as the original map. Each voxel value in the map precisely encodes the annualized volume change rate at that anatomical location. Table 1 shows the validation results of the Jacobian determinant values ​​for different tissue types.

[0053] Table 1. Validation data on the rate of change in microscopic volume of brain tissue regions. ; As shown in Table 1, the values ​​after logarithmic transformation can clearly distinguish between normal aging and pathological changes, providing a quantitative basis for subsequent automated diagnosis.

[0054] Please see Figure 1 and Figure 5S4: Map the standard brain anatomical structure template to the first three-dimensional neural image matrix to delineate the region of interest, extract the logarithmic Jacobian determinant value of the region of interest in the microscopic volume change rate map, calculate the arithmetic mean of the logarithmic Jacobian determinant value, compare the arithmetic mean with the lesion threshold, and determine the volume change state of the region of interest.

[0055] S41: Obtain a standard brain anatomical structure template including preset anatomical partition labels, calculate the affine transformation parameters and nonlinear deformation parameters of the standard brain anatomical structure template relative to the first three-dimensional neural image matrix, and use the parameters to inversely deform the standard brain anatomical structure template and resample it to the coordinate space where the first three-dimensional neural image matrix is ​​located.

[0056] S42: Generate a binary mask matrix based on the anatomical partition labels in the mapped standard brain anatomy template. Perform Boolean AND operation on the binary mask matrix and the microscopic volume change rate map to delineate the region of interest and extract the logarithmic Jacobian determinant values ​​of all effective voxels within the region.

[0057] S43: Calculate the logarithmic Jacobian determinant values ​​of all voxels within the region of interest and their arithmetic mean. Compare the arithmetic mean with the preset lesion threshold. If the mean exceeds the threshold range, generate an abnormal status indicator to determine the volume change status of the region of interest.

[0058] The process of determining the volume change state of the region of interest includes: Obtain the target atrophy detection threshold for the hippocampus region and the target expansion detection threshold for the ventricle region. Calculate the mean of the log-Jacobi determinant values ​​of all voxels in the hippocampus region and the mean of the log-Jacobi determinant values ​​of all voxels in the ventricle region.

[0059] Determine if the average number of hippocampal regions is less than the atrophy detection threshold; if so, generate diagnostic information.

[0060] Determine whether the average number of ventricular regions is greater than the expansion detection threshold. If so, generate diagnostic prompts and generate a report on the volume change status of the region of interest based on the above judgment results.

[0061] A standardized anatomical template was used to perform region-specific analysis of microscopic volume change rate maps to achieve automated assisted diagnosis of lesions. First, the system retrieved a pre-set MNI152 standard brain anatomical structure template, which included AAL anatomical region labels. The hippocampus region was labeled with IDs 37 and 38, and the lateral ventricle region with IDs 7901 and 7902. To accurately project these anatomical labels onto the subject's individual space, the standard template was calculated relative to a first three-dimensional neuroimaging matrix. The mapping parameters are obtained. This process involves affine transformation and subsequent nonlinear SyN deformation registration. After obtaining the inverse transform field, the label matrix of the AAL template is inversely resampled to... In the original coordinate space, a personalized anatomical partition mask is generated for the subject.

[0062] The aforementioned MNI152 standard brain anatomical template refers to the average brain template constructed by the Montreal Neurological Institute based on MRI scan data of 152 normal subjects. It is a commonly used spatially standardized reference system in neuroimaging research.

[0063] Using the generated binary mask matrix, voxel data from all regions of the hippocampus and lateral ventricle in the microscopic volume change rate map were extracted. For the approximately 3000 voxels covered by the hippocampal mask, their corresponding logarithmic Jacobian determinant values ​​were read. And calculate the arithmetic mean. Similarly, calculate the mean value of the ventricular region. This section introduces a lesion assessment threshold, which is derived from a database established by our laboratory containing 500 healthy controls and 300 confirmed AD patients. The specific threshold is set as follows: Hippocampal atrophy detection threshold. The corresponding annualized atrophy rate is approximately 4%; the threshold for detecting ventricular dilation. This corresponds to an annualized expansion rate of approximately 5%.

[0064] In a practical example, using the data from S3, the system calculates the voxels of all voxels in the subject's left hippocampus region. The mean value was -0.065, in the bilateral ventricle region. The mean value is 0.072. The system logically compares these calculation results with a preset threshold: first, it determines the hippocampal region, because... If the conditions are met, it is determined to be "significant atrophy"; secondly, the ventricular region is examined because... If the condition is met, it is judged as "significant dilation". Combining the above two judgments, the system automatically generates an abnormal state identifier and outputs a structured diagnostic report: "The subject's hippocampus exhibits pathological atrophy characteristics, accompanied by compensatory ventricular dilation, consistent with the imaging manifestations of neurodegenerative diseases." This process is fully automated, eliminating the subjective errors caused by manual delineation of the ROI. Experimental verification shows that... Under the specified threshold, the method achieved a sensitivity of 89% and a specificity of 92% in the diagnosis of early Alzheimer's disease, verifying the effectiveness of the technology in clinical auxiliary diagnosis.

[0065] A multi-time-point neural image comparison and analysis system is provided. This system is used to execute the aforementioned multi-time-point neural image comparison and analysis method. The system includes: The image acquisition and rigid registration module is used to acquire a first three-dimensional neural image matrix and a second three-dimensional neural image matrix using a magnetic resonance imaging instrument. It calculates the rigid transformation matrix of the second three-dimensional neural image matrix relative to the first three-dimensional neural image matrix through geometric center alignment and iterative optimization algorithms, and performs spatial resampling and registration operations on the second three-dimensional neural image matrix based on the rigid transformation matrix to generate a third three-dimensional neural image matrix.

[0066] The elastic mesh and vector field generation module is used to construct a B-spline-based three-dimensional elastic mesh model using the first and third three-dimensional neural image matrices as input data. Based on the maximum mutual information criterion, the node coordinates of the three-dimensional elastic mesh model are iteratively adjusted to match the local features of the image. The displacement vector of the third three-dimensional neural image matrix relative to the reference image is calculated by interpolation to generate a three-dimensional displacement vector field.

[0067] The micro-volume change analysis module is used to calculate the axial partial derivatives of the three-dimensional displacement vector field using the central difference algorithm. Based on the axial partial derivatives, a Jacobian matrix describing the local deformation gradient is constructed. The determinant of the Jacobian matrix is ​​calculated and subjected to a natural logarithmic transformation to generate a micro-volume change rate map.

[0068] The region status assessment and determination module is used to map the standard brain anatomical structure template to the first three-dimensional neural image matrix space through nonlinear transformation to delineate the region of interest. It uses masking technology to extract the logarithmic Jacobian determinant value of the corresponding region in the microscopic volume change rate map, calculates its arithmetic mean and compares it with the preset lesion threshold, and outputs the volume change status of the region of interest.

[0069] The above embodiments illustrate preferred embodiments of the present invention. Any equivalent adjustments to the technical solution based on software engineering methods are within the scope of protection, including but not limited to: implementing algorithm logic using different programming languages, refactoring functional modules into services, adjusting data interaction protocols, and optimizing resource scheduling strategies. Any implementation scheme derived from reasonable modifications to the data processing flow, service call chain, or system architecture layer without departing from the core technology of the present invention should be considered within the protection scope defined by the technical solution of the present invention.

Claims

1. A method for comparative analysis of neural images at multiple time points, characterized in that, Includes the following steps: S1: Use a magnetic resonance imaging (MRI) instrument to acquire a first three-dimensional neural image matrix and a second three-dimensional neural image matrix, calculate the rigid transformation matrix of the second three-dimensional neural image matrix relative to the first three-dimensional neural image matrix, register the second three-dimensional neural image matrix based on the rigid transformation matrix, and generate a third three-dimensional neural image matrix. S2: Construct a three-dimensional elastic mesh model using the first three-dimensional neural image matrix and the third three-dimensional neural image matrix, adjust the node coordinates of the three-dimensional elastic mesh model based on the mutual information criterion, calculate the displacement vector of the third three-dimensional neural image matrix, and generate a three-dimensional displacement vector field. S3: Calculate the axial partial derivative of the three-dimensional displacement vector field, construct the Jacobian matrix based on the axial partial derivative, calculate the determinant of the Jacobian matrix, perform a natural logarithmic transformation on the determinant, and generate a microscopic volume change rate map. S4: Map the standard brain anatomical structure template onto the first three-dimensional neural image matrix to delineate the region of interest, extract the logarithmic Jacobian determinant value of the region of interest in the microscopic volume change rate map, calculate the arithmetic mean of the logarithmic Jacobian determinant value, compare the arithmetic mean with the lesion threshold, and determine the volume change state of the region of interest.

2. The multi-timepoint neural image comparison and analysis method according to claim 1, characterized in that, The specific steps of S1 are as follows: S11: At the baseline time point, the original image data of the subject's brain is acquired using a magnetic resonance imaging device, and the original image data is processed by skull dissection and grayscale non-uniformity correction to generate a first three-dimensional neural image matrix. At the follow-up time point, brain image data is acquired using the same imaging parameters and the same preprocessing is performed to generate a second three-dimensional neural image matrix. S12: Extract the geometric center coordinates of the first three-dimensional neural image matrix and the second three-dimensional neural image matrix, calculate the translation vector and rotation matrix that can make the two geometric centers coincide and the principal axis directions aligned based on the least squares method, and combine the translation vector and rotation matrix to construct a rigid transformation matrix. S13: Based on the spatial mapping relationship defined by the rigid transformation matrix, perform inverse coordinate transformation on each voxel coordinate in the second three-dimensional neural image matrix, and use the trilinear interpolation algorithm to calculate the gray value of the transformed coordinate position to generate the third three-dimensional neural image matrix.

3. The multi-time-point neural image comparison and analysis method according to claim 2, characterized in that, The specific steps of S2 are as follows: S21: Define a set of uniformly distributed grid control points in the image space of the first three-dimensional neural image matrix, construct a three-dimensional elastic grid model using B-spline basis functions, and map the third three-dimensional neural image matrix as a floating image to the coordinate system defined by the model; S22: Calculate the joint gray-level histogram between the first three-dimensional neural image matrix and the third three-dimensional neural image matrix affected by mesh deformation, calculate the similarity measure between the two images according to the mutual information criterion, and iteratively adjust the spatial coordinates of multiple control points in the three-dimensional elastic mesh model using the gradient descent optimization algorithm until the similarity measure reaches the convergence condition. S23: Obtain the offset of the control point coordinates after optimization and convergence relative to the initial grid position, and calculate the non-rigid deformation component of each voxel point in the image space using the B-spline tensor product formula to generate a three-dimensional displacement vector field.

4. The multi-timepoint neural image comparison and analysis method according to claim 3, characterized in that, The specific steps of S3 are as follows: S31: Traverse each voxel node in the three-dimensional displacement vector field, and use the central difference method to calculate the rate of change of the displacement vector at the node in the X-axis, Y-axis and Z-axis directions of the spatial rectangular coordinate system, and obtain nine axial partial derivatives; S32: Construct a Jacobian matrix using the nine axial partial derivatives, calculate the determinant value of the Jacobian matrix to quantify the local volume expansion or contraction factor of the voxel point between two time points, and obtain the determinant value. S33: Perform a logarithmic operation with the natural constant as the base on the determinant value, map the volume shrinkage change to a negative value and the volume expansion change to a positive value, and construct a microscopic volume change rate map based on the operation result.

5. The multi-timepoint neural image comparison and analysis method according to claim 4, characterized in that, The specific steps of S4 are as follows: S41: Obtain a standard brain anatomical structure template including preset anatomical partition labels, calculate the affine transformation parameters and nonlinear deformation parameters of the standard brain anatomical structure template relative to the first three-dimensional neural image matrix, and use the parameters to inversely deform the standard brain anatomical structure template and resample it to the coordinate space where the first three-dimensional neural image matrix is ​​located. S42: Generate a binary mask matrix based on the anatomical partition labels in the mapped standard brain anatomy template, perform a Boolean AND operation between the binary mask matrix and the microscopic volume change rate map, delineate the region of interest, and extract the logarithmic Jacobian determinant values ​​of all effective voxels within the region; S43: Calculate the logarithmic Jacobian determinant value of all voxels in the region of interest and its arithmetic mean. Compare the arithmetic mean with a preset lesion threshold. If the mean exceeds the threshold range, generate an abnormal state identifier to determine the volume change status of the region of interest.

6. The multi-timepoint neural image comparison and analysis method according to claim 5, characterized in that, The process of constructing the rigid transformation matrix includes: Obtain the voxel grayscale sets of the first three-dimensional neural image matrix and the second three-dimensional neural image matrix, construct a rigid registration objective function based on normalized mutual information, and use a multi-resolution pyramid strategy to perform hierarchical downsampling processing on the image; At each resolution level, the Powell optimization algorithm is used to search for the optimal six-degree-of-freedom rigid parameters, which include translation along the three coordinate axes and rotation angles around the three coordinate axes. The convergent parameter combination at the highest resolution level is converted into a four-by-four homogeneous transformation matrix to establish the rigid transformation matrix.

7. The multi-timepoint neural image comparison and analysis method according to claim 6, characterized in that, The node coordinate adjustment process of the three-dimensional elastic mesh model includes: Obtain the probability distribution of the reference image in the current iteration step. Floating image probability distribution and the joint probability distribution of the two under the current mesh deformation. According to the following formula: ； Calculate the mutual information objective function as a registration similarity measure; in, Representative reference image With floating images The mutual information values ​​between them Represents the gray level of the reference image. Represents the gray levels of the floating image. Represents the numerical value of the joint probability distribution. Represents the edge probability distribution of the reference image. Represents the edge probability distribution of the floating image; The negative gradient direction of the mutual information objective function relative to the coordinates of the grid control points is calculated. The coordinates of the control points are updated along the negative gradient direction with an adaptive step size. The histogram and gradient are repeatedly calculated until the rate of change of the objective function is lower than the preset convergence threshold, thereby generating a non-rigid deformation grid.

8. The multi-timepoint neural image comparison and analysis method according to claim 7, characterized in that, The numerical calculation process of the determinant of the Jacobian matrix includes: Obtain the three-dimensional displacement vector field in coordinates displacement components at the location , and Construct the local deformation gradient tensor matrix as shown below: ； in, The value representing the determinant of the Jacobian matrix is... Displacement vector Quantity in Axial partial derivative of the direction, Displacement vector Quantity in Axial partial derivative of the direction, Displacement vector Quantity in The axial partial derivatives of the direction, and the remaining elements in the matrix represent the cross partial derivatives; The determinant value of the matrix is ​​calculated to obtain the local volume change factor at the coordinate point. When the value is greater than 1, it represents local volume expansion. When the value is less than 1 but greater than 0, it represents local volume contraction. The determinant value is generated.

9. The multi-timepoint neural image comparison and analysis method according to claim 8, characterized in that, The process of determining the volume change state of the region of interest includes: Obtain the target atrophy detection threshold for the hippocampal region and the target expansion detection threshold for the ventricle region, and calculate the average of the logarithmic Jacobian determinant values ​​of all voxels in the hippocampal region and the average of the logarithmic Jacobian determinant values ​​of all voxels in the ventricle region. Determine if the average number of hippocampal regions is less than the atrophy detection threshold; if so, generate diagnostic prompts. Determine whether the average number of ventricular regions is greater than the expansion detection threshold. If so, generate diagnostic prompts and generate a report on the volume change status of the region of interest based on the above judgment results.

10. A multi-time-point neural image comparison and analysis system, characterized in that, The system is used to implement the multi-time-point neural image comparison and analysis method according to any one of claims 1-9, the system comprising: The image acquisition and rigid registration module is used to acquire a first three-dimensional neural image matrix and a second three-dimensional neural image matrix using a magnetic resonance imaging instrument, calculate the rigid transformation matrix of the second three-dimensional neural image matrix relative to the first three-dimensional neural image matrix through geometric center alignment and iterative optimization algorithms, and perform spatial resampling and registration operations on the second three-dimensional neural image matrix based on the rigid transformation matrix to generate a third three-dimensional neural image matrix. The elastic mesh and vector field generation module is used to construct a three-dimensional elastic mesh model based on B-splines using the first three-dimensional neural image matrix and the third three-dimensional neural image matrix as input data, iteratively adjust the node coordinates of the three-dimensional elastic mesh model to match local features of the image based on the maximum mutual information criterion, and calculate the displacement vector of the third three-dimensional neural image matrix relative to the reference image by interpolation to generate a three-dimensional displacement vector field. The micro-volume change analysis module is used to calculate the axial partial derivative of the three-dimensional displacement vector field using the central difference algorithm, construct a Jacobian matrix describing the local deformation gradient based on the axial partial derivative, calculate the determinant of the Jacobian matrix and perform a natural logarithmic transformation on it, thereby generating a micro-volume change rate map. The region status assessment and determination module is used to map the standard brain anatomical structure template to the first three-dimensional neural image matrix space through nonlinear transformation to delineate the region of interest. It uses masking technology to extract the logarithmic Jacobian determinant value of the corresponding region in the microscopic volume change rate map, calculates its arithmetic mean and compares it with the preset lesion threshold, and outputs the volume change status of the region of interest.

Images