Bone quality assessment method based on neural network enhancement and finite element analysis

Abstract

This invention provides a bone quality assessment method based on neural network enhancement and finite element analysis, comprising: acquiring ultra-high resolution CT images of living vertebrae; performing domain transformation processing on the CT images using a generative neural network; resampling and segmenting the enhanced CT images; constructing a trabecular microstructure model within the cancellous bone mask region; assembling RVE elements and a cortical bone model into a complete vertebral body finite element model; applying multi-condition loads to the RVE elements to calculate anisotropic material parameters; and performing axial compression finite element analysis on the complete vertebral body model to determine bone health status. This invention comprehensively considers information from both bone mass and bone quality dimensions, particularly the contribution of trabecular microstructure to bone strength, enabling earlier identification of patients with bone degeneration and high fracture risk, providing a scientific basis for early clinical intervention.

Description

Technical Field

[0001] This invention relates to the fields of medical image analysis and artificial intelligence technology, specifically to a bone quality assessment method based on neural network enhancement and finite element analysis. Background Technology

[0002] The diagnosis and risk assessment of osteoporosis and bone metabolic diseases require consideration of both bone mass and bone quality. Bone mass is typically quantified using bone mineral density (BMD), while bone quality is primarily reflected in the microstructural characteristics of trabecular bone. Numerous clinical studies have shown that trabecular microstructural parameters are core factors determining bone strength, and their degeneration often precedes the decline in BMD, serving as a key indicator of early disease progression. Therefore, accurately acquiring and analyzing trabecular microstructure is of significant clinical importance for the early diagnosis of osteoporosis and the prediction of fracture risk.

[0003] Currently, the most commonly used bone quality assessment technique in clinical practice is dual-energy X-ray absorptiometry (DEXA). This method calculates bone mineral density by measuring the differences in bone tissue absorption of X-rays of different energies, offering advantages such as ease of operation and low radiation dose. However, DEXA only provides two-dimensional projection bone mineral density information and cannot directly observe and quantify the three-dimensional microstructure of trabecular bone. Its measurement results are easily affected by factors such as equipment resolution limitations, image noise, changes in patient body fat percentage, and artifacts from implanted metal devices, resulting in unreliable assessment results for some samples. More importantly, DEXA has low sensitivity to trabecular bone microstructural degradation, often only detecting abnormalities after a significant decrease in bone mineral density has occurred, making it difficult to meet the needs of early diagnosis.

[0004] To directly observe the microstructure of bone trabeculae, researchers have developed micro-computed tomography (MCT). This technique enables three-dimensional imaging of bone tissue at micrometer-level resolution, clearly displaying microstructural parameters such as the spatial arrangement, connectivity, and thickness of bone trabeculae. MCT has achieved significant results in the study of ex vivo bone samples. However, applying MCT to in vivo detection faces considerable challenges. To obtain sufficiently high resolution to display bone trabeculae details, MCT requires high radiation doses and long scan times, which raises safety and practicality issues in in vivo applications. Existing in vivo MCT research primarily focuses on small animal experiments; for human bone tissue, current in vivo imaging techniques struggle to obtain sufficiently clear images of bone trabeculae within safe radiation dose ranges.

[0005] In recent years, the development of ultra-high resolution computed tomography (UHD) technology has brought new possibilities to in vivo trabecular bone imaging. However, even with the most advanced UHD equipment, the image quality obtained from in vivo scanning is still far lower than that from ex vivo high-resolution scanning. In vivo scanning is affected by physiological movements such as the patient's breathing and heartbeat, resulting in motion blur in the images. To control radiation dose within a safe range, scanning parameters need to be compromised, which further reduces the signal-to-noise ratio of the images. These factors collectively lead to blurred boundaries and loss of detail in trabecular bone structures in in vivo computed tomography images, making them difficult to use directly for accurate microstructural parameter extraction and biomechanical analysis.

[0006] In the mechanical analysis of trabecular bone microstructure, the finite element method (FEM) is the most commonly used computational tool. Traditional FEM analysis typically builds models based on simplified bone density distributions, treating bone tissue as a continuous homogeneous or heterogeneous material, which fails to accurately reflect the discrete network structure of trabecular bone. Although some studies have attempted to construct finite element models containing trabecular bone details based on high-resolution images, these studies primarily use micro-computed tomography data from ex vivo samples, making them unsuitable for in vivo patients. How to construct high-fidelity finite element models of trabecular bone under conditions of limited in vivo image quality remains an unsolved technical challenge.

[0007] In summary, existing bone assessment techniques have significant limitations. Dual-energy X-ray absorptiometry (DXA) cannot acquire information on the microstructure of bone trabeculae and is insensitive to early lesions. While microcomputed tomography (CT) can provide high-resolution imaging, it requires invasive sampling or administers radiation doses exceeding safe limits, limiting its widespread application in in vivo clinical testing. Although in vivo ultra-high resolution CT methods are non-invasive, the image quality is insufficient to support accurate trabecular structure analysis and biomechanical assessment. Therefore, there is an urgent need to develop a new technical solution that can accurately extract trabecular microstructure information from in vivo CT images safely and non-invasively, and perform reliable biomechanical analysis based on the actual microstructure, thereby achieving early and accurate assessment of osteoporosis and fracture risk. Summary of the Invention

