JP6436258B1 - Computer program, image processing apparatus, and image processing method - Google Patents

Abstract

A computer program, an image processing apparatus, an image processing method, and voxel data capable of suppressing moire are provided.
A computer program causes a computer to acquire image data of a plurality of xy plane images captured at predetermined intervals along the z-axis direction, and to calculate RGB values and pixel values for the pixel values of the xy plane image. A process for generating a provisional voxel structure composed of voxels having RGB values and opacity corresponding to each pixel of each of a plurality of xy plane images, and a voxel of the provisional voxel structure. Processing to generate a voxel structure by correcting the opacity of voxels based on the opacity of neighboring voxels located in the vicinity of the voxel, and processing to generate a rendering image by performing volume rendering processing on the voxel structure And execute.
[Selection] Figure 13

Description

The present disclosure relates to a computer program, an image processing apparatus and an image processing how.

  In the medical field, images are taken continuously at a predetermined slice interval by a medical image diagnostic device such as CT (Computed Tomography), MRI (Magnetic Resonance Imaging), PET (Positron Emission Tomography), and DICOM (Digital Imaging and Communication in Medicine). ) Based on a two-dimensional tomographic image stored in the form, a technique for visualizing an observation object such as an organ three-dimensionally and using it for image diagnosis support, surgery, or treatment support has been developed.

  As a typical technique for displaying a three-dimensional image of the internal tissue of an observation object, a volume rendering technique for drawing an image having a three-dimensional structure from three-dimensional digital data (volume data) of an object is known. As volume rendering methods, ray casting and texture-based rendering are known.

  In Patent Document 1, a three-dimensional voxel space is constructed by stacking a plurality of tomographic images, a ray (virtual ray) is irradiated from an arbitrary viewpoint direction to the three-dimensional voxel space, and at a sampling point on the line of sight Ray casting is disclosed in which color and luminance values are cumulatively calculated along the line of sight. 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 irradiated in the Z-axis direction to the three-dimensional voxel space after the coordinate transformation. The projection result is calculated as a pixel value on the projection surface (screen).

  Patent Document 2 discloses texture-based rendering that uses the polygon display function of graphics hardware by expressing the entire volume data as a stack of slice data using the texture mapping function and the alpha blending function. ing.

  In Patent Document 3, a conversion table (also referred to as a color map) in which a relationship between a voxel value of volume data, RGB color, and opacity (opacity) is defined by a user is applied to volume rendering. Alternatively, a technique is disclosed in which a lesion site is colored with an arbitrary color, or an organ or a lesion site that is hidden with a transparent organ in front is extracted and drawn.

Japanese Patent Laid-Open No. 11-175743 Japanese Patent No. 3589654 JP 2009-160306 A

  In a tomographic image, there is a pattern that does not have perceptible regularity that exists naturally. By performing coordinate transformation (rotation around the xyz axis) disclosed in Patent Document 1 and Patent Document 2 for a plurality of tomographic images, the above-described patterns existing in the tomographic image are shifted and synthesized. The patterns are connected in a line, and as a result, a curved stripe pattern is formed in the volume rendering image, resulting in moire. Therefore, smoothing processing can be considered in which the voxel values of the volume data change continuously by averaging the voxel values corresponding to each coordinate value of the volume data. However, if a color map is applied, it is calculated by the smoothing processing. With the RGB color and opacity converted based on the new voxel value, the change in RGB color and opacity becomes discontinuous, and moire such as an unintended curved striped pattern or a spot-like color shift occurs.

The present disclosure has been made in view of such circumstances, reveal a computer program it is possible to suppress moire, image processing apparatus and how the.

The present application includes a plurality of means for solving the above-mentioned problem. To give an example, the computer program stores image data of a plurality of xy plane images captured at predetermined intervals along the z-axis direction. Based on a color map that associates the relationship between the value to be acquired, the numerical value in the predetermined range, the RGB value, and the opacity, the RGB value and the opacity corresponding to the value of each pixel of each of the plurality of xy planar images. Processing to generate a provisional voxel structure composed of voxels, and for the voxels of the provisional voxel structure, the opacity of the voxels is corrected based on the opacity of neighboring voxels located in the vicinity of the voxel. Processing to generate a voxel structure, and volume rendering processing is performed on the voxel structure to generate a rendering image. To execute and management.

  According to the present disclosure, moire can be suppressed.

It is a block diagram which shows an example of a structure of the 3D texture mapping apparatus as an image processing apparatus of this Embodiment. It is a schematic diagram which shows an example of the tomographic image of a 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 process by the 3D texture mapping apparatus of this Embodiment. FIG. 4 is an example of the color map shown in FIG. 3 cut out. As a comparative example, it is a schematic diagram illustrating an example in which smoothing processing is performed on a voxel value of a DICOM image before passing through the above-described color map. It is a schematic diagram which shows an example which performed the smoothing process with respect to the opacity of the temporary voxel structure produced through the color map of this Embodiment. It is a schematic diagram which shows an example which performed the smoothing process with respect to the RGB value and opacity of the provisional voxel structure produced through the color map of this Embodiment. It is a schematic diagram which shows an example of the algorithm of the smoothing process of the voxel of this Embodiment. It is a schematic diagram which shows the concept of 3D texture mapping by the 3D texture mapping apparatus of this Embodiment. It is a flowchart which shows an example of the procedure of the 3D texture mapping process by the 3D texture mapping apparatus of this Embodiment. It is a flowchart which shows an example of the procedure of the 3D texture mapping process by the 3D texture mapping apparatus of this Embodiment. It is explanatory drawing which shows an example of the transparency process of the voxel outside a ROI clipping area | region by the 3D texture mapping apparatus of this Embodiment. It is explanatory drawing which shows an example of the dummy transparent image addition process by the 3D texture mapping apparatus of this Embodiment. It is explanatory drawing which shows an example of the transparency process of the voxel of the outer edge part of the tomographic image by the 3D texture mapping apparatus of this Embodiment. It is explanatory drawing which shows an example of the projection screen setting by the 3D texture mapping apparatus of this Embodiment. It is a flowchart which shows an example of the procedure of the color interpolation process of the transparent voxel by the 3D texture mapping apparatus of this Embodiment. It is a flowchart which shows an example of the procedure of the rendering process by the 3D texture mapping apparatus of this Embodiment. It is explanatory drawing which shows an example of the matching method with 3D texture mapping. It is a block diagram which shows the other example of a structure of the image process part of this Embodiment. It is a schematic diagram which shows an example of a mode that a curved moire occurs. It is explanatory drawing which shows the 1st example of the moire countermeasure result by the 3D texture mapping apparatus of this Embodiment. It is explanatory drawing which shows the 1st example of the moire countermeasure result by the 3D texture mapping apparatus of this Embodiment. It is explanatory drawing which shows the 2nd example of the moire countermeasure result by the 3D texture mapping apparatus of this Embodiment. It is explanatory drawing which shows the 2nd example of the moire countermeasure result by the 3D texture mapping apparatus of this Embodiment. It is a block diagram which shows an example of a structure of the ray casting apparatus as an image processing apparatus of this Embodiment. It is a schematic diagram which shows an example of the coordinate transformation of the voxel structure by the ray casting apparatus of this Embodiment. It is explanatory drawing which shows an example of the coordinate transformation process by the ray casting apparatus of this Embodiment. It is explanatory drawing which shows an example of the search control mask data generation process 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 explanatory drawing which shows the 1st example of the moire countermeasure result by the ray casting apparatus of this Embodiment. It is explanatory drawing which shows the 1st example of the moire countermeasure result by the ray casting apparatus of this Embodiment. It is explanatory drawing which shows the 2nd example of the moire countermeasure result by the ray casting apparatus of this Embodiment. It is explanatory drawing which shows the 2nd example of the moire countermeasure result by the ray casting apparatus of this Embodiment.

