JP7003635B2 - Computer program, image processing device and image processing method - Google Patents

Description

The present disclosure relates to computer programs, image processing devices and image processing methods .

In the medical field, it is continuously photographed at predetermined slice intervals by medical diagnostic imaging equipment such as CT (Computed Tomography), MRI (Magnetic Resonance Imaging), and PET (Positron Emission Tomography), and DICOM (Digital Imaging and COmmunications in Medicine). ) A technique has been developed for three-dimensionally visualizing an observation object such as an organ based on a two-dimensional slice image stored in the format to be useful for diagnostic imaging support, surgery or treatment support.

A volume rendering method is known as a typical method for displaying a three-dimensional image of the internal structure of an observation object. Patent Document 1 discloses a ray casting process as one of the volume rendering methods. In the ray casting process, a plurality of slice images are stacked to construct a three-dimensional voxel space, and a ray (virtual ray) is applied to the three-dimensional voxel space from an arbitrary viewpoint direction, and a projection plane orthogonal to the viewpoint direction ( The projection result is obtained as a pixel value on the screen). Specifically, coordinate transformation is performed so that the viewpoint direction is parallel to the Z axis of the three-dimensional coordinate axis, and a ray is applied to the three-dimensional boxel space after the coordinate transformation in the Z-axis direction on the XY plane. The projection result is calculated as a pixel value on the projection surface (screen).

Japanese Unexamined Patent Publication No. 11-175743

This coordinate conversion is generally performed by a real number operation, but the coordinate values after the coordinate conversion are converted into integers into discrete values called voxel coordinates, and rounding errors are accumulated in this process, and a periodic pattern is generated. This periodic pattern is a weak one that cannot be recognized by the naked eye with a single slice image, but when trying to synthesize an image in which multiple slice images are stacked as in the present application, even if each slice image is subjected to coordinate conversion. Even if the pattern is deviated, interference fringes are generated if the phase condition of the periodicity of the pattern is satisfied, and it is recognized as a visually remarkable moire. In addition, if irregular patterns similar to a plurality of slice images exist, each irregular pattern shifts due to rotation when coordinate conversion is performed on a three-dimensional voxel space composed of stacked slice images. By overlapping and connecting, a remarkable pattern with regularity appears, and it is recognized as moire in the three-dimensional image of the observation object. When such moiré occurs, the image quality deteriorates, it becomes difficult to distinguish the lesion, or there is a possibility that the moiré is erroneously identified as the lesion.

The present disclosure has been made in view of such circumstances, and clarifies a computer program, an image processing device, and an image processing method capable of suppressing moire.

The present application includes a plurality of means for solving the above problems. For example, a computer program is a computer program for causing a computer to perform a volume rendering process using a raycasting method, and is a computer. In addition, the process of acquiring image data of a plurality of xy plane images captured at predetermined intervals along the z direction, and the pixel values of arbitrary attention pixels of each of the plurality of xy plane images and pixels in the vicinity of the attention pixels. Based on this, a process of smoothing the plurality of xy plane images, and a three-dimensional voxel image composed of voxels corresponding to each pixel of the plurality of xy plane images whose RGB values and opacity are determined and smoothed. And the process of generating a converted voxel image by performing predetermined coordinate conversion on the three-dimensional voxel image, and when a virtual ray is emitted from a predetermined viewpoint toward the converted voxel image. The line of sight intersects the projection plane with the RGB value calculated based on the cumulative sum of the RGB values and opacity of the voxels and the cumulative brightness value obtained by multiplying the transmission intensity of the virtual light beam for each voxel on the line of sight. The process of generating a rendered image as the pixel value of the pixel to be performed is executed.

According to the present disclosure, moire can be suppressed.

It is a block diagram which shows an example of the structure of the ray casting apparatus as the image processing apparatus of this embodiment. It is a schematic diagram which shows an example of the tomographic image of DICOM format. It is explanatory drawing which shows the 1st example of the color map data of this embodiment. It is explanatory drawing which shows the 2nd example of the color map data of this embodiment. It is explanatory drawing which shows an example of the gradation compression processing by the ray casting apparatus of this embodiment. It is explanatory drawing which shows an example of the smoothing process by the ray casting apparatus of this embodiment. It is a schematic diagram which shows an example of the voxel structure of this embodiment. It is explanatory drawing which shows an example of the coordinate transformation processing by the ray casting apparatus of this embodiment. It is explanatory drawing which shows an example of the search control mask data generation processing by the ray casting apparatus of this embodiment. It is explanatory drawing which shows an example of the ray casting process in the case of the parallel projection by the ray casting apparatus of this embodiment. It is explanatory drawing which shows an example of the ray casting process in the case of the perspective projection by the ray casting apparatus of this embodiment. It is a flowchart which shows an example of the procedure of the ray casting process by the ray casting apparatus of this embodiment. It is a flowchart which shows an example of the procedure of the ray casting process by the ray casting apparatus of this embodiment. It is a flowchart which shows an example of the procedure of the coordinate conversion processing by the ray casting apparatus of this embodiment. It is a flowchart which shows an example of the procedure of the coordinate conversion processing by the ray casting apparatus of this embodiment. It is a schematic diagram which shows an example of the interpolation processing by the ray casting apparatus of this embodiment. It is a flowchart which shows an example of the procedure of the ray casting process by the ray casting apparatus of this embodiment. It is a flowchart which shows an example of the procedure of the ray casting process by the ray casting apparatus of this embodiment. It is a block diagram which shows the other example of the structure of the image processing part of this embodiment. It is a schematic diagram which shows an example of how the curved moire occurs. It is explanatory drawing which shows the 1st example of the moire measures result by the ray casting apparatus of this embodiment. It is explanatory drawing which shows the 2nd example of the moire measures result by the ray casting apparatus of this embodiment. It is a schematic diagram which shows an example of the appearance of the stripe / grid-like moire. It is explanatory drawing which shows an example of the appearance of the rounding error at the time of coordinate conversion. It is explanatory drawing which shows an example of the state of the rounding error dispersion by the ray casting apparatus of this embodiment. It is explanatory drawing which shows the 3rd example of the moire measures result by the ray casting apparatus of this embodiment. It is explanatory drawing which shows the 4th example of the moire measures result by the ray casting apparatus of this embodiment. It is explanatory drawing which shows the 5th example of the moire measures result by the ray casting apparatus of this embodiment. It is explanatory drawing which shows the 6th example of the moire measures result by the ray casting apparatus of this embodiment.

Hereinafter, embodiments of the present disclosure will be described with reference to the drawings. FIG. 1 is a block diagram showing an example of the configuration of the ray casting device 100 as the image processing device of the present embodiment. The ray casting device 100 of the present embodiment includes an input unit 10, a storage unit 11, an output unit 12, an image processing unit 20, and the like. The image processing unit 20 includes a gradation compression processing unit 21, a smoothing processing unit 22, a voxel structure generation unit 23, a coordinate conversion processing unit 24, a raycasting processing unit 25, a search control mask generation unit 26, and a coordinate value calculation unit 27. , A random number generation unit 28, an interpolation processing unit 29, a ROI (Region of Interest) clipping processing unit 30, and the like.

The input unit 10 has a function as an acquisition unit, and acquires image data of a plurality of xy plane images captured at predetermined intervals along the z direction. More specifically, the input unit 10 is a plurality of tomographic images taken at predetermined intervals with respect to an observation object (for example, various parts of the human body / animal including organs, brains, bones, blood vessels, etc.). Acquire image data in DICOM (Digital Imaging and COmmunications in Medicine) format of a two-dimensional tomographic image.

FIG. 2 is a schematic diagram showing an example of a DICOM format tomographic image. When the xyz axis is set as shown in FIG. 2, a plurality of slice images (two-dimensional tomographic images), which are xy plane images, are arranged along the z direction. Sx is the number of pixels in the x-axis direction, Sy is the number of pixels in the y-axis direction, and Sz is the number of sliced images. Note that Sz represents the number of pixels in the z-axis direction. Rxy is the resolution in the x-direction and the y-direction, and indicates the reciprocal of the pixel spacing, that is, the number of pixels per unit distance. The resolutions in the x and y directions are the same, but may be different. Rz is the resolution in the z direction, and represents the reciprocal of the spacing between adjacent slice images, that is, the number of slices per unit distance. The pixel value at the coordinates (x, y, z) of each slice image is represented by Do (x, y, z). Do (x, y, z) can take a value of -32768 or more and 32767 or less. For example, in the case of a CT image, the pixel value Do (x, y, z) shows 0 for water, a tissue having a density higher than water (for example, bone) has a positive value, and a tissue having a density lower than water. (For example, adipose tissue) has a negative value. In the coordinates (x, y, z) of the slice image, the x coordinate takes a value of 0 or more and Sx-1 or less, the y coordinate takes a value of 0 or more and Sy-1 or less, and the z coordinate takes a value of 0 or more and Sz. Takes a value less than or equal to -1.

Further, the input unit 10 acquires color map data in which RGB values and opacity are defined in association with pixel values.

FIG. 3 is an explanatory diagram showing a first example of the color map data of the present embodiment. The color map data of the first example is shown by Cmap (v, c). The variable v indicates a pixel value of Do (x, y, z), and represents a numerical value within the range of -23768 or more and 32767 or less, similarly to the range of Do (x, y, z). The variable c includes the values of 0, 1, 2, and 3, c = 0 indicates the R value of RGB, c = 1 indicates the G value, and c = 2 indicates the B value. c = 3 indicates the opacity α. Cmap (v, c) can take a numerical value of 0 or more and 255 or less. As shown in FIG. 3, the color map data defines the relationship between the pixel value Do (x, y, z) of the sliced image, the RGB value, and the opacity α. For example, when the pixel value Do (x, y, z) is 32767, the R value is R (32767), the G value is G (32767), the B value is B (32767), and the opacity is α (32767). ). Here, R (32767), G (32767), B (32767), and α (32767) indicate for convenience that they correspond to the pixel value 32767, and are actually appropriate values. Further, the color map data does not need to be defined for all the pixel values from -32768 to 32767, and in the case of a CT image, it is usually adjusted to take the range of -2048 to 2048. The RGB values and opacity in the range from 2048 to 2048 and from 2048 to 32767 may all be set to 0.