[0008] To overcome the shortcomings of existing technologies, this invention proposes a bone quality assessment method based on neural network enhancement and finite element analysis, comprising the following steps: S1: Acquire ultra-high resolution CT images of living vertebrae; S2: A generative neural network is used to perform domain transformation processing on the CT images, converting in vivo CT images with unclear trabecular bone structure into enhanced CT images with clear trabecular bone structure. The generative neural network achieves domain transformation by learning the features of isolated high-resolution trabecular bone images. S3: Resample and segment the enhanced CT images to obtain the cortical bone mask and the cancellous bone mask; S4: Within the cancellous bone mask region, a microstructure model of bone trabeculae is constructed using the representative volume element (RVE) method. Multiple candidate units are extracted by traversing the cancellous bone region, and effective RVE units are selected based on the statistical characteristics and connectivity of bone volume fraction. S5: Assemble the RVE unit with the cortical bone model into a complete vertebral finite element model, and apply multi-condition loads to the RVE unit to calculate the anisotropic material parameters; S6: Perform axial compression finite element analysis on the complete vertebral model, and determine the bone health status based on the relationship between the stress borne by the trabecular structure and the preset failure stress threshold.

[0009] Furthermore, in step S2, the generative neural network adopts a recurrent generative adversarial network architecture, including: The two generators respectively implement bidirectional conversion from the in vivo CT domain to the high-resolution domain and from the high-resolution domain to the in vivo CT domain, forming a cyclic consistency constraint; Two discriminators determine the quality of the generated images in the two domains, respectively. The generator adopts an encoder-decoder structure. The encoder extracts features through multiple downsampling, and the decoder restores the resolution through multiple upsampling. Skip connections are established between the downsampled features and the upsampled features for feature fusion.

[0010] Furthermore, the training of the generative neural network employs a comprehensive loss function: ; in: To combat the losses; For structural similarity loss; For three-dimensional gradient loss; , , These are the weighting coefficients.

[0011] Furthermore, the structural similarity loss is calculated using a simplified method: For the generated image X and the target image Y, apply a Gaussian kernel. Local statistical properties of convolution calculations include: Local mean: ; Mean product: ; variance: ; Covariance: ; The structural similarity loss is calculated as follows: ;

[0012] in: X: Generated image, a 3D matrix with dimensions H×W×D, where H is the height, W is the width, and D is the depth; Y: Target image, same size as X; Gaussian kernel, used for local weighted averaging, κ represents the standard deviation parameter of the kernel, typically 1.5; : Convolution operator, indicating convolution of the image with a Gaussian kernel; , : Local mean map of image X and Y after Gaussian kernel convolution, with the same size as the original image; : Pixel-wise product of two local mean maps; , The local variance map of image X and Y reflects the degree of dispersion of local pixel values; The local covariance plot of two images reflects the correlation of local pixel values ​​between the two images; The averaging operation calculates the average value across all pixel locations in the entire image. C: Stability constant, to prevent the denominator from being zero, typically ranging from 0.01 to 0.1.

[0013] Furthermore, the three-dimensional gradient loss is calculated based on the gradient differences of the image in three spatial directions: Calculate the gradients of the generated image and the target image in the x, y, and z directions respectively: ;

[0014] The gradient loss is expressed in the form of squared error: ;

[0015] in: : Slicing operation, which means to crop the image from the second element to the last element in the x x direction; the index starts from 0, so [1:][1:][1:] means from the second element to the last element; : This indicates the sequence from the first element to the second-to-last element, i.e., removing the last element; The gradient map of the image in the x-direction (height direction) is calculated by the difference between adjacent pixels and has a size of (H−1)×W×D(H-1). : The gradient map of the image in the y-direction (width direction), with dimensions H×(W−1)×D; : The gradient map of the image in the z-direction (depth direction), with dimensions H×W×(D−1); : The gradient of the generated image Xgen in the x-direction, where the superscript gen indicates the generated image; : The gradient of the target image Xtar in the x-direction, where the superscript tar indicates the target image; : The squared error per pixel of the two gradient maps; : 3D gradient loss value, the sum of gradient losses in three directions. The smaller the value, the closer the gradients are.

[0016] Furthermore, to eliminate edge artifacts caused by image block stitching after domain transformation, a three-dimensional weighted fusion method based on window functions is adopted: Construct a one-dimensional window function based on the patch size:

[0017] For a voxel position (px, py, pz) in three-dimensional space, its three-dimensional weight is calculated by normalizing the product of the weights in each direction:

[0018] in: The one-dimensional Hann window function is a smooth cosine window function that smoothly transitions to 0 at the boundaries. n: The index position of the window function, an integer, with a value range from 0 to N−1; N: Width of the overlapping area of ​​the patch boundary, in pixels, typically 10%-20% of the patch side length; : Voxel coordinates in three-dimensional space, representing the position indices in the x, y, and z directions respectively; : The window function weight value of the voxel at position px in the x-direction, with a value range of 0-1; : The window function weight value of the voxel at position py in the y direction, with a value range of 0-1; : The window function weight value of the voxel at position pz in the z-direction, with a value range of 0-1; : Pixel-wise product of the weights in the three directions; The maximum value among all voxel position weight products, used for normalization; The pixel values ​​of the overlapping regions are then weighted and averaged using this weight: , where k is the index of the overlapping patch.

