CN1711968A

CN1711968A - Rapid progressive three-dimensional reconstructing method of CT image from direct volume rendering

Info

Publication number: CN1711968A
Application number: CNA2005100427347A
Authority: CN
Inventors: 朱虹; 邓杰航; 马湘旺; 楚勇; 张波; 刘磊; 曹艳艳; 叶丽丽; 王子玉
Original assignee: Xian University of Technology
Current assignee: Xian University of Technology
Priority date: 2005-05-26
Filing date: 2005-05-26
Publication date: 2005-12-28
Anticipated expiration: 2025-05-26
Also published as: CN100348158C

Abstract

A high-speed progressive 3D reconstruction method for CT images includes such steps as dividing the data about the 2D laminagraphic images of an organ, interpolating to generate isotropic 3D data body, calculating the sizes of target data body after rotation, coordinate transplanting of said data body, initializing the projected data region and projection template, rotational assigning, projecting and reinforcing the projected result.

Description

Rapid progressive direct volume rendering three-dimensional reconstruction method of CT image

Technical Field

The invention belongs to the technical field of three-dimensional reconstruction visualization of CT images, relates to a three-dimensional reconstruction method of a CT image, and particularly relates to a three-dimensional reconstruction method of rapid progressive direct volume rendering of a CT image.

Background

For the three-dimensional reconstruction of the CT tomographic image, the direct volume rendering method can more truly represent the three-dimensional effect compared with the surface rendering method, the direct volume rendering takes the voxel as a basic unit, and the three-dimensional volume image is directly generated by the two-dimensional tomographic image data set.

Disclosure of Invention

In order to solve the defects of large volume rendering calculation amount and low speed of the prior art, the invention aims to provide a fast progressive direct volume rendering method, which greatly reduces the calculation amount and obviously improves the speed when the method is used for carrying out three-dimensional reconstruction on a CT image.

The invention adopts the technical scheme that a CT image fast progressive direct volume rendering three-dimensional reconstruction method comprises the steps of firstly collecting two-dimensional sectional image data of an organ by CT equipment, carrying out image segmentation in the two-dimensional sectional image to obtain the two-dimensional sectional CT image data of the organ needing three-dimensional reconstruction and the organ associated with the organ, then carrying out three-dimensional data interpolation on the two-dimensional sectional images of different organs to form an isotropic three-dimensional data volume, and further carrying out the following steps,

1) calculating the size of the rotated target data volume

Calculating the size of the rotated target data volume according to the isotropic three-dimensional data volume value;

2) data volume coordinate translation

Translating the original data volume and the rotated coordinate origin of the target data volume to respective geometric centers;

3) initialization of x-axis and y-axis coordinates

For the x-axis, determining the first projection point on the x-axis, and setting other projection points with the interval of m, namely x (k) is 0, x (k) is 8, x (k) is 16, …, m is the power of 2,

setting a first pair of reconstructed images on a y axis, wherein when the number of total tomograms is an odd number, y0 is n, y1 is-n, and when the number of total tomograms is an even number, y0 is n-1, and y1 is-n;

4) initializing a projection data field and a projection template

Initializing the projection data area to a minimum CT value, and initializing the projection template to '0';

5) performing rotation assignment on data body

Traversing the data body of a first quadrant for the data body of original data at the same time, two pieces of two-dimensional fault data containing eight symmetrical points in the rotated data body, a projection plane data body and a projection template, calculating corresponding coordinates of the data of the rotated data body in the original data body, assigning the data to the voxel value of the corresponding coordinate of the rotated data body if the calculated coordinates do not exceed the original data space and the voxel value of the corresponding coordinate is not empty, correspondingly completing the traversal of the data of other seven quadrants after the voxel of the first quadrant is traversed, and setting the corresponding position of the projection template to be '1' when a non-empty voxel value appears in the Z-axis direction, or not changing;

6) projection processing

Performing ray tracing on the two obtained data along the Z axis, when the template value of the current point is 0, not processing, and when the template value is 1, performing ray tracing, and filling the value obtained by ray tracing to a corresponding pixel point in the projected image;

