A fast window optimization method for calculating spatial orientation of fibrous structures based on a Gaussian mixture model

Abstract

This invention discloses a fast window optimization method for calculating the spatial orientation of fiber structures based on a Gaussian mixture model. The method includes setting a grayscale threshold, binarizing the original fiber image (where pixels larger than the threshold represent fibers with a value of 1, and pixels smaller than the threshold represent background with a value of 0), and calculating the diameter of the fiber image in the binary image. This method achieves high-resolution resolution of fiber orientation at the voxel level. Traditional methods typically only reflect the overall orientation trend of an image, with limited ability to capture details of local structures. This method calculates orientation independently pixel by pixel, achieving fine characterization of individual fibrous tissues. This supports more accurate analysis of fiber orientation, curvature variations, and microscopic morphological heterogeneity, providing a more reliable data foundation for quantitative research on complex fiber systems. While maintaining high accuracy, this method can completely preserve and quantitatively represent detailed information in three-dimensional structures, possessing higher spatial resolution and computational precision.

Description

Technical Field

[0001] This invention belongs to the field of quantitative characterization technology of biological tissues, specifically involving a fast window optimization method for calculating the spatial orientation of fiber structures based on a Gaussian mixture model. Background Technology

[0002] The orientation of fibrous structures in biological tissues within three-dimensional space is one of the core parameters characterizing the microstructure of biological tissues. It not only determines the mechanical properties, material transport efficiency, and signal transduction function of tissues, but also reflects physiological processes and disease evolution through dynamic changes in orientation, serving as a crucial link between the "microstructure" and "macrofunction" of biological tissues. The spatial orientation of biological fibers plays an irreplaceable role in basic research and clinical diagnosis. For example, collagen fibers in tendons are arranged parallel to the direction of force, ensuring the tendon can withstand high-intensity tension; the fibrous skeleton of vascular endothelial cells is arranged along the direction of blood flow, reducing blood flow resistance; the directional alignment of neuronal axons ensures rapid transmission of nerve signals, and abnormal orientation can lead to cognitive impairment. The orientation of biological fibers is also a "peripheral biomarker" for disease diagnosis. Detecting the orderliness of collagen fiber arrangement can differentiate between benign and malignant tumor tissues; excessive collagen fiber deposition and abnormal orientation in pulmonary fibrosis can quantify the degree of fibrosis through orientation analysis; and the breakage and disordered orientation of elastic fibers can serve as early warning indicators of unstable atherosclerotic plaques.

[0003] Existing methods for calculating the spatial orientation of fibrous structures are limited and have certain limitations. First, most current methods are limited to two-dimensional image calculations, and orientation calculation in three-dimensional images remains challenging. Second, most methods can only calculate the global alignment trend of fibers, lacking detailed representation. Third, while existing weighted vector summation algorithms can achieve voxel-level fiber orientation representation, their accuracy is severely reduced for complex systems with large diameter variations because they use a fixed window size for all fibers in the image. The optimal window size is 2-4 times the fiber diameter.

[0004] Therefore, an adaptive method is needed to achieve fast and accurate fiber orientation calculation for complex systems with large differences in fiber diameter. Summary of the Invention

[0005] To address the shortcomings of the aforementioned technical solutions, this invention provides a fast window optimization method for calculating the spatial orientation of fiber structures based on a Gaussian mixture model. This method can achieve accurate and rapid calculation of the spatial orientation of three-dimensional fibers.

[0006] This method is implemented using the following technical solution:

[0007] This invention discloses a fast window optimization method for calculating the spatial orientation of fiber structures based on a Gaussian mixture model, comprising the following steps:

[0008] 1) Set a grayscale threshold to binarize the original fiber image, where pixels larger than the threshold are fibers with a value of 1, and pixels smaller than the threshold are background with a value of 0;

[0009] 2) Calculate the diameter of the fiber image in the binary image;

[0010] 3) Unfold the diameter spectrum of the fiber image into one-dimensional data, remove the voxels with a diameter of 0, and obtain a diameter-frequency histogram;

[0011] 4) Fit the diameter-frequency histogram using a Gaussian mixture model and calculate the diameter division threshold to divide the fibers into several groups according to their diameter.

[0012] 5) Apply the weighted vector summation algorithm to each group of fibers after division to calculate the spatial orientation, and obtain the fiber orientation map of each group;