[0019] Furthermore, in step S3: Large-size CT images are processed using a block resampling method. The data is divided into multiple blocks along the z-axis, and each block is independently interpolated in three dimensions. Overlapping regions are set between blocks and then fused. The resampled image is automatically segmented using a deep learning segmentation network to obtain cortical bone masks and cancellous bone masks.

[0020] Furthermore, the selection of RVE units in step S4 adopts the square root transformation criterion: Candidate cells were extracted by traversing the cancellous bone region using a three-dimensional sliding window, and the bone volume fraction (BV / TV) and connectivity index were calculated for each candidate cell. Calculate the variance and global mean of the bone volume fraction within the candidate unit, and filter out valid units whose variance is less than a preset threshold multiple of the mean. Select the optimal RVE cell from the available cells:

[0021] in : Optimal representative volume element, the most representative element selected from all valid candidate elements; i: Index number of the candidate unit, i=1,2,…; : The normalized connectivity index of the i-th candidate unit, with a value ranging from 0 to 1, reflects the connectivity of the trabeculae within the unit; : Connectivity deviation term, the difference between the connectivity of the measured unit and the ideal value of 1, the smaller the value, the better the connectivity; Bone volume refers to the volume occupied by bone tissue within a unit; Total volume refers to the total volume of the unit (including bone tissue and pores). : Bone volume fraction of the i-th candidate unit; : The arithmetic mean of the bone volume fraction of all candidate units, used to represent the average density of the overall cancellous bone; The deviation term of bone volume fraction after square root transformation; the smaller the value, the closer the unit is to the average state.

[0022] Furthermore, in step S5: Six basic load cases were applied to each RVE element, including three axial normal strain cases and three planar shear strain cases. The stress response is obtained through finite element analysis, and the anisotropic material parameters are extracted based on the stress-strain relationship, including the elastic modulus in three directions, the shear modulus in three planes, and Poisson's ratio.

[0023] Furthermore, in step S6: Set the lower end plate to a completely fixed boundary condition, and apply an axial compressive displacement load to the upper end plate. The analysis stops when the compression displacement reaches a preset proportion of the vertebral body height, and the maximum equivalent stress of the trabecular structure is extracted. The bone condition is determined by comparing the maximum stress with the preset failure stress threshold, and the safety factor is calculated:

[0024] in: Failure stress threshold The maximum equivalent stress is SF; SF>1 indicates safety, and SF≤1 indicates a risk of fracture.

[0025] Compared with the prior art, the beneficial effects of the present invention are: 1. This invention provides a bone quality assessment method based on neural network enhancement and finite element analysis, solving the problem that even high-resolution in vivo computed tomography images cannot clearly display fine structures. Through domain transformation technology using generative neural networks, blurred in vivo images are converted into enhanced images with clear trabecular bone structures, enabling image data acquired at safe radiation doses to be used for precise microstructural analysis.

[0026] 2. This invention provides a bone quality assessment method based on neural network enhancement and finite element analysis. The block resampling method effectively overcomes the drawbacks of traditional whole resampling methods, such as excessive memory consumption and long computation time when processing large-size 3D images. By processing images in blocks along the depth direction and setting overlapping regions for fusion, the computational resource requirements are significantly reduced while ensuring image quality, shortening the processing time and helping to improve the efficiency of medical diagnosis, making this technology more clinically practical.

[0027] 3. This invention provides a bone quality assessment method based on neural network enhancement and finite element analysis. Using the representative volume element method, representative trabecular microstructures are extracted from enhanced live images, constructing a finite element model containing the real trabecular network topology. This model accurately reflects the spatial arrangement, connectivity, and anisotropic mechanical properties of the trabeculae, meeting the urgent clinical need for precise assessment of bone status and prediction of fracture risk.

[0028] 4. This invention provides a bone quality assessment method based on neural network enhancement and finite element analysis. Doctors only need to input computed tomography (CT) image data to obtain bone health assessment results. The operation is simple and quick, significantly improving the efficiency of artificial intelligence-assisted diagnosis and treatment, lowering the technical threshold for use, and facilitating the promotion and application of this method in clinical practice. It accurately reflects the mechanical response characteristics of cancellous bone under different directions and loading modes. This microstructure-based material parameter extraction method is more accurate than traditional empirical formulas, and can capture the changes in mechanical properties caused by differences in trabecular bone structure between individuals, achieving personalized bone mechanical assessment. Attached Figure Description

[0029] To more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0030] Figure 1 This is a schematic diagram of the process of this invention; Figure 2 These are CT image comparison images; Figure 3 This is a diagram illustrating the comparison between the original and target images in the dataset. Figure 4 This is a diagram of a migration network; Figure 5 This is a schematic diagram of the loss curve; Figure 6 These are the original image from in vivo ultra-high resolution CT and the result image obtained from the network in this invention; Figure 7 The results obtained by the edge-clearing method in this invention (right) and the results of the original method (left) are shown. Figure 8 The finite element analysis results for the overall modeling of this invention (a total of four analysis steps). Detailed Implementation

[0031] The technical solution of the present invention will be more clearly and completely explained below with reference to the accompanying drawings and through the description of preferred embodiments of the present invention.

[0032] This invention provides a bone analysis method, such as... Figure 1 As shown, the present invention will be described in detail below with reference to specific embodiments.

