Brain MRA 3D Skeleton Extraction Based on Normal Plane Centroid Algorithm

INTRODUCTION: Analysis of magnetic resonance angiography image data is crucial for early detection and prevention of stroke patients. Extracting the 3D Skeleton of cerebral vessels is the focus and difficulty of analysis. OBJECTIVES: The objective is to remove other tissue components from the vascular tissue portion of the image with minimal loss by reading MRA image data and performing processing processes such as grayscale normalization, interpolation, breakpoint detection and repair, and image segmentation to facilitate 3D reconstruction of cerebral blood vessels and the reconstructed vascular tissues make extraction of the Skeleton easier. METHODS: Considering that most of the existing techniques for extracting the 3D vascular Skeleton are corrosion algorithms, machine learning algorithms require high hardware resources, a large number of learning and test cases, and the accuracy needs to be confirmed, an average plane center of mass computation method is proposed, which improves the average plane algorithm by combining the standard plane algorithm and the center of mass algorithm. RESULTS: Intersection points and skeleton breakpoints on the Skeleton are selected as critical points and manually labeled for experimental verification, and the algorithm has higher efficiency and accuracy than other algorithms in directly extracting the 3D Skeleton of blood vessels. CONCLUSION: The method has low hardware requirements, accurate and reliable image data, can be automatically modeled and calculated by Python program, and meets the needs of clinical applications under information technology conditions.


Introduction
With software and hardware technology development, medical imaging equipment has also been widely used in hospitals and rehabilitation institutions of different sizes.Among them, NMR [1] is favored by doctors and patients in the diagnosis of cerebrovascular diseases due to its unique imaging principle, no ionizing radiation, no skeletal artifacts, high resolution, multi-directional and multi-sequence imaging, and vascular imaging without the use of contrast agents."cerebrovascular accident" (CVA), also known as "cerebral stroke" (cerebral stroke), "stroke," is due to the blockage of blood vessels in the brain resulting in poor blood flow or because of the sudden rupture of blood vessels and brain tissue damage disease.The incidence and prevalence of stroke remain high in the world.The survey found that from November 2018 to December 2021, the incidence of stroke in China was about 28 million for men and 18 million for women [2], and the number was still increasing, among which Feng • Zhu • Li ischemic stroke accounted for 60%-80% of strokes, higher than the incidence of hemorrhagic stroke.Stroke has become the first cause of death in China and the leading cause of disability in adults [3][4].Therefore, most extensive and medium-sized hospitals in China have set up special stroke centers.For the diagnosis and treatment of stroke, due to the lack of effective treatment means, prevention is currently considered to be the best measure, and early detection and diagnosis can guide patients to take preventive measures in advance to save their health and lives.In recent years, imaging technology has made significant progress [5].The imaging clarity of Magnetic Resonance Angiography (MRA), CT angiography (CTA), and other images is getting higher and higher.The analysis effect of vascular diseases is also improving, especially the early analysis effect of stroke caused by cerebrovascular is more significant.By analyzing brain MRA, CTA, and other image data, it is of great significance to obtain the analysis results of cerebrovascular changes, search the lesion area of stroke, and guide the early clinical diagnosis of stroke in patients.Compared with CTA, MRA imaging does not require a contrast agent and has no damage to subjects.Brain MRA examination can be used as a physical screening item in middle-aged and older adults, so the importance of quantitative analysis of MRA imaging is self-evident.