[0013] 6) Merge the fiber orientation maps of each group to obtain the orientation distribution of the original fiber image;

[0014] 7) Perform pseudo-color encoding on the orientation of the fiber image to visualize the calculation results.

[0015] As a further improvement, step 2) of the present invention is specifically as follows: First, calculate the DT map of distance transformation. The DT value of the fiber pixel in the DT map is the shortest distance to the background pixel, and the DT value of the background pixel is 0. Then, take the maximum value of DT in the DT map and round it to get armd. Construct a cube window with a side length of 2armd+1. With each fiber pixel in the DT map as the center of the window, sort the non-center fiber pixels in the window according to the DT value from largest to smallest. Find the first pixel whose distance to the center of the window is less than its DT value. Take the DT value corresponding to this pixel as the DTT value of the center pixel of the window. Set the DTT of the non-fiber pixels to 0 to generate a DTT map. Construct a cube window with a size of 2*ArmdT+1 with each fiber pixel in the DTT map as the center. ArmdT is obtained by rounding the DTT value of the center pixel of the window. Calculate the average value of the DTT values ​​of all non-center fiber pixels in the window and multiply it by 2. The result is the diameter value of the center of the window, forming a diameter map.

[0016] As a further improvement, step 4) of the present invention specifically involves: weighting the diameter-frequency histogram so that its ordinate ranges from (-1, 1) to meet the input requirements of the Gaussian mixture model; inputting the weighted diameter-frequency histogram into the Gaussian mixture model; sorting the Gaussian distributions according to their mean values; and calculating the diameter value corresponding to the intersection of adjacent distributions as the diameter division threshold.

[0017] As a further improvement, in step 5) of this invention, the spatial orientation is calculated by applying a weighted vector summation algorithm to each group of fibers after division. Specifically, the average diameter of all fibers in the image is used as the half-window size, and the orientation of the center pixel of the window is calculated by multiplying the first and last two pixels symmetrical about the center of the window by a weight factor. and Then, sum them all together to obtain the vector orientation, which is the orientation of the center pixel of the window. The weighting factors are... Represented as:

[0018] ;

[0019] Where L is the length of the vector, and the weighting factor is... Represented as

[0020] ;

[0021] in These are the values ​​of the first, last, and center pixels of the vector. It is the average value of three pixels.

[0022] The beneficial effects of this invention are:

[0023] 1. Compared with existing fiber orientation calculation methods, this method can achieve high-resolution resolution of fiber orientation at the voxel level. Traditional methods can usually only reflect the overall orientation trend of an image, with limited ability to capture details of local structures. This method calculates orientation independently pixel by pixel, achieving fine characterization of individual fibrous tissues. This supports more accurate analysis of fiber orientation, curvature changes, and microscopic morphological heterogeneity, providing a more reliable data foundation for quantitative research on complex fiber systems.

[0024] 2. This invention is applicable to three-dimensional image processing, overcoming the limitations of most traditional techniques that are confined to two-dimensional analysis. While maintaining high accuracy, this method can achieve complete preservation and quantitative characterization of detailed information in three-dimensional structures, possessing higher spatial resolution and computational precision. This three-dimensional analysis capability significantly expands the applicability of the method, enabling it to meet the analysis needs of multi-scale biological images, from microscopic tissues to clinical imaging, and has significant scientific value for a deeper understanding of the three-dimensional structure-function relationship of biological tissues.

[0025] 3. Compared to the weighted vector algorithm using a fixed window size in the background technology, this method can achieve high-precision characterization of complex fiber systems with a wide range of diameter variations. Traditional methods, due to the use of a single-size calculation window, are difficult to adapt to the local characteristics of fibers with different diameters, resulting in large calculation deviations in complex systems. This invention significantly improves the computational accuracy and robustness of complex heterogeneous systems by calculating the fiber diameter at each voxel level and adaptively grouping fibers with different diameter ranges based on a Gaussian mixture model, thus providing a more reliable quantification tool for morphologically diverse biological fiber networks.

[0026] 4. Compared to traditional fiber orientation calculation methods that only output global statistical values, this invention uses pseudo-color coding technology to visualize fiber orientation based on voxel-level orientation characterization results. This visualization scheme can be used to display the three-dimensional spatial distribution and orientation changes of fibers in biological tissues, providing researchers with an intuitive and interactive analysis interface, and facilitating in-depth exploration of the relationship between biological tissue structure and function.