7) after the two images are processed, the plane is shifted by one unit along the Y axis, namely Y0 is Y0-1, Y1 is Y1+1, and the other two images in the data volume are taken to be processed in the 5 th image) and the 6 th image until Y0 is 0;

8) performing projection display

Refreshing a screen, and displaying a projection result;

9) then let x (k +1) ═ x (k) +1, return to above-mentioned 4) and repeat the process, until all points on the x axis have traversed, the projection of fast progressive three-dimensional reconstruction ends.

The present invention is also characterized in that the ray tracing along the Z-axis is performed as follows,

a, emitting rays from a projection surface, and traversing voxel points in the ray direction;

b, if the value of the coordinate point corresponding to the projection template is '0', setting the pixel as a background color, turning to a, and otherwise, turning to c;

c, recording a subscript for finding a first pixel point containing information;

d, calculating the opacity from the found subscript, and performing illumination mapping and depth modulation;

e, drawing the pixel of the point;

f, turning to a until all data processing is completed.

The method of the invention reduces the calculation amount of the data volume when rotating according to the inherent geometric symmetry of the three-dimensional data volume; and then, according to the fact that a large number of redundant voxels do not contribute to the projection value in the three-dimensional data volume, data which do not contribute to the projection are identified, and background pixel values are directly assigned to the voxels which do not contribute to the projection without processing during ray tracing iteration. And according to the ray tracing principle, the pixel value on the projection surface is only related to the voxel passed by the ray emitted from the point, and the voxel can be subjected to block iteration to realize progressive direct volume rendering. Therefore, the method greatly reduces the steps of three-dimensional reconstruction and improves the speed of three-dimensional reconstruction.

Drawings

FIG. 1 is a flow chart of the method of the present invention;

FIG. 2 is a schematic diagram of rotation, where a is the original data volume and b is the rotated data volume;

FIG. 3 is a three-dimensional data volume voxel value redundancy illustration;

fig. 4 is a fast progressive direct gray scale volume rendering effect diagram of the present invention, where a-e are the progressive process of lung imaging and f is the result of bone imaging.

Detailed Description

The invention is explained in further detail below by means of figures and examples.

The basic concept of the invention is that in the process of direct volume rendering, after obtaining the isotropic three-dimensional data volume after interpolation, the geometric symmetry of the three-dimensional data volume is utilized, and simple addition and subtraction method is used for replacing multiplication operation in rotation, so that the rotation calculation speed is accelerated; secondly, a large number of points with empty voxel values exist in the three-dimensional data volume, the points do not contribute to the final projection surface, the points can be omitted in the projection iteration process, and the iteration calculation amount is reduced; thirdly, according to the ray tracing principle, the pixel value on the projection surface is only related to the voxel passed by the ray emitted from the point, so that the voxel can be subjected to block iteration, progressive direct volume rendering is realized, and finally a three-dimensional projection image of the required organ is generated.

1. Geometric symmetry of three-dimensional data volume

Assuming that the resolution of each two-dimensional tomographic image is 512 × 512, a square is formed on a two-dimensional plane and is a centrosymmetric pattern, so that three-dimensional volume data forms a rectangular parallelepiped after three-dimensional reconstruction of the two-dimensional tomographic image. According to the solid geometry knowledge, the cuboid is an axisymmetric figure. As shown in fig. 2, in the rectangular parallelepiped ABCDEFGH, the point A, B, C, D is symmetric about the axes X and Y, and the points A, B, C, D and E, F, G, H are symmetric about the center O, respectively, with the geometric center O of the rectangular parallelepiped as the origin.

Assuming that the three-dimensional coordinates of point B are (X, Y, Z) and according to symmetry, A, C, D, E, F, G, H have coordinates (X, -Y, Z), (-X, Y, Z), (-X, -Y, Z), (X, -Y, -Z), (X, Y, -Z), (-X, Y, -Z), (-X, -Y, -Z), and (X, -Y, -Z), assuming that the rotation matrix from FIG. 2a to FIG. 2B is T, the volume data is rotated α degrees around the X axis, β around the Y axis, and γ degrees around the Z axis, and

T＝T_x·T_y·T_z

wherein,

namely, it is

Let A₁₁＝cosβcosγ，A₁₂＝-sinαsinβcosγ-cosαsinγ，