[0033] Three-dimensional CT images of living vertebrae were acquired using an ultra-high resolution CT scanner. During the scan, the patient was in a supine position, and the scan area covered the target vertebra and its adjacent vertebrae above and below it. The slice thickness and pixel size were set to 0.1 mm to 0.2 mm to ensure the capture of the microscopic structural features of the trabeculae. After the scan, the CT image data was exported in DICOM format and converted to a three-dimensional array format for subsequent processing.

[0034] Due to radiation dose limitations and motion artifacts in in vivo CT scans, the trabecular bone structure in directly acquired CT images is often insufficiently clear, making it difficult to use directly for microscopic structural analysis. Therefore, a generative neural network is needed to perform domain transformation processing on CT images. This neural network employs a recurrent generative adversarial network architecture, containing two generators and two discriminators. The first generator is responsible for converting the in vivo CT domain image into a high-resolution domain image, while the second generator performs the reverse transformation. This bidirectional transformation mechanism forms a cyclic consistency constraint, ensuring that the transformed image retains the anatomical information of the original image. The two discriminators are used to judge the authenticity of the generated images in the two domains, and the quality of the generated images is continuously improved through an adversarial training mechanism.

[0035] The generator employs an encoder-decoder architecture. The encoder consists of multiple convolutional and downsampling layers, progressively extracting multi-scale features from the image and reducing its spatial resolution. The decoder, on the other hand, progressively restores the image's spatial resolution through deconvolutional and upsampling layers. To better preserve detail, skip connections are established between the encoder's downsampling features and the decoder's upsampling features, directly transferring high-resolution features from shallower layers to deeper layers, thus achieving effective fusion of multi-scale features.

[0036] The training of the neural network employs a comprehensive loss function, which consists of three parts. The first part is the adversarial loss, used to measure the distributional difference between the generated image and the real high-resolution image, guiding the generator optimization through feedback signals from the discriminator. The second part is the structural similarity loss, used to maintain the similarity between the generated and target images in brightness, contrast, and structure. In calculating the structural similarity loss, a Gaussian kernel is first used to convolve the generated and target images, calculating the local mean at each location. Then, the pixel-wise product of the two local mean maps is calculated, along with their respective local variances and their local covariance. The value of the structural similarity loss is derived from these statistics; a smaller value indicates greater structural similarity between the two images. The third part is the three-dimensional gradient loss, used to constrain the edge and texture features of the generated image to remain consistent with the target image in three spatial directions. Specifically, the gradients of the generated and target images in the height, width, and depth directions are extracted respectively. Gradient maps are obtained by calculating the differences between adjacent pixels. Then, the squared error between the gradients of the generated and target images is calculated, and the sum of the errors in the three directions is the three-dimensional gradient loss. These three loss components are weighted and combined using weight coefficients to form the final comprehensive loss function that guides network training.

[0037] In the training data preparation phase, paired in vivo low-resolution CT images and ex vivo high-resolution trabecular bone images are collected. The ex vivo high-resolution images are acquired using a miniature CT scanner with scanning parameters set to a pixel size of ten to twenty micrometers, clearly displaying the three-dimensional network structure of the trabecular bone. These paired data are then input into a recurrent generative adversarial network (GAN) for training. During training, the generator continuously learns how to convert blurry in vivo images into clear, high-resolution images while maintaining the accuracy of anatomical structures. After thousands of iterations of training, the network effectively enhances the details of the trabecular bone structure in in vivo CT images.

[0038] When applying a trained generative neural network to new live CT images, a block-based processing strategy is required due to the large amount of data in the entire vertebral body. The 3D CT image is divided into multiple 3D blocks, each with a size of 128 to 256 voxels. Overlapping regions are established between blocks, with an overlap width typically ranging from 10% to 20% of the block's side length. After performing domain transformation on each block independently, the processed blocks need to be stitched back into a complete 3D image. To eliminate potential artifacts at the stitching boundaries, a window-based weighted fusion method is employed. First, a one-dimensional window function is constructed based on the block size and overlap width. This window function uses a Hanning window form, smoothly transitioning to zero at the boundaries and maintaining a higher weight in the central region. For each voxel location in 3D space, its fusion weight is obtained by multiplying the one-dimensional window function values ​​in three directions and then normalizing the result. In the overlapping region, pixel values ​​from different blocks are weighted and averaged according to their corresponding weights to achieve a smooth transition and seamless fusion.

[0039] After obtaining enhanced CT images with clear trabecular bone structure, the next step is image resampling and segmentation. Since the voxel sizes of the original CT images may be inconsistent in different directions, 3D interpolation is used to resample the images into an isotropic voxel grid. The resampled voxel size is uniformly set to 0.1 mm. For large CT images, a block resampling method is used to save memory and improve computational efficiency. The data is divided into multiple blocks along the z-axis, each containing approximately 200 to 300 slices, with 20 to 30 overlapping slices between blocks. 3D interpolation is performed independently on each block, and then weighted fusion is performed in the overlapping areas to obtain the complete resampled image.