FIG. 4 is an explanatory diagram showing a second example of the color map data of the present embodiment. The color map data of the second example is shown by Cmap8 (v, c). As will be described later, when the DICOM image is gradation-compressed from, for example, 16 bits to 8 bits, the pixel value can be in the range of 0 or more and 255 or less from the range of 32768 or more and 32767 or less. The variable v of Cmap8 (v, c) is a numerical value within the range of 0 or more and 255 or less. The variable c is the same as in the first example.

By using the color map data Cmap (v, c) and Cmap8 (v, c), for example, not only different colors can be added to each part of the human body, but also the opacity can be set for each part. It is possible to make it easier to recognize the observation target, such as making the located organ transparent and depicting the organ hidden in the back.

Further, the input unit 10 includes coordinate conversion parameters including the angle of rotation of the xyz axis, the offset value in the xyz axis direction, the enlargement or reduction magnification in the xyz axis direction, the scaling factor in the z axis direction, and the distance from the gazing point to the viewpoint. And the value of ROI (Region of Interest) in the xyz direction is acquired. Details of the coordinate conversion parameters and ROI values will be described later.

The storage unit 11 can store the data acquired by the input unit 10, the processing result of the image processing unit 20, and the like. The color map data, the coordinate conversion parameters, the ROI value, and the like may be stored in the storage unit 11 in advance.

The output unit 12 outputs the processing result of the image processing unit 20, for example, the image data of the rendered image generated by the raycasting processing unit 25 to a display device (not shown).

The gradation compression processing unit 21 converts a DICOM format tomographic image into, for example, an 8-bit gradation compression tomographic image.

FIG. 5 is an explanatory diagram showing an example of gradation compression processing by the ray casting apparatus 100 of the present embodiment. The gradation compression processing unit 21 specifies the minimum value and the maximum value of the pixels of a predetermined slice image among the plurality of slice images. For example, among the slice images from 1 to Sz, the minimum value Dmin and the maximum value Dmax of all the pixels in the Sz / second intermediate slice image can be specified. As the predetermined slice image, the method of specifying the minimum value and the maximum value of all the pixels in all the slice images is originally accurate. However, in that case, it becomes necessary to temporarily hold all the large-capacity 16-bit DICOM format tomographic images in the memory, and the gradation compression effect is halved. Further, the minimum value Dmin and the maximum value Dmax are only used as parameters for gradation compression, and do not directly affect the accuracy of the 8-bit gradation compression tomographic image. Therefore, a method of roughly specifying is taken by using the minimum value and the maximum value of the pixels of a single slice image. However, since the subject is often not properly reflected in the first slice image, a method of specifying the minimum value Dmin and the maximum value Dmax is adopted only in the intermediate slice image. As a result, it is possible to directly construct only an 8-bit gradation compressed tomographic image in the memory without holding a large-capacity 16-bit DICOM format tomographic image in the memory.

The 8-bit gradation compressed tomographic image D8 (x, y, z) is based on the formula D8 (x, y, z) = (Do (x, y, z) -Lmin) · 255 / (Lmax-Lmin). Can be generated. However, when D8 (x, y, z)> 255, D8 (x, y, z) = 255, and when D8 (x, y, z) <0, D8 (x, y, z). = 0. Here, Lmin = (Dmax-Dmin) · γ + Dmin, Lmax = (Dmax−Dmin) · (1-γ) + Dmin. γ is the contrast adjustment range of the gradation compressed image, and the closer it is to 0, the higher the contrast but the lower the brightness. Normally, it is set to γ = 0.1. Since the luminance contrast of the rendered image can be adjusted by various settings such as color map data, γ may be a fixed value. Further, 0 ≦ D8 (x, y, z) ≦ 255, 0 ≦ x ≦ Sx-1, 0 ≦ y ≦ Sy-1, and 0 ≦ z ≦ Sz-1. The resolution in the xy direction is Rxy, and the resolution in the z direction is Rz.

The color map Cmap8 (v, c) shown in FIG. 4 can be applied to the 8-bit gradation compressed tomographic image D8 (x, y, z).

As described above, the gradation compression processing unit 21 specifies the minimum value Dmin and the maximum value Dmax of the pixels of the predetermined xy plane image among the plurality of xy plane images, and the upper limit value Lmax smaller than the specified maximum value. And the lower limit Lmin larger than the specified minimum value is calculated. The gradation compression processing unit 21 compresses within the range of the upper limit value Lmax and the lower limit value Lmin of the pixel value of each pixel of the plurality of xy plane images. That is, the gradation compression processing unit 21 compresses the range of the upper limit value Lmax and the lower limit value Lmin into, for example, 256 steps, sets the upper limit value Lmax or more to 255, and sets the lower limit value Lmin or less to 0.

The smoothing processing unit 22 has a function as a smoothing unit, and smoothes a plurality of xy plane images based on the pixel values of an arbitrary pixel of interest of each of the plurality of xy plane images and pixels in the vicinity of the pixel of interest. do. The xy plane image may be a DICOM image Do (x, y, z) or a gradation compression tomographic image D8 (x, y, z) that is gradation-compressed by the gradation compression processing unit 21.

FIG. 6 is an explanatory diagram showing an example of smoothing processing by the ray casting apparatus 100 of the present embodiment. In the example of FIG. 6, the gradation compressed tomographic image D8 (x, y, z) is shown, but it may be a DICOM image Do (x, y, z). FIG. 6 shows the attention pixel D8 (x, y, z) of one slice image and the eight neighboring pixels D8 (x-1, y-1, z) and D8 (x,) existing in the vicinity of the attention pixel. y-1, z), D8 (x + 1, y-1, z), D8 (x-1, y, z), D8 (x + 1, y, z), D8 (x-1, y + 1, z), D8 (X, y + 1, z), D8 (x + 1, y + 1, z) are shown. Assuming that the smoothed image of the xyz coordinates is D (x, y, z), D (x, y, z) is the average value of the pixel values of the neighboring pixels including the pixel of interest.

That is, more specifically, the smoothing processing unit 22 has a pixel value of an arbitrary pixel of interest on one xy plane image of each of the plurality of xy plane images and a pixel value of a neighboring pixel in the xy plane image of one of the attention pixels. Based on, one xy plane image is smoothed.

Anisotropic (meaning excluding the z direction) smoothing is performed using the mean value near 8 on the xy plane of the gradation compressed tomographic image D8 (x, y, z) or the DICOM image Do (x, y, z). As a result, it is possible to suppress moire generated by overlapping and connecting irregular patterns that can be perceived in a plurality of slice images by shifting them due to rotation.

The voxel structure generation unit 23 has a function as a voxel image generation unit, and is composed of voxels corresponding to each pixel of each of a plurality of smoothed xy plane images in which RGB values and opacity α are determined. Generate a 3D voxel image. More specifically, the voxel structure generation unit 23 associates the RGB values and opacity of the color map data corresponding to the pixel values of each pixel of each of the smoothed plurality of xy plane images with the RGB values and opacity. Set α to generate a 3D voxel image. Hereinafter, the three-dimensional voxel image is also referred to as a voxel structure.

FIG. 7 is a schematic diagram showing an example of the voxel structure of the present embodiment. FIG. 7A shows a voxel structure before coordinate conversion, which will be described later, and FIG. 7B shows a voxel structure after coordinate conversion. First, the voxel structure before the coordinate conversion will be described. The voxel structure is a set of voxels in which volume data composed of RGB values and elements of opacity α are three-dimensionally packed and arranged in a grid pattern. One voxel corresponds to, for example, one pixel of the smoothed image D (x, y, z). When the voxel structure before coordinate conversion is represented by V (x, y, z, c), it can be represented by V (x, y, z, c) = Cmap8 (D (x, y, z), c). can. Here, 0 ≦ x ≦ Sx-1, 0 ≦ y ≦ Sy-1, 0 ≦ z ≦ Sz-1, and 0 ≦ c ≦ 3. The above equation means that the voxel structure V (x, y, z, c) is applied to the pixel values of the smoothed image D (x, y, z) by applying the color map Cmap8 (v, c). ) Is generated.

The data of the three-dimensional voxel image of the present embodiment is composed of voxels corresponding to each pixel of a plurality of xy plane images captured at predetermined intervals along the z direction, and the voxel data is composed of a plurality of xy planes. A plurality of xy plane images are smoothed based on the pixel values of an arbitrary pixel of interest in each image and pixels in the vicinity of the pixel of interest, and RGB values and opacity are determined in association with the color map data. .. In other words, the data of the three-dimensional voxel image is the voxel data of the voxels corresponding to each pixel of the plurality of xy plane images captured at predetermined intervals along the z direction, and is a plurality of xy plane images. A plurality of smoothed data structures in which a plurality of xy plane images are smoothed based on the pixel values of each arbitrary pixel of interest and pixels in the vicinity of the pixel of interest, and RGB values and opacity are determined. It is used to execute a process of generating a three-dimensional voxel image composed of voxels corresponding to each pixel of the xy plane image.

The coordinate conversion processing unit 24 has a function as a post-conversion voxel image generation unit, and converts the generated three-dimensional voxel image by performing predetermined coordinate conversion using the coordinate conversion parameters acquired by the input unit 10. Generate a post-voxel image. That is, the coordinate conversion processing unit 24 performs a predetermined coordinate conversion process on the voxel structure V (x, y, z, c) to generate the converted voxel structure. The converted voxel structure is represented as V'(x, y, z', c). The converted voxel structure V'(x, y, z', c) is V'(x, y, z', c) = Matrix (4 × 4) · V (x, y, z, c). Can be represented by. Here, 0 ≦ x ≦ Sx-1, 0 ≦ y ≦ Sy-1, 0 ≦ z ′ ≦ Sz ′-1, and 0 ≦ c ≦ 3. Note that z'is corrected so that the resolution matches the xy axis with respect to the z axis, and z'= z · Rxy / Rz and Sz ′ = Sz · Rxy / Rz. Matrix (4 × 4) is a transformation matrix for performing coordinate transformation, and a specific transformation formula will be described later.