A₁₃＝-cosαsinβcosγ+sinαsinγ，B₁₁＝cosβsinγ，

B₁₂＝-sinαsinβsinγ+cosαcosγ，B₁₃＝-cosαsinβsinγ-sinαcosγ，

C₁₁＝sinβ，C₁₂＝sinαcosβ，C₁₃＝cosαcosβ，

Then

T = [\begin{matrix} A_{11} & B_{11} & C_{11} \\ A_{12} & B_{12} & C_{12} \\ A_{13} & B_{13} & C_{13} \end{matrix}] - - - (3)

The reconstructed point A 'transformed by the rotation matrix T is denoted by A', and the coordinates of A 'are denoted by (x'_A，y′_A，z′_A) Then, then

Because of [ x_A y_A z_A]＝[x -y z]So the system of equations can be obtained

Similarly, point B, C, D is rotated to be B '(x'_B，y′_B，z′_B)、C′(x′_C，y′_C，z′_C)、D′(x′_D，y′_D，z′_D) The system of equations, have

Obviously, point a 'does not have a linear relationship with point B', point C ', and point D', respectively.

And because this is four points on the three-dimensional space, these four points are linearly related, and of the points a ', B', C ', D', the coordinates of one point can be linearly represented by three more points, k₁，k₂，k₃Are respectively real numbers other than zero, then

[x_A′，y_A′，z_A′]＝k₁·[x_B′，y_B′，z_B′]+k₂·[x_C′，y_C′，z_C′]+k₃·[x_D′，y_D′，z_D′](9) Substituting the above formulae (5), (6), (7) and (8) can give:

\{\begin{matrix} k_{1} - k_{2} - k_{3} = 1 \\ k_{1} + k_{2} - k_{3} = - 1 \\ k_{1} + k_{2} + k_{3} = 1 \end{matrix} - - - (10)

solving the equation set to k₁＝1，k₂＝-1，k ₃1. The formula (9) is abbreviated as

A′＝B′-C′+D′ (11)

In addition, since the points A, B, C, D are respectively symmetrical to the points G, H, E, F about the center O, the

\{\begin{matrix} A = - G \\ B = - H \\ C = - E \\ D = - F \end{matrix} - - - (12)

Rotational transformation of E

E′＝E·T (13)

By substituting formula (12) for formula (13)

E′＝E·T＝-C·T＝-C′ (14)

By the same way, obtain

F′＝F·T＝-D·T＝-D′ (15)

G′＝G·T＝-A·T＝-A′ (16)

H′＝H·T＝-B·T＝-B′ (17)

Therefore, the eight points A, B, C, D, E, F, G, H of the data volume are subjected to rotation transformation to obtain the points a ', B', C ', D', E ', F', G ', H', without eight rotation transformations, only three rotation transformations and three additions are needed, and four negations are needed.

2. Redundancy of three-dimensional data volumes

The three-dimensional data volume has great redundancy and a plurality of empty voxels with the voxel value of 0 according to the iterative formula of ray tracing

sum_n＝α_n·v_n(x)+(1-α_n)·sum_n-1 (18)

Wherein, Sum_nFor the iteration value (projection value) of the first n sampling points, alpha_nOpacity of voxel value, v, for the nth sample point_n(x) For the depth-modulated voxel value of the current voxel, Sum₀＝v₀(x)。

As shown in fig. 3, in the two-dimensional tomographic image, except for the region of interest (lung region), there are null pixels, and these null pixels are necessarily null pixels even after the three-dimensional data volume is formed by interpolation of the three-dimensional data, and if these pixels are also subjected to iterative processing without resolution, the computation time is wasted. It can be seen that the voxel values of the three-dimensional data volume have a large degree of redundancy.

3. According to the relation of pixel value and voxel

According to the ray tracing principle, the pixel value on the projection surface is only related to the voxel passed by the ray emitted from the point, so that the voxel can be subjected to block iteration to realize progressive direct volume rendering.

As shown in fig. 1, firstly, two-dimensional tomographic image data of a certain organ is acquired by a CT device, image segmentation is performed in the two-dimensional tomographic image to obtain two-dimensional tomographic CT image data of the certain organ and organs related to the certain organ which need to be three-dimensionally reconstructed, for example, data of parts such as skin, muscle, bone, lung, mediastinum and the like are respectively obtained for a chest film, then three-dimensional data interpolation is performed on the two-dimensional tomographic images of different organs to form an isotropic three-dimensional data body, the x-axis direction of the data body is a horizontal direction, the y-axis direction is a position direction of the two-dimensional tomographic image, and the z-axis direction is a direction far away from a projection plane, and then the following steps are performed according to the method of the present invention:

1) calculating the size of the rotated target data volume

Calculating a rotated target data volume from an isotropic three-dimensional data volume, wherein the length, width and height of the original data volume are M, N, H respectively, the rotation matrix is T, T is required to be a reversible matrix, the inverse matrix of T is T-1, and the length, width and height of the rotated target data volume are M ', N ' and H ', so that

[M′，N′，H′]＝[M，N，H]·abs(T)

Where abs (T) represents the absolute value of each element within matrix T.

2) Data volume coordinate translation