[0027] This invention first calculates the fiber diameter at the location of each fiber pixel in the binary fiber image through distance transformation and a series of operations, thereby obtaining the voxel-level fiber diameter. Attached Figure Description

[0028] Figure 1 This is a schematic diagram of the process of the present invention;

[0029] Figure 2 The result of fitting a diameter-frequency histogram to a Gaussian mixture model and calculating the diameter threshold is shown in the figure.

[0030] Figure 3 This is a schematic diagram comparing the effects of the present invention with those of a fixed window. Detailed Implementation

[0031] To describe the present invention in more detail, the technical methods of the present invention will be described in detail below with reference to the accompanying drawings and specific embodiments.

[0032] This invention presents a fast window optimization method based on weighted vector summation for calculating fiber three-dimensional orientation. The specific steps are as follows: Figure 1 As shown.

[0033] Taking a three-dimensional simulated fiber image with a fiber diameter ranging from 2 to 18 pixels as an example. The fiber image is a three-dimensional grayscale image. First, the grayscale values ​​of the image are normalized, and then filtered according to a set threshold to obtain a binarized image, in which the fiber part is identified as a value of 1, and the background part is identified as a value of 0.

[0034] Based on the binarized image, the fiber diameter is calculated. First, the distance transform (DT) map is calculated. In the DT map, the DT value of each fiber pixel is the shortest distance to the background pixel, and the DT value of the background pixel is 0. Next, the maximum DT value of the DT map is taken and rounded to obtain armd. A cube window with a side length of 2armd+1 is constructed. With each fiber pixel in the DT map as the center of the window, the non-center fiber pixels within the window are sorted in descending order of DT value. The first pixel whose distance to the center of the window is less than its DT value is found, and its corresponding DT value is used as the DTT value of the center pixel of the window. The DTT of non-fiber pixels is set to 0, generating a DTT map. A cube window with a size of 2*ArmdT+1 is constructed with each fiber pixel in the DDT map as the center, where ArmdT is obtained by rounding down the DTT value of the center pixel of the window. The average DTT value of all non-center fiber pixels within the window is calculated and multiplied by 2. The result is the diameter value of the fiber at that pixel, forming the diameter map.

[0035] The diameter map is expanded into one-dimensional data, and terms with zeros are removed to obtain a diameter-frequency histogram. This histogram is then fitted with a Gaussian mixture model, with the number of Gaussian distributions determined empirically. The Gaussian distributions are sorted according to their mean values, and the intersection points of adjacent distributions correspond to diameter values ​​used as thresholds for grouping fibers of different diameters.

[0036] Figure 2 The results of fitting the diameter-frequency histogram using a Gaussian mixture model and calculating the diameter threshold are presented. , , These are the means of three Gaussian distributions. , The diameter classification threshold is determined by the intersection of adjacent Gaussian distributions. The figure shows that the diameter classification thresholds for this example are 10.6 voxels and 14.7 voxels, respectively. Therefore, the fibers in this example are divided into three groups according to the diameter ranges of 0~10.6, 10.6~14.7, and greater than 14.7.

[0037] Next, the spatial orientation is calculated for each group of fibers. The mean of the fiber diameter in the grouped fiber image is rounded to obtain d. A cube window with a side length of 2*d+1 is constructed with each fiber pixel as the center. All two pixels symmetrical about the center of the window are multiplied by a weight factor as the first and last vectors. and Then, all the vectors are added together, and the resulting vector orientation is used as the orientation of the center pixel of the window, thus obtaining the spatial orientation map of each group of fiber images. and The expressions are respectively

[0038] and

[0039] Where L is the length of the vector. These are the values ​​of the first, last, and center pixels of the vector. It is the average value of three pixels.

[0040] Finally, the spatial orientation maps of each group of fibers are superimposed to obtain the spatial orientation map of the original fiber image. The obtained spatial orientation map is then pseudo-color encoded, and a MATLAB function is used to convert the spatial orientation map into a color image. Finally, HSV's colormap is used to convert the color spatial orientation map into a new pseudo-color encoded image, thus visualizing the result.