As shown in FIG. 7A, if the orientation of the projection plane on which the rendered image is generated and the z-axis direction of the voxel structure V (x, y, z, c) do not match, the raycasting process described later is complicated. become. Therefore, as shown in FIG. 7B, the xyz axis of the voxel structure V (x, y, z, c) is rotated so that, for example, the z'axis after coordinate conversion is in the direction perpendicular to the projection plane. Therefore, the ray casting process can be simplified.

FIG. 8 is an explanatory diagram showing an example of coordinate conversion processing by the ray casting apparatus 100 of the present embodiment. As shown in FIG. 8, the converted voxel structure V'(x, y, z', c) is a voxel structure before conversion according to ROI with respect to V (x, y, z, c). It can be obtained by clipping, setting the region of interest, and then performing coordinate transformation. In clipping by ROI, Xs-Xe is set as the ROI in the x-axis direction, Ys-Ye is set as the ROI in the y-axis direction, and Zs-Ze is set as the ROI in the z-axis direction. Here, 0 ≦ Xs <Xe ≦ Sx-1, 0 ≦ Ys <Ye ≦ Sy-1, and 0 ≦ Zs <Ze ≦ Sz-1.

The coordinate transformation can be performed in the order of scaling, z-direction scaling processing, offset, rotation, and perspective transformation, for example, but the order of each processing is as long as it can be defined by a 4 × 4 transformation matrix. It is not limited to the example of FIG. The order of each process in the inverse transformation described later is the reverse order of the order of each process in the coordinate transformation, and the inverse matrix of the transformation matrix is applied.

Scaling is, for example, a process of expanding an area of interest when an area other than the area of interest is excluded by the ROI process. In scaling, the same scaling factor Scale is used in the xyz axis direction.

The z-direction scaling process is a process for adjusting the physical spacing of pixels in the z-direction to the physical spacing of pixels in the xy direction. In the z-direction scaling process, the z-direction scaling factor Scz (= Rxy / Rz) is used.

The offset is a process for translating in the x-axis direction, the y-axis direction, and the z-axis direction. For the offset, the offset Xoff in the x-axis direction, the offset Yoff in the y-axis direction, and the offset Zoff in the z-axis direction are used.

The rotation is a process for rotating around the x-axis, around the y-axis, and around the z-axis. In rotation, the rotation angle Rx at the center of the x-axis, the rotation angle Ry at the center of the y-axis, and the rotation angle Rz at the center of the z-axis (all angles are in radians) are used.

Perspective transformation is a process for displaying an image in a so-called endoscope mode. In the perspective transformation, the distance Dist from the gaze point (origin, if there is an offset, the origin after the offset) to the viewpoint in the z-axis direction is used. In the case of parallel projection, Dist = 0 is set.

By the above-mentioned clipping and coordinate transformation, the value of the converted voxel structure V'(x, y, z', c) is determined. In addition, each processing of scaling, z-direction scaling processing, offset, rotation, and perspective transformation can be summarized in a transformation matrix Matrix (4 × 4).

The search control mask generation unit 26 is a voxel of the converted voxel structure V'(x, y, z', c) for each xy coordinate of the converted voxel structure V'(x, y, z', c). Specify the minimum and maximum values on the line of sight where the opacity α of is not 0 (transparent).

FIG. 9 is an explanatory diagram showing an example of the search control mask data generation process by the ray casting apparatus 100 of the present embodiment. FIG. 9 shows a case where the converted voxel structure V'(x, y, z', c) is viewed from the xz'plane. The projection plane is an xy plane. The search control mask generation unit 26 generates the search control mask data M (x, y, s). Here, s = 0 or 1.

As shown in FIG. 9, the range of effective voxels to be searched is set in the z-axis direction for each two-dimensional coordinate (x, y). A valid voxel is a voxel having a non-zero opacity α. When the opacity α is 0, it means that the voxel is transparent, and the voxel is excluded from the calculation target and is not displayed in the rendered image.

For each coordinate (x, y) of 0 ≦ x ≦ Sx-1 and 0 ≦ y ≦ Sy-1, increase z ′ by 1 from z ′ = 0 toward z ′ = Sz ′ -1, and first. Search for z1 that satisfies V'(x, y, z1,3)> 0. Further, for the same coordinates (x, y), z'is reduced one by one from z'= Sz'-1 toward z'= 0, and V'(x, y, z2,3)> 0 first. Search for z2 that satisfies. That is, when z'is at least z1 and z2, V'(x, y, z', 3)> 0 is satisfied, and the voxel opacity α is not 0, but when z'is between z1 and z2, V ′ (X, y, z ′, 3) = 0 may be obtained. That is, z1 is the minimum value of z'in which a valid voxel exists, and z2 is the maximum value. The search control mask data M (x, y, s) is M (x, y, 0) = z1 and M (x, y, 1) = z2.

If z1 = z2 or z1 and z2 are not found, M (x, y, 0) = M (x, y, 1) = -1.

By generating the search control mask data M (x, y, s), invalid voxels can be excluded and subsequent ray casting processing can be performed, so that the load of ray casting processing can be reduced. The volume rendering process can be speeded up.

The ray casting processing unit 25 has a function as a rendering image generation unit, and converts a virtual light beam having a predetermined transmission intensity from a predetermined viewpoint into a converted voxel structure V'(x, y, z', c). When irradiating toward the z'direction, the intensity after passing through the voxel is calculated for each voxel on the z'axis while considering the absorption / attenuation of virtual light rays based on the opacity α of the voxel, and the transmission intensity is determined. Determine the voxel at the deepest point that is below the predetermined lower limit. Then, while tracing each voxel in which the virtual light beam has passed from the deepest voxel in the reverse direction, the RGB value is cumulatively added based on the RGB value and the opacity α of each voxel, and the line of sight is the projection plane. A rendered image is generated as the pixel value of the intersecting pixels. In the process of irradiating the virtual light beam toward the converted voxel structure V'(x, y, z', c) in the z'direction and determining the voxel at the deepest point, the RGB value and opacity α of each voxel In addition, a method of simultaneously calculating the pixel values of the pixels intersecting the projection surface can be obtained by cumulatively adding the products of the transmission intensities of the virtual rays. In the present application, the latter method will be described below.

FIG. 10 is an explanatory diagram showing an example of a ray casting process in the case of parallel projection by the ray casting apparatus 100 of the present embodiment. In FIG. 10, the converted voxel structure V'(x, y, z', c) is illustrated as a cube filled with voxel data in three dimensions. The projected plane (rendered image) parallel to the xy axis of the converted voxel structure V'(x, y, z', c) is the converted voxel structure V'(x, y, z', c). It is assumed that they are arranged at a distance from the required distance. As shown in FIG. 10, in parallel projection generally used when the viewpoint is located at the point at infinity (originally, Dist = ∞, but in this application, Dist = 0 is set), and the inside of the voxel structure is set from the viewpoint. All lines of sight towards are parallel to the z'axis. In the ray casting process (volume rendering), the inside of the converted voxel structure V'(x, y, z', c) is visualized by performing semi-transparent display (α blending) based on the opacity α. .. That is, the visualization is realized by the voxel values (RGB values and opacity α) possessed by the voxels of the converted voxel structure V'(x, y, z', c).

FIG. 11 is an explanatory diagram showing an example of a ray casting process in the case of perspective projection by the ray casting apparatus 100 of the present embodiment. As shown in FIG. 11, when the viewpoint is relatively close to the projection surface (Dist> 0 in the present application), the perspective projection is performed, and in the case of endoscopic display, the projection surface and the viewpoint are the voxel structure V'(x, It may be immersed inside y, z', c). Similar to the case of FIG. 10, in the ray casting process (volume rendering), the converted voxel structure V'(x, y, z') is displayed by translucent display (α blending) based on the opacity α. , C) Visualize the inside. That is, the visualization is realized by the voxel values (RGB values and opacity α) possessed by the voxels of the converted voxel structure V'(x, y, z', c).

As shown in FIGS. 10 and 11, the line of sight is extended from the viewpoint to the inside of the converted voxel structure V'(x, y, z', c), and the voxel is located at the sampling point (voxel) on the line of sight. The RGB brightness value is cumulatively calculated along the line of sight based on the RGB value and the opacity α, and the cumulatively calculated value is used as a pixel (pixel value of the pixel) at the position where the line of sight intersects the projection plane to generate a rendered image. .. At this time, the cumulative calculation is performed for z'between the search control masks M (x, y, 0) and M (x, y, 1) described above. The details of the ray casting process will be described later.

Next, the operation of the ray casting device 100 of the present embodiment will be described.

12 and 13 are flowcharts showing an example of the procedure of the ray casting process by the ray casting apparatus 100 of the present embodiment. Hereinafter, for convenience, the subject of processing will be described as the image processing unit 20. The image processing unit 20 acquires the image data of the DICOM image (S11) and acquires the color map data (S12). The color map data is illustrated in FIGS. 3 and 4.

The image processing unit 20 acquires the coordinate conversion parameter (S13) and acquires the ROI value in the xyz direction (S14). The coordinate conversion parameters and ROI values are exemplified in FIG. The image processing unit 20 performs gradation compression processing of the DICOM image (S15), and defines a color map for the gradation compression tomographic image D8 (S16). The color map definition defines the RGB value and the opacity α with respect to the pixel value of the gradation compressed tomographic image D8. The process of step S15 is not essential, and if the process of step S15 is not performed, the color map may be defined for the DICOM image. However, in the case of a CT image, the color map can be generally defined for each part to be viewed and is not defined depending on the given DICOM image. Therefore, the process of step S16 is omitted. A predefined color map can also be used.

The image processing unit 20 performs smoothing processing of the gradation compressed tomographic image D8 (S17) to generate a voxel structure V (x, y, z, c) (S18). The image processing unit 20 performs coordinate conversion processing on the voxel structure V (x, y, z, c) (S19), and generates a converted voxel structure V'(x, y, z', c). (S20). The details of the coordinate conversion process will be described later.

The image processing unit 20 generates search control mask data M (x, y, s) (S21), performs ray casting processing (S22), and generates a rendered image Image (x, y, c) (S23). , End the process. Here, 0 ≦ Image (x, y, c) ≦ 255, c = 0 (R), c = 1 (G), c = 2 (B). The details of the ray casting process will be described later.