The origin of coordinates of the original data volume and the rotated target data volume are translated to the respective geometric centers, and are represented as coordinates of the data volume as shown in fig. 2.

3) Initialization of x-axis and y-axis coordinates

In order to increase the display speed, the invention adopts a mode of interval display every time, namely, projection points with the interval of m (m is a power of 2, as shown in fig. 2, m is 8) are arranged on an x axis, only one point is displayed in each iteration, and after eight times of iteration refreshing display, the complete projection display of the target object can be completed. For this purpose, during initialization, first, the first number of projection points on the x-axis is determined, i.e., x (k) is 0, x (k) is 8, and x (k) is 16, …;

in the y-axis direction, after the coordinate origin translation, if the total number of tomographic images after the interpolation is an odd number (set to 2n +1), the value of y is an integer between [ -n, n ]. If the total tomographic image number is an even number (set to 2n), y takes an integer value between [ -n, n-1 ]. In order to reduce the consumption of the memory and improve the processing speed, the invention only projects at most two images at each time during the three-dimensional reconstruction. For this reason, when initializing, the first reconstructed image is set to y0 ═ n, y1 ═ n (when the total number of tomograms is an odd number), or y0 ═ n-1, and y1 ═ n (when the total number of tomograms is an even number).

4) Initializing a projection data field and a projection template

The projection data area is used for storing the projected image data and initializing the projected image data to a minimum CT value (which is a minimum value provided by a CT equipment manufacturer, and is-2048 in general).

The projection template is used for storing '0' and '1' values, wherein '1' indicates that the ray emitted from the corresponding position of the projection surface meets a voxel point with a volume value which is not empty after passing through a rotated voxel, and otherwise, the ray is '0', and the projection template is initialized to be '0';

5) carrying out rotation assignment and projection processing on data volume

In consideration of the problem of processing speed, the data volume is converted from the rotated data volume to the data volume before rotation to realize rotation operation, and the assignment of the data volume is realized by adopting a nearest neighbor interpolation method.

As shown in fig. 2, in order to save memory, only two pieces of two-dimensional tomographic data including eight symmetric points in the original data volume and the rotated data volume, i.e., the plane a 'B' C 'D' and the plane E 'F' G 'H', as well as the projection plane data volume and the projection template, exist in the memory at the same time. Assuming that the coordinates of the first quadrant point B 'are (x, y, z), the coordinates of the second and third quadrants C', D 'are (-x, y, z), (-x, -y, z), and assuming that the coordinates of the rotated data volume at points a', B ', C', D ', E', F ', G', H 'correspond to the coordinates of the original data volume at points a, B, C, D, E, F, G, H, the following processes are performed for points B', C ', D':

B＝B′·T^-1，C＝C′·T^-1，D＝D′·T^-1

since the rotation is reversible, it follows from the above formula