[0041] Figure 3 This diagram illustrates a comparison of the effects of the present invention and the initial method on a three-dimensional simulated fiber with a diameter range of 2-18 pixels. RWO represents the effect of the present invention, and window=21, window=41, and window=61 represent the calculation of the fiber polar angle using cubic windows with side lengths of 21, 41, and 61 pixels, respectively, on the fiber image. and The results. The calculations of this invention. and The angle error is significantly smaller than that of the fixed window, indicating that this method can improve the computational accuracy of complex systems and realize the three-dimensional orientation visualization of fibrous structures, which is of great significance for the study of biological microstructures and disease analysis.

[0042] The above description of the examples is provided to enable those skilled in the art to understand and apply the present invention. It will be apparent to those skilled in the art that various modifications can be made to the above examples, and the general principles described herein can be applied to other embodiments without inventive effort. Therefore, the present invention is not limited to the above embodiments, and any improvements and modifications made to the present invention by those skilled in the art based on the disclosure thereof should be within the scope of protection of the present invention.

Claims

1. A fast window optimization method for calculating the spatial orientation of fiber structures based on a Gaussian mixture model, characterized in that... Includes the following steps: 1) Set a grayscale threshold to binarize the original fiber image, where pixels larger than the threshold are fibers with a value of 1, and pixels smaller than the threshold are background with a value of 0; 2) Calculate the diameter of the fiber image in the binary image; 3) Unfold the diameter spectrum of the fiber image into one-dimensional data, remove the voxels with a diameter of 0, and obtain a diameter-frequency histogram; 4) Fit the diameter-frequency histogram using a Gaussian mixture model and calculate the diameter division threshold to divide the fibers into several groups according to their diameter. 5) Apply the weighted vector summation algorithm to each group of fibers after division to calculate the spatial orientation, and obtain the fiber orientation map of each group; 6) Merge the fiber orientation maps of each group to obtain the orientation distribution of the original fiber image; 7) Perform pseudo-color encoding on the orientation of the fiber image to visualize the calculation results.

2. The fast window optimization method for calculating the spatial orientation of fiber structures based on a Gaussian mixture model as described in claim 1, characterized in that, Step 2) specifically involves: First, calculating the distance transformation (DT) map. The DT value of the fiber pixel in the DT map is the shortest distance to the background pixel, and the DT value of the background pixel is 0. Then, taking the maximum DT value of the DT map and rounding it to obtain armd, constructing a cube window with a side length of 2armd+1. Using each fiber pixel in the DT map as the center of the window, sorting the non-center fiber pixels in the window according to their DT values ​​from largest to smallest, finding the first pixel whose distance to the center of the window is less than its DT value, and using the DT value corresponding to this pixel as the DTT value of the center pixel of the window. The DTT of non-fiber pixels is set to 0, generating a DTT map. Using each fiber pixel in the DTT map as the center, constructing a cube window with a size of 2*ArmdT+1, where ArmdT is obtained by rounding the DTT value of the center pixel of the window. Calculating the average value of the DTT values ​​of all non-center fiber pixels in the window and multiplying it by 2, the result is the diameter value of the center of the window, forming a diameter map.

3. The fast window optimization method for calculating the spatial orientation of fiber structures based on a Gaussian mixture model as described in claim 1 or 2, characterized in that, Step 4) specifically involves: weighting the diameter-frequency histogram so that its ordinate ranges from (-1, 1) to meet the input requirements of the Gaussian mixture model; inputting the weighted diameter-frequency histogram into the Gaussian mixture model; sorting the Gaussian distributions according to their mean values; and calculating the diameter value corresponding to the intersection of adjacent distributions as the diameter division threshold.

4. The fast window optimization method for calculating the spatial orientation of fiber structures based on a Gaussian mixture model as described in claim 3, characterized in that, In step 5), the spatial orientation of each group of fibers after division is calculated using a weighted vector summation algorithm. Specifically, the average diameter of all fibers in the image is used as the half-window size. The orientation of the center pixel of the window is calculated by multiplying the first and last two pixels symmetrical about the center of the window by a weight factor. and Then, sum them all together to obtain the vector orientation, which is the orientation of the center pixel of the window. The weighting factors are... Represented as: ； Where L is the length of the vector, and the weighting factor is... Represented as ； in These are the values ​​of the first, last, and center pixels of the vector. It is the average value of three pixels.

Images