Fast sub-voxel accuracy estimation of surface area of region of interest in voxel model

Abstract

The application discloses a kind of sub-voxel precision fast estimation method of surface area of region of interest in voxel model, belong to image processing technical field, including the division of ROI region;Extract the surface voxel of ROI region, search the voxel 6 near neighbor of ROI region voxel, when there is non-ROI region voxel, current voxel is marked as surface voxel, initialization defines the balance center position of this surface voxel as absolute center position;Surface voxel sub-voxel precision balance center position is defined using iterative method;After obtaining the surface voxel balance center position after convergence, the surface area value corresponding to surface voxel is calculated;Surface area gray value of surface voxel is processed, the resolution of surface voxel is reduced until resolution is 1, and the surface area of ROI region is obtained.The application uses the above method, can adapt to complex structure and reduce the dependence on preset model, provide reliable support for high-precision three-dimensional analysis, and improve work efficiency.

Description

Technical Field

[0001] This invention relates to the field of image processing technology, and in particular to a fast sub-voxel-accurate estimation method for the surface area of ​​regions of interest in a voxel model. Background Technology

[0002] In fields such as 3D modeling and medical imaging, industrial inspection, and virtual reality, surface area estimation using voxel models plays a crucial role in quantitative analysis (e.g., tumor morphology assessment, material surface property calculation). However, existing methods still have significant limitations in terms of accuracy, efficiency, and adaptability to complex structures. Traditional mainstream technologies mainly include:

[0003] Triangulation: This method first converts the voxel model into a triangular mesh, then calculates the sum of the areas of the triangular meshes. A representative algorithm is Marching Cubes, which determines the triangular face type using voxel values ​​within a 2x2x2 area and ultimately generates the model. While this method can improve surface smoothness, it requires traversing all voxels and interpolating to generate triangular facets, resulting in high computational complexity, large memory consumption, and low efficiency. It is difficult to meet the needs of real-time computing.

[0004] Voxel counting method: The surface area is estimated by accumulating the number of exposed voxels. Although it is computationally efficient, it leads to overestimation error due to the staircase effect and cannot capture the geometric details at the sub-voxel scale.

[0005] Geometric fitting method: It relies on preset regular geometric primitives (such as ellipsoids and cylinders) to estimate the area. It has poor adaptability to complex topological structures (such as porous media and fractal surfaces) and requires manual intervention in model selection, resulting in low automation.

[0006] The common bottleneck of the aforementioned methods lies in the inherent contradiction between accuracy and efficiency: high-precision algorithms (such as MarchingCubes) are too computationally expensive to meet real-time requirements, while efficient methods (such as voxel counting) are limited by the discretization characteristics of voxels and cannot break through the technical ceiling of sub-voxel accuracy. Furthermore, existing technologies lack the ability to track surface areas of non-uniformly distributed, small features, and dynamically deformed models, further limiting their application in detailed scenarios. Therefore, there is an urgent need for a surface area estimation method that balances sub-voxel accuracy and computational efficiency, adapts to complex structures, and reduces reliance on pre-defined models, thereby providing reliable support for high-precision 3D analysis. Summary of the Invention

[0007] The purpose of this invention is to provide a fast sub-voxel-accurate estimation method for the surface area of ​​regions of interest in voxel models, in order to solve the problems existing in the background art.

[0008] To achieve the above objectives, this invention provides a rapid sub-voxel accuracy estimation method for the surface area of ​​a region of interest in a voxel model, comprising the following steps:

[0009] S1. Divide the ROI region and set the grayscale threshold to determine the voxels of the ROI region;

[0010] S2. Extract the surface voxels of the ROI region, and search the six nearest neighbor voxels of the voxels in the ROI region (up, down, left, right, front, and back). When there are voxels in the non-ROI region, the current voxel is marked as a surface voxel, and the equilibrium center position of the surface voxel is initialized and defined as the absolute center position of the surface voxel.

[0011] S3. Use an iterative method to define the equilibrium center position of the surface voxel subvoxel precision;

[0012] S4. After obtaining the position of the equilibrium center of the surface voxel after convergence, calculate the surface area value of the corresponding surface voxel.

