JPH04363775A - Method for three-dimensional display at high speed - Google Patents

Abstract

PURPOSE:To display a three-dimensional picture with quality level being proportionate to a linear interpolation within processing time being proportionate to be most adjacent and approximate. CONSTITUTION:A three-dimensional picture shaded in such a way that a density gradient vector is calculated by peripheral density changing quantity concerning the individual voxcel of three-dimensional picture data, a luminance value in the individual vox is obtained by the density gradient vector and the incident light vector to the voxcel and the luminance value is projection-processed in a two-dimensional plane is generated in a three-dimensional picture displaying method. In the method, density changing quantity in the respective reference axis directions in the coordinate system of pre-rotation is obtained at every voxcel, the density changing quantity in the respective reference axis directions in the coordinate system of post-rotation is estimated and the vector where the respective density changing quantity becomes a component is used as the density gradient vector at the time of displaying the three-dimensional picture from the optional direction.

Description

[Detailed description of the invention]

0001

[Industrial Application Field] The present invention relates to a high-speed three-dimensional display method, and in particular, to display three-dimensional images useful for diagnosis from human body three-dimensional data obtained by tomography equipment such as X-ray CT and MRI. Regarding the method.

0002

2. Description of the Related Art In medical images, since it is difficult to model the shape of an organ, shading processing is usually performed on a voxel (pixel) basis. In shading processing, pixel values on the projection plane are determined as the inner product of the incident light vector and the normal vector calculated for each voxel (pixel) on the object surface. There are two typical three-dimensional display methods:

[0003] Z-buffer shading calculates the normal vector of the surface based on the anteroposterior relationship of voxels (pixels) existing on the object surface, and performs three-dimensional display. This method is fast 3
Although it has the advantage of being able to obtain a dimensional image, normal vectors are quantized depending on the size of the voxel itself, so normal data (resolution of 128 to 256 pixels/column) cannot provide satisfactory image quality.

On the other hand, gray level gradient shading calculates the normal vector of the object surface from the density difference between adjacent pixels. In this method, the normal vector is obtained as a continuous quantity, so a high-quality three-dimensional image with excellent density resolution can be obtained.