Results
Vascular skeleton extraction is a necessary step in the quantitative analysis of MRA image data.Because most of the existing skeleton extraction algorithms are twodimensional (also known as binary image thinning), 3D image skeleton extraction algorithms are relatively few and are mainly classified as corrosion algorithms, including the iterative thinning method used by MASSEY et al. [6].The minimum path method proposed by JIN et al. [7] and other machine learning methods emerged with the development of artificial intelligence technologyhowever, the precision of corrosion algorithm refinement in this method is needed.The deep learning algorithm proposed by WAGNER et al. [8] requires many case models, high requirements on hardware resources, expensive Graphics Processing Unit (GPU) cards, and large resource consumption.Under the same hardware and software conditions, the improved average plane centroid algorithm has apparent advantages in efficiency and accuracy compared with the corrosion algorithm and machine learning algorithm to extract the three-dimensional Skeleton of MRA blood vessels.This method's defect is that the early stage preprocessing process will cause specific errors.With the development of artificial intelligence technology, the machine processing speed increases.In addition, in the case of using a large number of samples for deep learning, the machine learning algorithm has more room for improvement.The next step that must be verified is improving the hardware configuration.The effect will be better if the machine learning algorithm is combined with the average plane centroid algorithm to extract the Skeleton of three-dimensional images.Through the use of three-dimensional vascular Skeleton, the cross-sectional area of each section of blood vessel can be reshaped, the changes of vascular cross section can be more intuitively understood, the narrow and dilated weak areas of blood vessels can be accurately located, and the relationship between brain region and changes in cerebrovascular blood supply can be analyzed in combination with the brain structure data of MRIT1, which is of great significance for guiding the prevention of cerebrovascular diseases in patients, especially in the middle-aged and elderly.

Materials and methods
To extract vascular Skeleton from brain MRA data, image data must be removed first, the coordinate system must be transformed, image preprocessing must be carried out, noise, bone, and other tissue interference factors can be eliminated, pure blood vessels and single background images can be obtained, and then the vascular Skeleton can be extracted.

3.1Image transformation 3.1.1Reading data
Image data transformation refers to the original image data after processing it into a standard image file and adjusting the accuracy of the three dimensions to facilitate the following calculation process.Image transformation includes normalization (normalization of the check value to the 0-255 image space), interpolation processing (to make the accuracy of the sagittal, coronal, and axial directions of the image data consistent), breakpoint detection and repair (to supplement some pixel abnormalities or missing caused by noise, etc.) and other processing processes.The above process is also called image preprocessing in image analysis.MRA original image data (dicom format file) includes header and image data.As shown in Figure 1, image processing needs to read out part of the image.The Python pydicom package or SimpleITK reads the dicom file and extracts part of the image data.

Data normalization
Data normalization usually converts data in different ranges into the same data range to facilitate subsequent analysis and processing.It is an essential data preprocessing method.There are many methods for data normalization, including maximum and minimum value methods, Z-score methods, decimal scaling methods, logarithmic function methods, and deviation standardization methods [9,10].Most normalization methods are based on the need to normalize the value to between 0 and 1 or -1 and 1.The value of image data is usually between -1000 and 1000.Figure 2 shows the selected window width and position data during image inspection.It can be seen from the figure that the maximum pixel value of MRA image data with the adjusted window width and window position is 658, and the minimum pixel value is 0. Here, the improved maxmin method is used, as shown in equation ( 1), to convert the pixel values of 0~658 to the gray level of 0~255, that is, 256 color levels, to become the image data that genuinely meets the image standard, as shown in Figure 3.

Figure 2 MRA image data with adjusted window width and window level
Select a window width of 447 and a window center of 447.After changing the window width and window level, the maximum pixel value of the MRA image data is 658, the minimum pixel value is 0, and the scanning accuracy voxel value is 0.390625\0.390625.Use the onis2.5 software to open and display the header file data of the case and take a partial display.
The normalized pixel value is new x , the original pixel value is x , the maximum value of the actual pixel value is max x , and the minimum value of the initial pixel value is min x .

Interpolation processing
The original MRA data is three-dimensional, and the dimensionality reduction is displayed for convenience.The image needs to be interpolated to obtain the consistency of three-dimensional accuracy.Here, the following convolutional three-interpolation algorithm is used: 1 2abs abs , 0 abs 1 4 8abs 5abs abs ,1 abs 2 0 , abs 2 Using the gray values of 16 points around the sampling point for cubic interpolation, taking into account the influence of gray values of 4 directly adjacent points and the power of gray value change rate to each neighboring point, the BiCubic function is ( )