14 and 15 are flowcharts showing an example of the procedure of the coordinate conversion process by the ray casting apparatus 100 of the present embodiment. The image processing unit 20 multiplies the voxel structure V (x, y, z, c) before conversion by the transformation matrix Matrix (4 × 4), and the voxel structure V ′ (x, y, z ′) after conversion. , C) is generated (S101). That is, V'(x, y, z', c) = Matrix (4 × 4) · V (x, y, z, c). Here, the coordinate values (x, y, z) of the voxel structure V (x, y, z, c) before conversion and the voxel structure V'(x, y, z', c) after conversion. The coordinate value (x, y, z') is an integer value. Therefore, the voxel structure V (x, y, The coordinate values (xr, yr, zr) of the real numbers of z, c) are calculated as follows.

The image processing unit 20 converts the coordinate values (x, y, z') (integer values) of the converted voxel structure V'(x, y, z', c) into real values (xx, yy, zz). Convert (S102). The conversion from an integer value (x, y, z') to a real value (xx, yy, zz) can be performed by the equations (1), (2) and (3).

The image processing unit 20 multiplies the boxel structure V'(xx, yy, zz, c) obtained by converting the coordinate values into real values by the inverse matrix of the transformation matrix Matrix (4 × 4), and the boxel structure before conversion. A real value (xr, yr, zr) of a coordinate value corresponding to V (x, y, z, c) is calculated (S103).

The processing of steps S102 and S103 can be performed by the coordinate value calculation unit 27. That is, the coordinate value calculation unit 27 calculates the real value of the xyz coordinates of the voxel by performing the inverse conversion of the predetermined coordinate transformation with respect to the xyz coordinates of the voxel of the converted voxel structure. The process of step S103 can be performed as follows.

First, when the distance Dist> 0 of the perspective transformation, the perspective transformation is performed as in the equations (4), (5), and (6), and the converted (xx', yy', zz') is converted to (xx, yy). , Zz).

The meanings of the equations (4), (5), and (6) are that the smaller the distance Dist, the larger the near and the smaller the distance, that is, the so-called perspective is emphasized, and the larger the distance Dist, the more the perspective. The state in which the distance Dist is set to ∞ (Dist = 0 in the present application) is parallel projection.

Next, rotation in the reverse direction is sequentially performed in the order of rotation at the center of the z-axis, rotation at the center of the y-axis, and rotation at the center of the x-axis (the order opposite to the order in the coordinate transformation). That is, the rotation according to the equations (7) and (8) is performed, and (xx', yy') is replaced with (xx, yy). Next, rotation is performed according to the equations (9) and (10), and (zz', xx') is replaced with (zz, xx). Then, the rotation according to the equations (11) and (12) is performed, and (yy', zz') is replaced with (yy, zz).

Scaling, z-direction scaling, and offset (simultaneous in xyz direction) are performed at the same time, and equations (13), (14), and (15) correspond to the voxel structure V (x, y, z, c) before conversion. Obtain the real values (xr, yr, zr) of the coordinate values to be used.

The image processing unit 20 adds a random number to the calculated real value (xr, yr, zr) and updates the coordinate value (S104). This is to suppress the moire that occurs due to the accumulation of rounding errors in the process of the coordinate transformation described above and the generation of a periodic pattern. By adding random numbers, the periodicity of the pattern collapses and moire occurs. Can be suppressed. It should be noted that the process of adding a random number to the calculated real value (xr, yr, zr) to update the coordinate value is not essential.

More specifically, the random number generation unit 28 generates a random number within a predetermined range, for example, a random number Rnd from 0 to 32767. The interpolation processing unit 29 adds a random number within a predetermined range to the real value of the xyz coordinate of the voxel. That is, according to the equations (16), (17), and (18), random numbers in the range of -0.5 to +0.5 are added to the real-valued coordinates (xr, yr, zr) and updated. The coordinate values (xr', yr', zr') are obtained. Let the updated coordinate values (xr', yr', zr') be (xr, yr, zr).

The image processing unit 20 decomposes the updated coordinate values (xr, yr, zr) into integer values (xi, yi, zi) and fractions after the decimal point (wx, wy, wz) (S105). More specifically, the interpolation processing unit 29 specifies the difference (fraction) between the real value of the xyz coordinate of the voxel and the integer value obtained by rounding down the decimal point of the real value. That is, the updated coordinate values (xr, yr, zr) are decomposed as the equations (19), (20), and (21). Here, 0 ≦ wx, wy, and wz ≦ 1.

The image processing unit 20 adds a random number to the ROI value to update the clipping region (S106).

More specifically, the ROI clipping processing unit 30 sets the xyz coordinates that define the target area of the volume rendering processing. The random number generation unit 28 generates a random number within a predetermined range, for example, a random number Rnd from 0 to 32767. The ROI clipping processing unit 30 adds a random number within a predetermined range to the set xyz coordinates. That is, according to the equations (22), (23), (24), (25), (26), and (27), for example, with respect to the six coordinate values Xs, Xe, Ys, Ye, Zs, and Ze of the ROI. , The six coordinate values Xs', Xe', Ys', Ye', Zs', and Ze'of the ROI updated by adding random numbers in the range of -0.1 to +1.1 are obtained.

When Xs'<0, Xs'= 0, and when Xe'> Sx-1, Xe'= Sx-1. If Ys'<0, then Ys'= 0, and if Ye'> Sy-1, then Ye'= Sy-1. When Zs'<0, Zs'= 0, and when Ze'> Sz-1, Ze'= Sz-1.

The image processing unit 20 determines whether or not the integer value (xi, yi, zi) of the coordinate values corresponding to the voxel structure V (x, y, z, c) before conversion is in the updated clipping region. Judgment (S107). Here, as the condition 1, in the case of xi <Xs'or xi> Xe'or yi <Ys' or y> Ye'or zi <Zs' or zi> Ze', the integer value of the coordinate value (xi, yi). , Zi) is determined not to be within the updated clipping region (NO in S107), the image processing unit 20 deletes the voxel value (S108), and performs the process of step S113 described later. More specifically, the image processing unit 20 sets the converted voxel structure V'(x, y, z', c) = 0. That is, the RGB value and the opacity α are set to 0.

When the integer value (xi, yi, zi) of the coordinate value is within the updated clipping region (YES in S107), that is, when the above-mentioned condition 1 is not satisfied, the image processing unit 20 uses the voxel structure before conversion. It is determined whether or not the integer values (xi + 1, yi + 1, zi + 1) of the coordinate values corresponding to are within the updated clipping region (S109). As condition 2, in the case of xi + 1> Xe'or yi + 1> Ye'or zi + 1> Ze', it is determined that any of the integer values (xi + 1, yi + 1, zi + 1) of the coordinate values is not within the updated clipping area. (NO in S109), the image processing unit 20 performs the processing of step S113 described later without interpolating (S110). Here, "no interpolation" means that the converted voxel structure V'(x, y, z', c) = V (xi, yi, zi, c) and is an integer value of the coordinate values (xi, yi, zi). ) Is adopted as it is.

When the integer values of the coordinate values (xi + 1, yi + 1, zi + 1) are within the updated clipping region (YES in S109), that is, the integer values of the coordinate values (xi, yi, zi) and (xi + 1, yi + 1, zi + 1) are When it is in the updated clipping area, the image processing unit 20 sets 8 sets of the integer value of the coordinate value (xi, yi, zi) and the integer value of the coordinate value (xi, yi, zi) plus 1. The boxel value of the boxel structure before conversion with the numerical value as the coordinate value is multiplied by the weighting coefficient, the boxel value multiplied by the weighting coefficient is added to interpolate the boxel value (S111), and the interpolated boxel value is converted. The boxel value of the boxel structure V'(x, y, z', c) is set to (S112), and the process of step S113 described later is performed. As a result, the value of the converted voxel structure V'(x, y, z', c) is obtained.

That is, the interpolation processing unit 29 identifies a plurality of voxels whose xyz coordinates are integer values in the vicinity of the real value of the calculated xyz coordinates of the voxels, and after conversion based on the RGB values and opacity of the specified plurality of voxels. The voxel RGB value and opacity of the voxel structure are calculated. Further, the interpolation processing unit 29 can also specify a plurality of voxels having an integer value in the vicinity of the real value of the xyz coordinate to which a random number is added as the xyz coordinate.

The interpolation processing unit 29 multiplies the RGB values and opacity of the specified plurality of voxels by a weighting coefficient based on the difference (fraction) between the real value of the specified coordinate value and the neighboring integer value, and the converted voxel structure. Calculate the RGB value and opacity of the voxel of. Further, when the xyz coordinates of the voxels of the converted voxel structure are within the target area defined by the xyz coordinates obtained by adding random numbers, the interpolation processing unit 29 is based on the RGB values and opacity of the specified plurality of voxels. It is also possible to calculate the RGB value and opacity of the voxel of the converted voxel structure.

The processing of steps S111 and S112 can be expressed by the equation (28).

FIG. 16 is a schematic diagram showing an example of interpolation processing by the ray casting apparatus 100 of the present embodiment. In FIG. 16, the coordinates (xr, yr, zr) are random numbers within a predetermined range of the real value coordinates (xr, yr, zr) corresponding to the voxel structure V (x, y, z, c) before conversion. Is the coordinate value of the real value updated by adding and giving the fluctuation. In FIG. 16, it can be seen that the real numerical value (xr, yr, zr) is decomposed into an integer value (xi, yi, zi) and a fraction (wx, wy, wz). Then, around the coordinate values (xr, yr, zr), eight sets of integer values (xi, yi, zi), (xi + 1, yi, zi), (xi, yi + 1, zi), (xi + 1, yi + 1,) zi), (xi, yi, zi + 1), (xi + 1, yi, zi + 1), (xi, yi + 1, zi + 1), (xi + 1, yi + 1, zi + 1) exist.

As can be seen from Equation (28) and FIG. 16, in the interpolation of the boxel value, the smaller the fraction between the real value coordinate value and the integer value coordinate value, the more the boxel value of the boxel structure at the integer value coordinate value becomes. A large weighting factor is multiplied.