[0005] Volume rendering (I.E.E.
ECG & A May 1988 issue, pages 29-37 (IEEE Computer
Graphics & Applications) also uses the concept of gray level gradient shading.

[0006] A gray level gradient is a vector representing the gradient direction of density and the amount of gradient. This is the voxel (
It is obtained by taking the difference between the density values of pixels (pixels). In gray level gradient shading, this gray level gradient is defined as a normal vector on a microplane virtually set for each voxel (pixel), and the inner product of this normal vector and the incident light vector is calculated for each surface. Calculate the brightness.

[0007] Therefore, in order to improve the image quality and fidelity of a three-dimensional image by gray level gradient shading, it is necessary to accurately calculate the gray level gradient.

[0008] When the rotation process is performed, the alignment of the rotated lattice and the original lattice is shifted. Therefore, it is necessary to interpolate and find the density values of the shifted grid points. When nearest neighbor interpolation is used for interpolation processing, the value of the gray level gradient periodically increases or decreases due to variations in the distance between voxels (pixels) from which density differences are taken. This appears as striped artifacts on the three-dimensional image, significantly reducing the image quality of the three-dimensional image. Therefore,
Previously I.E.E.C.G.&.
May 1988 issue, pp. 29-37 (IEEE
Computer Graphics & Applications
ations), when performing three-dimensional display from any direction, linear interpolation (tri-linear
interpolation) was used. but,
Linear interpolation requires much more access to data and a much larger amount of calculations than nearest-neighbor interpolation, which has been a bottleneck in speeding up display processing.

[0009] Here, an outline of the three-dimensional display processing when using the conventional interpolation method using nearest neighbor approximation and the factors that cause artifacts to occur will be described with reference to the reference diagram of FIG.

[0010] The gray level gradient is calculated from the density difference between points adjacent to the point (reference point) for which the gray level gradient is to be determined. Note that for the sake of simplicity, the description will be limited to one direction (vertical direction) component of the gray level gradient.

FIG. 8(a) shows a case where rotation processing is not performed. At this time, since adjacent points are always on the grid, the distance between adjacent points is always equal. Therefore, the amount of concentration gradient (=density difference/adjacent point distance) can always be accurately determined.

When the rotation is performed, as shown in FIG. 8(b),
The alignment of the grid before and after rotation is misaligned. If the adjacent points are determined by nearest neighbor approximation, the distance and direction between the adjacent points will vary as shown in FIG. 8(c), so that the concentration gradient cannot be determined accurately. Moreover, since this is caused by interference between equally spaced grids, the distance between adjacent points increases and decreases periodically, and the result appears as a striped pattern on the image.

Linear interpolation refers to voxel (pixel) values around the point of interest and performs weighted addition to these values to determine the value of the point of interest. Therefore, no stripe artifacts occur. However, in linear interpolation calculations, real-time processing is difficult because each point requires data access to eight surrounding points (in the case of three dimensions) and weighting calculations.

As described above, when the nearest neighbor approximation is used, striped artifacts occur and the image quality is significantly degraded. On the other hand, when linear interpolation is used, there is no problem in terms of image quality, but the time required to create a three-dimensional image becomes a problem.

[0015]

SUMMARY OF THE INVENTION An object of the present invention is to make it possible to display a three-dimensional image having an image quality comparable to linear interpolation in a processing time comparable to nearest neighbor approximation.

[0016]

[Means for Solving the Problems] The present invention calculates a density gradient vector for each voxel of three-dimensional image data from the amount of surrounding density change, and calculates each voxel from the density gradient vector and the incident light vector to the voxel. In a 3D image display method that creates a shaded 3D image by determining the brightness value at a voxel and projecting the brightness value onto a 2D plane, when displaying a 3D image from any direction. , find the amount of density change in each reference axis direction in the coordinate system before rotation for each voxel, estimate the amount of density change in each reference axis direction in the coordinate system after rotation from this amount of density change,
The present invention is characterized in that a vector having each density change amount as a component is used as the density gradient vector.

[0017]

[Operation] In the present invention, nearest-neighbor interpolation is used as a basic method for rotation calculation. Thereby, the time required to create a three-dimensional image is no longer a problem. The difference between the present invention and the conventional method is that the conventional method directly obtains the value of the point used to calculate the amount of density change by nearest neighbor approximation, whereas in the present invention, the reference point of the gray level gradient is first determined. Obtained by nearest neighbor approximation, and from the values of points adjacent to this reference point,
This is where you find the amount of change in concentration. At this time, since the adjacent points become points on the grid of the coordinate system before rotation, it is possible to accurately calculate the gray level gradient in the coordinate system before rotation. The gray level gradient regarding the rotated coordinate system is obtained by subjecting the gray level gradient obtained here to linear transformation processing using the same transformation matrix used for the rotation of the three-dimensional data. By using the method of the present invention, a high-quality three-dimensional display image without striped artifacts can be obtained.

The concept of the method of the present invention will be described again using the reference diagram of FIG. As shown in FIG. 8(d), the reference point for the concentration gradient is placed at the point (black square) on the grid in the coordinate system before rotation that is closest to the original reference point (shown by the black circle). By resetting the reference point to a point on the grid in this way, the adjacent points also become points (white squares) on the grid, as shown in FIG. 8(e). As a result, the distance between adjacent points is always constant, making it possible to accurately determine the amount of concentration gradient.

However, the amount of concentration gradient determined here is determined with respect to the reference axis direction of the coordinate system before rotation, and not with respect to the coordinate system after rotation, which should be determined. Therefore, as shown in FIG. 8(e), the concentration gradient in the direction to be determined is calculated from the concentration gradient determined in the vertical and horizontal directions and the vector in the rotational direction. The concentration gradient vector (x, y, z) after rotation is obtained from the matrix A representing coordinate rotation and the concentration gradient vector (dx, dy, dz) before rotation [Equation 2]
It can be calculated using the formula (see formula 2 in FIG. 9).

The calculation point of the rotational density gradient vector thus obtained is the point represented by the white circle in FIG. 8(f), and the distance and direction between the two points are always equal.

[0021] Such a gray level gradient calculation method makes it possible to obtain a three-dimensional image in which striped artifacts are suppressed at approximately the same speed as the conventional method using nearest neighbor approximation. That is, a three-dimensional image having an image quality comparable to linear interpolation can be created faster than when linear interpolation is used.

[0022]

EXAMPLES The present invention will be explained based on the following examples. In this embodiment, a three-dimensional display of a skin region is performed from a plurality of head tomographic data captured using MRI.

FIG. 1 shows a flowchart relating to a gray level gradient calculation method, which is the essential part of the present invention. <Step 001> The reference point A of the gray level gradient is set to the nearest point A' of the original grid by nearest neighbor interpolation, and A' is hereinafter used as the reference point of the gray level gradient. <Step 002> Adjacent points B', C', D', to the reference point A' on the grid
E' is also a point on the grid. The density change amount dy' in the original coordinate system is determined as the difference between the density values of these adjacent points B' and C', and the density change amount dx in the original coordinate system is calculated as the difference between the density values of these adjacent points D' and E'. 'Seek.'<Step003> Since the amount of concentration change obtained in step 002 is with respect to the coordinate system of the original lattice before rotation, the inclination angle θ of the lattice after rotation with respect to the original lattice and the amount of concentration change dx', dy'
Therefore, the amount of concentration change in the lattice after rotation dx", dy"
seek.

FIG. 2 is a system configuration diagram of an apparatus to which the present invention is applied. This device includes an I/O device 10, a three-dimensional memory 11, an image memory 12, a CPU 13, and a CRT 14.
, console 15. First, the role of each component will be explained.

Head tomographic image data obtained by MRI is transferred to a three-dimensional memory 11 through an I/O device 10. The sent multiple tomographic image data are stored in 3D memory 1.
1 while maintaining their mutual positional relationship. As a result, the fault line data is configured as three-dimensional data on the three-dimensional memory 11.

The CPU 13 refers to the three-dimensional data on the three-dimensional memory 11 and performs processing such as rotation and three-dimensional image creation. This process is performed by the console 15 connected to the CPU 13.
This is done based on instructions from. The processing results are stored on the image memory 12 and displayed on the CRT 14.

FIG. 3 shows a flowchart of the processing procedure for three-dimensional image creation in the CPU 13 of FIG. 2. 2nd below
An embodiment will be described using the apparatus configuration diagram shown in the figure and the flowchart shown in FIG. <Step 100> I/O multiple tomographic data captured by MRI
The data is stored on the three-dimensional memory 11 through the device 10. Each tomographic data is configured as three-dimensional data on the three-dimensional memory 11. <Step 101> The operator inputs the rotation direction of the three-dimensional image from the console 15. Examples of this input method include a method of numerically inputting the rotation angle regarding the X, Y, and Z axes, and a method of intuitively specifying the rotation angle using three-dimensional graphics and a mouse. <Step 102> In order to efficiently obtain the coordinates after rotational transformation, a reference vector in the rotated coordinate system is calculated in this step. By doing this, each coordinate after rotation can be determined as a linear sum of the reference vectors. When rotating the rotation axes in the order of the X-axis, y-axis, and Z-axis, each reference vector after rotation is expressed by the formula
Refer to equation 3. ). <Step 103> A three-dimensional image is created on the image memory 12 while referring to the data in the three-dimensional memory 11. In the three-dimensional display method to which the present invention is applied, the pixel value of each pixel on the image memory 12 is calculated. The calculation method for each pixel value will be described in detail later, dividing it into surface rendering and volume rendering. For an example application of the present invention to surface rendering, see Step 2.
00-205, and an example of application to volume rendering will be described later in steps 300-308. <Step 104> After the calculation of pixel values for all points on the image memory 12 is completed, filtering and LUT processing are performed on the image on the image memory 12 as necessary, and the results are displayed on the display 14. . <Step 105> If you wish to observe from a different angle, return to step 101 and repeat the rotation angle input and rotation/display processing.

<<Example of Application to Surface Rendering>> An example of application of the present invention to surface rendering will be explained using the flowchart of FIG. 4 and FIG. 7.

FIG. 7 shows the positional relationship between the image memory and the three-dimensional memory. In the figure, 10 is a screen for creating a three-dimensional image, and corresponds to an image memory. 1
1 indicates a three-dimensional data space and corresponds to a three-dimensional memory. The three-dimensional data space 11 in this figure shows the one after rotation processing. First, the brightness at point 12 on the screen is calculated. The brightness at this point 12 is determined by the calculation described below from the value of the pixel column 14 in the three-dimensional memory space that intersects with the line of sight passing through this point 12, that is, the normal line 13 of the screen 10. <Step 200> In FIG. 7, normal line 13 passing through pixel 12 on the screen
The pixel rows 14 of the three-dimensional data existing above are scanned in order from the pixels closest to the screen. <Step 201> If scanning is performed to the end of the pixel row 14, it is assumed that the target object does not exist on the normal line 13 thereof. At this time, the value of pixel 12 on the screen is set to 0, and the calculation moves on to the next pixel on the screen. <Step 202> A pixel taken out from the pixel row 14 is compared with a preset threshold value to determine whether the pixel represents the surface of an object. If it is a front pixel, the brightness of the pixel 12 on the screen is calculated in step 203. If it is not the surface, the process returns to step 200 and continues scanning backwards in the pixel column. <Step 203> For the pixel representing the object surface determined in step 202, a concentration gradient vector on the surface is determined from the values of its neighboring pixels. The concentration gradient vector calculation process starts from step 500.
504 will be explained in detail. <Step 204> The brightness of the point 12 on the screen 10 is calculated from the concentration gradient vector, line-of-sight direction vector, and light ray vector obtained in step 203. <Step 205> The brightness of the point 12 on the screen 10 calculated in step 204 is stored in the image memory.

<<Example of Application to Volume Rendering>> An example of application of the present invention to volume rendering will be explained using the flowchart of FIG. 5 and FIG. 7. <Step 300> An initial value of the amount of incident light is set. (Initial value = 1) Sets the brightness value on the screen. (Initial value = 0) <Step 3
01> In FIG. 7, the pixel row 14 in the three-dimensional data 11 existing on the normal line 13 passing through one point 12 on the screen 10 is
Scan pixels starting from those closest to the screen. <Step 302> If all the pixels in the pixel row 14 on the normal line 13 have been scanned, or if the amount of light incident on each pixel has been attenuated to 0, the pixel row 14 calculated so far is All the values of each pixel above are integrated and stored in the image memory as the brightness of point 12 on the screen. <Step 303> If the density value of the scanned pixel is lower than the value preset as the nozzle level, the following calculation for that pixel is skipped. <Step 304> A density gradient vector at the scanned pixel is determined from the values of neighboring pixels of that pixel. The method for calculating the concentration gradient vector will be described in steps 500-502. <Step 305> The luminance of each pixel on the pixel row 14 is calculated from the density gradient vector, line of sight direction vector, and light ray vector obtained in step 304. <Step 306> Calculate the amount of incident light for the next pixel. Returning to step 301, the above calculation is repeated for the next pixel in the depth direction.

In the gray level gradient method to which the present invention is applied, the brightness of each pixel is basically determined from the angle formed between the incident light and the normal to a virtual plane considered for each pixel. The normal direction of the virtual plane is defined as the direction of maximum concentration gradient. The method for calculating the concentration gradient vector will be explained below using the flowchart shown in FIG. <Step 500> Calculate coordinates after rotational transformation. By performing rotation processing on the reference vector before rotation and obtaining the reference vector after rotation in advance, the coordinates after rotation can be obtained as a linear sum of these vectors. <Step 501> The point on the grid before rotation that is closest to the coordinates determined in the previous step is set as the reference point for the concentration gradient to be determined. <Step 502> The concentration gradient amount in each reference axis direction is calculated from the density value of the grid point adjacent to the reference point. Let the coordinates of the reference point be P(x0, y0,
z0), the amount of concentration gradient in each reference axis direction is given by [Equation 4]
] formula (see formula [4] in [Fig. 10]). Note that f(x, y, z) is the density value at the point (x, y, z). <Step 503> From the concentration gradient amounts dx, dy, and dz in each reference axis direction and the three-dimensional rotation matrices Az, Ay, and Az around each axis, the concentration gradient vector at the rotated coordinates is calculated using the formula [Equation 5] ([Figure 10] (See equation 5.).

[0032]

According to the present invention, it is possible to obtain a three-dimensional image in which striped artifacts are suppressed at almost the same speed as the conventional method using nearest neighbor approximation, and it is possible to obtain a three-dimensional image in which striped artifacts are suppressed. In comparison, a high-quality three-dimensional display image with minimal deterioration in image quality can be obtained in a much shorter time. Furthermore, since the method of the present invention is simpler in processing than the conventional method, it can be easily installed in hardware.

[Brief explanation of the drawing]

FIG. 1 is a flowchart of a gray level gradient calculation method proposed by the present invention.

FIG. 2 is a configuration diagram of an apparatus implementing the present invention.

FIG. 3 is a main flowchart of an embodiment of the present invention.

FIG. 4 is a detailed explanation of a part of FIG. 2, and is a flowchart of an embodiment in which the present invention is applied to surface rendering.

FIG. 5 is a flowchart of an embodiment similarly applied to volume rendering.

FIG. 6 is a flowchart of a concentration gradient vector calculation method in an embodiment.

FIG. 7 shows the positional relationship between a screen for creating a three-dimensional image and rotated three-dimensional data used to calculate values on this screen.

FIG. 8 is a reference diagram showing the concept of a conventional method and a method of the present invention.

FIG. 9 shows formulas [Math. 2] and [Math. 3].

FIG. 10 shows equations [4] and [5].

Claims (6)

[Claims] Claim 1: Calculate a density gradient vector for each voxel of three-dimensional image data from surrounding density changes,
The brightness value in each voxel is determined from this density gradient vector and the incident light vector to the voxel, and the brightness value is
3 shaded by projection processing onto a dimensional plane
In a three-dimensional image display method for creating a dimensional image, when displaying a three-dimensional image from any direction, the amount of density change in each reference axis direction in the coordinate system before rotation is determined for each voxel,
A high-speed three-dimensional display method characterized in that the density change amount in each reference axis direction in the rotated coordinate system is estimated from this density change amount, and a vector having each density change amount as a component is used as the concentration gradient vector. . 2. The method according to claim 1, wherein the projection process of the three-dimensional image includes a process of calculating the rotated coordinates, and the pixel range for calculating the rotated coordinates is set to the pixels necessary for the projection process. A high-speed three-dimensional display method characterized by limiting to a range. 3. A high-speed three-dimensional display method according to claim 1, characterized in that the values of the individual coordinates after rotation are determined using nearest-neighbor interpolation. 4. In the method according to claim 1, the concentration gradient vector has a reference position on a grid point of a coordinate system before rotation whose coordinates after rotation are obtained by nearest neighbor approximation. A high-speed three-dimensional display method characterized by: 5. In the method according to claim 4, the concentration gradient vector is obtained by first determining a vector whose components are the amount of change in concentration in each reference axis direction determined from the values of grid points adjacent to the reference position; A high-speed three-dimensional display method characterized in that the vector is then calculated by linearly transforming it using a transformation matrix used when rotating the three-dimensional data. 6. In the method according to claim 4, the concentration gradient vector is defined as x for three-dimensional data, with dx, dy, and dz being the difference amounts in each axis direction of grid points adjacent to the reference position.
A high-speed three-dimensional display method characterized in that rotation matrices around the axes, y-axes, and z-axes are determined by the following equation [Equation 1] when Ax, Ay, and Az are rotation matrices. [Formula 1] Note that the arrangement of Ax, Ay, and Az changes depending on the rotation order. θz, θy, and θx are rotation angles around the z, y, and x axes.