[0040] After resampling, a deep learning segmentation network was used to automatically segment the image, dividing the vertebral body into two regions: cortical bone and cancellous bone. The segmentation network adopted a 3D U-Net architecture, which was trained on a large amount of labeled data and could accurately identify the boundaries between cortical and cancellous bone. The segmentation results were output in the form of binary masks. In the cortical bone mask, the voxel value of the cortical bone region was one, and the voxel value of the other regions was zero; in the cancellous bone mask, the voxel value of the cancellous bone region was one, and the voxel value of the other regions was zero. To improve segmentation accuracy, morphological post-processing was performed after segmentation, including hole filling, small region removal, and boundary smoothing.

[0041] Within the cancellous bone mask region, a microstructural model of the trabecular bone was constructed using the representative volume element method. First, a three-dimensional sliding window was used to traverse the cancellous bone region. The window size was set according to the typical trabecular spacing of cancellous bone, typically a cube of five to ten millimeters. The sliding step size was set to one-quarter to one-half of the window size to ensure full coverage of the entire cancellous bone region. For each window location, a three-dimensional image patch at that location was extracted as a candidate cell, and the bone volume fraction and connectivity index of that cell were calculated. The bone volume fraction was defined as the ratio of bone tissue volume within the cell to the total volume, obtained by statistically analyzing the ratio of the number of bone voxels within the cell to the total number of voxels. The connectivity index was calculated through three-dimensional connected domain analysis. First, all independent trabecular connected regions within the cell were identified, and the proportion of the volume of the largest connected region to the total bone volume was calculated. This proportion, after normalization, was used as the connectivity index.

[0042] To select truly representative volumetric units, statistical analysis is required for all candidate units. The global mean and variance of bone volume fraction for all candidate units are calculated, and the homogeneity of each candidate unit is evaluated using a square root transformation criterion. Specifically, the variance of bone volume fraction in sub-regions within a candidate unit is calculated and compared to the square root of the global mean. If the variance is less than a preset multiple of the square root of the mean, the unit is considered to have a relatively uniform bone density distribution and can be considered a valid unit. The preset multiple is typically set to 0.5 to 1.0. Among the selected valid units, the optimal representative volumetric unit is further selected. The selection criteria consider two factors: the first is the deviation of the connectivity index from the ideal value of one; a smaller deviation indicates a more complete trabecular bone network. The second factor is the deviation of the unit's bone volume fraction from the global average bone volume fraction; the deviation after square root transformation better reflects differences in bone density. These two factors are added together to obtain a comprehensive score, and the unit with the lowest score is selected as the optimal representative volumetric unit.

[0043] After obtaining representative volumetric elements, they need to be assembled with the cortical bone model to form a complete vertebral finite element model. First, a three-dimensional geometric model of the cortical bone is generated based on the cortical bone mask, and the triangular mesh of the cortical bone surface is extracted using the traveling cube algorithm. Then, the cortical bone model is meshed using finite element methods, generating tetrahedral or hexahedral element meshes. In the cancellous bone region, the selected representative volumetric elements are copied and arrayed to fill the entire cancellous bone space. During the filling process, the scaling ratio of the representative volumetric elements is adjusted according to changes in local bone density to match the bone volume fraction with the actual measured value. Node merging and mesh transition processing are performed at the boundary between the cortical bone mesh and the cancellous bone mesh to ensure good connection between the two meshes.

[0044] Before performing mechanical analysis, it is necessary to determine the anisotropic material parameters of the representative volume element. Six basic load cases are applied to a single representative volume element, including normal strain cases along the three coordinate axes and shear strain cases in the three coordinate planes. Under each load case, displacement constraints are applied to one boundary surface of the element, uniform displacement loads are applied to the opposite boundary surface, and the remaining boundary surfaces are set as free boundaries or periodic boundary conditions. The stress response under each load case is obtained through finite element calculation, and material parameters are extracted based on the stress-strain relationship. Specifically, the elastic modulus in the three principal directions is extracted from the three normal strain cases, the shear modulus in the three planes is extracted from the three shear strain cases, and Poisson's ratio is extracted from the transverse strain response of the normal strain cases. These parameters together constitute the material constant matrix describing the anisotropic mechanical behavior of cancellous bone. Cortical bone is modeled as an isotropic material with an elastic modulus set to 10 to 20 GPa and a Poisson's ratio set to 0.3.

[0045] After assigning material parameters, an axial compression analysis was performed on the complete vertebral finite element model to assess bone health. Completely fixed boundary conditions were applied to the lower endplate of the model, constraining three translational and three rotational degrees of freedom for all nodes on this surface. An axially downward displacement load was applied to the upper endplate, with the displacement set to one to two percent of the vertebral height, simulating compressive deformation under physiological loading conditions. Static analysis was performed using an implicit finite element solver, considering geometric and material nonlinearities during the calculation. After the analysis, the stress distribution of the trabecular structure was extracted, the equivalent stress of each element was calculated, and the stress state was evaluated using the von Mises criterion.

[0046] Bone health is determined based on the failure stress threshold of trabecular bone materials. The failure stress threshold is determined using experimental test data, typically set at 60-80 MPa for normal bone and potentially lowered to 30-50 MPa for osteoporotic patients. The calculated maximum equivalent stress is compared to the failure stress threshold to calculate a safety factor. The safety factor is defined as the ratio of the failure stress threshold to the maximum equivalent stress. A safety factor greater than one indicates that the trabecular structure can safely withstand the load, and the bone condition is good. A safety factor less than or equal to one indicates that local trabecular stress exceeds the failure threshold, posing a risk of fracture and requiring further clinical evaluation and intervention.