The image processing unit 20 determines whether or not the voxel values of all voxel structures V'(x, y, z', c) have been determined (S113), and all voxel structures V'(x, y). , Z', c) If the voxel values have not been determined (NO in S113), the processing after step S105 is continued, and the voxel values of all voxel structures V'(x, y, z', c) are set. If it is determined (YES in S113), the process ends.

The data of the three-dimensional voxel image of the present embodiment has a structure including voxel data of voxels corresponding to each pixel of a plurality of xy plane images captured at predetermined intervals along the z direction, and voxel data. Is data that has undergone a predetermined coordinate conversion, and the actual value of the voxel's xyz coordinate is calculated and calculated by performing the inverse conversion of the predetermined coordinate conversion on the voxel's xyz coordinate after the coordinate conversion. A plurality of voxels having an integer value in the vicinity of the real value of the xyz coordinate of the voxel as the xyz coordinate are identified, and the RGB of the voxel after the predetermined coordinate conversion is based on the RGB values and opacity of the identified voxels. It is the data for which the value and opacity are calculated, and when a virtual ray is irradiated from a predetermined viewpoint toward a predetermined coordinate-converted three-dimensional voxel image, the voxel RGB of the coordinate-converted voxel for each voxel on the line of sight. Processing to generate a rendered image as the pixel value of the pixel whose line of sight intersects the projection plane, using the RGB value calculated based on the cumulative sum of the cumulative brightness value based on the product of the value, opacity, and transmission intensity of the virtual light beam. Used to perform.

17 and 18 are flowcharts showing an example of the procedure of the ray casting process by the ray casting apparatus 100 of the present embodiment. The image processing unit 20 sets the light source Light (c) for color calculation (S121). Here, 0 ≦ Light (c) ≦ 255, c = 0 (R), 1 (G), and 2 (B). The image processing unit 20 sets the transmission intensity Trans = 1.0 of the virtual light beam and sets the cumulative luminance value Energy (c) = 0.0 (S122).

The image processing unit 20 sets z'= M (x, y, 1) (S123). That is, the value of z'is set to the maximum value z2 of the search control mask data. The image processing unit 20 determines whether or not the value of z'is smaller than 0 (S124), and if z'<0 (NO in S124), the opacity α = V'(x, y, z', 3) / 255 (S125). According to this equation, the opacity α is 0 ≦ α ≦ 1. When z ′ <0 (YES in S124), the image processing unit 20 performs the processing of step S135 described later. The image processing unit 20 determines whether or not the opacity α is 0 (S126), and when the opacity α = 0 (YES in S126), the image processing unit 20 performs the processing of step S133 described later.

When the opacity α is not 0 (NO in S126), that is, when the opacity α> 0, the image processing unit 20 determines whether or not the shadow processing is performed (S127). When there is a shadow process (YES in S127), the image processing unit 20 performs a shadow Shade calculation process (S128) and performs the process of step S130 described later. The shadow Shade calculation can be performed, for example, as follows.

The light source vector (Lx, Ly, Lz) can be specified as a unit vector with a parallel light source, for example, (Lx, Ly, Lz) = (0.57735, 0.57735, 0.57735). Further, the ambient light component is Ab. Here, 0 ≦ Ab ≦ 1. For example, Ab = 0.2 can be set.

The raycasting processing unit 25 calculates the difference between the opacity of the voxel of interest in the converted voxel structure and the opacity of each of the plurality of neighboring voxels located in the vicinity of the voxel of interest. The raycasting processing unit 25 calculates the gradient vector of the voxel of interest based on the logical distance between the voxel of interest and the neighboring voxel and the difference in opacity of each of the voxel of interest and the neighboring voxel.

The gradient vector (Gx, Gy, Gz) of the voxel (x, y, z') in the range of 0≤x≤Sx-1, 0≤y≤Sy-1, 0≤z'≤Sz'-1 is an equation. It can be calculated by (29), (30), and (31).

In the case of x = 0 or x = Sx-1 or y = 0 or y = Sy-1 or z'= 0 or z'= Sz'-1 (that is, the voxel at the boundary), the voxel (x + i, y + j). , Z'+ k) may not exist. In that case, V'(x + i, y + j, z'+ k, 3) = 0.

The raycasting processing unit 25 calculates the voxel shadow value Shade based on the calculated gradient vector and a predetermined light source vector. G = {Gx ² + Gy ² + Gz ² } ^1/2 (G is the square root of the sum of the square of Gx, the square of Gy, and the square of Gz). When G ≧ 1, the shadow value Shadow is given by the equation (32) when calculating the diffuse reflection component. Since the gradient vector also changes following the change of the line of sight, the specular reflection component in consideration of the line of sight vector can be calculated and added.

The raycasting processing unit 25 can correct the brightness value Shadow of the voxel based on the opacity α of the voxel. For example, the shadow value Shade can be corrected by multiplying the opacity α of the voxel. That is, the shadow value Shade can be calculated by the equation (33). When G <1, the shadow value Shade = 0.

When there is no shading process (NO in S127), the image processing unit 20 sets Shade = α (S129). The image processing unit 20 updates the cumulative luminance value (S130). That is, the cumulative luminance value of the three components 0 ≦ c ≦ 2 is updated by the formula Energy (c) = Energy (c) + Trans · α · Shadow · V ′ (x, y, z ′, c) / 255. .. The image processing unit 20 updates the transmission intensity (S131). That is, the transmission intensity is updated according to the equation Trans = Trans · (1-α).

The image processing unit 20 determines whether or not the opacity α = 1.0 or whether or not the transmission intensity Transs <predetermined value (S132). Here, the predetermined value may be 0 or a value close to 0 (for example, 0.001). When the opacity α = 1.0 and the transmission intensity Trans <predetermined value (NO in S132), the image processing unit 20 sets z ′ = z ′ -1 (S133), and the value of z ′ is M (NO). It is determined whether or not it is smaller than x, y, 0) (S134). When z'<M (x, y, 0) (YES in S134), the image processing unit 20 determines the RGB value in the two-dimensional coordinates (x, y) (S135). The RGB values in the two-dimensional coordinates (x, y) of the projection plane can be calculated by Image (x, y, c) = k ・ Energy (c) ・ Light (c). Here, k is an intensity magnification, and the RGB value can be corrected. The value of k can be, for example, 1.0, but is not limited to this. Light (c) is an RGB value of the light source.

When the opacity α = 1.0 or the transmission intensity Trans <predetermined value (YES in S132), the image processing unit 20 performs the process of step S135 without performing the processes of steps S133 and S134. .. If z'<M (x, y, 0) is not satisfied (NO in S134), the image processing unit 20 performs the processing after step S125.

The image processing unit 20 determines whether or not the RGB values of all the two-dimensional coordinates have been determined (S136), and if the RGB values of all the two-dimensional coordinates have not been determined (NO in S136), step S122 and subsequent steps. Process. When the RGB values of all the two-dimensional coordinates are determined (YES in S136), the image processing unit 20 ends the processing.

As described above, the ray casting processing unit 25 does not have the RGB value (V'(x, y, z', c)) of the voxel for each voxel within the range of the minimum value z1 and the maximum value z2 on the line of sight. The cumulative brightness value Energy (c) was cumulatively added (updated) by the product of the transparency α (V'(x, y, z', 3)), the shadow value Shadow calculated in the voxel, and the transmission intensity Trans of the virtual light beam. A rendered image can be generated using the RGB values calculated based on the values as the pixel values of the pixels whose line of sight intersects the projection plane.

The ray casting processing unit 25 determines that the box cell on the line of sight is outside the range of the minimum value z1 and the maximum value z2, the minimum value z1 and the maximum value z2 are the same value, or the minimum value z1 or the maximum value z2. If is not found, a predetermined pixel value can be given to the pixel whose line of sight intersects the projection plane to generate a rendered image. The predetermined pixel value may be any value as long as it can be given without performing ray casting processing, and may be, for example, a background color (Rb, Gb, Bb). As a result, the load related to the ray casting process can be reduced. Specifically, when the background color component is added to the equation of S134 and Back (0) = Rb, Back (1) = Gb, Back (2) = Bb, Image (x, y, c) = K. Energy (c) / Light (c) + Back (c) may be used.

Further, when performing the shadow calculation, the raycasting processing unit 25 sets and calculates a parallel light source in the same direction as the virtual light source when the virtual light ray is irradiated toward the voxel structure after conversion from a predetermined viewpoint. You can also. In that case, the light source vector (Lx, Ly, Lz) may be set to (Lx, Ly, Lz) = (0,0,1), and the mathematical formula (32) may be executed in the same manner. Then, for each voxel on the line of sight, the RGB value (V'(x, y, z', c)) of the voxel, the opacity α (V'(x, y, z', 3)) and the voxel are calculated. Cumulative brightness value by the product of the shadow value {Shade · V'(x, y, z', 3) / 255} corrected by multiplying the shade value Shadow by the opacity α and the transmission intensity Trans (c) of the virtual light beam. A rendered image can be generated using the RGB values calculated based on the cumulative addition (updated) of the Energy (c) as the pixel values of the pixels whose line of sight intersects the projection plane.

FIG. 19 is a block diagram showing another example of the configuration of the image processing unit 20 of the present embodiment. As shown in FIG. 19, the image processing unit 20 can be composed of a CPU 201, a ROM 202, a RAM 203, a recording medium reading unit 204, and the like. The computer program recorded on the recording medium 1 can be read by the recording medium reading unit 204 and stored in the RAM 203. By executing the computer program stored in the RAM 203 by the CPU 201, the image processing unit 20 has a gradation compression processing unit 21, a smoothing processing unit 22, a voxel structure generation unit 23, a coordinate conversion processing unit 24, and a ray casting process. The processing performed by the unit 25, the search control mask generation unit 26, the coordinate value calculation unit 27, the random number generation unit 28, the interpolation processing unit 29, and the ROI clipping processing unit 30 can be executed. The CPU 201 is a general term for a processor, a processor device, a processing unit, and the like, and may be configured by one or a plurality. Further, the computer program can be downloaded via a network such as the Internet instead of the configuration of reading by the recording medium reading unit 204.