(First embodiment)
Hereinafter, embodiments of the present disclosure will be described with reference to the drawings. FIG. 1 is a block diagram illustrating an example of a configuration of a 3D texture mapping apparatus 100 as an image processing apparatus according to the present embodiment. The 3D texture mapping apparatus 100 according to 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 provisional voxel structure generation unit 22, a voxel structure generation unit 23, a conversion processing unit 24, a rendering processing unit 25, a 3D texture image generation unit 26, and a layered rectangle setting unit 27. A color interpolation unit 28, a ROI (Region of Interest) clipping processing unit 29, 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-axis direction. More specifically, the input unit 10 includes a plurality of tomographic images taken at predetermined intervals with respect to an observation target (for example, various parts of a human body or animal including an organ, brain, bone, blood vessel, or the like). Obtain image data in DICOM (Digital Imaging and Communication in Medicine) format of a two-dimensional tomographic image.

  FIG. 2 is a schematic diagram showing an example of a tomographic image in DICOM format. When the xyz axis is set as shown in FIG. 2, a plurality of two-dimensional tomographic images that are xy plane images are arranged along the z-axis 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 tomographic images. In the three-dimensional voxel space, Sx and Sy represent the number of voxels in the x-axis direction and the y-axis direction, respectively, and Sz corresponds to the number of voxels in the z-axis direction. Rxy is the resolution in the x and y directions, and indicates the reciprocal of the pixel interval, that is, the number of pixels per unit distance. The resolutions in the x-axis direction and the y-axis direction are the same, but they may be different. Rz is the resolution in the z-axis direction and represents the reciprocal of the interval between adjacent tomographic images, that is, the number of images per unit distance. The pixel value at the coordinates (x, y, z) of each tomographic 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) indicates that water is 0, a tissue having a higher density than water (for example, a bone) has a positive value, and a tissue having a lower density than water. (For example, adipose tissue) has a negative value. In the coordinates (x, y, z) of the voxel 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, the z coordinate takes a value of 0 or more, Sz It takes a value less than -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 color map data according to the present embodiment. The color map data of the first example is indicated by Cmap (v, c). The variable v indicates a pixel value of Do (x, y, z), and represents a numerical value within a range of −32768 or more and 32767 or less similarly to the range of Do (x, y, z). The variable c includes values of 0, 1, 2, and 3. c = 0 indicates an RGB R value, c = 1 indicates a G value, and c = 2 indicates a 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 slice 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 that they correspond to the pixel value 32767 for convenience, and are actually appropriate numerical values. The color map data does not need to be defined for all pixel values from −32768 to 32767. In the case of a CT image, the color map data is normally adjusted to take a range of −2048 to 2048. The RGB value and opacity in the range from -2048 and in the range from 2048 to 32767 may all be set to zero.

  FIG. 4 is an explanatory diagram showing a second example of color map data according to the present embodiment. The color map data of the second example is indicated by Cmap8 (v, c). As will be described later, when a DICOM image is gradation-compressed from, for example, 16 bits to 8 bits, the pixel value can be set in a range from −32768 to 32767, and from 0 to 255. 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), Cmap 8 (v, c), for example, not only can different colors be given to different parts of the human body, but also opacity can be set for each part. It is possible to make the observation target easier to recognize, for example, by rendering the organ that is located transparent and depicting the organ hidden behind.

  Further, the input unit 10 includes a rotation angle about the x-axis, a rotation angle about the y-axis, a rotation angle about the z-axis, an offset value in the xyz-axis direction, an enlargement or reduction magnification in the xyz-axis direction, and a scaling factor in the z-axis direction. , And a coordinate conversion parameter including the distance from the gazing point to the viewpoint, a value of ROI (Region of Interest) in the xyz direction, a perspective conversion parameter, and the like. Details of the coordinate transformation parameters, ROI values, and perspective transformation parameters will be described later.

  The storage unit 11 can store data acquired by the input unit 10, processing results by the image processing unit 20, and the like. The storage unit 11 may be configured to store color map data, coordinate conversion parameters, ROI values, perspective conversion parameters, and the like in advance.

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

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

  FIG. 5 is an explanatory diagram showing an example of gradation compression processing by the 3D texture mapping apparatus 100 according to the present embodiment. The gradation compression processing unit 21 specifies the minimum value and the maximum value of pixels of a predetermined tomographic image among a plurality of tomographic images. For example, among the tomographic images from 1 to Sz, the minimum value Dmin and the maximum value Dmax of all the pixels in the Sz / 2-th intermediate tomographic image can be specified. As the predetermined slice image, a method for specifying the minimum value and the maximum value of all the pixels in all the slice images is accurate. However, in that case, it is necessary to once hold all large-capacity 16-bit DICOM format tomographic images in the memory, and the gradation compression effect is halved. The minimum value Dmin and the maximum value Dmax are only used as gradation compression parameters, and do not directly affect the accuracy of the 8-bit gradation compressed tomographic image. Therefore, a method of roughly specifying using the minimum value and the maximum value of the pixels of a single slice image is adopted. However, since there are many cases where the subject is not properly reflected in the first slice image, a method of specifying the minimum value Dmin and the maximum value Dmax using only the intermediate slice image is used. Thus, only an 8-bit gradation compressed tomographic image can be directly constructed 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 expressed by the following equation: 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, and Lmax = (Dmax−Dmin) · (1−γ) + Dmin. γ is the contrast adjustment width of the gradation-compressed image. The closer to 0, the greater the contrast but the lower the luminance. Normally, γ = 0.1 is set. Since the luminance contrast of the rendered image can be adjusted by various settings such as color map data, γ may be a fixed value. Also, 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.

  Note that 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 that is 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 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 in, for example, 256 levels, sets the upper limit value Lmax to 255, and sets the lower limit value Lmin to 0.

  Based on the definition information that defines the relationship between a plurality of discrete numerical values, RGB values, and opacity, the image processing unit 20 is a relationship between continuous numerical values (voxel values) in a predetermined range, RGB values, and opacity. Generate a color map that associates. The definition information can be defined by the user, for example.

  The provisional voxel structure generation unit 22 is configured with voxels in which RGB values and opacity are determined with respect to the pixel values of the xy plane image, and the RGB values and opacity correspond to each pixel of the plurality of xy plane images. A provisional voxel structure is generated. That is, the provisional voxel structure generation unit 22 generates a provisional voxel structure composed of voxels having RGB values and opacity corresponding to the values of the respective pixels of the plurality of xy plane images based on the generated color map. Generate. The xy plane image may be a DICOM image Do (x, y, z) or a gradation compressed tomographic image D8 (x, y, z) compressed by the gradation compression processing unit 21. In this specification, the provisional voxel structure is represented by Vo (x, y, z, c). Here, 0 ≦ Vo (x, y, z, c) ≦ 255, c = 0 (R), 1 (G), 2 (B), 3 (α).

  The provisional voxel structure is a set of voxels in which volume data in which vector values composed of RGB values and opacity elements are three-dimensionally packed are discretized and arranged in a grid pattern. For example, one pixel of the DICOM image Do (x, y, z) corresponds to one voxel. The provisional voxel structure Vo (x, y, z, c) can be expressed as Vo (x, y, z, c) = Cmap (Do (x, y, z), c). Here, 0 ≦ x ≦ Sx−1, 0 ≦ y ≦ Sy−1, 0 ≦ z ≦ Sz−1, and 0 ≦ c ≦ 3. The above equation means that the provisional voxel structure Vo (x, y, z, c) is applied by applying the color map Cmap (v, c) to the pixel values of the DICOM image Do (x, y, z). It is to generate.

  The provisional voxel structure Vo (x, y, z, c) can also be generated from the gradation compressed tomographic image D8 (x, y, z) subjected to gradation compression by the gradation compression processing unit 21, and Vo ( x, y, z, c) = Cmap8 (D8 (x, y, z), c). This expression means that the provisional voxel structure Vo (x, y, z, c) is applied by applying the color map Cmap8 (v, c) to the pixel values of the gradation compressed tomographic image D8 (x, y, z). ) Is generated.

  For each voxel of the provisional voxel structure Vo (x, y, z, c), the voxel structure generation unit 23 determines the opacity of the voxel based on the opacity of the neighboring voxel located in the vicinity of the voxel. To generate a voxel structure. In this specification, the voxel structure is represented by V (x, y, z, c). Here, 0 ≦ V (x, y, z, c) ≦ 255, c = 0 (R), 1 (G), 2 (B), 3 (α). In the present specification, a process for generating a voxel structure from a provisional voxel structure is also referred to as a voxel smoothing process. The voxel smoothing process includes an opacity smoothing process, an RGB value and opacity smoothing process, and an RGB value smoothing process.

  Next, the voxel smoothing process will be specifically described. The smoothing processing is not performed on the RGB values and the individual opacity values of the provisional voxel structure created through the color map described in this specification. The method performed on the voxel values of the image is more efficient. However, taking the latter method causes the following problems.

  FIG. 6 shows an example of the color map shown in FIG. As shown in FIG. 6, the color map associates the relationship between voxel values, RGB values, and opacity. For convenience, in FIG. 6, the voxel values are 296 to 304, and the RGB value and α (opacity) are determined corresponding to each voxel value.

  FIG. 7 is a schematic diagram showing an example in which smoothing processing is performed on the voxel values of the DICOM image before passing through the above-described color map as a comparative example. As shown in FIG. 7, for example, it is assumed that the voxel values of three voxels arranged in the x-axis direction are 296, 298, and 304. For convenience, the smoothing process will be described as an average of three voxels adjacent in the x-axis direction. However, in practice, it is desirable to take an average of 27 voxels including the y-axis direction and the z-axis direction. In the smoothing of the voxel value, the voxel value of the center voxel of the three voxels is smoothed by the voxel values (296, 304) of the voxels on both sides. Thereby, the voxel value of the center voxel is changed from 298 to 300. When the color map illustrated in FIG. 6 is applied to the smoothed voxel value (300), RGBα of the central voxel becomes 0, 0, 0, and 80, respectively. Compared to RGBα of both adjacent voxels, a relatively large value change occurs after two identical values of R value, G value, and α value overlap between adjacent voxels, resulting in discontinuous stripes Appears and appears as moire.

  FIG. 8 is a schematic diagram showing an example in which the smoothing process is performed on the opacity of the temporary voxel structure created through the color map of the present embodiment. As in FIG. 7, it is assumed that the voxel values of three voxels arranged in the x-axis direction are 296, 298, and 304. When the color map illustrated in FIG. 6 is applied, RGBα of the three voxels are (40, 0, 0, 80), (20, 0, 0, 80), and (0, 40, 0, 150), respectively. . In the opacity smoothing process, the opacity of the center voxel of the three voxels is smoothed by the opacity (80, 150) of the adjacent voxels. This changes the opacity of the central voxel from 80 to 115. Focusing on the opacity before and after the opacity smoothing process, it can be seen that the opacity between adjacent voxels is smoothly changed by the opacity smoothing process.

  FIG. 9 is a schematic diagram illustrating an example in which smoothing processing is performed on the RGB values and opacity of the provisional voxel structure created through the color map of the present embodiment. As in FIG. 8, it is assumed that the voxel values of three voxels arranged in the x-axis direction are 296, 298, and 304. When the color map illustrated in FIG. 6 is applied, RGBα of the three voxels are (40, 0, 0, 80), (20, 0, 0, 80), and (0, 40, 0, 150), respectively. . In the smoothing process of RGB value and opacity, the RGB value and opacity of the center voxel of the three voxels are smoothed with the RGB value and opacity of the adjacent voxels. As a result, the G value of the central voxel is changed from 0 to 20, and the opacity is changed from 80 to 115. It can be seen that the RGB value (G value in the example of FIG. 9) and opacity between adjacent voxels smoothly change due to the smoothing processing of the RGB value and the opacity. Although not shown, the same applies when the smoothing process is performed only on the RGB values. Further, in the examples of FIGS. 7 to 9, for convenience, the description has been given using the neighboring voxels in the vicinity of 2 in the x-axis direction with respect to the central voxel, but the voxel smoothing processing of the present embodiment is limited to 2 neighborhoods. Is not to be done. Actually, it is desirable to perform the smoothing process using RGBα of 26 neighboring voxels with respect to the xyz axis.

  For each voxel Vo (x, y, z, c) in the provisional voxel structure, Vo (x-1, y-1, z), which are the voxel and 26 neighboring voxels existing in the vicinity of the voxel. −1, c), Vo (x, y−1, z−1, c), Vo (x + 1, y−1, z−1, c), Vo (x−1, y, z−1, c) , Vo (x, y, z-1, c), Vo (x + 1, y, z-1, c), Vo (x-1, y + 1, z-1, c), Vo (x, y + 1, z-) 1, c), Vo (x + 1, y + 1, z-1, c), Vo (x-1, y-1, z, c), Vo (x, y-1, z, c), Vo (x + 1, c) y-1, z, c), Vo (x-1, y, z, c), Vo (x, y, z, c), Vo (x + 1, y, z, c), Vo (x-1, c) y + 1, z, c), Vo (x, y + 1, z, c), o (x + 1, y + 1, z, c), Vo (x-1, y-1, z + 1, c), Vo (x, y-1, z + 1, c), Vo (x + 1, y-1, z + 1, c) ), Vo (x-1, y, z + 1, c), Vo (x, y, z + 1, c), Vo (x + 1, y, z + 1, c), Vo (x-1, y + 1, z + 1, c), The value of 0 ≦ c ≦ 3 is replaced with the average value of Vo (x, y + 1, z + 1, c) and Vo (x + 1, y + 1, z + 1, c). That is, if a smoothed voxel with respect to the voxel Vo (x, y, z, c) is V (x, y, z, c), the smoothed voxel V (x, y, z, c) is obtained. Can be calculated by equation (1). Smoothing is synonymous with voxel smoothing processing.

  In Expression (1), smoothing is not performed on the RGB values (c = 0, 1, 2), but only the opacity (c = 3) value is smoothed (opacity smoothing process). However, four values of 0 ≦ c ≦ 3 can be smoothed (RGB value and opacity smoothing process). Further, only RGB values can be smoothed (RGB value smoothing processing).

  In this way, each voxel Vo (x, y, z, c) of the provisional voxel structure is subjected to isotropic smoothing with an average value in the vicinity of 26 in the xyz direction to obtain a voxel structure V (x, By generating y, z, c), a smoothing process is performed on the voxel values of the DICOM image described in FIG. 7, and then a method of generating a voxel structure through the color map is used. It is possible to suppress moiré that occurs due to overlapping of RGB values and opacity, and steps. In this specification, the voxel structure before the smoothing process is referred to as a temporary voxel structure in order to distinguish and represent the voxel structure before and after performing the voxel smoothing process.

  As described above, in the voxel data of the voxel structure according to the present embodiment, RGB values and opacity are determined with respect to pixel values of a plurality of xy plane images captured at predetermined intervals along the z-axis direction. , Having a structure including voxel data of a voxel structure composed of voxels having RGB values and opacity corresponding to each pixel of a plurality of xy plane images, and the voxel data is in relation to the voxels of the voxel structure Based on the opacity of neighboring voxels located in the vicinity of the voxel, the data including the data in which the opacity of the voxel is corrected is performed, and the volume rendering process is performed on the voxel structure to generate a rendering image. Used for.

  In addition, the voxel data of the voxel structure according to the present embodiment is data obtained by correcting the voxel RGB values based on the RGB values of the neighboring voxels located in the vicinity of the voxel with respect to the voxels of the voxel structure. including.

  Next, an example of a method for speeding up the voxel smoothing process according to the present embodiment will be described.

  In FIG. 7, compared with the case where the smoothing process is performed on the voxel values of the DICOM image and the gradation compressed tomographic image given as the comparative example, the calculation load increases in the smoothing process of the temporary voxel structure according to the present embodiment. . For example, if the smoothing process is performed on all of RGBα, the calculation load increases four times as compared with the case where the smoothing process is performed on the voxel values of the DICOM image. Even when OpenGL / GPU is used, since the smoothing process is executed by software, the waiting time becomes longer until the initial rendering image display is activated. Hereinafter, this point will be described.

  First, the case of smoothing processing for voxel values of a DICOM image or a gradation compressed tomographic image as a comparative example will be described. If the 16 / 8-bit monochrome voxel of the comparative example is Do (x, y, z), a new buffer D (x, y, z) of the smoothed voxel is secured, and x = 1,. , Sx-2, y = 1,..., Sy-2, z = 1,..., Sz-2, the average calculation of 26 neighboring voxels is (Sx-1) × (Sy− 1) It is necessary to repeat x (Sz-1) times.