Sx
x the distance from the interpolation point to the control point.Does any point on the image represent a tensor function , and i 、 j 、 u 、 v are coordinates (the third direction is a constant value or is measured by the modulus of the third direction tensor),

 
B is a tensor function of the sampling point,   A and   C are a tensor function of the affected neighborhood.After changes, the above algorithm adjusts the three directional value arrays.Then, the array values are written into the image data, and the updated image data is obtained through GetArrayFromImage [11][12][13][14].Because it is necessary to compare the experimental machinegenerated value with the manual annotated value to reduce the comparison error between the two, 0 interpolation is required, and standard interpolation is needed when only machine processing is done.

Breakpoint detection repair
Gaussian filters can detect and smooth the image in medical images to remove noise and repair minor breakpoints [15,16].It filters all pixels in the embodiment to be processed through the specified template (convolution, mask) and replaces the scanned pixel value with the weighted average gray matter of the pixel in the neighborhood obtained by calculation.Its functional form is as follows: ( ) x 、 y 、 z represents the pixel values at three corresponding positions in three-dimensional space, e represents the constant of nature,  and represents PI.
Here, the gaussian_filter function of Python is used to write a three-dimensional filtering program and perform three-dimensional Gaussian filtering.Selecting a sigma value below 0.4 is best because significant parameters easily blur small blood vessels and cause unclear boundaries.The two-dimensional display of Gaussian filtering images using sigma 0.3 is shown in Figure 4.