Next, the pattern of moire occurrence will be described. The reason why moire occurs in a volume rendered image is that tomographic images are laminated and combined. As a type of moire, there is a curved pattern as the first type. Further, as the second pattern, there is a stripe / grid pattern. It should be noted that the curved pattern and the striped / grid-like pattern are not strictly separated, and in the volume rendering image, both are present in a complex mixture. Hereinafter, they will be described separately for convenience of explanation. First, the first curved pattern will be described.

FIG. 20 is a schematic diagram showing an example of how curved moire is generated. The cause of occurrence is a pattern that does not have perceptible regularity that originally exists naturally in a tomographic image (slice image). By performing coordinate transformation (rotation around the xyz axis) of a plurality of tomographic images having such a pattern, the above-mentioned non-regular patterns existing in each tomographic image are deviated and synthesized. The patterns are linearly connected, resulting in a curved striped pattern in the volume rendered image, which is perceived as a perceptible moire.

FIG. 21 is an explanatory diagram showing a first example of a moire countermeasure result by the ray casting device 100 of the present embodiment. FIG. 21A is a comparative example and shows a case where the moire countermeasure of the present embodiment is not implemented. FIG. 21B shows a case where each tomographic image is smoothed in a 2D space (xy plane) as described in FIG. 6 by the smoothing processing unit 22 of the raycasting apparatus 100 of the present embodiment. .. In FIG. 21A, a curved striped pattern is generated, but in FIG. 21B, it can be seen that the striped pattern has disappeared.

FIG. 22 is an explanatory diagram showing a second example of the moire countermeasure result by the ray casting device 100 of the present embodiment. FIG. 22A is a comparative example and shows a case where the moire countermeasure of the present embodiment is not implemented. FIG. 22B shows a case where interpolation processing at the time of coordinate conversion is performed by the ray casting apparatus 100 of the present embodiment. The interpolation process at the time of coordinate conversion includes the processes such as steps S111 and S112 in FIG. 15 described above. In FIG. 22A, a curved striped pattern is generated, but in FIG. 22B, it can be seen that the striped pattern has disappeared.

Next, a stripe / grid pattern, which is the second pattern of moiré, will be described.

FIG. 23 is a schematic diagram showing an example of how stripe-lattice moire is generated. The characteristic of the striped / grid pattern is that unlike the first pattern, the moire source is not in the tomographic image itself, but in the artificial pattern generated by calculation in the process of ray casting. By performing coordinate conversion (rotation around the xyz axis) of each tomographic image, the fraction is rounded down to an integer in the process of real number operation called coordinate conversion, and rounding error is performed in the process of referring to the boxel value before conversion with the xyz coordinate of the integer value. A periodic pattern is generated by the accumulation of, and it is added to the tomographic image after the coordinate conversion. The periodic pattern added to each tomographic image is so weak that it cannot be perceived. By strengthening each other, so-called interference fringes are generated, and the level becomes perceptible as stripes and grid-like moire.

FIG. 24 is an explanatory diagram showing an example of the occurrence of rounding error during coordinate conversion. In FIG. 24, voxels are represented in two dimensions for convenience. The voxels before coordinate conversion are represented by V11, V12, ..., V49 with the horizontal direction as the x-axis direction and the vertical direction as the z-axis direction. Due to the coordinate transformation that is rotated around the z-axis, the rounding error that accompanies the integerization of the x-coordinate of the real value is accumulated in the x-coordinate value to be referenced. For example, the calculated x-coordinate value of the real number is 1.3, 2.6, 3.9, 5.2, 6.5, 7.8 from left to right.

For example, when the calculated x-coordinate value of the real number is 1.3, the voxel after the coordinate conversion refers to V11 and V12 before the coordinate conversion by V11 × 0.7 + V12 × 0.3, and 0.3. It can be obtained by interpolating while multiplying the weights based on the fraction. Further, when the calculated x-coordinate value of the real number is 2.6, the voxel after the coordinate conversion refers to V12 and V13 before the coordinate conversion by V12 × 0.4 + V13 × 0.6, and is 0.6. It can be obtained by interpolating while multiplying the weights based on the fraction. Further, when the calculated x-coordinate value of the real number is 3.9, the voxel after the coordinate conversion refers to V13 and V14 by V13 × 0.1 + V14 × 0.9 before the coordinate conversion, and 0.9. It can be obtained by interpolating while multiplying the weights based on the fraction. Further, when the calculated x-coordinate value of the real number is 5.2, the voxel after the coordinate conversion refers to V15 and V16 before the coordinate conversion by V15 × 0.8 + V16 × 0.2, and 0.2. It can be obtained by interpolating while multiplying the weights based on the fraction.

In this way, the fourth voxel from the left after the coordinate conversion is not interpolated by the voxels V14 and V15 before the coordinate conversion, but is interpolated by V15 and V16, and the third voxel and 4 from the left after the coordinate conversion. There is a step in the voxel value with the third voxel. However, this step is not a perceptible level in a single layer (tomographic image) formed by V11 to V19. However, as shown in FIG. 24, when similar steps occur at the same location in the three layers formed by V21 to V49, as a result, they are perceived as stripes and grid-like moire along the z-axis direction. It will be possible.

As a countermeasure to the above-mentioned problem, we propose a method to give fluctuation to the rounding error. FIG. 25 is an explanatory diagram showing an example of a state of rounding error dispersion by the ray casting apparatus 100 of the present embodiment. The rounding error variance includes the process of step S104 in FIG. 14 or the process of step S106 in FIG. 15 described above. As shown in FIG. 25, the voxels before the coordinate conversion are the same as those in FIG. 24. Since the x-coordinate value to be referred to has a fluctuation added to the voxel coordinate value by a random number, or a fluctuation is added to the voxel coordinate value and the ROI set coordinate by a random number, the x-coordinate value to be referred to is for each z-coordinate value. Different to. For example, the second x-coordinate value from the top is 2.1, 3.2, 4.9, 6.1, 7.1, 8.0 from left to right. The third x-coordinate value from the top is 1.2, 2.4, 3.5, 4.8, 5.9, and 6.9 from left to right. The fourth x-coordinate value from the top is 2.3, 3.1, 4.4, 5.9, 6.5, and 7.8 from left to right. Then, in the first layer and the second layer, a step remains between the third and fourth voxels from the left as in FIG. 23, but in the third layer and the fourth layer, the step disappears and these overlap. Then, the stripe / grid-like moire generated along the z-axis direction cannot be perceived.

As described above, according to the present embodiment, the rounding error that occurs when the coordinate value of the real value is converted into an integer is dispersed (actually, it is dispersed in the xyz direction), and the periodic unevenness is eliminated. Since there is no difference in stripes or grids in the converted boxel values, stripes and grids of moire are reduced.

FIG. 26 is an explanatory diagram showing a third example of the moire countermeasure result by the ray casting device 100 of the present embodiment. FIG. 26A is a comparative example and shows a case where the moire countermeasure of the present embodiment is not implemented. FIG. 26B shows a case where fluctuations are given to the coordinate values of real numbers of voxels by random numbers by the ray casting apparatus 100 of the present embodiment. The process of giving fluctuations to the real coordinate values of voxels by random numbers includes the processes of step S104 and the like in FIG. 14 described above, the calculations by the equations (16) to (18), and the like. In FIG. 26A, a grid-like striped pattern is generated due to a rounding error when converting the coordinate value of the real number calculated by the coordinate conversion process into an integer, but in FIG. 26B, it can be seen that the striped pattern disappears.

FIG. 27 is an explanatory diagram showing a fourth example of the moire countermeasure result by the ray casting device 100 of the present embodiment. FIG. 27A is a comparative example and shows a case where the moire countermeasure of the present embodiment is not implemented. FIG. 27B shows a case where the ray casting device 100 of the present embodiment gives fluctuations to the coordinate values of real numbers of voxels by random numbers. The process of giving fluctuations to the coordinate values of the real numbers of voxels by random numbers includes the processes of step S104 and the like in FIG. 14 described above, the calculations by the equations (16) to (18), and the like. In FIG. 27A, a grid-like striped pattern is generated due to a rounding error when converting the coordinate value of the real number calculated by the coordinate conversion process into an integer, but in FIG. 27B, it can be seen that the striped pattern disappears.

FIG. 28 is an explanatory diagram showing a fifth example of the moire countermeasure result by the ray casting device 100 of the present embodiment. FIG. 28A is a comparative example and shows a case where the moire countermeasure of the present embodiment is not implemented. FIG. 28B shows a case where the ray casting apparatus 100 of the present embodiment gives fluctuations to the clipping coordinate values by random numbers during the ROI clipping process. The process of giving fluctuation to the clipping coordinate value by a random number includes the process of step S106 and the like in FIG. 15 described above, the calculation by the equations (22) to (27), and the like. In FIG. 28A, a vertical striped pattern is generated on the ROI clipping cut surface due to a rounding error when the coordinate values associated with the coordinate conversion are converted into integers, but in FIG. 28B, it can be seen that the striped pattern disappears.

FIG. 29 is an explanatory diagram showing a sixth example of the moire countermeasure result by the ray casting device 100 of the present embodiment. FIG. 29A is a comparative example and shows a case where the moire countermeasure of the present embodiment is not implemented. FIG. 29B shows a case where the ray casting apparatus 100 of the present embodiment performs a process of correcting the shadow value using the opacity α. The process of correcting the shadow value using the opacity α includes the process of step S128 and the like in FIG. 17 described above, the calculation by the equation (33), and the like. In FIG. 29A, since the contrast between light and dark of the shadow is insufficient, a grid-like striped pattern is generated and moire can be perceived. However, in FIG. 29B, the calculated shadow value Shadow is multiplied by the opacity α. You can see that the striped pattern has disappeared.

As described above, according to the present embodiment, it is possible to suppress the generation of various moire (curve-shaped moire, stripe-shaped moire) while suppressing deterioration of the image quality of the rendered image.