D (x, y, z) = {Σ _{i = −1, 1} Σ _{j = −1, 1} Σ _{k = −1, 1} Do (x + i, y + j, z + k)} / 27

  Next, the case of the smoothing process for the provisional voxel structure according to the present embodiment will be described. Assuming that the 32-bit provisional voxel structure of RGBα format color after conversion by the color map is Vo (x, y, z, c), the buffer V (32-bit voxel structure of RGBα format after smoothing processing ( x, y, z, c) are newly secured, x = 1,..., Sx-2, y = 1,..., Sy-2, z = 1,. , 3, the average calculation of 26 neighboring voxels needs to be repeated (Sx−1) × (Sy−1) × (Sz−1) × 4 times as in the following equation. Thus, it can be seen that the smoothing process of the provisional voxel structure increases the calculation cost four times as compared with the smoothing process for the voxel values of the DICOM image.

V (x, y, z, c) = {Σ _{i = −1, 1} Σ _{j = −1, 1} Σ _{k = −1, 1} Vo (x + i, y + j, z + k, c)} / 27

  The above formula is a formal rewrite of formula (1) and is substantially the same as formula (1). Hereinafter, a method of performing the smoothing process without using a working voxel buffer and at a low calculation cost will be described.

  The average calculation of 26 neighborhoods in Equation (1) can be divided into an average calculation of 2 neighborhoods in the x-axis direction, an average calculation of 2 neighborhoods in the y-axis direction, and an average calculation of 2 neighborhoods in the z-axis direction.

  In other words, an average calculation of two neighborhoods in the x-axis direction is performed on the voxel Vo (x, y, z, c) before the voxel smoothing processing based on the equations (2) to (10).

  Next, for the voxel Vo ′ (x, y, z, c) updated by Expression (2) to Expression (10), two neighborhoods in the y-axis direction based on Expression (11) to Expression (13) Calculate the average of.

  Finally, for the voxel Vo ″ (x, y, z, c) updated by Expression (11) to Expression (13), the average calculation of two neighbors in the z-axis direction is performed based on Expression (14). Do. Thereby, the same result as the average calculation of 26 neighborhoods of Formula (1) can be obtained.

  When the above-described method is applied to all voxels of the temporary voxel structure, the result is as follows. That is, if a 32-bit voxel in RGBα format is Vo (x, y, z, c), x = 1,..., Sx-2, y = 0,..., Sy-1, z = 0,. −1, c = 0,..., 3, the average calculation of 2 neighborhoods in the x-axis direction is repeated (Sx−2) × Sy × Sz × 4 times, and the updated voxel is Vo ′ (x, y, If z, c), Vo ′ (x, y, z, c) is expressed by the following equation.

Vo ′ (x, y, z, c) = {Σ _{i = −1, 1} Vo (x + i, y, z, c)} / 3

  Using the updated voxel Vo ′ (x, y, z, c), x = 1,..., Sx-2, y = 1,..., Sy-2, z = 0,. = 0,..., 3, the average calculation of 2 neighborhoods in the y-axis direction is repeated (Sx−2) × (Sy−2) × Sz × 4 times, and the updated voxel is Vo ″ (x, y , Z, c), Vo ″ (x, y, z, c) is expressed by the following equation.

Vo ″ (x, y, z, c) = {Σ _{j = −1, 1} Vo ′ (x, y + j, z, c)} / 3

  Using the updated voxel Vo ″ (x, y, z, c), x = 1,..., Sx-2, y = 1,..., Sy-2, z = 1,. For c = 0,..., the average calculation of 2 neighborhoods in the z-axis direction is repeated (Sx−2) × (Sy−2) × (Sz−2) × 4 times, and the updated voxel is represented by V ( If x, y, z, c), V (x, y, z, c) is represented by the following equation.

V (x, y, z, c) = {Σk _{= −1, 1} Vo ″ (x, y, z + k, c)} / 3

  The voxel V (x, y, z, c) updated in the final stage becomes the same as the smoothing processing result averaged by the neighboring voxels in the vicinity of 26. That is, a voxel addition operation involving 27 memory accesses is performed on 27 voxels obtained by summing the voxels to be smoothed and the neighboring voxels in the vicinity of 26 of the voxel (Sx−2) × (Sy−2) × (Sz−). 2) The operation to be performed four times can be realized by a voxel addition operation involving three memory accesses to three voxels obtained by summing up the voxel to be smoothed and neighboring voxels in the vicinity of the voxel, and memory access The number of times can be reduced. Thereby, the smoothing process of the voxel can be speeded up.

  As described above, the voxel structure generation unit 23 calculates, for each voxel of the provisional voxel structure, the voxel of the voxel based on the opacity or the RGB value of 26 neighboring voxels located in the vicinity of the voxel in the xyz direction. The voxel structure can be generated by correcting opacity or RGB values.

  In addition, the voxel structure generation unit 23 determines, for each voxel of the provisional voxel structure, the opacity of two neighboring voxels located in the vicinity of the first axial direction among the three axial directions of the voxel in the xyz direction or Based on the RGB value, the opacity or RGB value of the voxel is corrected and the provisional voxel structure is updated. Next, the voxel structure generation unit 23 applies the second of the three axial directions to the voxel of the temporary voxel structure in which the opacity or RGB value of the voxel is corrected by two neighboring voxels in the first axial direction. Based on the opacity or RGB value of two neighboring voxels located in the vicinity of the axial direction, the opacity or RGB value of the voxel is corrected and the provisional voxel structure is updated. Next, the voxel structure generation unit 23 applies the third of the three axial directions to the voxel of the temporary voxel structure in which the opacity or RGB value of the voxel is corrected by the two neighboring voxels in the second axial direction. A voxel structure can be generated by correcting the opacity or RGB value of a voxel based on the opacity or RGB value of two neighboring voxels located in the axial vicinity.

  As described above, in the voxel data of the present embodiment, RGB values and opacity are determined with respect to pixel values of a plurality of xy plane images captured at predetermined intervals along the z-axis direction, and a plurality of xy The voxel data has a structure including voxel data of a provisional voxel structure composed of voxels having RGB values and opacity corresponding to each pixel of the planar image, and the voxel of the provisional voxel structure is in the xyz direction of the voxel. Processing for correcting the opacity or RGB value of a voxel based on the opacity or RGB value of two neighboring voxels located in the vicinity of the first axial direction among the three axial directions, and two neighboring voxels in the first axial direction 2 neighboring voxels located in the vicinity of the second axial direction of the three axial directions with respect to the voxels of the temporary voxel structure whose opacity or RGB value is corrected by Voxel of the provisional voxel structure in which the opacity or RGB value of the voxel is corrected by the process of correcting the opacity or RGB value of the voxel based on the opacity or RGB value of the image and the two neighboring voxels in the second axis direction On the other hand, a process for correcting the opacity or RGB value of a voxel based on the opacity or RGB value of two neighboring voxels located in the vicinity of the third axis direction of the three axial directions is performed. Used.

  Next, a method that does not use a working voxel buffer will be described. As described above, the smoothing process of the voxel values of the DICOM image (comparative example shown in FIG. 7) requires a voxel buffer for work having the same capacity in order to perform the averaging process of the voxels in the vicinity of the three-dimensional 26. On the other hand, the smoothing process of the provisional voxel structure according to the present embodiment can be divided into three one-dimensional 2-neighbor averaging processes, so that data for all voxels in the voxel structure can be stored. A work voxel buffer is not required, and only a buffer for two voxels necessary for the one-dimensional two-neighbor averaging process is required. This will be specifically described below.

  FIG. 10 is a schematic diagram showing an example of a voxel smoothing algorithm according to the present embodiment. As shown in FIG. 10, a buffer capable of storing opacity or RGB values for two voxels is defined as a first buffer V1 (c) and a second buffer V2 (c). As shown in FIG. 10, the voxel smoothing processing of the present embodiment can be divided into smoothing processing in the x-axis direction, y-axis direction, and z-axis direction (averaging process in the vicinity of 2). The lower part of FIG. 10 illustrates details of the smoothing process in the x-axis direction, and 5 voxels continuous in the x-axis direction are represented by Vo (x−1, y, z, c), Vo (x, y, z, c). ), Vo (x + 1, y, z, c), Vo (x + 2, y, z, c), Vo (x + 3, y, z, c),. Actually, these voxels are followed by voxels in the x-axis direction, but the leftmost Vo (x-1, y, z, c) is regarded as the first voxel in the x-axis direction, A description will be given below.

  When the average value of the two neighborhoods of Vo (x, y, z, c) is calculated, the average value of the two neighborhoods of Vo (x, y, z, c) is stored in the first buffer V1 (c). . At this time, the first buffer V1 (c) must be empty, but is empty in the initial state. Next, when the average value of the two neighborhoods of Vo (x + 1, y, z, c) is calculated, the average value of the two neighborhoods of Vo (x + 1, y, z, c) is also the first buffer V1 (c However, the first buffer V1 (c) is not empty at this stage. Therefore, the average value in the vicinity of Vo (x, y, z, c) stored in the first buffer V1 (c) is moved to the second buffer V2 (c), and the first buffer V1 (c) is moved. Free. Then, an average value in the vicinity of Vo (x + 1, y, z, c) is stored in the first buffer V1 (c).

  Next, when the average value of the two neighborhoods of Vo (x + 2, y, z, c) is calculated, the average value of the two neighborhoods of Vo (x + 2, y, z, c) is also the first buffer V1 (c ), But the first buffer V1 (c) is not empty. For this reason, an average value in the vicinity of Vo (x + 1, y, z, c) stored in the first buffer V1 (c) is moved to the second buffer V2 (c), but the second buffer V2 (C) is also not available. Therefore, the average value of two neighbors of Vo (x, y, z, c) stored in the second buffer V2 (c) is used as the voxel Vo (x, y, z, c) of the original provisional voxel structure. Are overwritten (updated) to free the second buffer V2 (c). Subsequently, an average value in the vicinity of Vo (x + 1, y, z, c) stored in the first buffer V1 (c) is moved to the second buffer V2 (c), and the first buffer V1 (c). Free up. Then, an average value of two neighborhoods of Vo (x + 2, y, z, c) is stored in the first buffer V1 (c). To summarize the above processing, the following three operations may be performed every time an average value of two neighborhoods is calculated. (1) If there is data in the second buffer V2 (c), the data is transferred to a predetermined address of the temporary voxel structure and overwritten. (2) When there is data in the first buffer V1 (c), the data in the first buffer V1 (c) is transferred to the second buffer V2 (c). (3) The average value in the vicinity of 2 is stored in the first buffer V1 (c). However, since (1) and (2) are not executed at the start of processing and (3) is executed first, in this embodiment, the order of operations is changed to (3), (1), and (2) in the following order. explain. (A) The average value in the vicinity of 2 is stored in the first buffer V1 (c). (A) If there is data in the second buffer V2 (c), the data is transferred to a predetermined address of the temporary voxel structure and overwritten. (C) The data in the first buffer V1 (c) is transferred to the second buffer V2 (c). When these three operations are repeated until the end is completed, data remains in the second buffer V2 (c), so the data remaining in the second buffer V2 (c) is transferred to a predetermined address of the temporary voxel structure and overwritten. Then, a series of processing ends.

  That is, for each of y = 0,..., Sy-1, z = 0,..., Sz-1, (y, z), for x = 1,. ) V1 (c) = {Vo (x-1, y, z, c) + Vo (x, y, z, c) + Vo (x + 1, y, z, c)} / 3 (c = 0,..., 3 ), (B) When x> 1, Vo (x-1, y, z, c) = V2 (c), (c) V2 (c) = V1 (c) (c = 0,..., 3) The above three operations (a) to (c) are repeated. Finally, Vo (Sx−2, y, z, c) = V2 (c) is set, and voxel Vo (x, y, z, c) ( Update one column (x = 1,..., Sx−2) in the x-axis direction in y, z).

  Next, in the updated voxel Vo (x, y, z, c), x = 0,..., Sx-1, z = 0,. , Y = 1,..., Sy-2, (a) V1 (c) = {Vo (x, y-1, z, c) + Vo (x, y, z, c) + Vo (x, y + 1) , Z, c)} / 3 (c = 0,..., 3), (b) When y> 1, Vo (x, y-1, z, c) = V2 (c), (c) V2 ( c) = V1 (c) (c = 0,..., 3) (a) to (c) are repeated, and finally, Vo (x, Sy-2, z, c) = V2 (c ) And one column (y = 1,..., Sy−2) in the y-axis direction in (x, z) of voxel Vo (x, y, z, c) is updated.

  Next, in the updated voxel Vo (x, y, z, c), x = 0,..., Sx-1, y = 0,. Z = 1,..., Sz-2, (a) V1 (c) = {Vo (x, y, z-1, c) + Vo (x, y, z, c) + Vo (x, y) , Z + 1, c)} / 3 (c = 0,..., 3), (a) When z> 1, Vo (x, y, z-1, c) = V2 (c), (c) V2 ( c) = V1 (c) (c = 0,..., 3) (a) to (c) are repeated, and finally, Vo (x, y, Sz−2, c) = V2 (c ), And updates one column (z = 1,..., Sz−2) in the z-axis direction in (x, y) of voxel Vo (x, y, z, c).

  As described above, the voxel structure generation unit 23 stores, in the first buffer, the correction value obtained by correcting the opacity or the RGB value of each voxel by using two neighboring voxels of the voxel of the provisional voxel structure. The immediately preceding voxel of the provisional voxel structure is corrected by the immediately preceding correction value obtained by correcting the opacity or RGB value of the immediately preceding voxel by two neighboring voxels of the immediately preceding voxel processed immediately before each of the voxels of the provisional voxel structure stored in The RGB value or opacity of each voxel of the provisional voxel structure is updated by repeating the three operations of updating the opacity or RGB value of the image and transferring the correction value stored in the first buffer to the second buffer. Thus, a voxel structure can be generated.

  As described above, the voxel data according to the present embodiment is stored in the first buffer with a correction value obtained by correcting the opacity or RGB value of each voxel by using two neighboring voxels of each voxel of the provisional voxel structure. Temporary voxel with the previous correction value obtained by correcting the opacity or RGB value of the previous voxel by processing and two neighboring voxels of the previous voxel processed immediately before each voxel of the temporary voxel structure stored in the second buffer This is used to execute a second process for updating the opacity or RGB value of the immediately preceding voxel of the structure and a third process for transferring the correction value stored in the first buffer to the second buffer.

  As described above, the required buffer can be reduced to two voxels, and the memory capacity can be greatly reduced. In addition, the average value of two neighbors calculated two voxels before the one-dimensional two-near voxel averaging process can be advanced while overwriting the original voxel structure, and all voxels of the voxel structure can be processed. Smoothing processing can be performed at high speed without using a working memory.

  In the following, a process for the voxel structure after the voxel smoothing process is performed will be described. That is, in the following, it is assumed that the voxel structure has been subjected to a voxel smoothing process.

  The voxel structure is represented by V (x, y, z, c). Here, 0 ≦ V (x, y, z, c) ≦ 255, c = 0 (R), 1 (G), 2 (B), 3 (α), 0 ≦ x ≦ Sx−1, 0 ≦ It is y <= Sy-1 and 0 <= z <= Sz-1.

  Further, the voxel structure generation unit 23 can generate a voxel structure including voxels whose RGB values and opacity are predetermined values. As the predetermined values, the RGB values are R = G = B = 0 (that is, colorless or black) and the opacity α = 0 (that is, transparent), but these do not contribute to the calculation of the volume rendering image. It means that the value is a value and does not contribute to the calculation of the volume rendering image, it is not necessary to limit to R = G = B = 0 and opacity α = 0.

  The color interpolation unit 28 performs color interpolation on the target voxel having the RGB value and the opacity of the voxel structure having predetermined values based on the RGB values of neighboring voxels located in the vicinity of the target voxel, and Correct the RGB values. Note that color interpolation is also referred to as color correction. Details of the color interpolation will be described later.

  The 3D texture image generation unit 26 generates a 3D texture image composed of voxels having a voxel structure in which RGB values are color-interpolated. That is, the 3D texture image generation unit 26 generates a 3D texture image based on the voxels generated by the voxel structure generation unit 23 (including voxels whose RGB values and opacity are predetermined values). The data of the 3D texture image is also called a 3D texture map, and the substance of the data is the same as the voxel structure. The 3D texture image generation unit 26 registers the volume of the voxel structure as a 3D texture map of a graphic API (for example, OpenGL). More specifically, the volume of the voxel structure is transferred from the CPU main memory to the GPU video memory via the graphic API, converted into the graphic API texture format, and held in the video memory. In the present application, the contents of the process will be described below for OpenGL as a graphic API for convenience of explanation, but the present invention is not limited to this.

  A 3D texture map is represented by Tex (xt, yt, zt, c). Here, 0 ≦ Tex (xt, yt, zt, c) ≦ 255, c = 0 (R), 1 (G), 2 (B), 3 (α), 0 ≦ xt ≦ Sx−1, 0 ≦ It is yt <= Sy-1 and 0 <= zt <= Sz-1.

  Note that the coordinate system of the voxel structure is defined by the OpenGL world coordinate system, expressed by XYZ coordinates (-1 ≦ x ≦ 1, −1 ≦ y ≦ 1, −1 ≦ z ≦ 1), and the coordinates of the 3D texture map. The system is defined by a texture coordinate system, and OpenGL has a custom of representing a texture coordinate system by STR coordinates (0 ≦ s ≦ 1, 0 ≦ t ≦ 1, 0 ≦ r ≦ 1). The coordinate s corresponds to the coordinate x, the coordinate t corresponds to the coordinate y, and the coordinate r corresponds to the coordinate z. In the present embodiment, for convenience, XYZ in the world coordinate system is represented as XYZ coordinates, and STR coordinates in the texture coordinate system are represented as xyz coordinates.

  The stacked quadrangle setting unit 27 sets a stacked quadrangle in which a plurality of quadrilaterals defined on the XY coordinate plane are arranged with a predetermined interval in the Z-axis direction in the world coordinate system. A stacked quadrangle is a quadrangle defined by four vertices on a three-dimensional space.

  The conversion processing unit 24 performs a predetermined conversion in the texture coordinate system on the 3D texture image composed of a plurality of slice images defined in the texture coordinate system, and generates a converted 3D texture image. The predetermined conversion includes a perspective conversion process, a coordinate conversion process, and the like.

  FIG. 11 is a schematic diagram showing the concept of 3D texture mapping by the 3D texture mapping apparatus 100 of the present embodiment. FIG. 11A shows a state in which a 3D texture map before coordinate conversion is mapped to the layered quadrangle, and FIG. 11B shows a state in which a 3D texture map after coordinate conversion is mapped to the layered quadrangle. In the world coordinates (X, Y, Z), a plurality of quadrangles on the XY coordinate plane are arranged along the Z-axis direction. In FIG. 11, the viewpoint exists above the projection plane (the plane on which the rendering image is generated) (that is, in the positive direction of the Z axis). Since the stacked rectangle defined in the world coordinate system is fixed, each rectangle is always perpendicular to the line of sight. In FIG. 11, the oval shape illustrated inside the stacked rectangle indicates a mapped 3D texture map (voxel data). The 3D texture map is appropriately converted in the texture coordinate system, but since the correspondence (mapping) between the world coordinate system and the texture coordinate system is fixed, the 3D texture map is mapped to the stacked rectangle after the conversion of the 3D texture map. The content of the existing 3D texture map is automatically changed.

  As shown in FIG. 11A, in 3D texture mapping, voxel data is defined as a 3D texture map in the texture coordinate system, and the 3D texture map is pasted on each quadrangle of the stacked rectangle defined in the world coordinate system. That is, by mapping each cross section (slice) of the 3D texture map to a quadrangle, the 3D texture map is reconstructed as a collection of cross sections perpendicular to the line of sight.

  As shown in FIG. 11B, the viewpoint can be changed by pasting the rectangle defined in the world coordinates (X, Y, Z) and the rectangle without changing the correspondence between the world coordinate system and the texture coordinate system. The attached 3D texture map is rotated in the texture coordinate system. As a result, the four vertexes of each quadrangle and the coordinate position of the 3D texture map (represented by xyz coordinates in FIG. 7) are changed by rotation, and an image with a changed viewpoint can be obtained.

  The rendering processing unit 25 performs a volume rendering process on the voxel structure to generate a rendering image.

  Specifically, the rendering processing unit 25 converts the RGB values of voxels of the converted 3D texture image corresponding to the XY coordinates of the rectangle on the line of sight parallel to the Z-axis direction from a predetermined viewpoint in the order of the rectangles farther from the viewpoint. A rendering image is generated by using the RGB value obtained by alpha blending based on the alpha value of the pixel as the pixel value of a pixel whose line of sight intersects the projection plane. In this case, the rendering processing unit 25 calculates the intersection coordinates (XYZ coordinates) of a quadrangle that intersects the XY coordinates of the rendering image in the line-of-sight direction, and calculates the xyz coordinates of the converted 3D texture image corresponding to the calculated intersection coordinates. Then, alpha blending is performed based on the RGB value and opacity of the voxel corresponding to the calculated xyz coordinates. That is, volume rendering is performed by alpha-blending the quadrangle with the 3D texture map pasted in order from the viewpoint. By instructing these procedures using a graphic API such as OpenGL, high-speed processing using a GPU function is possible. By using the GPU function, it is possible to generate a rendering image at high speed while suppressing waiting time. In the present embodiment, even a GPU mounted on a video card mounted on a general-purpose computer can perform rendering processing within 0.1 seconds, and can be interactively (or assured of interactivity). ) The rendered image can be displayed in real time.

  In the present embodiment, the 3D texture mapping method is used. However, the general 2D texture mapping method has the following problems. That is, in the 2D texture mapping method, the outline shape of each slice constituting the voxel data is represented by a rectangle in the world coordinate system, as in the present application, but each slice constituting the voxel data in the two-dimensional texture coordinate system. Are defined as a plurality of 2D texture maps (in this embodiment, all slice images are defined as a single 3D texture map in a three-dimensional texture coordinate system). Then, any 2D texture map corresponding to each square is pasted on each square (in this embodiment, four coordinates in the 3D texture map corresponding to the four vertices of each square are pasted in association with each other. ). Since the 2D texture map is defined in a two-dimensional texture coordinate system, rotation can be performed only two-dimensionally. For this reason, in the world coordinate system, a method of rotating the stacked rectangles three-dimensionally must be taken. Then, the depth relationship is broken depending on the rotation angle, or the 2D texture map is reversed, and when the stacked rectangle is rotated 90 degrees about the Y axis, the rectangle (thickness is 0) is not drawn and the rendered image becomes black.

  In addition, as the angle between the slice and the line of sight decreases, the projected area decreases, the number of sampled data decreases, the significance of volume rendering is lost, and a practical rendered image may not be obtained. . That is, an appropriate rendering image cannot be obtained unless the correspondence between the layered rectangle in the world coordinate system and the 2D texture map in the texture coordinate system is corrected by the rotation angle.

  According to the present embodiment, the problem in the 2D texture mapping method or the problem that occurs when the world coordinates are rotated can be solved.

  Next, the operation of the 3D texture mapping apparatus 100 according to the present embodiment will be described.

  12 and 13 are flowcharts illustrating an example of a 3D texture mapping process performed by the 3D texture mapping apparatus 100 according to the present embodiment. Hereinafter, for the sake of convenience, the subject of processing will be described as the image processing unit 20. The image processing unit 20 acquires image data of the DICOM image (S11), and acquires color map data (S12). The color map data is illustrated in FIGS. 3 and 4.

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

  The image processing unit 20 generates the provisional voxel structure Vo (x, y, z, c) by applying the color map (S17). The image processing unit 20 performs voxel smoothing processing on the provisional voxel structure Vo (x, y, z, c) to generate the voxel structure V (x, y, z, c) (S18). . The image processing unit 20 performs transparency processing for voxels outside the ROI clipping region (S19).

  FIG. 14 is an explanatory diagram illustrating an example of the voxel transparency processing outside the ROI clipping region by the 3D texture mapping apparatus 100 according to the present embodiment. FIG. 14 illustrates an example of the xy plane image located at the zth position. 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. Although not shown in FIG. 14, as ROI in the z-axis direction, Set Zs-Ze. Here, 0 ≦ Xs, Xe ≦ Sx−1, 0 ≦ Ys, Ye ≦ Sy−1, 0 ≦ Zs, Ze ≦ Sz−1. That is, the ROI clipping processing unit 29 sets xyz coordinates that define the target area for volume rendering processing. The voxel structure generation unit 23 performs a transparency process for voxels outside the ROI clipping region.

  In the voxel transparency processing, when x <Xs or x> Xe or y <Ys or y> Ye or z <Zs or z> Ze, V <x, y, z, c ) = 0, that is, R = G = B = α = 0 (that is, the pixel color is colorless, that is, black and transparent). In the example of FIG. 14, the color of the pixel at coordinates (x, y) satisfying x <Xs, x> Xe, y <Ys, and y> Ye is colorless (black), and the opacity is α = 0 (transparent ).

  It should be noted that the color interpolating unit 28 determines a neighboring voxel (opaque voxel) located in the vicinity of the target voxel with respect to a colorless / transparent voxel including a voxel that has been subjected to transparency processing such as outside the ROI clipping region in step S24 described later. Based on the RGB values, color interpolation is performed on the RGB values of the target voxel. However, even if the voxel is transparent, it is not transparentized, and color interpolation is not performed when one or more of RGB values is set (such as lung field air).

  As described above, in the voxel data of the voxel structure according to the present embodiment, the voxel structure is outside the target area of the volume rendering process, and the RGB values and the opacity are predetermined values. Based on the RGB values of the neighboring voxels located in the vicinity of the target voxel, the data including the color interpolation of the RGB values of the target voxel is included.

  Not only the ROI clipping plane, but an opaque area (area where the opacity has a value of 1 or more) that is the target area of volume rendering processing and a transparent area (opacity is 0) that is outside the target area of volume rendering processing. The boundary surface with the region having the value of) becomes the source of moire. (Details will be described later, but this is because a step-like pattern is generated on the boundary surface by coordinate transformation.) In the first place, the tomographic image has a background portion other than the subject such as a human body and a photographing jig (such as a bed). In general, since the background portion is colorless and transparent with air, a boundary surface with an opaque subject originally exists and can be a source of moire. It is necessary to perform color interpolation together.

  Not only the ROI clipping section but also the first slice and / or the last slice of the tomographic image may be the boundary surface between the opaque region and the transparent region as it is, and the boundary surface may be the source of moire. There is no image data before and after the first slice of the tomographic image that is CT-scanned for the chest and abdomen that are part of the human body, and volume rendering processing is performed as a colorless and transparent space. Therefore, it has the same characteristics as the ROI clipping section. That is, the first slice and the last slice of the tomographic image inevitably become a boundary surface between the opaque region and the transparent region, and become a moire generation source. Therefore, the generation of moire can be suppressed at the same time by adding voxels to the area before the first slice and behind the last slice to make it colorless and transparent and performing color interpolation as the subsequent processing. Details of the color interpolation will be described later.

  The image processing unit 20 newly adds dummy transparent images with RGB values and opacity set to predetermined values (here, colorless and transparent) at both ends in the z direction (S20), and the dummy transparent image is added to the voxel structure. A voxel corresponding to is added (S21). A voxel structure to which a voxel corresponding to the dummy transparent image is added is represented as V (x, y, z, c). Here, 0 ≦ x ≦ Sx−1, 0 ≦ y ≦ Sy−1, and 0 ≦ z ≦ Sz + 1.

  FIG. 15 is an explanatory diagram showing an example of dummy transparent image addition processing by the 3D texture mapping apparatus 100 of the present embodiment. Dummy transparent images in which the RGB values and opacity of all the pixels are set to predetermined values at both ends in the z-axis direction of the original tomographic image arranged at the slice interval in the z-axis direction (here, set to colorless and transparent) Add That is, the voxel structure generation unit 23 corresponds to each pixel of the dummy xy plane image arranged at both ends in the z-axis direction with a plurality of xy plane images in between, and the RGB value and the opacity are set to predetermined values. Generate a voxel structure with dummy voxels added. Specifically, the voxel structure generator 23 adds voxels corresponding to the dummy transparent image to both ends in the z-axis direction of the tomographic image. As the predetermined values, the RGB value can be set to R = G = B = 0 and the opacity can be set to α = 0 (that is, colorless and transparent) so that the voxel does not contribute to the generation of the volume rendering image. As long as the value does not contribute to the generation of the volume rendering image, the RGB value need not be limited to R = G = B = 0, and the opacity need not be limited to α = 0.

  Note that the color interpolation unit 28 determines the RGB values of the neighboring voxels located in the vicinity of the target voxel with respect to the target voxel having a predetermined value for the RGB value and opacity of the voxel structure generated in step S24 described later. On the basis of the color interpolation of the RGB values of the target voxel. However, even if it is a transparent voxel, if a value other than a predetermined value (here, a value of 1 or more) is set in any of the RGB values (for example, air in the lung field), color interpolation is performed. Not performed.

  As described above, in the voxel data of the voxel structure according to the present embodiment, the voxel structure corresponds to each pixel of the dummy xy plane image arranged at both ends in the z-axis direction with a plurality of xy plane images in between. And a voxel data including a dummy voxel having a predetermined value for the RGB value and the opacity, and the voxel data for the target voxel having the RGB value and the opacity for the predetermined value. , Including data obtained by color-interpolating the RGB values of the target voxel based on the RGB values of the neighboring voxels located in the vicinity of the target voxel.

  Not only the ROI clipping section but also the peripheral part of the tomographic image may be the boundary surface between the opaque region and the transparent region as it is and may be a source of moire. For example, depending on the CT scan magnification setting, an image may be captured with the peripheral portion of the tomographic plane of the human body cut off. In this case, x = 0, y in FIG. = 0, x = Sx-1, and y = Sy-1 are the outer regions of the tomographic image of x <0, y <0, x> Sx-1, y> Sy-1 Since no image data exists and volume rendering processing is performed as a colorless and transparent space, it has the same characteristics as the ROI clipping cut surface. That is, a plane satisfying any of x = 0, y = 0, x = Sx-1, and y = Sy-1 inevitably becomes a boundary surface between the opaque region and the transparent region, and becomes a moire generation source. Therefore, for all voxels at the outer edge portion of 1 voxel width satisfying any of x = 0, y = 0, x = Sx-1, and y = Sy-1 regardless of the presence or absence of the ROI clipping instruction, As in the shaded area in FIG. 14, the color is made colorless and transparent, and color interpolation which is subsequent processing is also performed. In other words, the image processing unit 20 performs the voxel transparency processing corresponding to the pixels at the outer edge (also referred to as border) of the gradation compressed tomographic image (S22). Note that the process of step S22 is performed regardless of whether or not the process of step S23 described later is performed. The outer edge means an area in the image, not the outer side of the image, and a portion corresponding to the edge of the image.

  FIG. 16 is an explanatory diagram showing an example of the voxel transparency processing of the outer edge portion of the tomographic image by the 3D texture mapping apparatus 100 of the present embodiment. The color of the voxel corresponding to the pixel at the outer edge of each tomographic image aligned at the slice interval in the z-axis direction is set to R = G = B = 0 (colorless), and the opacity α is set to 0 (transparent). That is, the voxel structure generation unit 23 makes voxels corresponding to pixels at the outer edge of the tomographic image colorless and transparent. Here, the outer edge portion is, for example, a region in the tomographic image and can be one pixel from the edge. The voxel structure V (x, y, z, c) has V (x, y, z, c) = 0 (x = 0 or x = Sx-1 or y = 0 or y = Sy-1 or z = 0 or z = Sz + 1, 0 ≦ c ≦ 3). That is, the voxel structure generation unit 23 generates a voxel structure including voxels corresponding to each pixel of the outer edge portion of the plurality of xy planar images and having RGB values and opacity values set to predetermined values. Specifically, the voxel structure generation unit 23 makes voxels corresponding to pixels at the outer edge of the tomographic image colorless and transparent. Note that the outer edge portion need not be limited to one pixel. As the predetermined value, the RGB value can be R = G = B = 0 and the opacity α = 0 (that is, the color of the pixel is colorless, that is, black and transparent), but the voxel that does not contribute to the volume rendering can be used. If it is a value, the RGB value may not be limited to R = G = B = 0, and the opacity may not be limited to α = 0.

  Note that, in step S24 described later, the color interpolation unit 28 is based on the RGB values of the neighboring voxels positioned in the vicinity of the target voxel with respect to the target voxel of the voxel corresponding to the colorless and transparent pixel at the outer edge portion of the tomographic image. Thus, the RGB value of the target voxel is color-interpolated. However, even if it is a transparent voxel, if a value other than a predetermined value (here, a value of 1 or more) is set in any of the RGB values (for example, air in the lung field), color interpolation is performed. Not performed.

  As described above, in the voxel data of the voxel structure according to the present embodiment, the voxel structure corresponds to each pixel of the outer edge portion of the plurality of xy plane images, and the voxel in which the RGB value and the opacity are predetermined values. The voxel data includes the RGB values of the target voxel based on the RGB values of the neighboring voxels located in the vicinity of the target voxel with respect to the target voxel having the predetermined RGB value and opacity of the voxel structure including the voxel. Contains data with color interpolated.

  The image processing unit 20 calculates the shadow value of the voxel (S23). However, the process of step S23 is not essential, and if the process of step S23 is not performed, the process may proceed to step S24.

  The voxel shadow value can be calculated as follows. That is, the voxel structure generation unit 23 calculates a gradient vector related to the opacity of the voxel of the voxel structure, and calculates a shadow value of the voxel based on the calculated gradient vector and the predetermined light source vector.

  More specifically, the voxel structure generation unit 23 calculates a difference between the opacity of the voxel of the voxel structure and the opacity of each of the plurality of neighboring voxels located in the vicinity of the voxel. The voxel structure generator 23 calculates a voxel gradient vector based on the distance between the voxel and the neighboring voxel and the difference in opacity between the voxel and the neighboring voxel. In this case, considering that the resolution in the XY direction is different from that in the Z direction when using the distance to the neighboring voxel, the voxel structure generation unit 23 converts the z-axis direction component of the gradient vector to the scaling factor Rxy in the z-axis direction. It is necessary to correct with / Rz.

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

  The gradient vector (Gx, Gy, Gz) of the voxel (x, y, z) in the range of 0 ≦ x ≦ Sx−2, 0 ≦ y ≦ Sy−2 and 0 ≦ z ≦ Sz is expressed by the equations (15), ( 16) and (17). In the present embodiment, the gradient vector can be calculated with high accuracy by the difference between the voxel values in the vicinity of 26 and the voxel values of the center voxel. When calculating the gradient vector, since the voxel has not been subjected to z-direction scaling before coordinate conversion, it should be noted that the resolution differs between the xy direction and the z-direction, and the change in the z-direction component Gz of the gradient vector is performed. Double correction (multiplication of Rxy / Rz) is performed.

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 S (x, y, z) is calculated by only the diffuse reflection component, and is given by the equation (18). Thereby, when the line-of-sight vector is changed and the voxel is rotated, a stereoscopic effect can be expressed. Even if the line of sight is changed, it is before the coordinate rotation process and the gradient vector does not change, so no specular reflection component is added. The voxel structure generator 23 corrects the RGB value of the voxel based on the calculated luminance value. When the voxel after the shadow calculation is represented by V ′ (x, y, z, c), the correction of the RGB value of the voxel can be performed by applying Expression (19) in the range of 0 ≦ c ≦ 2. .

  In the 3D texture mapping method using OpenGL, processing related to voxel gradient vector and shadow calculation is not supported in hardware rendering, and therefore an α value is always given as a shadow value. However, in the present embodiment, pseudo shadow calculation can be performed. In other words, in order to accurately perform the shadow calculation, it is necessary to perform the voxel coordinate conversion. However, since the coordinate conversion is performed together in hardware rendering, the voxel data to be input to hardware rendering before coordinate conversion is previously stored. By incorporating the shadow calculation result, it is possible to obtain a rendering image with a pseudo shadow.

  G <1 may occur in voxels at the outer edge or background of the tomographic image, and in this case, the calculation is uniformly colorless and transparent (R = G = B = α = 0) without performing shadow calculation. That is, when the voxel after the shadow calculation is represented by V ′ (x, y, z, c), V ′ (x, y, z, c) = 0 (c = 0, 1, 2, 3). This is a process equivalent to step S22 described above.

  The image processing unit 20 performs transparent voxel color interpolation processing (S24). Details of the color interpolation process will be described later.

  The image processing unit 20 registers the color-interpolated voxel structure as a 3D texture map (S25), and sets a projection screen (S26).

  FIG. 17 is an explanatory diagram showing an example of a projection screen setting by the 3D texture mapping apparatus 100 according to the present embodiment. The conversion processing unit 24 sets the screen size (number of vertical and horizontal pixels, vertical and horizontal aspect ratio) of the rendered image. The conversion processing unit 24 sets either parallel projection (normal appearance rendering) or perspective projection (endoscopic mode). Then, when the perspective projection is set, the conversion processing unit 24 sets the perspective projection parameter. The perspective projection parameters are, for example, a camera viewing angle (focal length), a viewpoint position (a viewpoint and a gazing point, the former being the position of the eye and the latter being the position on the object being viewed, both on the Z axis. In general, the gazing point is fixed at the origin of the world coordinate system), a clipping position (distance on the Z-axis from the viewpoint, two near and far points), and the like. Note that the clipping position may be only near. As shown in FIG. 17, in the case of parallel projection, it is assumed that the line of sight from the viewpoint is all parallel to the Z axis, and the viewpoint is virtually located at an infinite distance in the left direction. All the stacked rectangles defined in the Z-axis direction are to be rendered. On the other hand, in the case of perspective projection, the line of sight spreads according to the viewing angle from the viewpoint, and the rendering target is also limited by near clipping.

  The image processing unit 20 performs rendering processing (S27), generates a rendering image on the frame memory (S28), and ends the processing. Details of the rendering process will be described later.

  FIG. 18 is a flowchart illustrating an example of a procedure of transparent voxel color interpolation processing by the 3D texture mapping apparatus 100 according to the present embodiment. The image processing unit 20 extracts colorless and transparent voxels from all voxels of 0 ≦ x ≦ Sx−1, 0 ≦ y ≦ Sy−1, and 0 ≦ z ≦ Sz + 1 in the voxel structure (S101). That is, voxels satisfying the equation (20) are extracted.

  The image processing unit 20 identifies a voxel that is an opaque voxel near the target voxel among the extracted voxels (S102). Assuming that the xyz coordinates of the target voxel are (x, y, z), from 26 neighboring voxels in the vicinity of the target voxel in the range of −1 ≦ i ≦ 1, −1 ≦ j ≦ 1, and −1 ≦ k ≦ 1. An opaque voxel satisfying V (x + i, y + j, z + k, 3)> 0 is specified, and the number thereof is set to C according to the definition of Expression (21). That is, the voxel structure generation unit 23 specifies neighboring voxels having an opacity greater than a predetermined value among neighboring voxels related to the target voxel. The predetermined value may be opacity α = 0 that does not contribute to volume rendering (that is, transparent), but the opacity may not be limited to α = 0 if it does not contribute to volume rendering.

  The image processing unit 20 determines whether C> 0, that is, whether there is an opaque voxel (S103). If C> 0, that is, if there is an opaque voxel (YES in S103), the average of the RGB values of the opaque voxels is determined. The RGB value of the target voxel is updated with the value (S104), and the process of step S106 described later is performed. Assuming that the average value of RGB values is Ra, Ga, Ba, the average value is the R value of each opaque voxel satisfying V (x + i, y + j, z + k, 3)> 0 according to equations (22), (23), and (24). V (x + i, y + j, z + k, 0), G value V (x + i, y + j, z + k, 1), and B value V (x + i, y + j, z + k, 2) can be calculated. In addition, the RGB value of the colorless and transparent target voxel can be updated by the equations (25), (26), and (27).

  When C = 0, that is, when there is no opaque voxel (NO in S103), the image processing unit 20 sets the RGB value of the target voxel to 0 (that is, Ra = Ga = Ba = 0) (S105) and extracts it. It is determined whether or not the color interpolation of all target voxels of colorless and transparent voxels has been completed (S106). When the color interpolation of all target voxels has not been completed (NO in S106), the image processing unit 20 continues the processing from step S102 onward, and when the color interpolation of all target voxels has been completed (YES in S106), The process ends.

  The color interpolation process is performed by the color interpolation unit 28 of the image processing unit 20. More specifically, the color interpolation unit 28 performs color interpolation on the RGB values of the target voxel based on the statistical values of the RGB values of the specified neighboring voxels. As the statistical value, an average value can be used, but the statistical value is not limited to the average value. For example, a median value may be used.

  As described above, the voxel data of the voxel structure according to the present embodiment includes RGB of neighboring voxels having an opacity greater than a predetermined value among neighboring voxels related to the target voxel (that is, voxels having an opacity α of 1 or more). Based on the statistical value, RGB data of the target voxel is included in the color interpolated data.

  FIG. 19 is a flowchart illustrating an example of a rendering process performed by the 3D texture mapping apparatus 100 according to the present embodiment. The image processing unit 20 performs scaling and z-direction scaling processing on the 3D texture map (S111), and performs rotation processing on the 3D texture map (S112).

  The image processing unit 20 performs offset processing on the 3D texture map (S113), and sets a stacked quadrangle composed of a plurality of quadrangles (S114). The image processing unit 20 associates each quadrangle constituting the stacked quadrangle with the 3D texture map (S115).

  FIG. 20 is an explanatory diagram illustrating an example of a method of associating with 3D texture mapping. As shown in FIG. 20A, the DICOM tomographic image is converted into the RGBα format and registered as a 3D texture map via the graphic API. Then, as shown in FIG. 20B, a stacked quadrangle is defined and 3D texture mapping is performed. In this case, the three-dimensional world coordinates of the four vertices of each rectangle in the world coordinate system are associated with the three-dimensional texture coordinates of the 3D texture map in the texture coordinate system. As shown in FIG. 20, a quadrangle composed of four vertices (-1, -1, z), (-1, 1, z), (1, -1, z), and (1, 1, z). On the other hand, the world coordinates (-1, -1, z) are associated with the texture coordinates (0, 0, r).

  When 3D texture mapping is set, the layered rectangle setting unit 27 sets a layered rectangle in which the same number of rectangles as the total number of the plurality of xy plane images and the number of dummy xy plane images are arranged in the Z-axis direction. .

  The number of quadrangles in the layered quadrangle can be arbitrarily set, but when the number of quadrangles and the number of tomographic images do not match, interpolation processing works in the Z-axis direction, and rounding errors of coordinate values due to coordinate transformation accumulate, Striped / lattice moire occurs.

  Therefore, the number of quadrangles in the layered quadrangle and the number of 2D tomographic images registered in the 3D texture map are set to the number of original tomographic images in DICOM format + 2 (dummy transparent image added in the z-axis direction). When constructing a 3D voxel structure, it is easier to perform post-processing if the XYZ three-dimensional resolution is unified, so that the number of slices in the 2D tomographic image is interpolated to increase the number of tomographic images in the original DICOM format. Although processing is generally performed, the present embodiment is more preferable from the viewpoint of suppressing moire.

  The image processing unit 20 performs scan conversion while pasting the 3D texture map on the layered rectangle, performs alpha blending processing in the order of the rectangle far from the viewpoint (S116), and ends the processing.

  More specifically, all the RGB values of the rendered image are initialized in advance with a background color (for example, R = G = B = 0). The rendering processing unit 25 pastes the 3D texture map on the quadrilateral farthest from the viewpoint, performs parallel projection or perspective projection on the rendered image defined in FIG. 17, and (-1) of the XYZ coordinates in the world coordinate system of the projected quadrangle. , -1, -1) to the X-axis direction and the Y-axis direction, the interval u corresponding to the pixel interval in the projected rendered image (the interval of world coordinates corresponding to one pixel of the rendered image is u, for example, u. = 0.002), scan conversion is performed. Each world coordinate (−1 + iu, −1 + ju, −1) (i, j is a coordinate value of the rendered image, for example, an integer of 0 ≦ i, j ≦ 511) subjected to scan conversion, and corresponding texture coordinates (2iu) , 2ju, 0), RGB values and opacity are obtained by referring to the 3D texture map based on the calculated texture coordinates, and the RGB already recorded at the coordinates (i, j) of the corresponding rendered image The value and the alpha blending process are performed, and the RGB value of the coordinate (i, j) of the corresponding rendering image is updated. In this way, when the alpha blending process is completed for the single rectangle farthest from the viewpoint, the Z coordinate is v (the interval between world coordinates where the rectangle is arranged in the Z-axis direction is v. For example, v = 2. / Sz) is increased, and the same processing is repeated for the quadrangle (−1 + iu, −1 + ju, −1 + v) closest to the viewpoint direction to update the rendered image. Rendered image when scan conversion is performed for all squares and alpha blending processing is completed for world coordinates (−1 + iu, −1 + ju, −1 + kv) (k is a square number and 0 ≦ k ≦ Sz−1) Is completed.

  In the alpha blending process, equations (28), (29), and (30) can be used.

  Here, R ′, G ′, and B ′ are RGB values of the rendered image updated on the projection plane. R, G, and B are RGB values in the voxel (2iu, 2ju, 2kv) of the 3D texture map corresponding to the square world coordinates (-1 + iu, -1 + ju, -1 + kv), and correspond to the values of the colors to be superimposed. α is also an α value in the voxel (2iu, 2ju, 2kv) of the 3D texture map corresponding to the rectangular world coordinates (-1 + iu, -1 + ju, -1 + kv), and controls the overlapping ratio. Rb, Gb, and Bb are RGB values that are already recorded at the coordinates (i, j) of the rendered image corresponding to the world coordinates (-1 + iu, -1 + ju, -1 + kv), and the viewpoint with respect to the rectangle Is the farthest from the viewpoint, which matches the values of R ′, G ′, and B ′ calculated based on the equations (28), (29), and (30) with respect to the quadrangle located on the opposite side of When calculating with respect to a square, a background color (for example, Rb = Gb = Bb = 0) is given.

  FIG. 21 is a block diagram illustrating another example of the configuration of the image processing unit 20 according to the present embodiment. As shown in FIG. 21, the image processing unit 20 can be composed of a CPU 201, a ROM 202, a RAM 203, a GPU 204, a video memory 205, a recording medium reading unit 206, and the like. A computer program or data recorded on the recording medium 1 (for example, an optically readable disk storage medium such as a CD-ROM) can be read by the recording medium reading unit 206 (for example, an optical disk drive) and stored in the RAM 203. It may be stored in a hard disk (not shown) and stored in the RAM 203 when the computer program is executed. Computer-readable storage medium such as recording medium 1, ROM 202, RAM 203, hard disk, etc., any non-transitory tangible storage that provides recorded instructions (computer program) and / or data for execution and / or processing on a computer Means medium. By causing the CPU 201 to execute the computer program stored in the RAM 203, the gradation compression processing unit 21, the provisional voxel structure generation unit 22, the voxel structure generation unit 23, the conversion processing unit 24, the rendering processing unit 25, and the 3D texture image Processing performed by the generation unit 26, the layered rectangle setting unit 27, the color interpolation unit 28, and the ROI clipping processing unit 29 can be executed. Note that the CPU 201 is an abbreviation for a central processing unit, and is configured by one or a plurality of cores (processors) and operates in cooperation with the RAM 203. On the other hand, the GPU 204 is a graphics processing unit, and many graphics. It is composed of a processor specialized for processing and operates in cooperation with the video memory 205. Further, the computer program can be downloaded from another computer or a network device via a network such as the Internet, instead of being read by the recording medium reading unit 206.

  Further, the CPU 201 outputs a command such as rendering processing to the GPU 204. The GPU 204 interprets the command and executes necessary processing such as rendering processing. In general, a graphics API (for example, OpenGL) can be used for programming the GPU 204. For example, the computer program can cause the CPU 201 to execute a process of outputting a command for executing the conversion process and the rendering process to the GPU 204.

  Next, a moire generation pattern will be described. The reason why moire occurs in the volume rendering image is that the tomographic images are laminated and synthesized. As a type of moire, there is a curved pattern.

  FIG. 22 is a schematic diagram illustrating an example of a state in which a curved moire is generated. The cause of the occurrence is a pattern having no perceptible regularity that naturally exists in the tomographic image (slice image). By performing coordinate transformation (rotation around the xyz axis) of a plurality of tomographic images having such a pattern, the aforementioned non-regular patterns present in the respective tomographic images are shifted and synthesized. The patterns are connected in a line, and as a result, a curved striped pattern is formed in the volume rendering image, which is recognized as a perceptible moire.

  23 and 24 are explanatory diagrams illustrating a first example of a moire countermeasure result by the 3D texture mapping apparatus 100 according to the present embodiment. FIG. 23A is a comparative example and shows a case where a color map is applied after smoothing processing is performed on voxel values of a DICOM image. FIG. 23B, FIG. 24A, and FIG. 24B show the case where the smoothing process is performed on the provisional voxel structure after applying the color map of this embodiment, and FIG. 23B performs the smoothing process only on the opacity α of the voxel. FIG. 24A shows a case where the smoothing process of the voxel RGB value and the opacity α is performed, and FIG. 24B shows a case where only the voxel RGB value is smoothed. In FIG. 23 and FIG. 24, the above-described color interpolation of the transparent pixels is performed by applying a color map without shading in the OpenGL / GPU version based on the 3D texture mapping method.

  As shown in FIG. 23A, a curved striped pattern is generated in the comparative example. It can be seen in FIG. 23B that the striped pattern has disappeared. Also, the striped pattern disappears in FIG. 24A, but the image appears slightly blurred as compared with FIG. 23B. However, either FIG. 23B or FIG. 24A can be used depending on the user's preference. In FIG. 24B, some striped patterns remain. Therefore, it is preferable to perform smoothing only for the opacity α or all smoothing processing for the RGB value and the opacity α, rather than performing the smoothing processing only for the RGB values of the voxels.

  25 and 26 are explanatory diagrams illustrating a second example of the moire countermeasure result by the 3D texture mapping apparatus 100 according to the present embodiment. FIG. 25A is a comparative example and shows a case where a color map is applied after smoothing processing is performed on voxel values of a DICOM image. FIG. 25B, FIG. 26A, and FIG. 26B show the case where the smoothing process is performed on the provisional voxel structure after the color map of the present embodiment is applied, and FIG. 25B performs the smoothing process only on the opacity α of the voxel. FIG. 26A shows a case where the smoothing process of the voxel RGB value and the opacity α is performed, and FIG. 26B shows a case where only the voxel RGB value is smoothed. In FIG. 25 and FIG. 26, a color map is applied with a shadow in the OpenGL / GPU version based on the 3D texture mapping method, and the above-described color interpolation of the transparent pixels is performed.

  As shown in FIG. 25A, a curved striped pattern is generated in the comparative example. In FIG. 25B, it can be seen that the striped pattern has disappeared. In addition, the striped pattern disappears in FIG. 26A, but the image appears slightly blurred as compared with FIG. 25B. However, either FIG. 25B or FIG. 26A can be used according to the user's preference. In FIG. 26B, some striped patterns remain. Therefore, it is preferable to perform smoothing only for the opacity α or all smoothing processing for the RGB value and the opacity α, rather than performing the smoothing processing only for the RGB values of the voxels.