[0047] Furthermore, fracture-prone sites can be identified by analyzing spatial patterns of stress distribution. High-stress concentration areas are marked in the vertebral model; these areas often correspond to weak or poorly connected trabecular networks and are key areas of focus for clinical treatment. By comparing the analysis results at different time points, trends in bone changes can be assessed, providing quantitative evidence for adjusting treatment plans.

[0048] The method of this invention enhances the details of the trabecular bone structure in in vivo CT images through generative neural networks, and constructs a high-fidelity finite element model of the vertebral body by combining the representative volume element method. This enables bone biomechanical analysis based on the real trabecular bone microstructure, which can more accurately assess fracture risk and provides an effective technical means for the early diagnosis and personalized treatment of osteoporosis.

[0049] As a preferred embodiment, the present invention focuses on a method for bone analysis of the microstructure of trabecular bone in intact vertebral bodies in living organisms, the main steps of which are as follows: Generative networks accomplish style transfer This invention is based on the analysis of ultra-high resolution CT images (3072x3072x1090), such as... Figure 2 As shown, the left image is an ultra-high resolution image of a living vertebra, and the right image is a MicroCT image of a vertebra. In comparison, although direct modeling of the trabecular structure of a living vertebra remains difficult at high resolution, faint traces still exist. These traces can be transferred to the living CT image using a generative network to enhance the display of the trabecular structure.

[0050] Dataset Configuration: The dataset for this method comes from the Department of Radiology, First Affiliated Hospital of Soochow University, and includes 1600 ultra-high resolution image slices of living vertebrae and 2151 MicroCT image slices of vertebrae. The ultra-high resolution images are the original images requiring style transfer, and the MicroCT images are used as label images. Since the original and label images are not directly corresponding, no additional data pairing or registration is required. Figure 3 The image shown is a partial representation of the dataset used in this invention.

[0051] Reference Figure 4 This invention proposes an improved 3D style transfer network based on CycleGAN, which is suitable for inputs such as ultra-high resolution CT images. At that time, the generated graph is first obtained through generator C. The generated graph is then obtained through generator F. When the input is a target domain image (MicroCT image) At that time, the generated graph is first obtained through generator F. The generated graph is obtained through generator C. The first and second generated images are then judged by discriminators T and D, respectively. This invention replaces the generator network with SegresNet, introduces a residual structure and a deep network, downsamples the original image domain three times to obtain low-level features, and then upsamples these low-level features three times back to the original image size. During this process, downsampled and upsampled features are continuously fused to capture multi-scale feature information, ultimately completing the image generation. The generated image is then fed into the discriminator and compared with a threshold to determine whether the generated image quality is acceptable.

[0052] Specific model configuration: For an input image of size... The original image slice data first enters the generator network, referring to... Figure 1 The generator structure in the code, the basic block structure of the generator network is normalization, ReLU activation operation, Kernel-sized convolution operations, normalization operations, ReLU activation operations, Convolution operations with kernel size . In the diagram, the basic blocks marked with "*" will undergo downsampling with a kernel stride of 2, and the basic blocks marked with "=" will undergo upsampling with a kernel stride of 1. The size of the lowest-level feature map is... After three upsampling operations, the image is generated, and the size of the generated image is the same as the original image. The discriminator structure is as follows: Figure 1 As shown, the model input is The structure includes the core size. Convolution operations, ReLU activation operations, kernel size Convolution operations, normalization, ReLU activation operations, The kernel-sized convolution operation is followed by a sigmoid function to convert the features into probability values ​​for output.

[0053] Parameter configuration: The learning rate during model training is 0.0001, with a learning rate decay of 0.1 every 50 epochs, for a total of 400 training epochs. All the dataset samples mentioned above will be used, with the training and test sets split in a 9:1 ratio. The model training process is as follows: Figure 5 As shown, the left side of the figure shows the loss curves of the forward generated image and the original image, while the right side shows the loss curves of the reverse generated image and the original image. When these two curves no longer show a significant decrease, the model is considered to have completed training.

[0054] In the loss function section, this invention incorporates Structural Similarity Loss (SSIM Loss) and Gradient Loss to evaluate the structural similarity and gradient differences of the generated graphs. For SSIM Loss, a Gaussian kernel is assumed to exist. Then generate the graph And the original image mean and Square of the mean and Cross product This can be expressed by the following formula:

[0055] The mean of the squared pixel values ​​of the image and It can be expressed by the following formula:

[0056] Finally, the mean of the product of the two images is calculated. for:

[0057] The final loss function is:

[0058] In the formula, To calculate the mean, These are non-zero constants. For Gradient Loss, the gradient information along the three axes of the image is first calculated. , , :

[0059] Then gradient loss It can be expressed as the following formula:

[0060] In the formula, , , To generate the three-dimensional gradient values ​​of the image, , , To generate the three-dimensional gradient values ​​of the image, the comprehensive loss function is... for:

[0061] In the formula, This is the original loss function. Furthermore, to address image artifacts resulting from segmenting images into 3D patches and then stitching them back together, this invention proposes an edge-sharpening method in the post-processing section. For a set of 3D CT image sequences, firstly, one-dimensional windows in three axial directions are determined based on the patch size, using the following formula:

[0062] In the formula, n represents the nth element (n = 0, 1, ..., N-1). For the magnitude of the weight, Let be the number of elements. Then, at any position within the 3D window... 3D weight size It can be expressed by the following formula:

[0063] In the formula, , and These are the weight values ​​for the indices corresponding to the three directions. This is for retrieving the maximum value.

[0064] Reference Figure 6 The image on the left is the original CT image, while the image on the right is a CT image generated by a model that has learned clear features, showing a distinct trabecular structure. This is helpful for further bone microstructural biomechanical analysis of in vivo data. (See reference) Figure 7 The image shown is the result of the edge-clearing method described in this paper. This method effectively eliminates artifacts at the edges.

[0065] Image resampling and segmentation First, the image is resampled from the original 0.165mm to 1.0mm to ensure consistent voxel spacing in all three directions. This invention first resamples the original CT image data to a slice thickness of 1.0mm, thus determining the target size. It can be expressed by the following formula:

[0066] In the formula, The spatial size of the original image voxels. This represents the original image size. The scaling factor is...

[0067] To avoid excessive memory usage, a block resampling method is used. Using the axis as a reference, first set the maximum allowed memory value to... and the number of blocks to be divided Then the number of slices contained in each block is:

[0068] Next, three-dimensional linear interpolation was used to obtain resampled data for each block of CT image content. Finally, the results from different blocks were stitched together to obtain the final downsampled result.

[0069] The resampled image was segmented, and deep learning methods were used to segment the cortical bone and cancellous bone to obtain a cortical bone mask. and cancellous bone mask You can choose the skellytour method to achieve this step.

[0070] RVE Unit Construction Method To preserve both macroscopic and microscopic characteristics during overall modeling, this invention uses the RVE cell construction method to reconstruct the trabecular bone structure. After obtaining the segmented regions, the masked regions are resampled and then restored to their original voxel side lengths. mm.

[0071] Select If the individual element is used as the side length of the RVE unit, then the length of one side is... This can be represented by the following formula:

[0072] The volume of the element can then be expressed as:

[0073] For cancellous bone region Use a 3D sliding window for traversal, with a step size selected as [value]. Then, a cell at a certain moment can be represented as:

[0074] Calculate the bone volume fraction in this region:

[0075] Calculate the variance and mean of the bone volume fraction within the region. A cell is considered valid if the variance is less than 0.05 times the mean. Finally, select the RVE cell with the highest tangential connectivity that is closest to the average bone volume, satisfying the following formula:

[0076] In the formula It is the best RVE unit in a small area. To obtain the RVE unit operation that minimizes the right side of the equation, This represents the mean bone volume. For the first The connectivity index of each unit. Similarly, perform the same operation on all regions within the cancellous bone to obtain all RVE units, and save them as trabecular inp files.

[0077] Superior anisotropic materials Read the trabecular bone .inp file and set boundary constraints on the three axes based on the basic parameters. Generate six basic working cases based on the coordinate axes and the plane: first, generate them along the three axes respectively. , , Normal strain condition, and three shear strain conditions generated in the three xy, xz, and yz planes respectively. , , Then, the corresponding working condition calculations are performed to obtain the normal stress and shear stress, and finally the elastic modulus, shear modulus and Poisson's ratio are calculated respectively.

[0078] Finite element analysis like Figure 8 As shown, the trabecular bone region and cancellous bone, after modeling, are assembled to obtain a composite model. Then, a point is selected in the middle region of the endplate on the model. As a reference point, the lower endplate of the model is set as a fully fixed boundary condition. Axial compression is applied to the upper endplate, and displacement is applied at a rate of 0.1 mm / min. When the vertebral body produces a displacement of 5% of the vertebral body height, it is determined whether the stress on the trabecular structure is greater than 138 MPa. If it is less than 138 MPa, it indicates good bone quality. If it is greater than 138 MPa, it indicates a certain risk of low bone quality.

[0079] The above-described specific embodiments are merely preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Various modifications, substitutions, and improvements made by those skilled in the art to the technical solutions of the present invention based on the provided textual description and drawings, without departing from the design concept and spirit of the present invention, should all fall within the scope of protection of the present invention. The scope of protection of the present invention is determined by the claims.

Claims