The computer program of the present embodiment is a computer program for causing a computer to perform a volume rendering process using a raycasting method, and is a plurality of xy planes imaged by the computer at predetermined intervals along the z direction. A process of acquiring image data of an image, a process of smoothing the plurality of xy plane images based on the pixel values of arbitrary attention pixels of each of the plurality of xy plane images and pixels in the vicinity of the attention pixels, and RGB. A process of generating a three-dimensional voxel image composed of voxels corresponding to each pixel of each of the plurality of xy plane images whose values and opacity are determined and smoothed, and predetermined coordinates with respect to the three-dimensional voxel image. A process of performing conversion to generate a converted voxel image, and when a virtual ray is irradiated toward the converted voxel image from a predetermined viewpoint, the RGB value and opacity of the voxel and the above-mentioned for each voxel on the line of sight. The processing of generating a rendered image is executed by using the RGB value calculated based on the cumulative sum of the cumulative brightness values obtained by the product of the transmission intensities of the virtual rays as the pixel values of the pixels whose line of sight intersects the projection plane.

The computer program of the present embodiment tells the computer the pixel values of any pixel of interest on one xy plane image of each of the plurality of xy plane images and neighboring pixels of the pixel of interest in the one xy plane image. Based on this, the process of smoothing the one xy plane image is executed.

The computer program of the present embodiment performs a process of calculating the real value of the xyz coordinate of the voxel by performing the inverse conversion of the predetermined coordinate conversion on the xyz coordinate of the voxel of the converted voxel image on the computer. The process of specifying a plurality of voxels whose xyz coordinates are integer values in the vicinity of the real values of the calculated xyz coordinates of the voxels, and the voxels of the converted voxel image based on the RGB values and opacity of the specified multiple voxels. The process of calculating the RGB value and the opacity is executed.

In the computer program of the present embodiment, a process of adding a random number within a predetermined range to the real value of the xyz coordinate of the boxel and an integer value in the vicinity of the real value of the xyz coordinate obtained by adding the random number are used as the xyz coordinate. To execute the process of identifying multiple box cells.

In the computer program of the present embodiment, the computer is subjected to a process of specifying the difference between the real value of the xyz coordinate of the voxel and the integer value in the vicinity of the real value, and the RGB values and non-RGB values of the specified plurality of voxels. The process of multiplying the transparency by a weighting coefficient based on the specified difference to calculate the RGB value and opacity of the voxel of the converted voxel image is executed.

The computer program of the present embodiment has a process of setting xyz coordinates that define a target area of volume rendering processing in a computer, a process of adding a random number within a predetermined range to the set xyz coordinates, and a voxel image after conversion. When the voxel xyz coordinates of the voxel are within the target area defined by the xyz coordinates obtained by adding random numbers, the voxel RGB values of the converted voxel image and the voxel RGB values of the converted voxel image are based on the RGB values and opacity of the specified plurality of voxels. The process of calculating the opacity is executed.

The computer program of the present embodiment specifies to the computer the minimum value and the maximum value on the line of sight where the voxel opacity of the converted voxel image is larger than a predetermined value for each xy coordinate of the converted voxel image. The line of sight is the projection plane of the processing and the value obtained by cumulatively adding the product of the RGB value and opacity of the voxel and the transmission intensity of the virtual light beam for each voxel within the range of the minimum value and the maximum value on the line of sight. The process of generating a rendered image as the pixel value of the pixel intersecting with is executed.

In the computer program of the present embodiment, when the voxels on the line of sight are outside the range of the minimum value and the maximum value, a predetermined pixel value is given to the pixels whose line of sight intersects the projection plane. Execute the process to generate the rendered image.

In the computer program of the present embodiment, a process of specifying to the computer the minimum value and the maximum value of the pixel of the predetermined xy plane image among the plurality of xy plane images, and the upper limit value smaller than the specified maximum value and the upper limit value and the maximum value are specified. A process of calculating a lower limit value larger than the specified minimum value, a process of compressing the pixel value of each pixel of the plurality of xy plane images within the range of the upper limit value and the lower limit value, and a process of compressing the plurality of xy planes. The process of smoothing the plurality of xy plane images is executed based on the pixel values of any pixel of interest in each image and pixels in the vicinity of the pixel of interest.

The computer program of the present embodiment has a process of acquiring color map data in which RGB values and opacity are defined in association with pixel values, and a process of acquiring each pixel of each of the smoothed plurality of xy plane images. The process of generating a three-dimensional voxel image is executed by setting the RGB value and the opacity in association with the RGB value and the opacity of the color map data corresponding to the pixel value.

The computer program of the present embodiment includes, in the computer, the rotation angle of the xyz axis, the offset value in the xyz axis direction, the enlargement or reduction magnification in the xyz axis direction, the scaling factor in the z axis direction, and the distance from the gazing point to the viewpoint. The process of acquiring the parameters of the predetermined coordinate conversion and the process of performing the predetermined coordinate conversion using the acquired patternometer on the generated three-dimensional boxel image to generate the converted boxel image are executed. ..

The computer program of the present embodiment causes a computer to calculate a gradient vector relating to voxel opacity of the converted voxel image, and calculates a voxel shading value based on the calculated gradient vector and a predetermined light source vector. Processing, processing to correct the shadow value of the voxel based on the opacity of the voxel, and processing of irradiating a virtual light beam from a predetermined viewpoint toward the converted voxel image, the voxel on the line of sight Execution of processing to generate a rendered image as the pixel value of the pixel whose line of sight intersects the projection plane, with the cumulative sum of the RGB value of the voxel, the opacity and the corrected shadow value, and the transmission intensity of the virtual light beam. Let me.

The computer program of the present embodiment causes the computer to calculate the difference between the opacity of the voxel of interest in the converted voxel image and the opacity of each of the plurality of neighboring voxels located in the vicinity of the voxel of interest. The process of calculating the gradient vector of the voxel of interest is executed based on the distance between the voxel of interest and the voxel in the vicinity and the difference in opacity between the voxel of interest and the voxel in the vicinity.

The image processing device of the present embodiment is an image processing device that performs volume rendering processing using a ray casting method, and acquires image data of a plurality of xy plane images captured at predetermined intervals along the z direction. An acquisition unit, a smoothing unit that smoothes the plurality of xy plane images based on the pixel values of an arbitrary pixel of interest in each of the plurality of xy plane images and pixels in the vicinity of the pixel of interest, an RGB value, and opacity. Is defined, and a boxel image generation unit that generates a three-dimensional boxel image composed of a boxel corresponding to each pixel of each of the plurality of xy plane images smoothed by the smoothing unit, and a boxel image generation unit that generates the image. When the converted boxel image generator that performs predetermined coordinate conversion on the three-dimensional boxel image to generate the converted boxel image and the virtual light beam from a predetermined viewpoint are irradiated toward the converted boxel image. The line of sight intersects the projection plane with an RGB value calculated based on a value obtained by cumulatively adding the cumulative brightness value by the product of the RGB value and opacity of the boxel and the transmission intensity of the virtual light beam for each boxel on the line of sight. It includes a rendered image generation unit that generates a rendered image as a pixel value of a pixel.

The image processing method of the present embodiment is an image processing method using an image processing device that performs volume rendering processing using a raycasting method, and is an image of a plurality of xy plane images captured at predetermined intervals along the z direction. Data is acquired, the plurality of xy plane images are smoothed based on the pixel values of an arbitrary pixel of interest in each of the plurality of xy plane images and pixels in the vicinity of the pixels of interest, and RGB values and opacity are determined. A three-dimensional voxel image composed of voxels corresponding to each pixel of each of the plurality of smoothed xy plane images is generated, and a predetermined coordinate conversion is performed on the generated three-dimensional voxel image to perform a converted voxel. When an image is generated and a virtual ray is irradiated from a predetermined viewpoint toward the converted voxel image, it is accumulated by the product of the RGB value and opacity of the voxel and the transmission intensity of the virtual ray for each voxel on the line of sight. A rendered image is generated using the RGB value calculated based on the cumulative sum of the brightness values as the pixel value of the pixel whose line of sight intersects the projection surface.

The image data of the present embodiment is data of a three-dimensional voxel image, which is voxel data of voxels corresponding to each pixel of a plurality of xy plane images captured at predetermined intervals along the z direction. , The plurality of xy plane images have a smoothed data structure based on the pixel values of any of the pixels of interest in each of the plurality of xy plane images and the pixels in the vicinity of the pixels of interest, and the RGB values and opacity are high. It is used to perform a process of generating a three-dimensional voxel image composed of voxels corresponding to each pixel of each of the plurality of defined and smoothed xy plane images.

The image data of the present embodiment is data of a three-dimensional voxel image, which is predetermined with respect to voxel data of voxels corresponding to each pixel of a plurality of xy plane images captured at predetermined intervals along the z direction. The voxel that has undergone the coordinate conversion of A plurality of voxels having xyz coordinates as integer values in the vicinity of the real value of the xyz coordinates of Has a calculated data structure, and when a virtual ray is emitted from a predetermined viewpoint toward the three-dimensional voxel image after the predetermined coordinate conversion, the RGB of the voxel after the coordinate conversion is applied to each voxel on the line of sight. A rendered image is generated using the RGB value calculated based on the cumulative sum of the value, the opacity, and the transmission intensity of the virtual light beam as the pixel value of the pixel whose line of sight intersects the projection plane. Used to perform processing.

10 Input unit 11 Storage unit 12 Output unit 20 Image processing unit 21 Gradation compression processing unit 22 Smoothing processing unit 23 Voxel structure generation unit 24 Coordinate conversion processing unit 25 Raycasting processing unit 26 Search control mask generation unit 27 Coordinate value calculation Part 28 Random number generation part 29 Interpolation processing part 30 ROI clipping processing part 100 Raycasting device 201 CPU
202 ROM
203 RAM
204 Recording medium reader

Claims (13)