[0013] S5. Process the surface area grayscale value of the surface voxels to reduce the resolution of the surface voxels until the resolution is 1, and obtain the surface area of ​​the ROI region.

[0014] Preferably, in S1, a threshold gray value g is set, and voxels are traversed. Voxels with a gray value greater than g are defined as ROI region voxels.

[0015] Preferably, the surface voxels in S2 are denoted as 1, and the rest are 0, to obtain a binary voxel model of the surface voxels.

[0016] Preferably, the content of S3 is as follows:

[0017] For a surface voxel, query its surrounding 3x3x3 nearest neighbor voxels. If a surface voxel exists, extract the equilibrium center position C of that voxel. i The obtained nearest-neighbor voxel equilibrium positions are weighted and averaged, and the equilibrium center position of the current voxel is updated with the value. The equilibrium center position is constrained to be no more than 1 voxel unit different from the absolute center position. The specific formula is as follows:

[0018]

[0019] C(t+1) = clamp(C(t), C 0 +(-1,-1,-1),C 0 +(1,1,1));

[0020] Among them, C i C represents the subvoxel equilibrium center location of the nearest voxel. 0 Let w be the absolute center position of the current voxel, t be the iteration step, and w be the position of the absolute center position of the voxel. i The weight value for this voxel is the weight of the surface voxel (if it is a surface voxel, then w). i =1, otherwise w i=0), where i is the nearest voxel. For a voxel, there are 26 nearest voxels, and i ranges from 0 to 25.

[0021] If the iterative process of this surface voxel contains unstable details, then each iteration can use only the average of the old equilibrium center position and the new equilibrium center position as the new equilibrium center position, i.e.

[0022] C(t+1)=0.5*(C(t)+clamp(C(t),C 0 +(-1,-1,-

[0023] 1), C 0 +(1,1,1)));

[0024] The above is the calculation process for one equilibrium center. Repeating this step multiple times will gradually converge the equilibrium center to a stable value.

[0025] Preferably, for complex curved surfaces, there are three strategies for obtaining the position of the equilibrium center:

[0026] a. Limit the number of iterations to 32 steps. Generally, the voxel model corresponding to a sphere will not collapse within 32 iterations.

[0027] b. After each iteration, c(t+1) is decayed, as shown in the following formula:

[0028] c(t+1)=C(t+1)*α

[0029] Where α is the attenuation coefficient, which is generally controlled at around 0.9 for spherical surfaces;

[0030] c. Calculate the offset of the equilibrium center position of the nearest neighbor surface voxels, and adjust the equilibrium center position of the current voxel accordingly. The formula is:

[0031]

[0032] C(t+1) = C(t+1) - D*β;

[0033] Where β is the offset coefficient, which is generally controlled at around 0.1 for spherical surfaces.

[0034] Preferably, the content of S4 is as follows:

[0035] Within a 3x3x3 area around the voxel, the surface normal is calculated based on the position of the equilibrium center of the nearest voxel. Then, the surface area occupied by the current voxel is calculated based on the direction of the normal and the position of the equilibrium center of the voxel. The surface area of ​​the surface voxel is stored in the current voxel model as a gray value.

[0036] The normal is calculated not using the absolute center position of the voxel, but using the iteratively derived equilibrium center position, with the Z direction as the reference. The specific formula is as follows:

[0037]

[0038]

[0039] Where Cx, Cy, and Cz are the components of vector C in the three directions; i, j, and k are the voxel indices in the X, Y, and Z directions, taking values ​​of -1, 0, and 1 within a 3x3x3 range; the above formula is based on the z-direction and is not applicable to surfaces parallel to the z-direction. This can be easily transformed into recalculating the gradient based on either the x-axis or the y-axis.

[0040] Given the normal N, find the minimum angle θ between the normal and the X, Y, and Z planes. The surface area corresponding to this voxel is denoted as 1 / cos(θ), and the formula is as follows:

[0041] θ=arcos(max(abs(Nx),abs(Ny),abs(Nz)));

[0042] A=1 / cos(θ)=1 / max(abs(Nx),abs(Ny),abs(Nz));

[0043] Where abs is the absolute value and A is the surface area corresponding to the voxel.