Image segmentation
MRA image segmentation must remove other tissues as much as possible, mainly bone tissue.Other soft tissues have little impact on vascular imaging, and the vascular part is preserved to the maximum extent.There are many image segmentation technologies, and the segmentation algorithms include the combination of traditional segmentation algorithms and a variety of segmentation algorithms, such as threshold-based segmentation algorithm, region-based segmentation algorithm (region growth algorithm, split and merge algorithm, watershed algorithm, etc.), and edge detection-based segmentation algorithm.There are also algorithms based on specific tools (segmentation algorithm based on clustering, segmentation algorithm based on wavelet transform, segmentation algorithm based on active contour model, genetic algorithm, etc.) [17][18][19][20] and algorithms based on deep learning [8,[21][22][23], etc. Deep learning algorithms are also subdivided into various algorithms.Each of these segmentation algorithms has its advantages and disadvantages.To improve the computing speed and reduce the requirement on hardware computing power without losing accuracy, the traditional watershed algorithm and threshold algorithm are combined to perform MRA image segmentation, as shown in Figure 5.The watershed algorithm is a mathematical segmentation algorithm based on topological theory.It obtains the edge information of the image by calculating the gradient.The algorithm is as follows: Suppose that for a point ( ) ,, P x y z on the image, in the range PD  , D is the image space that conforms to the function ( , , ) u f x y z = , then there is the gradient value at the end P : i ， j k are unit vectors in three directions.
Then, threshold processing is performed on the gradient image to eliminate excessive segmentation caused by small changes in the grayscale.After threshold processing, the segmentation function becomes: ( ) g represents the threshold variable.
The calculation of threshold variables usually includes the maximum inter-class variance method (otsu algorithm), median method, mean method, maximum entropy method (entropy method), iterative method, etc.Here, the complete inter-class variance method is used, and its threshold algorithm is as follows: ( ) The interclass variance is m

Three-dimensional reconstruction
At present, image 3D reconstruction is divided into two directions: traditional 3D reconstruction based on visual geometry [24][25][26][27] and 3D reconstruction based on deep learning [28,29].Since voxel values of 3D images are already available here, the SIFT detection algorithm [24] is only needed to obtain feature points in the above photos.Then, the FLANN matching algorithm [25] matches the feature points to different images.Finally, the traditional 3D reconstruction algorithm -SFM algorithm [26] is used to reconstruct the above-processed photos, which can be divided into incremental, global, hybrid, hierarchical, and semantic-based methods of the SFM algorithm [30][31][32].A hierarchical algorithm is used for image 3D reconstruction, as shown in Figure 6.
Figure 6 3D Reconstruction Effect of MAR Preprocessed Image Although the image was segmented using a watershed algorithm and threshold correction, minor omissions may occur due to noise, and skull bones may be difficult to remove due to high pixel values and poor continuity.
To obtain the simple vascular part image, the connectivity analysis of 6 neighborhood areas was performed on the three-dimensional image, and the number, location, surface area, and volume of each part of the connectivity were calculated according to the voxel number.The statistical results are shown in Table 1.According to the obtained data and 3D images, it can be seen that there are still some tiny blood vessels (less than 600 in volume), bone fragments (less than 1700 in volume), and interference noise (less than 100 in book as recommended by clinical experts) that are not connected with large blood vessels.A Gaussian filter must filter out the bone fragments and interference noise.The rest is the blood supply vessels outside the brain that can also be filtered out.After Gaussian filtration, pure blood vessels supplying blood to the brain were reconstructed in three dimensions, as shown in Figure 7.  From the perspective of the 3D display, only the vascular part has been preserved.

Skeleton Extraction
The most commonly used 3D skeleton extraction algorithm is the corrosion algorithm [6,33,34], which corrodes the outermost outline from the outer layer and iterates continuously until the innermost Skeleton of the image remains.The centroid algorithm [7,[35][36][37] is also used in the skeleton extraction of two-dimensional images.In this paper, the centroid algorithm is improved, and the cross-section of blood vessels is located and reshaped by adding average plane calculation.The procedure of the standard plane centroid algorithm is first to regard the blood vessel as a multi-segment curve.The middle plane of the curve and the cross-section of the blood vessel form many small irregular elliptical crosssections D, and the centroids of these small elliptical cross-sections are connected to form the Skeleton of the blood vessel.Let a small segment of the blood vessel curve be derivable in this range, and the derivative ()  is not all 0, then the ordinary plane equation at the point ( ) For each normal plane and the cross-section of blood vessels, according to the plane sheet centroid theory, there is a centroid coincidence: The particle distributed on the cross-section is i m the total mass of the particle system, the static moment of the particle system concerning the Y-axis, and the static moment of the particle system relating to the X-axis.
Our sections here all meet the conditions of homogeneous flakes, so the above equation ( 12) can be converted to: Where A is the area of section D and  is the microregion that meets the requirements of the function.Python's NumPy library calculates the normal plane and cross-section points determined by the points on the periphery of blood vessels with the threshold value in the image.Then, the facts on the Skeleton are calculated according to the homogenized slice centroid equation.The positions and voxel values of these points are extracted by changing the voxel values of this point.Make it have a uniform and distinguishable voxel value with other centroid points) to form a new threedimensional model [37,38], as shown in Figure 8.
Figure 8 3D skeleton rendering From the horizontal, coronal, and sagittal positions, it can be seen that the vascular Skeleton extracted through the average plane centroid algorithm has good uniformity and continuity and is complete.

Data source
Moreover, from a total of 113 normal subjects, the data from the brain magnetic resonance (NMR) public data sets in the https://public.kitware.com/Wiki/TubeTK/Datafor 100 cases, with the approval of the patient collected 13 subjects of MRA image data in the brain, Acquisition equipment Siemens magnetic resonance instrument model MAGNETOM skyra3.0T,scanning sequence 3D-fl3d1r_t70, scanning accuracy 0.39mm, layer thickness 0.6mm, image size 512*464.

Experimental environment
The hardware environment used here is similar to the current mainstream computer, and the software uses a lower open-source stable version, which is convenient for practical application units such as hospitals and medical institutions.Hardware environment: The computer is configured with i7-8520U CPU, 16GDDRIII memory, and 500GSSD hard disk.Software environment: The operating system environment is CentOS6.5,Python version 3.6, torch version 1.7, and other dependency packages and libraries.

Verification method 4.3.1 Efficiency Comparison Experiment
The efficiency of 3D modeling and skeleton extraction has a great relationship with the software and hardware configuration.The size, accuracy, and thickness of the sweeping layer of the image all affect the operation efficiency.Under the above software and hardware configuration environment and sample data, the maximum CPU usage of the 3D modeling machine is 10%, the minimum 0.5%, the maximum memory consumption is 4%, the minimum 1.5% and the average time is 48ms; the maximum CPU usage of the skeleton extraction is 8%, the minimum memory consumption is 1.6%, the maximum memory consumption is 3%, the minimum memory consumption is 0.7%, and the average time is 34ms.The skeleton extraction efficiency pairs under similar conditions are shown in Table 2.  [39]   3D erosion algorithm 21.1mins 2Jin et al.,2016 [7] 3D erosion algorithm 82sec

Accuracy analysis
By taking the key points in the three positions of coronal position, disfigures-position, and axial position of the image, the machine automatically selects the endpoint specified by the Skeleton, as shown in Figure 9; the manual selection is verified through communication with clinical experts, as shown in Figure 10; and then the original data position measurement through the subject data and the corresponding position measurement after the formation of the vascular Skeleton is verified.Based on the structure of the blood vessel image, the branch intersection points and the endpoints of each segment of blood vessels (that is, the breakpoint) are selected as the key points.The critical point locations measured in the original image (Table 3) are compared with the corresponding points of the blood vessel skeleton (Table 4).If the fitting probability of the position points exceeds 97% of the experimental design, it is considered that the experimental requirements are met.The expected probability is calculated as follows:  That is to say, from the original image coordinate to the converted image coordinate also occurs rotation, with the Y axis as the rotation axis, the vertical plane rotates 90º clockwise, the actual Z axis becomes the present X axis, and then with the current X axis as the rotation axis, rotate the X-axis normal plane 90º counterclockwise, the Y axis becomes the present Z axis.The X axis becomes the present Y axis.Then, according to the above equation ( 14), the P values of image skeleton line fitting for all subjects were calculated to be the maximum 99.72% and the minimum 99.35%, respectively.The expected probability [40,41] E P =99.61% is obtained by calculating the average value.The 3D skeleton extraction accuracy pairs under similar conditions are shown in Table 3.

Figure. 1
Figure. 1 Raw imaging data of brain MRA Randomly select a subject case, use the onis2.5 software to open the header file data and 2D image of the original image file that displays the case, and capture some photos for display.

Figure 3
Figure 3 Comparison of images before and after normalization processing Through image normalization, the voxel values of the image data, in this case, changed from 0-658 to 0-255, indicating a significant increase in edge clarity from the image.min max min

Figure. 4
Figure. 4 Comparison of images before and after Gaussian filtering Through Gaussian filtering of the image, the edge of the image, in this case, is enhanced, and the clarity is more evident in the image comparison display.

u 0 w 1 w
the ratio of foreground voxel number to population prime number, the ratio of background voxel number to population prime number, 0 u the average gray value of voxel foreground, 1 u and the average gray value of voxel background.

Figure. 5
Figure. 5 Comparison of images before and after watershed and threshold correction processing By performing watershed algorithm segmentation and threshold correction on the filtered image, only the vascular images were retained, while other tissues were removed.

Figure 7 3D
Figure 7 3D Reconstruction Effect of Filtered Image Calculating the volume and position of each discontinuous part of the 3D reconstructed image can eliminate the bone tissue within the skull region and position deviation or noise data suggested by experts.From the perspective of the 3D display, only the vascular part has been preserved.

14 )P
represents the probability of fitting the attempted image skeleton, i p the probability of fitting the critical point, n the selected vital point, i A the difference vector of the vital point, and i B the actual vector value of the critical point (the skeleton position vector is used here).

Figure 9
Figure 9 Selection of Key Points for 3D SkeletonAll key points that need to be verified are selected at the top and fork positions of the Skeleton.The advantage of this selection method is that the work is easy to determine, manual annotation is more accessible, and it is convenient for later comparison and verification.

Table 1
Statistical Table of Connectivity of 3D Modeling of Images after MAR Preprocessing

Table 1 -
4Statistics of surface area and volume of main connected blood vessels ID Center-X

Table 3
Comparison of accuracy of 3D image skeleton extraction ID references methods accuracy