(Second Embodiment)
In the first embodiment described above, an example of 3D texture mapping has been described as a volume rendering method. However, the voxel smoothing processing of the present embodiment can also be applied to ray casting. Hereinafter, ray casting will be described.

  FIG. 27 is a block diagram showing an example of the configuration of the ray casting apparatus 120 as the image processing apparatus of the present embodiment. The ray casting apparatus 120 according to the present embodiment includes an input unit 10, a storage unit 11, an output unit 12, an image processing unit 30, and the like. The image processing unit 30 includes a gradation compression processing unit 31, a provisional voxel structure generation unit 32, a voxel structure generation unit 33, a coordinate conversion processing unit 34, a ray casting processing unit 35, a search control mask generation unit 36, and an ROI clipping. A processing unit 37 and the like are provided.

  The input unit 10, the storage unit 11, the output unit 12, the gradation compression processing unit 31, the provisional voxel structure generation unit 32, the voxel structure generation unit 33, and the ROI clipping processing unit 37 are the same as those in the first embodiment. Since this is the same, the description is omitted.

FIG. 28 is a schematic diagram showing an example of coordinate conversion of a voxel structure by the ray casting apparatus 120 of the present embodiment. FIG. 28A shows a voxel structure before coordinate conversion described later, and FIG. 28B shows a voxel structure after coordinate conversion. First, the voxel structure before coordinate conversion is a voxel structure generated by the voxel structure generation unit 33. The voxel structure V (x, y, z, c) is V (x, y, z, c) = {Σ _{i = −1} , where the temporary voxel structure is Vo (x, y, z, c). _{, 1} Σ _{j = −1, 1} Σ _{k = −1, 1} Vo (x + i, y + j, z + k, c)} / 27. Here, 0 ≦ x ≦ Sx−1, 0 ≦ y ≦ Sy−1, 0 ≦ z ≦ Sz−1, and 0 ≦ c ≦ 3. Hereinafter, the voxel structure V (x, y, z, c) is also referred to as a three-dimensional voxel image.

  The coordinate conversion processing unit 34 performs predetermined coordinate conversion on the generated three-dimensional voxel image using the coordinate conversion parameter acquired by the input unit 10 to generate a converted voxel image. That is, the coordinate conversion processing unit 34 performs a predetermined coordinate conversion process on the voxel structure V (x, y, z, c) to generate a 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 expressed as Here, 0 ≦ x ≦ Sx−1, 0 ≦ y ≦ Sy−1, 0 ≦ z ′ ≦ Sz′−1, and 0 ≦ c ≦ 3. Note that z ′ is obtained by correcting the z axis so that the resolution matches that of the xy 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. 28A, if the direction of the projection plane where the rendering image is generated and the z-axis direction of the voxel structure V (x, y, z, c) do not match, the later-described ray casting process is complicated. become. Therefore, as shown in FIG. 28B, the xyz axis of the voxel structure V (x, y, z, c) is rotated so that, for example, the z ′ axis after coordinate transformation is perpendicular to the projection plane. Thus, the ray casting process can be simplified.

  FIG. 29 is an explanatory diagram showing an example of coordinate conversion processing by the ray casting apparatus 120 according to the present embodiment. As shown in FIG. 29, the converted voxel structure V ′ (x, y, z ′, c) is obtained by ROI with respect to the voxel structure before conversion with respect to V (x, y, z, c). It can be obtained by performing clipping, setting a 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. In the voxel structure before conversion, the voxel RGB values and opacity other than the region of interest are all considered to be 0, and coordinate conversion processing is performed.

  For example, coordinate transformation can be performed in the order of scaling, z-direction scaling, offset, rotation, and perspective transformation. If the range can be defined by a 4 × 4 transformation matrix, the order of each processing is as follows: It is not limited to the example of FIG.

  Scaling uses the same scale factor Scale in the xyz-axis direction.

  The z-direction scaling process is a process for adjusting the physical interval between pixels in the z direction to the physical interval between pixels in the xy direction. In the z-direction scaling process, z-direction scaling ratio 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. In the offset, an offset Xoff in the x-axis direction, an offset Yoff in the y-axis direction, and an offset Zoff in the z-axis direction are used.

  The rotation is processing for rotating around the x axis, the y axis, and the z axis. In the rotation, a rotation angle Rx about the x axis, a rotation angle Ry about the y axis, and a rotation angle Rz about the z axis (all angle units are radians) are used.

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

  The value of the voxel structure V ′ (x, y, z ′, c) after conversion is determined by the above clipping and coordinate conversion. Note that each process of scaling, z-direction scaling, offset, rotation, and perspective transformation can be summarized in a transformation matrix Matrix (4 × 4).

  The search control mask generation unit 36 performs 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). The minimum value and the maximum value on the line of sight where the opacity α is not 0 (transparent) are specified.

  FIG. 30 is an explanatory diagram illustrating an example of search control mask data generation processing by the ray casting apparatus 120 according to the present embodiment. FIG. 30 shows a case where the voxel structure V ′ (x, y, z ′, c) after conversion is viewed from the xz ′ plane. The projection plane is an xy plane. The search control mask generation unit 36 generates search control mask data M (x, y, s). Here, s = 0 (the start point of a valid voxel range) or 1 (the end point of a valid voxel range).

  As shown in FIG. 30, for each two-dimensional coordinate (x, y), an effective voxel range to be searched is set in the z-axis direction. An effective voxel is a voxel having an opacity α that is not zero. When the opacity α is 0, the opacity α 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, z ′ is incremented by one from z ′ = 0 to z ′ = Sz′-1, Search for z1 satisfying V ′ (x, y, z1, 3)> 0. For the same coordinate (x, y), z ′ is decreased one by one from z ′ = Sz′−1 to 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) may be zero. That is, z1 is the minimum value of z ′ where there is a valid voxel, 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), it is possible to exclude the voxels that are not valid and perform the subsequent ray casting process, thereby reducing the load of the ray casting process. The volume rendering process can be speeded up.

  The ray casting processing unit 35 performs a volume rendering process on the voxel structure (converted voxel structure) to generate a rendering image.

  Specifically, the ray casting processing unit 35 in the z ′ direction toward the voxel structure V ′ (x, y, z ′, c) after conversion of the virtual ray having a predetermined transmission intensity from a predetermined viewpoint. When irradiated, the intensity after passing through the voxel is calculated for each voxel on the z ′ axis while taking into account the absorption / attenuation of the virtual ray based on the opacity α of the voxel, and the transmission intensity falls below a predetermined lower limit. Determine the deepest voxel. Then, while tracing back each voxel through which the virtual ray has passed in the direction of the viewpoint from the voxel at the deepest point, a value obtained by accumulating the RGB values based on the RGB value and opacity α of each voxel is used as the projection plane. A rendering image is generated as pixel values of intersecting pixels. In the process of irradiating a virtual ray toward the voxel structure V ′ (x, y, z ′, c) after conversion in the z ′ direction and determining the deepest voxel, the RGB value and opacity α of each voxel are determined. In addition, a method of simultaneously calculating pixel values of pixels intersecting with the projection plane by accumulatively adding products of transmission intensities of virtual rays. The present application will be described below based on the latter method.

  FIG. 31 is an explanatory diagram showing an example of ray casting processing in the case of parallel projection by the ray casting apparatus 120 of the present embodiment. In FIG. 31, the voxel structure V ′ (x, y, z ′, c) after conversion is illustrated as a cube filled with voxel data three-dimensionally. A projection plane (rendered image) parallel to the xy axis of the voxel structure V ′ (x, y, z ′, c) after conversion is converted into a voxel structure V ′ (x, y, z ′, c) after conversion. It is assumed that they are arranged at a required distance from each other. As shown in FIG. 31, the case where the viewpoint is located at an infinite point is a parallel projection that is generally used (originally Dist = ∞, but set to Dist = 0 in this application), and the inside of the voxel structure from the viewpoint. All the lines of sight toward are parallel to the z ′ axis. In the ray casting process (volume rendering), the inside of the voxel structure V ′ (x, y, z ′, c) after conversion is visualized by performing translucent display (α blending) based on the opacity α. . That is, visualization is realized by the voxel values (RGB value and opacity α) of the voxels of the voxel structure V ′ (x, y, z ′, c) after conversion.

  FIG. 32 is an explanatory diagram showing an example of ray casting processing in the case of perspective projection by the ray casting apparatus 120 of the present embodiment. As shown in FIG. 32, when the viewpoint is relatively close to the projection plane (Dist> 0 in the present application) is perspective projection, the projection plane and viewpoint are voxel structures V ′ (x, It may be immersive inside y, z ', c). As in the case of FIG. 31, in the ray casting process (volume rendering), the voxel structure V ′ (x, y, z ′) after conversion is obtained by performing translucent display (α blending) based on the opacity α. , C) is visualized. That is, visualization is realized by the voxel values (RGB value and opacity α) of the voxels of the voxel structure V ′ (x, y, z ′, c) after conversion.

  As shown in FIGS. 31 and 32, the line of sight is extended from the viewpoint to the inside of the converted voxel structure V ′ (x, y, z ′, c), and at the sampling point (voxel) on the line of sight, Based on the RGB value and the opacity α, the RGB luminance value is cumulatively calculated along the line of sight, and a rendering image is generated using the cumulatively calculated value as a pixel (pixel value of the pixel) at a position where the line of sight intersects the projection plane. . At this time, cumulative calculation is performed for z ′ between the above-described search control masks M (x, y, 0) and M (x, y, 1). Details of the ray casting process will be described later.

  Next, the operation of the ray casting apparatus 120 of the present embodiment will be described.

  33 and 34 are flowcharts showing an example of a procedure of ray casting processing by the ray casting apparatus 120 according to the present embodiment. Hereinafter, for the sake of convenience, the subject of processing will be described as the image processing unit 30. The image processing unit 30 acquires image data of the DICOM image (S41), and acquires color map data (S42). The color map data is illustrated in FIGS. 3 and 4.

  The image processing unit 30 acquires coordinate conversion parameters (S43), and acquires an ROI value in the xyz direction (S44). The coordinate conversion parameters and the ROI values are those illustrated in FIG. The image processing unit 30 performs gradation compression processing of the DICOM image (S45), and performs color map definition for the gradation compressed tomographic image D8 (S46). The color map definition defines RGB values and opacity α for the pixel values of the gradation compressed tomographic image D8. Note that the process of step S45 is not essential, and if the process of step S45 is not performed, color map definition may be performed for the DICOM image.

  The image processing unit 30 generates the provisional voxel structure Vo (x, y, z, c) by applying the color map (S47). The image processing unit 30 performs a voxel smoothing process on the provisional voxel structure Vo (x, y, z, c) to generate the voxel structure V (x, y, z, c) (S48). . The voxel smoothing process may be an opacity only smoothing process or an RGB value and opacity smoothing process.

  The image processing unit 30 determines whether or not to perform the shadow calculation (S49). When performing the shadow calculation (YES in S49), the image processing unit 30 performs the shadow calculation (S50) and performs the process of step S51 described later. When the shadow calculation is not performed (NO in S49), the image processing unit 30 performs a coordinate conversion process on the voxel structure V (x, y, z, c) (S51), and the converted voxel structure V ′. (X, y, z ', c) is generated (S52).

  The image processing unit 30 generates search control mask data M (x, y, s) (S53), performs ray casting processing (S54), and generates a rendering image Image (x, y, c) (S55). The process is terminated. Here, 0 ≦ Image (x, y, c) ≦ 255, c = 0 (R), c = 1 (G), and c = 2 (B).

  The image processing unit 30 of the present embodiment can be configured as illustrated in FIG. 21, similarly to the image processing unit 20 of the first embodiment. However, the GPU 204 and the video memory 205 may not be provided.

  FIG. 35 and FIG. 36 are explanatory diagrams showing a first example of the moire countermeasure result by the ray casting apparatus 120 of the present embodiment. FIG. 35A is a comparative example and shows a case where a color map is applied after smoothing processing is performed on voxel values of a DICOM image. FIG. 35B, FIG. 36A, and FIG. 36B show the case where the smoothing process is performed on the provisional voxel structure after applying the color map of the present embodiment, and FIG. 35B shows the case where the smoothing process is performed only for the opacity α of the voxel. FIG. 36A shows a case where the smoothing process of the voxel RGB value and the opacity α is performed, and FIG. 36B shows a case where only the voxel RGB value is smoothed. In FIGS. 35 and 36, a color map is applied without shading in the CPU version based on the ray casting method.

  As shown in FIG. 35A, a spot-like color shift occurs in the comparative example. In FIG. 35B, it can be seen that the spot-like color shift disappears. Also, the spot-like color shift disappears in FIG. 36A, but the image looks slightly blurred as compared with FIG. 35B. However, either FIG. 35B or FIG. 36A can be used according to the user's preference. In FIG. 36B, some spot-like color misregistration remains. Therefore, it is preferable to perform smoothing only for the opacity α or all smoothing processing for the RGB value and the opacity α, rather than performing the smoothing processing only for the RGB values of the voxels.

  FIGS. 37 and 38 are explanatory diagrams showing a second example of the moire countermeasure result by the ray casting apparatus 120 of the present embodiment. FIG. 37A is a comparative example and shows a case where a color map is applied after performing a smoothing process on the voxel values of the DICOM image. FIG. 37B, FIG. 38A, and FIG. 38B show the case where the smoothing process is performed on the provisional voxel structure after the color map of the present embodiment is applied, and FIG. 37B performs the smoothing process only on the opacity α of the voxel. FIG. 38A shows a case where the smoothing process of the voxel RGB value and opacity α is performed, and FIG. 38B shows a case where only the voxel RGB value is smoothed. In FIG. 37 and FIG. 38, a color map is applied with shading in the CPU version based on the ray casting method.

  As shown in FIG. 37A, a spot-like color shift occurs in the comparative example. In FIG. 37B, it can be seen that the spot-like color shift disappears. Also, the spot-like color shift disappears in FIG. 38A, but the image appears slightly blurred as compared with FIG. 37B. However, either FIG. 37B or FIG. 38A can be used according to the user's preference. In FIG. 38B, some spot-like color misregistration remains. Therefore, it is preferable to perform smoothing only for the opacity α or all smoothing processing for the RGB value and the opacity α, rather than performing the smoothing processing only for the RGB values of the voxels.

  As described above, according to the first and second embodiments, it is possible to suppress the occurrence of moiré when applying a color map while maintaining the same real-time property and processing speed as in the past.

  The computer program according to the present embodiment performs a process of acquiring image data of a plurality of xy plane images captured at predetermined intervals along the z-axis direction, and performs RGB processing on pixel values of the xy plane image. Processing for generating a provisional voxel structure composed of voxels having RGB values and opacity corresponding to each pixel of each of the plurality of xy plane images, and a voxel of the provisional voxel structure On the other hand, processing for generating a voxel structure by correcting the opacity of the voxel based on the opacity of neighboring voxels located in the vicinity of the voxel, and performing volume rendering processing on the voxel structure And processing for generating a rendered image.

  The computer program according to the present embodiment corrects the RGB value of the voxel based on the RGB value of the neighboring voxel located in the vicinity of the voxel with respect to the voxel of the provisional voxel structure. A process of generating a body and a process of generating a rendering image by performing volume rendering processing on the voxel structure.

  The computer program according to the present embodiment causes the voxel to be detected based on the opacity or RGB value of 26 neighboring voxels located in the xyz direction of the voxel with respect to the voxels of the provisional voxel structure. A process of correcting the transparency or the RGB value is executed.

  The computer program according to the present embodiment causes a computer to execute the determination of two neighboring voxels positioned in the vicinity of the first axial direction among the three axial directions of the voxel in the xyz direction with respect to the voxels of the provisional voxel structure. Based on the transparency or RGB value, the processing for correcting the opacity or RGB value of the voxel and the voxel of the provisional voxel structure in which the opacity or RGB value of the voxel is corrected by the two neighboring voxels in the first axis direction On the other hand, the process of correcting the opacity or RGB value of the voxel based on the opacity or RGB value of two neighboring voxels located in the vicinity of the second axial direction of the three axial directions; For the voxels of the temporary voxel structure in which the opacity or RGB value of the voxel is corrected by two neighboring voxels in the axial direction, the three axes Based on the opacity or RGB value of 2 neighboring voxels located in the vicinity of the third axis direction of the direction, to execute a process for correcting the opacity or RGB values of the voxels.

  The computer program according to the present embodiment stores, in the first buffer, a correction value obtained by correcting the opacity or RGB value of each voxel by using two neighboring voxels of each voxel of the provisional voxel structure. The provisional voxel structure according to the immediately preceding correction value obtained by correcting the opacity or RGB value of the immediately preceding voxel by two neighboring voxels of the immediately preceding voxel processed immediately before each voxel of the provisional voxel structure stored in two buffers The process of updating the opacity or RGB value of the immediately preceding voxel and the process of transferring the correction value stored in the first buffer to the second buffer are executed.

  The computer program according to the present embodiment allows a computer to determine a continuous range of numerical values, RGB values, and opacity based on definition information that defines a relationship between a plurality of discrete numerical values, RGB values, and opacity. A process for generating a color map in which the relationship is associated, and a provisional image composed of voxels having RGB values and opacity corresponding to the values of each pixel of each of the plurality of xy plane images based on the generated color map And a process for generating a voxel structure.

  The computer program according to the present embodiment causes the computer to execute the target voxel based on the RGB values of neighboring voxels located in the vicinity of the target voxel with respect to the target voxel having the opacity of the voxel structure being a predetermined value. Processing for correcting the RGB value, processing for generating a 3D texture image composed of voxels having a corrected voxel structure, and performing a predetermined conversion on the 3D texture image to convert the 3D texture A process of generating an image, a process of setting a stacked quadrangle in which quadrilaterals on an XY coordinate plane of a three-dimensional space are arranged in the Z-axis direction, and XY coordinates of the quadrangle on a line of sight parallel to the Z-axis direction from a predetermined viewpoint The RGB values of voxels in the converted 3D texture image corresponding to are in the order of rectangles farther from the viewpoint based on the opacity of the voxels. RGB values obtained by Alpha blending, the line of sight to execute a process of generating a rendered image as a pixel value of the pixel intersects the projection plane.

  The computer program according to the present embodiment corresponds to each pixel of the dummy xy plane image arranged at both ends in the z-axis direction with the plurality of xy plane images in between, and the opacity is set to the predetermined value. Processing for generating a voxel structure to which the dummy voxel is added, and RGB of the neighboring voxels located in the vicinity of the target voxel with respect to the target voxel having the opacity of the predetermined value among the voxels of the generated voxel structure And a process of correcting the RGB value of the target voxel based on the value.

  The computer program according to the present embodiment performs processing for setting a stacked quadrangle in which the same number of the quadrilaterals as the total number of the plurality of xy plane images and the number of the dummy xy plane images are arranged in the Z-axis direction. Let it run.

  The computer program according to the present embodiment generates a voxel structure that includes a voxel corresponding to each pixel at the outer edge of the plurality of xy planar images and having an opacity set to the predetermined value, and Processing for correcting the RGB value of the target voxel based on the RGB value of a neighboring voxel located in the vicinity of the target voxel with respect to the target voxel having an opacity of the predetermined value among the voxels of the voxel structure Is executed.

  The computer program according to the present embodiment is a computer program for setting xyz coordinates that define a target area for volume rendering processing, and among the voxels of the voxel structure, the xyz coordinates are for voxels outside the target area. A process of setting the predetermined value for the opacity of the voxel and a target voxel having an opacity of the predetermined value among the voxels of the voxel structure in which the predetermined value is set, is positioned in the vicinity of the target voxel And a process of correcting the RGB value of the target voxel based on the RGB value of the neighboring voxel to be executed.

  The computer program according to the present embodiment causes the computer to calculate the xyz coordinates of the 3D texture image corresponding to the rectangular XY coordinates on the line of sight, and the converted 3D texture image corresponding to the calculated xyz coordinates. And alpha blending based on the RGB value and opacity of the voxel corresponding to the xyz coordinates.

  The computer program according to the present embodiment causes the computer to calculate a gradient vector related to the opacity of the voxel of the voxel structure, and to calculate the shadow value of the voxel based on the calculated gradient vector and a predetermined light source vector. A process and a process of correcting the RGB value of the voxel of the voxel structure based on the calculated shadow value are executed.

  The computer program according to the present embodiment causes the computer to calculate a difference between the opacity of the voxel of the voxel structure and the opacity of each neighboring voxel located in the vicinity of the voxel, and the voxel and the neighborhood A process of calculating a gradient vector of the voxel based on a distance from the voxel and a difference in opacity between the voxel and the neighboring voxel is executed.

  The computer program according to the present embodiment, for a voxel that cannot calculate a gradient vector related to the opacity of the voxel of the voxel structure in the computer, sets the predetermined value to the opacity of the voxel; For the target voxel having the predetermined value of opacity among the voxels of the voxel structure in which the predetermined value is set, the RGB of the target voxel is determined based on the RGB values of neighboring voxels located in the vicinity of the target voxel. And a process of correcting the value.

  The computer program according to the present embodiment causes the computer to rotate the x-axis center, the y-axis center, and the z-axis center, the viewing angle, the viewpoint position, the clipping position, the offset value in the xyz axis direction, and the enlargement or reduction in the xyz axis direction. Processing for acquiring the predetermined conversion parameters including the magnification and the scaling factor in the z-axis direction, and performing predetermined conversion using the acquired parameters on the 3D texture image to generate a converted 3D texture image Process.

  The computer program according to the present embodiment causes a computer to execute a process of generating a post-conversion 3D texture image and a command for executing a process of generating the rendering image to the GPU.

  The computer program according to the present embodiment performs a process of generating a voxel image by performing predetermined coordinate transformation on the voxel structure and irradiating a virtual ray from a predetermined viewpoint toward the voxel image. In addition, for each voxel on the line of sight, the RGB value calculated based on a value obtained by accumulating the accumulated luminance value by the product of the RGB value and opacity of the voxel and the transmission intensity of the virtual ray, the line of sight is the projection plane. And generating a rendering image as the pixel value of the intersecting pixel.

  The computer program according to the present embodiment causes a computer to specify a minimum value and a maximum value of pixels of a predetermined xy plane image among the plurality of xy plane images, an upper limit value smaller than the specified maximum value, and A process of calculating a lower limit value greater than the specified minimum value, a process of compressing pixel values of each pixel of the plurality of xy plane images within the range of the upper limit value and the lower limit value, and the compressed plurality of xy planes And a process of generating a provisional voxel structure composed of voxels having RGB values and opacity corresponding to each pixel of each image.

  The image processing apparatus according to the present embodiment is an image processing apparatus that performs volume rendering processing, and an acquisition unit that acquires image data of a plurality of xy plane images captured at predetermined intervals along the z-axis direction; RGB values and opacity are defined for the pixel values of the xy plane image, and a provisional voxel structure composed of voxels having RGB values and opacity corresponding to each pixel of the plurality of xy plane images is generated. And generating a voxel structure for the voxel of the temporary voxel structure by correcting the opacity of the voxel based on the opacity of a neighboring voxel located in the vicinity of the voxel A voxel structure generation unit for performing rendering processing for generating a rendering image by performing volume rendering processing on the voxel structure And an image generation unit.

  The image processing method according to the present embodiment is an image processing method that performs volume rendering processing, acquires image data of a plurality of xy plane images captured at predetermined intervals along the z-axis direction, and the xy plane image RGB values and opacity are determined with respect to the pixel values of the plurality of xy plane images, and a provisional voxel structure including voxels having RGB values and opacity corresponding to each pixel of the plurality of xy plane images is generated, For the voxel of the temporary voxel structure, the voxel structure is generated by correcting the opacity of the voxel based on the opacity of the neighboring voxel located in the vicinity of the voxel, and the volume for the voxel structure. A rendering process is performed to generate a rendering image.

  The voxel data of the present embodiment is voxel data of a voxel structure, and RGB values and opacity are determined with respect to pixel values of a plurality of xy plane images captured at predetermined intervals along the z-axis direction. The voxel structure includes voxel data composed of voxels having RGB values and opacity corresponding to the respective pixels of the plurality of xy plane images, and the voxel data includes the voxel data of the voxel structure. The voxel includes data in which the opacity of the voxel is corrected based on the opacity of a neighboring voxel located in the vicinity of the voxel, and performs a volume rendering process on the voxel structure to generate a rendered image. Used to execute the process to generate.

  In the voxel data according to the present embodiment, the voxel data is data obtained by correcting the RGB values of the voxels based on the RGB values of neighboring voxels located in the vicinity of the voxel with respect to the voxels of the voxel structure. including.

  The voxel data of the present embodiment is voxel data of a voxel structure, and RGB values and opacity are determined with respect to pixel values of a plurality of xy plane images captured at predetermined intervals along the z-axis direction. A voxel data of a provisional voxel structure composed of voxels having RGB values and opacity corresponding to each pixel of each of the plurality of xy plane images, the voxel of the provisional voxel structure A process for correcting the opacity or RGB value of the voxel based on the opacity or RGB value of two neighboring voxels located in the vicinity of the first axial direction of the three axial directions of the voxel in the xyz direction; , The three axes for the voxel of the temporary voxel structure in which the opacity or RGB value of the voxel is corrected by two neighboring voxels in the first axis direction Processing for correcting the opacity or RGB value of the voxel based on the opacity or RGB value of two neighboring voxels located in the vicinity of the second axis direction in the direction, and two neighboring voxels in the second axis direction Based on the opacity or RGB value of two neighboring voxels located in the vicinity of the third axial direction of the three axial directions with respect to the voxel of the provisional voxel structure in which the opacity or RGB value of the voxel is corrected. , And a process for correcting the opacity or RGB value of the voxel.

  The voxel data of the present embodiment includes a process of storing a correction value obtained by correcting the opacity or RGB value of each voxel with two neighboring voxels of each voxel of the provisional voxel structure in the first buffer, and the second buffer. Immediately before the provisional voxel structure by the immediately preceding correction value obtained by correcting the opacity or RGB value of the immediately preceding voxel by two neighboring voxels of the immediately preceding voxel processed immediately before each voxel of the provisional voxel structure stored. It is used to execute a process of updating the opacity or RGB value of the voxel and a process of transferring the correction value stored in the first buffer to the second buffer.

DESCRIPTION OF SYMBOLS 10 Input part 11 Storage part 12 Output part 20, 30 Image processing part 21, 31 Gradation compression processing part 22, 32 Temporary voxel structure generation part 23, 33 Voxel structure generation part 24 Conversion processing part 25 Rendering process part 26 3D Texture image generation unit 27 Laminated rectangle setting unit 28 Color interpolation unit 29, 37 ROI clipping processing unit 34 Coordinate transformation processing unit 35 Ray casting processing unit 36 Search control mask generation unit 100 3D texture mapping device 120 Ray casting device 201 CPU
202 ROM
203 RAM
204 GPU
205 Video memory 206 Recording medium reader

Claims (21)

On the computer,
processing for acquiring image data of a plurality of xy plane images captured at predetermined intervals along the z-axis direction;
Based on a color map that correlates the relationship between a numerical value in a predetermined range, RGB values, and opacity, the voxel having RGB values and opacity corresponding to the values of each pixel of each of the plurality of xy plane images. Processing to generate a provisional voxel structure;
A process for generating a voxel structure by correcting the opacity of the voxel based on the opacity of a neighboring voxel located in the vicinity of the voxel for the voxel of the temporary voxel structure;
A computer program for executing a volume rendering process on the voxel structure to generate a rendered image. On the computer,
For the voxels of the temporary voxel structure, based on the RGB values of neighboring voxels located in the vicinity of the voxel, processing to generate a voxel structure by correcting the RGB values of the voxels;
The computer program according to claim 1, wherein a volume rendering process is performed on the voxel structure to generate a rendered image. On the computer,
A process for correcting the opacity or RGB value of the voxel based on the opacity or RGB value of 26 neighboring voxels located in the xyz direction vicinity of the voxel for the voxel of the provisional voxel structure. Item 3. A computer program according to item 1 or 2. On the computer,
processing for acquiring image data of a plurality of xy plane images captured at predetermined intervals along the z-axis direction;
A provisional voxel structure composed of voxels having RGB values and opacity corresponding to respective pixels of each of the plurality of xy plane images, wherein RGB values and opacity are determined with respect to the pixel values of the xy plane image. Process to generate,
A process for generating a voxel structure by correcting the opacity of the voxel based on the opacity of a neighboring voxel located in the vicinity of the voxel for the voxel of the temporary voxel structure;
Processing to generate a rendering image by performing volume rendering processing on the voxel structure;
Based on the opacity or RGB value of two neighboring voxels located in the vicinity of the first axial direction of the three axial directions of the voxel in the xyz direction with respect to the voxels of the provisional voxel structure. Processing to correct transparency or RGB values;
Two neighbors located in the vicinity of the second axial direction of the three axial directions with respect to the voxel of the temporary voxel structure in which the opacity or RGB value of the voxel is corrected by the two neighboring voxels in the first axial direction Correcting the voxel opacity or RGB value based on the voxel opacity or RGB value;
Two neighbors located in the vicinity of the third axial direction of the three axial directions with respect to the voxel of the temporary voxel structure in which the opacity or RGB value of the voxel is corrected by the two neighboring voxels in the second axial direction based on the opacity or RGB values of the voxels, Turkey computer program to execute the processing for correcting opacity or RGB values of the voxels. On the computer,
Storing a correction value obtained by correcting the opacity or RGB value of each voxel by two neighboring voxels of each voxel of the provisional voxel structure in the first buffer;
The provisional voxel structure based on the immediately preceding correction value obtained by correcting the opacity or RGB value of the immediately preceding voxel by two neighboring voxels of the immediately preceding voxel processed immediately before each voxel of the provisional voxel structure stored in the second buffer Updating the opacity or RGB value of the previous voxel of the body;
The computer program according to claim 4, wherein the correction value stored in the first buffer is transferred to the second buffer. On the computer,
processing for acquiring image data of a plurality of xy plane images captured at predetermined intervals along the z-axis direction;
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;
Processing to calculate an upper limit value smaller than the identified maximum value and a lower limit value greater than the identified minimum value;
Processing 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;
Generating a provisional voxel structure composed of voxels having RGB values and opacity corresponding to each pixel of each of the plurality of compressed xy planar images ;
A process for generating a voxel structure by correcting the opacity of the voxel based on the opacity of a neighboring voxel located in the vicinity of the voxel for the voxel of the temporary voxel structure;
Turkey computer program by performing volume rendering processing to execute the processing for generating the rendered image to the voxel structures. On the computer,
Processing to perform a predetermined coordinate transformation on the voxel structure to generate a voxel image;
A value obtained by accumulatively adding a cumulative luminance value 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 when a virtual ray is irradiated toward the voxel image from a predetermined viewpoint. The computer program according to any one of claims 1 to 6, wherein a process of generating a rendering image using the RGB value calculated based on the image as a pixel value of a pixel whose line of sight intersects a projection plane is executed. . On the computer,
processing for acquiring image data of a plurality of xy plane images captured at predetermined intervals along the z-axis direction;
A provisional voxel structure composed of voxels having RGB values and opacity corresponding to respective pixels of each of the plurality of xy plane images, wherein RGB values and opacity are determined with respect to the pixel values of the xy plane image. Process to generate,
A process for generating a voxel structure by correcting the opacity of the voxel based on the opacity of a neighboring voxel located in the vicinity of the voxel for the voxel of the temporary voxel structure;
Processing for correcting the RGB value of the target voxel based on the RGB value of a neighboring voxel located in the vicinity of the target voxel for a target voxel in which the opacity of the voxel structure is a predetermined value;
Processing to generate a 3D texture image composed of voxels of a voxel structure with corrected RGB values;
A process of performing a predetermined conversion on the 3D texture image to generate a converted 3D texture image;
A process of setting a stacked quadrangle in which quadrangles on the XY coordinate plane of the three-dimensional space are arranged in the Z-axis direction;
Alpha blending of RGB values of voxels of the converted 3D texture image corresponding to the XY coordinates of the rectangle on the line of sight parallel to the Z-axis direction from a predetermined viewpoint in the order of the rectangles farthest from the viewpoint based on the opacity of the voxels the RGB values obtained by said visual axis is benzalkonium computer program to execute a process of generating a rendered image as a pixel value of the pixel intersects the projection plane. On the computer,
A voxel structure is generated in which dummy voxels corresponding to the pixels of the dummy xy plane image disposed at both ends in the z-axis direction with the plurality of xy plane images in between are added with dummy voxels having the predetermined opacity. Processing,
Processing for correcting the RGB value of the target voxel based on the RGB value of a neighboring voxel located in the vicinity of the target voxel with respect to the target voxel having an opacity of the predetermined value among the voxels of the generated voxel structure The computer program according to claim 8 , wherein: On the computer,
The computer program according to claim 8, wherein a process for setting a stacked quadrangle in which the same number of the quadrangle as the total number of the plurality of xy plane images and the number of the dummy xy plane images is arranged in the Z-axis direction is executed. On the computer,
A process of generating a voxel structure including voxels corresponding to each pixel of an outer edge portion of the plurality of xy plane images and having an opacity set to the predetermined value;
Processing for correcting the RGB value of the target voxel based on the RGB value of a neighboring voxel located in the vicinity of the target voxel with respect to the target voxel having an opacity of the predetermined value among the voxels of the generated voxel structure The computer program according to claim 8 , wherein: On the computer,
A process of setting xyz coordinates that define a target area of the volume rendering process;
A process of setting the predetermined value as the opacity of the voxel with respect to voxels whose xyz coordinates are out of the target area among the voxels of the voxel structure;
For the target voxel having the predetermined value of opacity among the voxels of the voxel structure in which the predetermined value is set, the RGB of the target voxel is determined based on the RGB values of neighboring voxels located in the vicinity of the target voxel. The computer program according to claim 8 , which executes a process of correcting a value. On the computer,
A process of calculating xyz coordinates of the 3D texture image corresponding to XY coordinates of the rectangle on the line of sight;
The computer program according to claim 8 , wherein alpha blending is executed based on an RGB value and opacity of a voxel corresponding to the xyz coordinate of the converted 3D texture image corresponding to the calculated xyz coordinate. On the computer,
Processing to calculate a gradient vector related to the opacity of the voxel of the voxel structure;
A process of calculating a shadow value of the voxel based on the calculated gradient vector and a predetermined light source vector;
The computer program according to claim 8 , wherein a process of correcting an RGB value of a voxel of the voxel structure based on the calculated shadow value is executed. On the computer,
Processing to calculate the difference between the opacity of the voxel of the voxel structure and the opacity of each neighboring voxel located in the vicinity of the voxel;
The computer program according to claim 8 , wherein a process of calculating a gradient vector of the voxel based on a distance between the voxel and the neighboring voxel and a difference in opacity between the voxel and the neighboring voxel is executed. On the computer,
Processing for setting the predetermined value to the opacity of the voxel for a voxel that cannot calculate a gradient vector related to the opacity of the voxel of the voxel structure;
For the target voxel having the predetermined value of opacity among the voxels of the voxel structure in which the predetermined value is set, the RGB of the target voxel is determined based on the RGB values of neighboring voxels located in the vicinity of the target voxel. The computer program according to claim 8, which executes a process of correcting a value. On the computer,
The predetermined including the rotation angle of the x-axis center, the y-axis center, and the z-axis center, the viewing angle, the viewpoint position, the clipping position, the offset value in the xyz axis direction, the enlargement or reduction magnification in the xyz axis direction, and the scaling factor in the z axis direction. Process to get the conversion parameters of
The computer program according to claim 8 , wherein a predetermined conversion using the acquired parameter is performed on the 3D texture image to generate a converted 3D texture image. On the computer,
The computer program according to claim 8 , wherein the computer program causes the GPU to execute a process for generating the converted 3D texture image and a process for generating the rendered image. On the computer,
A process for generating a color map that associates a relationship between a continuous numerical value in a predetermined range, an RGB value, and opacity based on definition information that defines a relationship between a plurality of discrete numerical values and RGB values and opacity When,
Based on the generated color map, according to claim 4 for executing the process of generating the composed provisional voxel structure voxels having RGB values and opacity corresponding to the value of each pixel of the plurality of the xy plane image, respectively The computer program according to claim 18 . An image processing apparatus that performs volume rendering processing,
an acquisition unit that acquires image data of a plurality of xy plane images captured at predetermined intervals along the z-axis direction;
Based on a color map that correlates the relationship between a numerical value in a predetermined range, RGB values, and opacity, the voxel having RGB values and opacity corresponding to the values of each pixel of each of the plurality of xy plane images. A provisional voxel structure generation unit for generating a provisional voxel structure;
A voxel structure generating unit that generates a voxel structure by correcting the opacity of the voxel based on the opacity of a neighboring voxel located in the vicinity of the voxel with respect to the voxel of the temporary voxel structure,
An image processing apparatus comprising: a rendering image generation unit configured to perform a volume rendering process on the voxel structure to generate a rendering image. An image processing method for performing volume rendering processing,
obtaining image data of a plurality of xy plane images captured at predetermined intervals along the z-axis direction;
Based on a color map that correlates the relationship between a numerical value in a predetermined range, RGB values, and opacity, the voxel having RGB values and opacity corresponding to the values of each pixel of each of the plurality of xy plane images. Create a provisional voxel structure
For the voxels of the provisional voxel structure, based on the opacity of neighboring voxels located in the vicinity of the voxel, correct the opacity of the voxels to generate a voxel structure,
An image processing method for generating a rendering image by performing volume rendering processing on the voxel structure.