\{\begin{matrix} A = B - C + D \\ E = - C \\ F = - D \\ G = - A \\ H = - B \end{matrix}

The coordinates of the corresponding points of the other seven quadrants can be obtained by only traversing the data volume of the first quadrant.

After the corresponding coordinates of the data of the rotated data volume in the original data volume are calculated, if the calculated coordinates do not exceed the original data space and the voxel value of the corresponding coordinates is not null, the point data is assigned to the voxel value of the corresponding coordinates of the rotated data volume, and the data is searched at intervals of m data on the X axis during traversal, wherein m is usually the integral power of 2.

When the voxel of the first quadrant of the plane A 'B' C 'D' is traversed, the data of other seven quadrants of the two corresponding pieces of data are also traversed, and in the process of traversing assignment, when a non-empty voxel value appears in the Z-axis direction, the corresponding (x, y) position of the projection template is set as '1', otherwise, the initial value is still unchanged as '0'.

6) Projection processing

And performing ray tracing on the two obtained data along the Z axis, and determining whether the ray tracing is performed or not according to the value of the projection template. If the template value of the current point is 0, it indicates that the light is not shielded by the target object, and no processing is performed, if the template value is 1, the light tracking is performed, and the value obtained by the light tracking is filled in the corresponding pixel point in the projected image.

Ray tracing iteration is carried out on the two pieces of data along the Z axis, and the specific steps are as follows:

a, emitting rays from a projection surface, and traversing voxel points in the ray direction (Z-axis direction);

e, drawing the pixel of the point;

f, turning to a until all data processing is completed.

7) After the two images are processed, the plane is moved by one unit along the Y axis (i.e. Y0 is Y0-1, Y1 is Y1+1), and the other two images in the data volume are taken for processing of 5 th and 6 th images until Y0 is 0;

8) performing projection display

After the processing of 4), 5), 6), 7), refreshing a screen, wherein the refreshing aims to overlap the current projection result and the previous projection result and then display the projection result;

9) then let x (k +1) ═ x (k) +1, return to 4) repeat the process until after all points on the x-axis have traversed, the fast progressive projection ends.

The following is a comparison of the results of the tests:

the test hardware is configured to: CPU, p42.4g; a memory 768M; a display memory, a shadowless tyrant 128M; hard disk, Xijie 60G, 7200rpm.

30 two-dimensional tomographic images (88 images after interpolation) are read in, and the space occupied by the two-dimensional tomographic images is 512 × 512 × 88, which is 230,686,762 bytes of data. The time spent in the projection of the rotation (X: 50 degrees, Y: 0 degrees, Z: 340 degrees) is:

the traditional algorithm is as follows: 41 seconds, the algorithm of the method: for 12 seconds.

Generally, volume rendering is performed on about 80 frames of two-dimensional tomographic images by the method, and under the above hardware conditions, about 15 seconds are required.

As can be seen from the comparative experiment result, the method of the invention can greatly reduce the time of rotation and projection, thereby accelerating the three-dimensional reconstruction of the CT image.

Fig. 4 shows an example of the processing by the method, which shows the process and effect of three-dimensional imaging of lung tissue in a chest film and the effect of three-dimensional imaging of bone tissue at the same position. In this example, the coordinate interval on the x-axis is set to m-8, and fig. a is the superimposed projection effect when x-3, that is, fig. a is a value 0, 1, 2, 3 on the x-axis; as a result of performing four iterative projections 0+8, 1+8, 2+8, 3+8, …, when x is 7, all eight iterative projections of 0, 1, 2, 3, 4, 5, 6,7, 0+8, 1+8, 2+8, 3+8, 4+8, 5+8, 6+8, 7+8, … are completed, because m is 8, and as shown in fig. e, a clear three-dimensional display of the lung organ tissue can be obtained.

The figure f is the result of three-dimensional imaging of bone tissues of the same group of CT tomograms by adopting the method of the invention. The time required from the reading of the raw data to the obtaining of the graph e is about 12 seconds, and the time required from the reading of the raw data to the obtaining of the graph f is about 13 seconds.

Claims

The fast progressive direct volume drawing three-dimensional reconstruction method of CT image includes the steps of firstly, collecting two-dimensional tomographic image data of organs by CT equipment, carrying out image segmentation in the two-dimensional tomographic image to obtain two-dimensional tomographic image data of organs needing three-dimensional reconstruction and organs related to the organs, and then carrying out three-dimensional data interpolation on the two-dimensional tomographic images of different organs to form isotropic three-dimensional data volume, and is characterized in that: the method is also carried out as follows,