[0044] Preferably, the content of S5 is as follows:

[0045] The grayscale values ​​of the surface voxels are processed to prepare a voxel mesh with a resolution half that of the current surface voxels. For the voxel with coordinates (i,j,k), the surface areas of the original voxel model with coordinates (2*i,2*j,2*k), (2*i+1,2*j,2*k), (2*i,2*j+1,2*k), (2*i+1,2*j+1,2*k), (2*i,2*j,2*k+1), (2*i+1,2*j,2*k+1), (2*i,2*j+1,2*k+1) are summed and recorded in the current voxel. At the same time, the resolution of the voxel model is reduced by half. This step is repeated until the voxel resolution is 1. The grayscale value of the stored surface area is the surface area of ​​the ROI region.

[0046] Therefore, the present invention employs a sub-voxel-precision fast estimation method for the surface area of ​​the region of interest in the above-mentioned voxel model, which has the following beneficial effects:

[0047] (1) By iterating the equilibrium center position of the surface voxel, sub-voxel precision information is obtained, and the surface normal corresponding to the surface voxel is calculated based on this information. The normal is smoother based on the sub-voxel precision information.

[0048] (2) It has good adaptability to various inclined planes or curves;

[0049] (3) It has high computational efficiency and the entire process can be processed in parallel, which improves processing efficiency and has a significant time advantage in the field of high-resolution image processing.

[0050] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0051] Figure 1 This is a flowchart illustrating the rapid sub-voxel precision estimation method for the surface area of ​​the region of interest in the voxel model of this invention.

[0052] Figure 2 This is a comparative diagram of the voxel equilibrium center positions of the method in this embodiment and the conventional method, wherein (a) is the conventional method and (b) is the method in this embodiment;

[0053] Figure 3 This is a comparative diagram of the voxel normal solution of the method in this embodiment and the conventional method, wherein (a) is the conventional method and (b) is the method in this embodiment;

[0054] Figure 4 This is a comparative diagram of the voxel-corresponding surface area calculation method of the present invention and the conventional method, wherein (a) is the conventional method and (b) is the method of the present invention. Detailed Implementation

[0055] The following detailed description of embodiments of the invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the invention without inventive effort are within the scope of protection of the invention.

[0056] Please see Figures 1-4 A fast sub-voxel precision estimation method for the surface area of ​​the region of interest in a voxel model includes the following steps:

[0057] S1. Divide the ROI region and set the grayscale threshold to determine the voxels of the ROI region.

[0058] The division of the ROI region can be based on setting a grayscale threshold according to actual conditions or needs. In this embodiment, a grayscale threshold method is used, and a threshold grayscale value g is set. The voxels are traversed, and voxels with grayscale values ​​greater than g are defined as ROI voxels, resulting in a binary voxel model of the ROI region, that is, the voxel value of the ROI region is 1, and the value of the other voxels is 0.

[0059] S2. Extract surface voxels from the ROI region. Search for the six nearest neighbor voxels (top, bottom, left, right, front, and back) of each voxel within the ROI region. If a voxel exists outside the ROI region, the current voxel is marked as a surface voxel and assigned a value of 1. The rest are assigned a value of 0, resulting in a binary voxel model of the surface. Initialize and define the equilibrium center position of this surface voxel as the absolute center position of the surface voxel.

[0060] S3. Use an iterative method to define the equilibrium center position of the surface voxel subvoxel precision.

[0061] For a surface voxel, query its surrounding 3x3x3 nearest neighbor voxels. If a surface voxel exists, extract the equilibrium center position C of that voxel. i The obtained nearest-neighbor voxel equilibrium positions are weighted and averaged, and the equilibrium center position of the current voxel is updated with the value. The equilibrium center position is constrained to be no more than 1 voxel unit different from the absolute center position. The specific formula is as follows:

[0062]

[0063] C(t+1) = clamp(C(t), C 0 +(-1,-1,-1),C 0 +(1,1,1));

[0064] Among them, C i C represents the subvoxel equilibrium center location of the nearest voxel. 0 Let w be the absolute center position of the current voxel, t be the iteration step, and w be the position of the absolute center position of the voxel. i The weight value for this voxel is the weight of the surface voxel (if it is a surface voxel, then w). i =1, otherwise w i =0), where i is the nearest voxel. For a voxel, there are 26 nearest voxels, and i ranges from 0 to 25.