1. A bone quality assessment method based on neural network enhancement and finite element analysis, characterized in that, Includes the following steps: S1: Acquire ultra-high resolution CT images of living vertebrae; S2: A generative neural network is used to perform domain transformation processing on the CT images, converting in vivo CT images with unclear trabecular bone structure into enhanced CT images with clear trabecular bone structure. The generative neural network achieves domain transformation by learning the features of ex vivo high-resolution trabecular bone images. The generative neural network in step S2 adopts a recurrent generative adversarial network architecture, including: The two generators respectively implement bidirectional conversion from the in vivo CT domain to the high-resolution domain and from the high-resolution domain to the in vivo CT domain, forming a cyclic consistency constraint; Two discriminators determine the quality of the generated images in the two domains, respectively. The generator adopts an encoder-decoder structure. The encoder extracts features through multiple downsampling, and the decoder restores the resolution through multiple upsampling. Skip connections are established between the downsampled features and the upsampled features to perform feature fusion. The generative neural network is trained using a comprehensive loss function: ； in: To combat the losses; For structural similarity loss; For three-dimensional gradient loss; , , These are the weighting coefficients; The three-dimensional gradient loss is calculated based on the gradient differences of the image in three spatial directions: Calculate the gradients of the generated image and the target image in the x, y, and z directions respectively: ； The gradient loss is expressed in the form of squared error: ； in: : Slicing operation, which means to cut the image in the x direction from the second element to the last element; : This indicates the sequence from the first element to the second-to-last element, i.e., removing the last element; : Gradient plot of the image in the x-direction; : Gradient map of the image in the y-direction; : Gradient map of the image in the z-direction; : The gradient of the generated image Xgen in the x-direction, where the superscript gen indicates the generated image; : The gradient of the target image Xtar in the x-direction, where the superscript tar indicates the target image; : The squared error per pixel of the two gradient maps; : Three-dimensional gradient loss value, the sum of gradient losses in three directions; To eliminate edge artifacts caused by image block stitching after domain transformation, a three-dimensional weighted fusion method based on window functions is adopted: Construct a one-dimensional window function based on the patch size: ； For a voxel position (px, py, pz) in three-dimensional space, its three-dimensional weight is calculated by normalizing the product of the weights in each direction: ； in: One-dimensional Hann window function; n: The index position of the window function; N: Width of the overlapping area of ​​the patch boundary; Voxel coordinates in three-dimensional space; : The window function weight value of the voxel at position px in the x-direction; : The window function weight of the voxel at position py in the y direction; : The window function weight value of the voxel at position pz in the z-direction; : Pixel-wise product of the weights in the three directions; The maximum value among all voxel position weight products, used for normalization; The pixel values ​​of the overlapping regions are then weighted and averaged using this weight: , where k is the index of the overlapping patch; S3: Resample and segment the enhanced CT images to obtain the cortical bone mask and the cancellous bone mask; S4: Within the cancellous bone mask region, a microstructure model of bone trabeculae is constructed using the representative volume element (RVE) method. Multiple candidate units are extracted by traversing the cancellous bone region, and effective RVE units are selected based on the statistical characteristics and connectivity of bone volume fraction. The selection of the RVE units adopts the square root transformation criterion: Candidate cells were extracted by traversing the cancellous bone region using a three-dimensional sliding window, and the bone volume fraction (BV / TV) and connectivity index were calculated for each candidate cell. Calculate the variance and global mean of the bone volume fraction within the candidate unit, and filter out valid units whose variance is less than a preset threshold multiple of the mean. Select the optimal RVE cell from the available cells: ； in : Optimal representative volume element, the most representative element selected from all valid candidate elements; i: Index number of the candidate unit, i=1,2,…; : The normalized connectivity index of the i-th candidate unit; : Connectivity deviation term; Bone volume; Total volume; : Bone volume fraction of the i-th candidate unit; : The arithmetic mean of the bone volume fraction of all candidate units; : Bone volume fraction deviation term after square root transformation; S5: Assemble the RVE unit with the cortical bone model into a complete vertebral finite element model, and apply multi-condition loads to the RVE unit to calculate the anisotropic material parameters; S6: Perform axial compression finite element analysis on the complete vertebral model, and determine the bone health status based on the relationship between the stress borne by the trabecular structure and the preset failure stress threshold.

2. The bone quality assessment method based on neural network enhancement and finite element analysis according to claim 1, characterized in that, The structural similarity loss is calculated using a simplified method: For the generated image X and the target image Y, apply a Gaussian kernel. Local statistical properties of convolution calculations include: Local mean: ; Mean product: ; variance: ; Covariance: ; The structural similarity loss is calculated as follows: ； in: X: Generated image, a 3D matrix with dimensions H×W×D, where H is the height, W is the width, and D is the depth; Y: Target image, same size as X; Gaussian kernel, used for local weighted averaging, κ represents the standard deviation parameter of the kernel; , : Local mean map of image X and Y after Gaussian kernel convolution; : Pixel-wise product of two local mean maps; , Local variance plots of images X and Y; Local covariance plots of two images; : Calculate the average value; C: Stability constant.

3. The bone quality assessment method based on neural network enhancement and finite element analysis according to claim 1, characterized in that, In step S3: Large-size CT images are processed using a block resampling method. The data is divided into multiple blocks along the z-axis, and each block is independently interpolated in three dimensions. Overlapping regions are set between blocks and then fused. The resampled image is automatically segmented using a deep learning segmentation network to obtain cortical bone masks and cancellous bone masks.

4. The bone quality assessment method based on neural network enhancement and finite element analysis according to claim 1, characterized in that, In step S5: Six basic load cases were applied to each RVE element, including three axial normal strain cases and three planar shear strain cases. The stress response is obtained through finite element analysis, and the anisotropic material parameters are extracted based on the stress-strain relationship, including the elastic modulus in three directions, the shear modulus in three planes, and Poisson's ratio.

5. The bone quality assessment method based on neural network enhancement and finite element analysis according to claim 1, characterized in that, In step S6: Set the lower end plate to a completely fixed boundary condition, and apply an axial compressive displacement load to the upper end plate. The analysis stops when the compression displacement reaches a preset proportion of the vertebral body height, and the maximum equivalent stress of the trabecular structure is extracted. The bone condition is determined by comparing the maximum stress with the preset failure stress threshold, and the safety factor is calculated: ； in: Failure stress threshold The maximum equivalent stress is SF; SF > 1 indicates safety, and SF ≤ 1 indicates a risk of fracture.

Images