1) calculating the size of the rotated target data volume

Calculating the size of the rotated target data volume according to the isotropic three-dimensional data volume value;

2) data volume coordinate translation

Translating the original data volume and the rotated coordinate origin of the target data volume to respective geometric centers;

3) initialization of x-axis and y-axis coordinates

For the x-axis, determining the first projection point on the x-axis, and setting other projection points with the interval of m, namely x (k) is 0, x (k) is 8, x (k) is 16, …, m is the power of 2,

setting a first pair of reconstructed images on a y axis, wherein when the number of total tomograms is an odd number, y0 is n, y1 is-n, and when the number of total tomograms is an even number, y0 is n-1, and y1 is-n;

4) initializing a projection data field and a projection template

Initializing the projection data area to a minimum CT value, and initializing the projection template to '0';

5) performing rotation assignment on data body

Traversing the data body of a first quadrant for the data body of original data at the same time, two pieces of two-dimensional fault data containing eight symmetrical points in the rotated data body, a projection plane data body and a projection template, calculating corresponding coordinates of the data of the rotated data body in the original data body, assigning the data to the voxel value of the corresponding coordinate of the rotated data body if the calculated coordinates do not exceed the original data space and the voxel value of the corresponding coordinate is not empty, correspondingly completing the traversal of the data of other seven quadrants after the voxel of the first quadrant is traversed, and setting the corresponding position of the projection template to be '1' when a non-empty voxel value appears in the Z-axis direction, or not changing;

6) projection processing

Performing ray tracing on the two obtained data along the Z axis, when the template value of the current point is 0, not processing, and when the template value is 1, performing ray tracing, and filling the value obtained by ray tracing to a corresponding pixel point in the projected image;

7) after the two images are processed, the plane is shifted by one unit along the Y axis, namely Y0 is equal to Y0-1, Y1 is equal to Y1+1, and the other two images in the data volume are taken to be processed in the 5 th image and the 6 th image until Y0 is equal to 0;

8) performing projection display

Refreshing a screen, and displaying a projection result;

9) then let x (k +1) ═ x (k) +1, return to above-mentioned 4) and repeat the process, until all points on the x axis have traversed, the projection of fast progressive three-dimensional reconstruction ends.
2. The three-dimensional reconstruction method as claimed in claim 1, wherein the size of the rotated target data volume is calculated by assuming that the length, width and height of the original data volume are M, N, H respectively, the rotation matrix is T, T is a reversible matrix, the inverse of T is T-1, and the length, width and height of the rotated data volume are M ', N ', H ' respectively, and then

[M′，N′，H′]＝[M，N，H]·abs(T)

Where abs (T) represents the absolute value of each element within matrix T.
3. A method of three-dimensional reconstruction as claimed in claim 1 wherein the data is searched on the X-axis every m data, m being an integer power of 2, during traversal.
4. The three-dimensional reconstruction method of claim 1, wherein the coordinates of the data of the rotated data volume in the original data volume are calculated, assuming that the coordinates of the first quadrant point B 'are (x, y, z), and the coordinates of the second and third quadrants C', D 'are (-x, y, z), -x, -y, z), and assuming that the coordinates of the rotated data volume points a', B ', C', D ', E', F ', G', H 'correspond to the coordinates of the original data volume points a, B, C, D, E, F, G, H, and the points B', C ', D' are processed as follows: b ═ B'. T^-1，C＝C′·T^-1，D＝D′·T^-1

Because the rotation is reversible, it is obtained

$\{\begin{matrix} A = B - C + D \\ E = - C \\ F = - D \\ G = - A \\ H = - B \end{matrix} .$
5. A three-dimensional reconstruction method as set forth in claim 1, wherein the step of ray tracing along the Z-axis is as follows,

a, emitting rays from a projection surface, and traversing voxel points in the ray direction;

b, if the value of the coordinate point corresponding to the projection template is '0', setting the pixel as a background color, turning to a, and otherwise, turning to c;

c, recording a subscript for finding a first pixel point containing information;

d, calculating the opacity from the found subscript, and performing illumination mapping and depth modulation;

e, drawing the pixel of the point;

f, turning to a until all data processing is completed.