[0065] If the iterative process of this surface voxel contains unstable details, then each iteration can use only the average of the old equilibrium center position and the new equilibrium center position as the new equilibrium center position, i.e.

[0066] C(t+1)=0.5*(C(t)+clamp(C(t),C 0 +(-1,-1,-1),C 0 +(1,1,1)));

[0067] The above is the calculation process for one equilibrium center. Repeating this step multiple times will gradually converge the equilibrium center to a stable value.

[0068] The above method is suitable for planes with various slopes; however, for complex surfaces, additional steps are required. Firstly, when the voxel model resolution is low, different complex surfaces may correspond to the same voxel model. In other words, the curvature information of the surface is lost during voxel model acquisition, and excessive iterations can lead to collapse for convex surfaces. Three solutions can be adopted to address this:

[0069] a. Limit the number of iterations to 32 steps. Generally, the voxel model corresponding to a sphere will not collapse within 32 iterations.

[0070] b. After each iteration, C(t+1) is decayed, as shown in the following formula:

[0071] C(t+1)=C(t+1)*α

[0072] Where α is the attenuation coefficient, which is generally controlled at around 0.9 for spherical surfaces;

[0073] c. Calculate the offset of the equilibrium center position of the nearest neighbor surface voxels, and adjust the equilibrium center position of the current voxel accordingly. The formula is:

[0074]

[0075] C(t+1) = C(t+1) - D*β;

[0076] Where β is the offset coefficient, which is generally controlled at around 0.1 for spherical surfaces.

[0077] All three schemes improve the adaptability of this embodiment to curved surfaces. The specific selection method is to draw the voxel equilibrium position normal obtained in the subsequent steps onto the voxel model and observe the smoothness of the normal to determine the required parameters and effects of the above three schemes.

[0078] S4. After obtaining the position of the equilibrium center of the surface voxel after convergence, calculate the surface area value of the corresponding surface voxel.

[0079] Within a 3x3x3 area around the voxel, the surface normal is calculated based on the position of the equilibrium center of the nearest voxel. Then, the surface area occupied by the current voxel is calculated based on the direction of the normal and the position of the equilibrium center of the voxel. The surface area of ​​the surface voxel is stored in the current voxel model as a gray value.

[0080] The normal is calculated using the traditional gradient formula, but instead of the absolute center position of the voxel, it is calculated based on the equilibrium center position obtained through iteration, with the Z direction as the reference. The specific formula is as follows:

[0081]

[0082] Where Cx, Cy, and Cz are the components of vector C in the three directions; i, j, and k are the voxel indices in the X, Y, and Z directions, taking values ​​of -1, 0, and 1 within a 3x3x3 range; the above formula is based on the Z direction and is not applicable to surfaces parallel to the Z direction. This can be easily transformed into recalculating the gradient based on the x-axis or y-axis.

[0083] Given the normal N, find the minimum angle θ between the normal and the X, Y, and Z planes. The surface area corresponding to this voxel is denoted as 1 / cos(θ), and the formula is as follows:

[0084] θ=arcos(max(abs(Nx),abs(Ny),abs(Nz)));

[0085] A=1 / cos(θ)=1 / max(abs(Nx),abs(Ny),abs(Nz));

[0086] Where abs is the absolute value and A is the surface area corresponding to the voxel.

[0087] S5. Process the surface area grayscale value of the surface voxels to reduce the resolution of the surface voxels until the resolution is 1, and obtain the surface area of ​​the ROI region.

[0088] The grayscale values ​​of the surface voxels are processed to prepare a voxel mesh with a resolution half that of the current binary voxel model. For the voxel with coordinates (i,j,k), the surface areas of the original voxel model coordinates (2*i,2*j,2*k), (2*i+1,2*j,2*k), (2*i,2*j+1,2*k), (2*i+1,2*j+1,2*k), (2*i,2*j,2*k+1), (2*i+1,2*j,2*k+1), (2*i,2*j+1,2*k+1) are summed and recorded in the current voxel. At the same time, the voxel model resolution is reduced by half. This step is repeated until the voxel resolution is 1. The grayscale value of the stored surface area is the surface area of ​​the ROI region.