A computer program for causing a computer to perform volume rendering processing using the ray casting method.
On the computer
Processing to acquire image data of a plurality of xy plane images taken at predetermined intervals along the z direction, and
A process of smoothing the plurality of xy plane images based on the pixel values of any pixel of interest in each of the plurality of xy plane images and pixels in the vicinity of the pixels of interest.
A process of generating a three-dimensional voxel image composed of voxels corresponding to each pixel of each of the plurality of xy plane images whose RGB values and opacity are determined and smoothed.
A process of performing predetermined coordinate transformation on the three-dimensional voxel image to generate a converted voxel image, and
When a virtual ray is irradiated toward the converted voxel image from a predetermined viewpoint, the cumulative brightness value is cumulatively added by the product of the RGB value and opacity of the voxel and the transmission intensity of the virtual ray for each voxel on the line of sight. A process of generating a rendered image using the RGB value calculated based on the calculated value as the pixel value of the pixel whose line of sight intersects the projection plane.
A process of calculating the real value of the xyz coordinate of the voxel by performing the inverse conversion of the predetermined coordinate transformation on the xyz coordinate of the voxel of the converted voxel image.
Processing to specify multiple voxels whose xyz coordinates are integer values in the vicinity of the real value of the calculated xyz coordinates of the voxels.
A computer program that calculates the RGB values and opacity of voxels in the converted voxel image based on the RGB values and opacity of a plurality of specified voxels. On the computer
The one xy plane image is smoothed based on the pixel values of any pixel of interest on one xy plane image of each of the plurality of xy plane images and neighboring pixels in the one xy plane image of the attention pixel. The computer program according to claim 1, wherein the processing is executed. On the computer
The process of adding a random number within a predetermined range to the real value of the xyz coordinates of the voxel, and
The computer program according to claim 1 or 2, wherein a process of specifying a plurality of voxels having an integer value in the vicinity of the real value of the xyz coordinate to which a random number is added as the xyz coordinate is executed. On the computer
Processing to specify the difference between the real value of the xyz coordinate of the voxel and the integer value in the vicinity of the real value.
From claim 1, the process of multiplying the RGB values and opacity of the specified plurality of voxels by the weighting coefficient based on the specified difference to calculate the RGB values and opacity of the voxels of the converted voxel image is executed. The computer program according to any one of claims 3. On the computer
The process of setting the xyz coordinates that define the target area of the volume rendering process, and
Processing to add random numbers within a predetermined range to the set xyz coordinates,
When the xyz coordinates of the voxels of the converted voxel image are within the target area defined by the xyz coordinates obtained by adding random numbers, the converted voxel image of the converted voxels image is based on the RGB values and opacity of the specified plurality of voxels. The computer program according to any one of claims 1 to 4, which executes a process of calculating voxel RGB values and opacity. A computer program for causing a computer to perform volume rendering processing using the ray casting method.
On the computer
Processing to acquire image data of a plurality of xy plane images taken at predetermined intervals along the z direction, and
A process of smoothing the plurality of xy plane images based on the pixel values of any pixel of interest in each of the plurality of xy plane images and pixels in the vicinity of the pixels of interest.
A process of generating a three-dimensional voxel image composed of voxels corresponding to each pixel of each of the plurality of xy plane images whose RGB values and opacity are determined and smoothed.
A process of performing predetermined coordinate transformation on the three-dimensional voxel image to generate a converted voxel image, and
When a virtual ray is irradiated toward the converted voxel image from a predetermined viewpoint, the cumulative brightness value is cumulatively added by the product of the RGB value and opacity of the voxel and the transmission intensity of the virtual ray for each voxel on the line of sight. A process of generating a rendered image using the RGB value calculated based on the calculated value as the pixel value of the pixel whose line of sight intersects the projection plane.
A process of specifying a minimum value and a maximum value on the line of sight in which the voxel opacity of the converted voxel image is larger than a predetermined value for each xy coordinate of the converted voxel image.
The line of sight intersects the projection plane with the cumulative sum of the RGB values and opacity of the voxels and the product of the transmission intensity of the virtual light beam for each voxel within the range of the minimum value and the maximum value on the line of sight. Processing to generate a rendered image as the pixel value of a pixel,
When the voxel on the line of sight is out of the range between the minimum value and the maximum value, a computer program that assigns a predetermined pixel value to a pixel whose line of sight intersects the projection plane and executes a process of generating a rendered image. A computer program for causing a computer to perform volume rendering processing using the ray casting method.
On the computer
Processing to acquire image data of a plurality of xy plane images taken at predetermined intervals along the z direction, and
A process of smoothing the plurality of xy plane images based on the pixel values of any pixel of interest in each of the plurality of xy plane images and pixels in the vicinity of the pixels of interest.
A process of generating a three-dimensional voxel image composed of voxels corresponding to each pixel of each of the plurality of xy plane images whose RGB values and opacity are determined and smoothed.
A process of performing predetermined coordinate transformation on the three-dimensional voxel image to generate a converted voxel image, and
When a virtual ray is irradiated toward the converted voxel image from a predetermined viewpoint, the cumulative brightness value is cumulatively added by the product of the RGB value and opacity of the voxel and the transmission intensity of the virtual ray for each voxel on the line of sight. A process of generating a rendered image using the RGB value calculated based on the calculated value as the pixel value of the pixel whose line of sight intersects the projection plane.
A process of specifying a minimum value and a maximum value of pixels of a predetermined xy plane image among the plurality of xy plane images, and a process of specifying the minimum value and the maximum value.
Processing to calculate the upper limit value smaller than the specified maximum value and the lower limit value larger than the specified minimum value,
A process of compressing the pixel value of each pixel of the plurality of xy plane images within the range of the upper limit value and the lower limit value, and
A computer program that executes a process of smoothing the plurality of xy plane images based on the pixel values of an arbitrary pixel of interest of each of the compressed plurality of xy plane images and pixels in the vicinity of the pixels of interest. On the computer
The process of acquiring color map data in which RGB values and opacity are defined in association with pixel values, and
A three-dimensional voxel image is generated by setting RGB values and opacity in association with the RGB values and opacity of the color map data corresponding to the pixel values of each pixel of each of the smoothed plurality of xy plane images. The computer program according to any one of claims 1 to 7, wherein the processing and execution are performed. On the computer
Processing to acquire the parameters of the predetermined coordinate conversion including the rotation angle of the xyz axis, the offset value in the xyz axis direction, the enlargement or reduction magnification in the xyz axis direction, the variable magnification in the z axis direction, and the distance from the gazing point to the viewpoint. ,
One of claims 1 to 8, wherein the generated three-dimensional voxel image is subjected to the predetermined coordinate transformation using the acquired parameters to generate a converted voxel image. The computer program described in section. On the computer
The process of calculating the gradient vector regarding the voxel opacity of the converted voxel image and
Processing to calculate the shadow value of the voxel based on the calculated gradient vector and a predetermined light source vector,
The process of correcting the shadow value of the voxel based on the opacity of the voxel, and
When a virtual ray is irradiated toward the converted voxel image from a predetermined viewpoint, the product of the RGB value, opacity and corrected shadow value of the voxel and the transmission intensity of the virtual ray is accumulated for each voxel on the line of sight. The computer program according to any one of claims 1 to 8, wherein the added value is used as a pixel value of a pixel whose line of sight intersects the projection surface to generate a rendered image. On the computer
The process of calculating the difference between the opacity of the voxel of interest in the converted voxel image and the opacity of each of the plurality of neighboring voxels located in the vicinity of the voxel of interest.
The computer program according to claim 10, wherein the process of calculating the gradient vector of the voxel of interest is executed based on the distance between the voxel of interest and the voxel in the vicinity and the difference in opacity between the voxel of interest and the voxel in the vicinity. An image processing device that performs volume rendering processing using the ray casting method.
An acquisition unit that acquires image data of a plurality of xy plane images captured at predetermined intervals along the z direction, and an acquisition unit.
A smoothing unit that smoothes the plurality of xy plane images based on the pixel values of any pixel of interest in each of the plurality of xy plane images and pixels in the vicinity of the pixels of interest.
A voxel image generation unit that generates a three-dimensional voxel image in which RGB values and opacity are determined and are composed of voxels corresponding to each pixel of each of the plurality of xy plane images smoothed by the smoothing unit.
A post-conversion voxel image generation unit that generates a post-conversion voxel image by performing predetermined coordinate conversion on the three-dimensional voxel image generated by the voxel image generation unit.
When a virtual ray is irradiated toward the converted voxel image from a predetermined viewpoint, the cumulative brightness value is cumulatively added by the product of the RGB value and opacity of the voxel and the transmission intensity of the virtual ray for each voxel on the line of sight. A rendered image generation unit that generates a rendered image as a pixel value of a pixel whose line of sight intersects the projection plane, using an RGB value calculated based on the calculated value.
A coordinate value calculation unit that calculates the real value of the xyz coordinates of the voxel by performing the inverse conversion of the predetermined coordinate transformation on the xyz coordinates of the voxel of the converted voxel image.
Multiple voxels whose xyz coordinates are integer values in the vicinity of the real value of the calculated xyz coordinates of the voxels are specified, and the RGB values of the voxels of the converted voxel image are based on the RGB values and opacity of the specified multiple voxels. An image processing device including an interpolation processing unit for calculating opacity. It is an image processing method by an image processing device that performs volume rendering processing using the ray casting method.
Image data of a plurality of xy plane images taken at predetermined intervals along the z direction are acquired, and the image data is acquired.
The plurality of xy plane images are smoothed based on the pixel values of any pixel of interest in each of the plurality of xy plane images and pixels in the vicinity of the pixels of interest.
A three-dimensional voxel image composed of voxels corresponding to each pixel of each of the plurality of xy plane images whose RGB values and opacity are determined and smoothed is generated.
A predetermined coordinate transformation is performed on the generated 3D voxel image to generate a converted voxel image.
When a virtual ray is irradiated toward the converted voxel image from a predetermined viewpoint, the cumulative brightness value is cumulatively added by the product of the RGB value and opacity of the voxel and the transmission intensity of the virtual ray for each voxel on the line of sight. A rendered image is generated by using the RGB value calculated based on the calculated value as the pixel value of the pixel whose line of sight intersects the projection plane.
The xyz coordinate of the voxel of the converted voxel image is subjected to the inverse conversion of the predetermined coordinate transformation to calculate the real value of the xyz coordinate of the voxel.
Identify multiple voxels whose xyz coordinates are integer values in the vicinity of the real value of the calculated xyz coordinates of the voxels.
An image processing method for calculating the RGB values and opacity of voxels in the converted voxel image based on the RGB values and opacity of a plurality of specified voxels.

Images