[0089] Therefore, the present invention employs a sub-voxel-precision rapid estimation method for the surface area of ​​the region of interest in the above-mentioned voxel model, which can adapt to complex structures and reduce dependence on preset models, thereby providing reliable support for high-precision three-dimensional analysis.

[0090] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit them. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical solutions of the present invention, and these modifications or equivalent substitutions cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A method for fast estimation of sub-voxel precision of surface area of a region of interest in a voxel model, characterized in that, The method comprises the following steps: S1, dividing the ROI region and setting a gray threshold to determine the ROI region voxel; S2, extracting the surface voxel of the ROI region, searching the 6 adjacent voxels of the ROI region voxel, when there is a non-ROI region voxel, the current voxel is marked as a surface voxel, and the balance center position of the surface voxel is initialized as the absolute center position of the surface voxel; S3, using the iterative method to define the balance center position of the surface voxel sub-voxel accuracy; the content is as follows: For the surface voxel, query the 3x3x3 adjacent voxels around it, if there is a surface voxel, extract the balance center position of the surface voxel, and perform a weighted average on the obtained adjacent voxel balance position, and update the balance center position of the current surface voxel, and the balance center position is constrained to be not more than 1 voxel unit different from the absolute center position, and the specific formula is as follows: ； ； wherein, is the absolute center position of the current surface voxel, t is the iteration step number, is the weight value that the voxel is a surface voxel, i is the index of the neighboring voxel; If the iterative process of the surface voxel has unstable details, the average value of the old balance center position and the new balance center position is used as the new balance center position each time: ； The above is the calculation process of the balance center position, and the balance center gradually converges to a stable value by repeating the step for multiple times; For complex curved surfaces, there are the following three strategies for obtaining the balance center position: a. Limit the number of iterations to 32 steps or less; b. After each iteration of step Decay is performed, with the formula as follows: wherein is the attenuation coefficient; c. Calculate the balance center position offset of the adjacent surface voxel, and adjust the balance center position of the current voxel in turn, and the formula is: ； ； wherein is an offset coefficient; S4, after obtaining the converged surface voxel balance center position, calculate the surface area value corresponding to the surface voxel; S5, processing the surface area gray value of the surface voxel, reducing the resolution of the surface voxel until the resolution is 1, and obtaining the surface area of the ROI region.

2. The method of claim 1, wherein: In S1, the threshold gray value g is set, the voxels are traversed, and the voxels with a gray value greater than g are defined as ROI region voxels.

3. The method of claim 2, wherein: In S2, the surface voxel is marked as 1, and the rest is 0, to obtain a binary voxel model of the surface voxel.

4. The method of claim 3, wherein, The content of S4 is as follows: Within the range of 3x3x3 around the voxel, the surface normal is calculated according to the balance center position of the adjacent voxel, and then the surface area occupied by the current voxel is calculated based on the normal direction and the balance center position of the voxel, and the surface area of the surface voxel is stored in the current voxel model as a gray value; The calculation of the normal does not use the absolute center position of the voxel, but uses the balance center position obtained by iteration, and the Z direction is taken as the reference, and the specific formula is as follows: ； ； ； wherein , , are the components of the C vector in the three directions; i, j, k are the voxel index variables for the three directions X, Y, Z, in the range -1, 0, 1 for 3x3x3. According to the normal N, the minimum angle between the normal and the X, Y, Z three planes is calculated , the surface area corresponding to the voxel is recorded , and the formula is as follows: ； ； wherein A is the surface area of the voxel.

5. The method of claim 4, wherein, The content of S5 is as follows: The grayscale values ​​of the binary voxel model of the surface voxels are processed to prepare a voxel mesh with half the resolution of the current binary voxel model. This mesh is then applied to coordinates... The voxels, taking the original voxel model coordinates as The surface areas are summed and recorded in the current voxel. At the same time, the voxel model resolution is reduced by half. This step is repeated until the voxel resolution is 1. The gray value of the stored surface area is the surface area of ​​the ROI region.

Images