JP7206846B2 - volume rendering device - Google Patents

Description

The present disclosure relates to volume rendering technology for three-dimensional visualization using a plurality of two-dimensional tomographic images.

Conventionally, multiple images are taken continuously at predetermined slice intervals by medical imaging diagnostic equipment (modality) such as CT, MRI, and PET, and stored in a medical information system such as PACS (Picture Archiving and Communication Systems) in DICOM format. Three-dimensional visualization of organs and the like is performed based on the tomographic image.

In the field of 3D computer graphics, a rendering image is generated with limited colors such as 256 colors, and the colors of the rendering image are changed in real time by interactively changing the color map (color lookup table) implemented in hardware. method is used (see Patent Document 1).

A method of rendering with a VIS volume (a voxel having three parameters of a signal value, a shadow luminance value, and an opacity intensity coefficient) that does not depend on a color map has also been proposed (see Patent Document 2).

JP-A-1-321577 Japanese Patent No. 2651787

However, the technique described in Patent Document 1 can only change the RGB values, and does not disclose controlling the α value expressing opacity, which is important in volume rendering. Moreover, in the technology described in Patent Document 2, suppression of the processing load is more limited than in the case of re-rendering by applying a color map to signal values.

Therefore, an object of the present disclosure is to provide a volume rendering apparatus capable of improving the speed of rendering image generation processing when generating a rendering image by applying a color map to a plurality of tomographic images. and

The present disclosure includes multiple means for solving the above problems, but if one example is given,
Rendering with reference to a color map defined by associating color component values and opacity with signal values, based on multiple 2D tomographic images of a given object taken at given intervals. A volume rendering device for generating an image, comprising:
calculating an opaque circumscribed region, which is a region that circumscribes in each three-dimensional axial direction the region in which pixels having non-zero opacity obtained by referring to the color map exist among the pixels of each of the tomographic images; opaque circumscribed area calculation means;
Among the tomographic images, a region corresponding to the opaque circumscribed region is defined as a voxelized region, and the signal value of each pixel in the voxelized region is converted by referring to the color map, and color components and non-color components are obtained as pixel values. voxel image creation means for creating a voxel image in which two-dimensional element images with defined transparency are superimposed at predetermined intervals;
and rendering means for generating a rendered image using the voxel image.

Also, in this disclosure:
A program for causing a computer to function as the volume rendering device is provided.

Also, in this disclosure:
A voxel structure composed of voxels determined by referring to a predefined color map corresponding to each pixel of a plurality of two-dimensional tomographic images of a predetermined object photographed at predetermined intervals. voxel data,
All the voxels on the outermost perimeter of the xyz axes have an opacity of 0, providing voxel data having at least one non-zero opacity voxel on a surface adjacent to the boundary.

According to the present disclosure, it is possible to apply a color map to a plurality of tomographic images and perform processing for generating rendering images at high speed.

FIG. 3 is a diagram showing a plurality of tomographic images and a color rendering image generated as a volume rendering image; 1 is a hardware configuration diagram of a volume rendering device 100 according to an embodiment of the present disclosure; FIG. 1 is a functional block diagram showing the configuration of a volume rendering device according to an embodiment; FIG. 4 is a flowchart showing processing operations of the volume rendering device according to the embodiment; It is a figure which shows the color map used by this embodiment. FIG. 5 is a flow chart showing a case where opaque circumscribed area calculation processing in step S30 shown in FIG. 4 is performed in parallel by a plurality of CPU cores. FIG. FIG. 7 is a flow chart showing calculation processing of an opaque circumscribed region for a group of tomographic images in step S330 of FIG. 6. FIG. FIG. 7 is a flow chart showing calculation processing of an opaque circumscribed region of the entire group of tomographic images in step S360 of FIG. 6. FIG. FIG. 11 compares ROIs and opaque circumscribing regions; 5 is a flowchart showing a case where voxel image creation processing and the like in step S40 shown in FIG. 4 are performed in parallel. FIG. 7 is a flowchart showing element image creation processing in steps S510 and S530 shown in FIG. 6; FIG. FIG. 4 is a diagram showing correspondence in 3D texture mapping; FIG. 4 is a diagram showing the configuration of the rendering means 70 when performing 3D texture mapping method processing; 4 is an explanatory diagram showing an example of projection screen setting by rendering means 70 of the present embodiment; FIG. 4 is a flowchart showing details of rendering processing by a 3D texture mapping method; FIG. 4 is a diagram showing parallel processing when the volume rendering device 100 is implemented by a computer; FIG. 4 is a diagram showing a software flow of parallel processing when the volume rendering device 100 is implemented by a computer.

Preferred embodiments of the present disclosure will be described in detail below with reference to the drawings.
The present disclosure generates a rendering image by referring to a color map defined by associating color component values and opacities with signal values based on a plurality of tomographic images. FIG. 1 is a diagram showing a plurality of tomographic images and a rendering image. FIG. 1(a) shows a plurality of tomographic images, FIG. 1(b) shows a single tomographic image, and FIG. 1(c) shows a rendering image. The tomographic images shown in FIGS. 1A and 1B are in DICOM format. FIG. 1(a) is a group of 370 tomographic images, and FIG. 1(b) is an enlarged view of one of them. A volume rendering apparatus according to an embodiment, which will be described later, executes rendering processing on a plurality of tomographic images as shown in FIG. Generate. The rendered image shown in FIG. 1(c) has 512×512 pixels, and each pixel is represented by a total of 24 bits of 8 bits for each color of RGB.

<1. Device configuration>
FIG. 2 is a hardware configuration diagram of the volume rendering device 100 according to an embodiment of the present disclosure. The volume rendering apparatus 100 according to this embodiment can be realized by a general-purpose computer, and as shown in FIG. , a large-capacity storage device 3 such as a hard disk, SSD (Solid State Drive), flash memory, etc. for storing programs and data executed by the CPU 1, and an instruction input I / F (interface) 4 such as a keyboard, mouse, etc. , a data input/output I/F (interface) 5 for data communication with an external device such as a data storage medium, a display unit 6 which is a display device such as a liquid crystal display, and an arithmetic processing unit specializing in graphics. A GPU (Graphics Processing Unit) 7 and a frame memory 8 for holding an image to be displayed on the display unit 6 are provided, and are connected to each other via a bus. Since the calculation result by the GPU 7 is written to the frame memory 8, the GPU 7 and the frame memory 8 are often mounted on a video card having an interface to the display unit 6 and attached to a general-purpose computer via a bus.

In this embodiment, the CPU 1 is a multi-core CPU, has a plurality of CPU cores, and is capable of parallel processing. Although only one RAM2 is shown in the example of FIG. 2, each CPU core of the CPU1 is configured to access one RAM2. Note that the number of CPUs 1 may be plural. A multi-core CPU may also be a CPU in which each CPU core logically accepts multiple threads of execution.

FIG. 3 is a functional block diagram showing the configuration of the volume rendering device according to this embodiment. 3, 10 is tomographic image reading means, 20 is color map reading means, 30 is ROI clipping setting means, 40 is mask setting means, 50 is opaque circumscribed area calculation means, 60 is voxel image creation means, and 61 is opacity. 62 is shading means; 70 is rendering means; and 80 is rendered image output means.

The tomographic image reading means 10 is means for reading a plurality of tomographic images from PACS collected and stored by medical image diagnostic equipment such as CT, MRI, and PET via a storage medium and an input unit. A tomographic image is obtained by photographing an object at predetermined intervals, and is a two-dimensional tomographic image in which a signal value is assigned to each pixel. The color map reading means 20 is a means for reading a color map created in advance by a color map creating means (not shown) from a storage medium or an input section. The ROI clipping setting means 30 is a means for setting an ROI, which is a region of interest, as a voxel image creation target among a plurality of read tomographic images. The mask setting means 40 is a means for setting a mask for an area not to be displayed among the plurality of read tomographic images. The opaque circumscribing region calculating means 50 refers to a color map defined by associating color component values and opacities with signal values, associates opacities with respective pixels of a plurality of tomographic images, and calculates the A means of calculating an opaque bounding area based on opacity. Although the details will be described later, the opaque circumscribing area is an area that circumscribes a range in which pixels whose opacity is not 0 exists in each three-dimensional axial direction.

The voxel image creating means 60 refers to a color map defined by associating color component values and opacities with signal values, and refers to each pixel of a plurality of tomographic images arranged three-dimensionally. This is a means of creating a voxel image by giving color component values and opacity to voxels. The voxel image creating means 60 further comprises opacity correcting means 61 and shadow adding means 62 . The opacity correction means 61 is means for correcting the opacity α of each pixel of the voxel image using a predetermined table. The shading means 62 calculates a gradient vector for the color component (R, G, B) of each voxel of the voxel image based on the opacity α of voxels in the vicinity of the voxel. Then, using the calculated gradient vector, the shade value is calculated, and each value of the color components (R, G, B) is multiplied by the shade value to replace the value.

The tomographic image reading means 10 and the color map reading means 20 are realized mainly in the data input/output I/F 5 with the assistance of the CPU 1 . The ROI clipping setting means 30 and the mask setting means 40 are realized mainly in the instruction input I/F 4 with the assistance of the CPU 1 . The opaque circumscribed area calculation means 50 , the voxel image creation means 60 , the opacity correction means 61 and the shadow addition means 62 are implemented by the CPU 1 executing programs stored in the storage device 3 . The rendering means 70 is implemented mainly by executing a program on the GPU 7 with the assistance of the CPU 1 . The rendered image output means 80 is realized mainly in the frame memory 8 and the display section 6 with the assistance of the CPU 1 .

Each component shown in FIG. 3 is actually implemented by installing a dedicated program in hardware such as a computer and its peripherals, as shown in FIG. That is, the computer executes the contents of each means according to a dedicated program. In this specification, a computer means a device that has an arithmetic processing unit such as a CPU or GPU and is capable of data processing. It also includes computers embedded in mobile terminals and various devices.

<2. Processing operation>
Next, processing operations of the volume rendering apparatus shown in FIGS. 2 and 3 will be described. FIG. 4 is a flowchart showing processing operations of the volume rendering apparatus according to this embodiment. First, the tomographic image reading means 10 reads a plurality of tomographic images. The plurality of tomographic images are in the form shown in FIG. 1(a). In this embodiment, a DICOM format tomographic image group Do (x, y, z) (-32768≤Do (x, y, z)≤32767, 0≤x≤Sx-1, 0≤y≤Sy-1, 0≦z≦Sz−1; resolution: Rxy, Rz). Since one pixel of the tomographic image group Do(x, y, z) is recorded with 16 bits, signal values Do(x, y, z) of -32768 to 32767 are recorded. The tomographic image group Do (x, y, z) is obtained by stacking Sz tomographic images each having an x-coordinate size (number of pixels) of Sx and a y-coordinate size (number of pixels) of Sy. In the following description, the number of tomographic images Sz is also expressed as M. The horizontal direction of one tomographic image in FIG. 1(b) corresponds to the x-axis, and the vertical direction corresponds to the y-axis. The lamination direction when superimposed at a predetermined slice interval (difference in Image Position tag information of adjacent tomographic images) corresponds to the z-axis.

After reading a plurality of DICOM-format tomographic images, the ROI clipping setting means 30 sets the ROI clipping (step S10). ROI is an abbreviation for Region of Interest, and in this specification means a region of interest in a three-dimensional space (scientifically, ROI refers to a two-dimensional region, and in the case of a three-dimensional region, VOI (Volume of Interest) , but in the medical imaging industry, ROI is often used even in a three-dimensional area). Here, the ranges for which rendering images and voxel images are to be created are shown. By clipping (cutting) and setting the ROI, it is possible to generate a volume rendering image in which the subject is cropped at any desired position three-dimensionally, and to depict the internal organs and organs hidden behind the body surface and bones. used for In the processing after step S40, the ROI is set in step S10 so that the voxel image creation target is limited to a predetermined range. This setting is performed by specifying each coordinate range in each of the xyz axis directions with numerical values. The numeric value may be specified directly, or may be specified by pointing a position on the screen with a mouse or the like and converting it into a numeric value via the instruction input I/F 4 .

Xso is the coordinate corresponding to the left end in the x-axis direction of each read tomographic image, Xeo is the coordinate corresponding to the right end in the x-axis direction, Yso is the coordinate corresponding to the upper end in the y-axis direction, and Yso is the coordinate corresponding to the lower end in the y-axis direction. Let Yeo be the coordinates, Zso be the coordinates corresponding to the tomographic image of the first slice (first: first), and Zeo be the coordinates corresponding to the last slice (end: Mth) tomographic image. A range defined by a rectangular parallelepiped of [Xs, Xe], [Ys, Ye], and [Zs, Ze] satisfying ≦Xeo, Yso≦Ys<Ye≦Yeo, Zso≦Zs<Ze≦Zeo is defined as an ROI and processed. Identify as a target. That is, six coordinates Xs, Xe, Ys, Ye, Zs, and Ze are defined that specify the voxel image area for which the rendering image is to be generated. These six coordinates define a rectangular parallelepiped with sides along the xyz axes, and the inside of the rectangular parallelepiped becomes the clipped ROI. In this embodiment, an opaque circumscribed area is calculated in step S30, which will be described later, and this opaque circumscribed area is also a rectangular parallelepiped area defined by a predetermined range on each of the xyz axes.

As a result of the ROI clipping setting in step S10, the number of pixels to be processed is (Xeo−Xso+1)×(Yeo−Yso+1)×(Zeo−Zso+1) pixels, and (Xe−Xs+1)×(Ye−Ys+1)×(Ze− Zs+1) pixels. By performing the ROI clipping setting in step S10 by the ROI clipping setting means 30, in addition to the effect of obtaining a volume rendering image in which the organ to be observed is exposed, the number of pixels to be processed is reduced, thereby suppressing the processing load and improving the responsiveness. Since there is a double merit that it can be improved, it is often set, but it is possible to omit it because it is not an essential process. At the stage of ROI clipping setting in step S10, only the range to which ROI clipping is applied is defined, and specific ROI clipping processing is executed in each processing after step S40. In each process, the processing load is reduced by performing calculations only within the ROI of the voxel image, and a clipped (cut) volume rendering image is generated in the rendering process of step S60.

After the ROI clipping setting is performed by the ROI clipping setting means 30, the mask setting means 40 next performs mask setting (step S20). A mask setting is a setting as to whether or not to display each pixel of a tomographic image. The mask setting method is not particularly limited, but in this embodiment, the group of tomographic images read in step S10 is displayed on the screen as an MPR image (three screens of axial, coronal, and sagittal), and the method described in this specification is displayed. Alternatively, a volume rendering image created by a known method is displayed, an image that is easiest to indicate is selected from a plurality of these two-dimensionally projected images, and a two-dimensional region of the selected image is indicated while masking. to specify the range that should not be displayed and the area that should be displayed without masking. With this mask setting, mask data in a three-dimensional space corresponding to the group of tomographic images is obtained. The mask data Mask (x, y, z) is set to either "0" or "1" corresponding to each pixel of the tomographic image group Do (x, y, z). . Although the process of creating such a three-dimensional mask is almost essential, it can be omitted. Although the above-mentioned ROI clipping is also available as a method of depicting an organ to be observed, the ROI clipping process is often used in combination because the operation is simpler than that of the three-dimensional mask and the processing load is reduced. A specific area to be displayed is set by the ROI clipping setting means 30 and the mask setting means 40. If both are set, the area set as the area to be displayed is the effective area. , and if either one is set, the set area becomes the valid area. When neither the ROI clipping setting means 30 nor the mask setting means 40 is used, the entire area of the loaded tomographic image group becomes the effective area.

Next, the opaque circumscribed area calculation means 50 performs opaque circumscribed area calculation processing (step S30). An opaque circumscribing area is an area that circumscribes a range in which pixels whose opacity is not 0 exists in each three-dimensional axial direction. The direction of the circumscribing side (line segment) is not particularly limited as long as the opaque circumscribing region is a region circumscribing the range in which pixels whose opacity is not 0 exist. This circumscribing side is not limited to a straight line and may be a curved line. However, it is preferable to set the region inside a rectangular parallelepiped defined by the directions of the two coordinate axes of the tomographic images and the side along the direction in which the tomographic images are stacked. As a specific process, using the color map read from the color map reading means 20, the opacity corresponding to the signal value of each pixel included in the ROI of the tomographic image after clipping setting is applied. Calculate the opaque bounding area based on the opacity. FIG. 5 is a diagram showing a color map used in this embodiment. In FIG. 5, the signal value is the signal value recorded in each pixel of the tomographic image, R, G, and B are the color components of red, green, and blue, respectively, and α is the degree of display of the voxel. Opacity to specify. FIG. 5 is a color map when signal values of a tomographic image are recorded in 16 bits and take values from -32768 to 32767. FIG. In the case of a CT image, it is often normalized to a range from -2048 to 2048 with the water component set to 0. As shown in FIG. 5, in the color map, the signal value of each pixel of the tomographic image is associated with each of the R, G, and B color components and the opacity α. This color map is also used in voxel image creation processing by the voxel image creation means 60 (described later: step S40). In the process of calculating the opaque circumscribed area in step S30, only the opacity α of the color map is used. By using the opacity α obtained by referring to the color map, it is possible to calculate, as an opaque circumscribed region, an area in which there are opaque pixels where α=0 in the tomographic image. A voxel image is created in step S40 using the area corresponding to this opaque circumscribed area as a voxelized area. As the "region corresponding to the opaque circumscribing region", the opaque circumscribing region can be used as it is as a voxelized region to be voxelized. The area that is covered by the rectangle is defined as the voxelized area corresponding to the opaque circumscribed area.

Next, the voxel image creation means 60 performs voxel image creation processing (step S40). A voxel image is a voxel structure in which each pixel of a tomographic image is arranged as a voxel in a three-dimensional space. As a specific process in step S40, a value corresponding to the signal value of each pixel included in the ROI range of the tomographic image after clipping setting is added using the color map read from the color map reading means 20. , to create a voxel image included in the voxelized area calculated in step S30. Using the color map as shown in FIG. 5, for each voxel of the voxel image corresponding to each pixel of the tomographic image, each color component of R, G, and B corresponding to the signal value of that pixel and the opacity α is given as a voxel value, a voxel image for generating a full-color volume rendering image can be created.

By using a color map such as that shown in FIG. 5, a portion having a predetermined signal value in a tomographic image can be expressed with an arbitrary color component and an arbitrary opacity. In a tomographic image, since signal values differ according to parts such as organs, specific organs can be artificially colored based on differences in signal values by defining a color map in a volume rendering image. or make certain organs transparent to reveal hidden organs behind them. If the opacity α of a specific part is set to the maximum value (α=255), only that part will be clearly displayed, and the parts located behind that part will not be rendered. Conversely, if the opacity α of a specific portion is set to the minimum value (α=0), that portion becomes invisible (transparent), and the portion located behind that portion is exposed. If the opacity α of a specific portion is set to a halftone (0<α<255), that portion is displayed semi-transparently, and the portion located behind it is displayed as a watermark. Moreover, as the color components R, G, and B, arbitrary values can be set so as to facilitate visual recognition of the target portion.

Therefore, a color map can be set for each site to be observed, and different volume rendering images can be obtained by using different color maps for the same group of tomographic images. That is, if the color map for lungs is used, a volume rendering image in which the lung portion is conspicuous is obtained, and if the color map for heart is used, a volume rendering image in which the heart portion is conspicuous is obtained. However, since the signal value distribution of each organ may overlap each other, it is often impossible to create a color map that renders only a specific organ. For example, in the color map for the heart, the ribs and sternum are also drawn at the same time, and the heart is hidden by these bones.

In this embodiment, it is possible to change the color map in use or load a different color map by operating via the instruction input I/F 4 . The volume rendering device then refers to the changed color map or the newly loaded color map to generate a rendered image to be presented as a volume rendered image. Since the volume rendering apparatus of this embodiment can generate a rendering image at high speed, it can update the volume rendering image in real time in response to changes in the color map.

In step S40, the voxel image generating means 60 performs voxelization of the plurality of tomographic images within the range of Xs<x<Xe, Ys<y<Ye, and Zs<z<Ze and calculated in step S30. Referencing the color map using the signal value of each pixel (x, y, z) included in the region, obtaining the values of the color components R, G, B and the opacity α corresponding to the signal value, Let that value be the value of the corresponding voxel. By performing this process on all pixels in a plurality of tomographic images, a three-dimensional voxel image corresponding to the plurality of tomographic images can be obtained. As a result, each voxel in the voxel image has four values of R, G, B, and α. Since each color component R, G, B, and α is recorded with 8 bits per pixel, each pixel is recorded with 32 bits in the voxel image.

In this embodiment, the opaque circumscribed area calculation process in step S30 and the voxel image creation process in step S40 are performed in parallel by a plurality of CPU cores. Here, the details of the calculation process of the opaque circumscribed area in step S30 will be described. FIG. 6 is a flow chart showing a case where the processing for calculating the opaque circumscribed area and the like in step S30 are performed in parallel. The opaque circumscribed area calculation means 50 divides the tomographic image group read in step S10 into N (an integer of N≧2), and performs parallel processing in each CPU core. Since the direction in which each tomographic image having (x, y) coordinates is stacked is the z-axis direction, it is divided into N pieces according to the value of the z-coordinate. Specifically, the coordinates of the leading tomographic image of each tomographic image group are set to Z1, and the coordinates of the trailing tomographic image are set to Z2, and the read tomographic images are assigned to the leading coordinate Z1 and the trailing coordinate Z2. For example, Z1=Zs and Z2=Zs+(Ze-Zs+1)/N-1 are set for the first CPU core. In each tomographic image group, tomographic images are arranged in order, with one end being the head and the other end being the tail.

Next, an opaque circumscribed region is calculated for the group of tomographic images to be processed (step S330). The processing in step S330 is performed in parallel by N CPU cores. Details of the processing in step S330 will be described later.

After the calculation of the opaque circumscribing region for the tomographic image group in step S330, the opaque circumscribing region is calculated for the entire tomographic image consisting of the N tomographic image groups (step S360). Details of the processing in step S360 will also be described later. In this embodiment, a plurality of divided tomographic image groups are processed in parallel by a plurality of CPU cores.

Details of the processing of step S330 described above will be described. FIG. 7 is a flow chart showing the process of calculating an opaque circumscribed region for a group of tomographic images in step S330. The processing of step S330 is performed on each of the divided N tomographic image groups. In step S330, the coordinates of the top tomographic image of each tomographic image group are set to Z1 and the coordinates of the last tomographic image are set to Z2, and the read tomographic images are assigned to the top coordinate Z1 and the bottom coordinate Z2. For example, Z1=Zs and Z2=Zs+(Ze-Zs+1)/N-1 are set for the first CPU core. In each tomographic image group, tomographic images are arranged in order, with one end being the head and the other end being the tail.

First, initial setting is performed (step S331). Specifically, z specifying a tomographic image, x and y as coordinate values of the tomographic image, variables Xmin(tid) and Xmax(tid) for determining the minimum and maximum values of x, respectively, Variables Ymin(tid) and Ymax(tid) for determining the minimum and maximum values of y, respectively, and variables Zmin(tid) and Zmax(tid) for determining the minimum and maximum values of z are initialized. give a value. (tid) is identification information for specifying a tomographic image group, and a value corresponding to each of the N tomographic image groups is given. This tid is also identification information that identifies the processing thread.

As x and y, when ROI clipping is set in step S10, the minimum value Xs and the minimum value Ys are given, respectively. In addition, as Xmin(tid), Xmax(tid), Ymin(tid), Ymax(tid), Zmin(tid), and Zmax(tid), the minimum values are set to The variable for determination is given the maximum value by the ROI clipping setting, and the variable for determining the maximum value is given the minimum value by the ROI clipping setting. Therefore, when ROI clipping is set, the initial values are set as shown in [Formula 1] below.

[Formula 1]
Xmin(tid) = Xe
Xmax(tid) = Xs
Y min (tid) = Ye
Ymax(tid)=Ys
Zmin(tid)=Ze
Zmax(tid)=Zs

Note that if the ROI clipping setting is omitted in step S10, [Xso , Xeo], [Yso, Yeo], and [Zso, Zeo] are set as initial values.

After the initial setting is performed, the color map shown in FIG. 5 is referenced by the signal value of the target pixel of the tomographic image to acquire the opacity α (step S332). Next, it is determined whether or not at least one of the acquired opacity α and mask data Mask (x, y, z) is "0" (step S333). If the mask setting in step S20 is omitted, it is not determined whether the mask data Mask (x, y, z) is "0", and only the opacity α is "0". It is determined whether or not there is

In step S333, if "No", that is, if both the opacity α and the mask data Mask (x, y, z) are not "0", the pixel is considered to be opaque, so the minimum value is determined. variables Xmin(tid), Ymin(tid), and Zmin(tid) for determining the maximum value, and variables Xmax(tid), Ymax(tid), and Zmax(tid) for determining the maximum value (step S334). Specifically, the value of the variable is changed by executing the processing according to [Formula 2] below.

[Formula 2]
If x<Xmin(tid) then Xmin(tid)=x
If x>Xmax(tid) then Xmax(tid)=x
If y<Ymin(tid) then Ymin(tid)=y
If y>Ymax(tid) then Ymax(tid)=y
If z<Zmin(tid) then Zmin(tid)=z
If z>Zmax(tid) then Zmax(tid)=z

By [Equation 2], each value of the minimum value or maximum value in each of the xyz axis directions defining the opaque circumscribed region in the target tomographic image group (tid) is the above-mentioned opacity α, the mask data Mask (x, y , z) are not "0".

In step S333, that is, when it is determined that the opacity α or the like is "0", or when it is determined not to be "0" and the variable is changed in step S334, the process for the target pixel is performed. Since this is the end, the target pixel is changed in the x-axis direction (step S335). Specifically, the x-coordinate value is incremented. As a result of changing the x-coordinate, if x≤Xe, the process returns to step S333 to process the next pixel. On the other hand, if x>Xe in step S335, the target pixel is changed in the y-axis direction (step S336). Specifically, the value of the y-coordinate is incremented. As a result of changing the x-coordinate, if y≤Ye, the process returns to step S333 to process the next pixel. When all the pixels (x, y) of the tomographic image specified by a certain z have been processed in this way, the target pixel is changed in the z-axis direction (step S337). That is, the target tomographic image is changed. If Z>z2 in step S337, it means that (Z2-Z1+1) tomographic images in one tomographic image group have been processed. finish.

By performing processing on all tomographic images and pixels in one tomographic image group, all tomographic images and pixels for which both the opacity α and the mask data Mask (x, y, z) are not "0" are processed in step S334. A process according to [Formula 2] in is executed. As a result, even if there is a pixel in which either the opacity α or the mask data Mask (x, y, z) is “0”, both the opacity α and the mask data Mask (x, y, z) are “0”. The outermost tomographic image that is not 0″, the position of the pixel can be acquired.

By performing the processing of step S330 shown in FIG. 6 from z=Z1 to z=Z2, an opaque circumscribed region in one divided tomographic image group is calculated. Then, in step S360, an opaque circumscribed region of the entire group of tomographic images obtained by collecting all the N tomographic image groups is calculated. Details of the processing of step S360 will be described. FIG. 8 is a flow chart showing the calculation processing of the opaque circumscribed region of the entire group of tomographic images in step S360. First, initialization is performed (step S361). Specifically, tid, which is identification information for specifying a group of tomographic images, a fixed value N indicating the number of processing threads, variables Xmin and Xmax for determining the minimum and maximum values of x, respectively, and the minimum value of y Initial values are given to variables Ymin and Ymax for determining the value and maximum value, respectively, and variables Zmin and Zmax for determining the minimum value and maximum value of z, respectively.

Next, using the minimum values Xmin(tid), Ymin(tid), Zmin(tid) and the maximum values Xmax(tid), Ymax(tid), Zmax(tid) calculated for each divided image group, the minimum value The variables Xmin, Ymin, and Zmin for determining , and the variables Xmax, Ymax, and Zmax for determining the maximum value are changed (step S362). Specifically, the value of the variable is changed by executing the processing according to [Formula 3] below.

[Formula 3]
If Xmin(tid)<Xmin then Xmin=Xmin(tid)
If Xmax(tid)>Xmax then Xmax=Xmax(tid)
If Ymin(tid)<Ymin then Ymin=Ymin(tid)
If Ymax(tid)>Ymax, then Ymax=Ymax(tid)
If Zmin(tid)<Zmin then Zmin=Zmin(tid)
If Zmax(tid)>Zmax then Zmax=Zmax(tid)

According to [Equation 3], when the minimum or maximum values in the xyz axis directions of the opaque circumscribed region in the target tomographic image group (tid) exceed the range defining the entire opaque circumscribed region, the tomographic The range defining the entire opaque circumscribing region is expanded and updated according to the minimum or maximum values in the xyz axis directions of the opaque circumscribing region in the image group (tid).

Next, the target tomographic image group is changed (step S363). Specifically, the value of tid is incremented. If tid is smaller than N as a result of incrementing the value of tid, it means that processing has not been completed for all tomographic image groups. Perform processing to change variables for determining values and maximum values. If tid is equal to or greater than N as a result of incrementing the value of tid, this means that processing has been completed for all tomographic image groups, and thus boundary processing of opaque circumscribed regions is performed (step S364). In step S364, as a measure against moire, a boundary plane (x=Xs' or x=Xe' or y=Ys' or y=Ye' or z=Zs' or z =Ze') is added to the edge of the circumscribing rectangular parallelepiped that is the outermost periphery of the opaque circumscribing area by one pixel width. Specifically, a process of expanding the calculated area by repeating step S362 by one pixel in each of the xyz axis directions is performed. More specifically, the value of the variable is changed by executing the processing according to [Formula 4] below.

[Formula 4]
When Xs'=Xmin-1 and Xs'<Xs, Xs'=Xs
If Xe′=Xmax+1 and Xe′>Xe, then Xe′=Xe
If Ys'=Ymin-1 and Ys'<Ys, then Ys'=Ys
If Ye'=Ymax+1 and Ye'>Ye, then Ye'=Ye
When Zs'=Zmin-1 and Zs'<Zs, Zs'=Zs
When Ze'=Zmax+1 and Ze'>Ze, Ze'=Ze

In step S364, the boundary pixels are moved outward by one pixel. The area expanded outward by one pixel in this manner is used as a voxelized area for which a voxel image is to be created. Since the outermost pixels of the voxelized area should have an opacity α of "0", the outermost pixels of the voxelized area, so-called pixels located on the boundary surface of the voxelized area, are all transparent pixels. . Therefore, the opaque circumscribed area calculation means 50 calculates the voxelized area so that the opacity corresponding to the outermost pixels in each axial direction is zero. However, regardless of whether the mask data Mask (x, y, z) is 0 or outside the range of the ROI, the color components (RGB values) of the pixels located on the boundary surface of the voxelized area are The RGB values defined by the color map are given as they are based on the signal values. As a result, the opacity is discontinuous inside and outside the voxelized area, but the continuity of the RGB values is maintained. It is possible to suppress moiré that occurs on the boundary surface of the area. A range defined by a rectangular parallelepiped of [Xs', Xe'], [Ys', Ye'], and [Zs', Ze'] obtained as a result of the processing in step S364 is defined as a voxelized region and specified as a processing target. . In step S364 of FIG. 8, by setting the area expanded by one pixel to the outside of the opaque circumscribed area as the voxelized area, the outermost periphery in each axial direction of the voxelized area when creating a voxel image to be described later. The opacity corresponding to each voxel (pixel) is set to 0. However, it is not limited to this, and the area defined by the opaque circumscribed area is used as it is as a voxelized area. By forcibly replacing the opacity of the voxel that becomes 0 with 0, the opacity corresponding to the outermost voxel may be set to 0.

As described above, the opaque circumscribed area is calculated by the process of step S30. Now, let us compare the ROI obtained by the ROI clipping process of step S10 and the opaque circumscribed area obtained by the process of step S30. FIG. 9 is a diagram comparing an ROI and an opaque circumscribed region. Of these, FIG. 9(a) shows the ROI, and FIG. 9(b) shows the opaque circumscribed area. In FIGS. 9A and 9B, elliptical solid lines indicate the distribution range of opaque pixels (α>0) in the space of multiple tomographic images. Here, opaque pixels are pixels whose opacity α is greater than 0 as a result of referring to the color map, and the distribution range of opaque pixels indicates the maximum range in which opaque pixels exist. . In other words, transparent pixels with an opacity α of 0 also exist within the distribution range of opaque pixels. Here, the distribution range of opaque pixels shown in FIGS. 9A and 9B is the same.

In FIG. 9A, solid-line rectangular parallelepipeds indicate ROIs, and broken-line rectangular parallelepipeds indicate opaque circumscribed regions. In FIG. 9B, solid-line rectangular parallelepipeds indicate opaque circumscribed regions. As described above, the ROI is performed by numerically specifying the respective coordinate ranges in the xyz axis directions in the ROI clipping setting in step S10. The opaque circumscribing area is calculated by calculating the opaque circumscribing area that circumscribes the distribution range of the opaque pixels in the calculation of the opaque circumscribing area in step S30. Normally, the ROI clipping setting in step S10 is set by the user so as to include the distribution range of opaque pixels, and in the calculation of the opaque circumscribing region in step S30, the minimum opaque circumscribing region including the distribution range of opaque pixels is Calculated. Therefore, as shown by solid-line rectangular parallelepipeds and broken-line rectangular parallelepipeds in FIG. 9A, the opaque circumscribing region is completely included in the ROI, and the opaque circumscribing region is smaller than the ROI.

Therefore, in step S40 described later, a voxel image is created in the voxelized area based on the opaque circumscribed area. However, since the ROI can be arbitrarily set in the ROI clipping setting in step S10, it is not necessarily set so as to include the entire distribution range of opaque pixels. In that case, a voxel image is created within the bounds of both the ROI and the opaque bounding region. Further, when a mask is defined in the mask setting process (S20), the outside of the range surrounded by an arbitrary three-dimensional shape such as an ellipsoid shown as the range of α>0 in FIG. 9 is masked by α = 0, but is defined independently of the distribution range of opaque pixels. That is, there are two types of ellipsoid-like distributions shown as the range of α>0 in FIG. is set, and the opaque circumscribing area becomes even smaller.

Now, the distribution of pixels inside the opaque circumscribing region will be described. FIG. 9(c) is obtained by adding doughnut-shaped mask data Mask (x, y, z) to FIG. 9(b). is the mask area where Mask(x, y, z)=0. In FIG. 9C, pixels where both the opacity α and the mask data Mask (x, y, z) are not “0” are shaded. z) is "0" is shown without shading. As shown in FIG. 9(c), even if the distribution of pixels where both the opacity α and the mask data Mask (x, y, z) are not “0” is doubled inside and outside, the opacity The opaque circumscribed area is set so as to include the outermost part of the distribution of pixels where both α and mask data Mask (x, y, z) are not "0". No opaque circumscribed area is set for the distribution of pixels in which both inner mask data Mask (x, y, z) are not "0".

In this way, even if volume rendering is performed by constructing a voxel image reduced to a range corresponding to the opaque circumscribed region, the resulting volume rendering image is similar to that obtained by constructing a voxel image for the entire conventional tomographic image. There is no difference from rendering. However, compared to the latter, the processing load for constructing the voxel image, the memory load for holding the voxel image, and the load for transferring the voxel image as a 3D texture image to the GPU side and executing it, which will be described later, are each suppressed, and the time is shorter than before. It is possible to obtain a volume rendering image with

Next, the details of the voxel image creation process in step S40, which enables the creation of a voxel image reduced to the range corresponding to the opaque circumscribed area at a higher speed, will be described. [Xs', Xe'], [Ys', Ye'], and [Zs', Ze'] obtained as a result of the process in step S364 are voxelized areas defined by rectangular parallelepipeds. Xe], [Ys, Ye], and [Zs, Ze]. FIG. 10 is a flow chart showing the process of performing the voxel image creation process in step S40 in parallel by each CPU core. The voxel image creating means 60 divides the group of Ze-Zs+1 tomographic images included in the voxelized area reduced by the opaque circumscribed area calculating means 50 into N pieces (an integer of N≧2), and each CPU core Do parallel processing. Specifically, the coordinates of the leading tomographic image of each tomographic image group are set to Z1, and the coordinates of the trailing tomographic image are set to Z2, and the read tomographic images are assigned to the leading coordinate Z1 and the trailing coordinate Z2. For example, Z1=Zs and Z2=Zs+(Ze-Zs+1)/N-1 are set for the first CPU core. In each tomographic image group, tomographic images are arranged in order, with one end being the head and the other end being the tail.

After the tomographic image group is divided into N pieces and each is assigned to each CPU core, each CPU core executes processing according to the flowchart shown in FIG. First, the outline configuration of the processing of the core portion from steps S520 to S550 in FIG. 10 will be described. In this process, for each tomographic image of the group of divided tomographic images, elemental images are basically created in order from the beginning Z1 to the end Z2 (step S530), and shading is applied to the created elemental images. (Step S540) is performed to update the element image. In the elemental image, each pixel value (signal value) of the pixels included in the range of [Xs, Xe] and [Ys, Ye] corresponding to the voxelized area of the tomographic image is used as the color value (R, G, B). An image converted to transparency (α) and reduced in number of pixels to (Xe−Xs+1)×(Ye−Ys+1) pixels. Department. That is, a voxel image is obtained by arranging two-dimensional elemental images so as to overlap each other at predetermined intervals. Therefore, the values of each voxel of the elemental image are four, R, G, B, and α, as described above. The elemental image of the coordinate z is created in one-to-one correspondence with each pixel in the range of [Xs, Xe] and [Ys, Ye] of the tomographic image of the coordinate z. In order to execute the shading process (step S540), it is necessary to refer to the element image at the coordinate z−1 and the element image at the coordinate z+1. . Therefore, for the previously created elemental image at the coordinate z−1, shadow addition processing is performed while referring to the previously created elemental image at the coordinate z−2 and the currently created elemental image at the coordinate z. to update the previously created element image at coordinates z−1. That is, the shading process (step S540) is always executed with a delay of one z coordinate. However, in the first two steps of z=Z1 and z=Z1+1, since there is no elemental image at the coordinate z−2 or the elemental image at the coordinate z−1 at the stage of executing the shading process (step S540), the shading The addition process (step S540) is skipped, the shadow addition process for z=Z1+1 is performed at the stage of z=Z1+2, and the shadow addition process for z=Z1 is performed at the final step S560.

In this way, if the processing from steps S520 to S550 is repeated in order from the beginning Z1 to the end Z2, all the element images from the beginning Z1 to the end Z2 are created (in the processing of FIG. , the element image of the end Z2 is not created in this process), and the element images of the beginning Z1 and the end Z2 are not subjected to shading. Therefore, the shading process for these two elemental images is executed in the final step S560. In order to execute the shading process on the element images at the beginning Z1 and the end Z2, it is necessary to refer to the element image at the coordinate Z1−1 and the element image at the coordinate Z2+1. created by other CPU cores that are Therefore, after completing the series of processes from steps S520 to S550, shadows are applied to the element images at the beginning Z1 and the end Z2 while referring to the element images created by other CPU cores arranged in the shared memory. Perform additional processing. This addition processing corresponds to step S560 in FIG.

Another problem arises here. When executing shading processing on the element image with the last Z2, the element image with coordinates Z2+1 created by another CPU core placed in the shared memory is referenced, but the element image is Since the element image corresponds to the element image at the beginning Z1 and is created first in order, the element image at the coordinate Z2+1 becomes There is no problem because there is a high possibility that it has already been created. However, when executing the shading process for the element image at the beginning Z1, the element image at the coordinate Z1-1 created by another CPU core is referred to, but the element image is the element image at the end Z2 in the other CPU core. Since the element image corresponds to the element image and is created last in order, the element image at the coordinate Z1-1 is created at the stage when the CPU core executes the shading process for the element image at the beginning Z1. may not be in time (more than half of the cases in the experiment were not in time). Therefore, each CPU core gives top priority to the creation of the element image of the end Z2 (step S510). After that, the processing from steps S520 to S550 is executed, but since the creation of the element image of the end Z2 has already been completed, the repetition range is shortened from the beginning Z1 to the end Z2-1 as shown in step S550. do. Then, in addition to the element images at the beginning Z1 and the end Z2, the element image at the coordinate Z2-1 is not subjected to the shading process. Therefore, as shown in step 560, only the shading process is additionally performed for the three element images of the head Z1, the tail Z2, and the coordinate Z2-1.

Based on the above, the processing contents will be described in detail according to the flowchart shown in FIG. First, an elemental image corresponding to the tomographic image that is the last z-coordinate is created (step S510). In step S510, processing for correcting the opacity of the created elemental image is also performed as an optional processing. Details of the processing in step S510 will be described later. After specifying the voxel values and correcting the opacity of the element image corresponding to the last tomographic image in step S510, the next tomographic image to be processed is shifted to the top z-coordinate (step S520). Specifically, the processing target is shifted to z=Z1 at the beginning. Z1 is an integer that satisfies Zs≦Z1<Z2≦Ze.

Next, an elemental image corresponding to the range of [Xs, Xe] and [Ys, Ye] of the z-coordinate tomographic image to be processed is created (step S530). Also in step S530, a process of correcting the opacity of the created elemental image is performed as a preferable but not essential process. The specific processing in step S530 is exactly the same as step S510, except that the element images to be processed are different. Therefore, the details of the processing in step S530 will also be described later.

After creating the elemental image in step S530, z is compared with (Z1+2), and if z≧(Z1+2), shading is added (step S540). Therefore, if z<(Z1+2), the shading process in step S540 is not performed. That is, the first element image (z=Z1) is not subjected to shading (the second element image (z=Z1+1) from the top is subjected to shading at the stage of z=Z1+2). The shading process in step S540 is performed on the elemental image one slice before. Details of the processing in step S540 will also be described later.

After finishing the shading process in step S540, z is incremented (z←z+1). If z<Z2, the process returns to step S530 and repeats the processes up to step S550. On the other hand, if z≧Z2, shadow Addition processing is performed (step S560). The specific processing in step S560 is exactly the same as step S540, except that the element images to be processed are different. Therefore, the details of the processing in step S560 will also be described later.

As a result, in the process shown in the flowchart of FIG. 10, first, the last (z=Z2) element image is created (step S510), and then the last (z=Z1) element image is sequentially created. Elemental images up to one (z=Z2-1) before (z=Z2-1) and shading processing are repeatedly performed (steps S530 to S550). However, the top element image and the last two element images are not shaded. Finally, shadow addition processing is performed on the three element images for which shadow addition processing has not been performed (step S560). That is, in steps S510 to S550, element images are created for all tomographic images from the top tomographic image to the last tomographic image in the group of tomographic images. Shading is added by correcting the color component values of the pixels of the elemental image corresponding to the image and the elemental images corresponding to the last two tomographic images. More specifically, in steps S510 to S550, after the voxel image creating means 60 creates element images for the last tomographic image of the tomographic image group, the voxel image creating unit 60 creates the tomographic image from the first tomographic image to the last tomographic image of the tomographic image group. In step S560, the shading means 62 generates an element image corresponding to the top tomographic image and the element image corresponding to the last tomographic image. and the element images corresponding to the tomographic images at the end are shaded by correcting the color component values.

In this embodiment, a group of divided tomographic images are processed in parallel by a plurality of CPU cores and finally integrated to create a voxel image. Further, shading processing for a given elemental image must be performed using adjacent elemental images. Therefore, shadow addition processing for the leading elemental image requires the trailing elemental image of another divided tomographic image group, and shadowing processing for the trailing elemental image requires the leading elemental image of the other divided tomographic image group. becomes. However, the leading element image corresponding to z=Zs and the trailing elemental image corresponding to z=Ze are not subjected to brightness value calculation in the shading process, and are given a brightness value of 1.0. Do not correct all the original values of the color components of the image.

Similarly, (Ye−Ys+1)×2 pixels in the Y-axis direction corresponding to x=Xs and x=Xe, which are both ends in the X-axis direction of each elemental image, and y=Ys, which is both ends in the Y-axis direction and (Xe−Xs+1)×2 pixels in the X-axis direction corresponding to y=Ye also have no neighboring pixels, so the luminance value is not calculated and 1.0 is given as the luminance value. , all the values of the color components of the elemental image are not corrected as they are. As described above, the opacity of pixels whose brightness values are not calculated is all set to 0, and the color component values are defined by the color map based on the signal values of the corresponding tomographic image pixels. By setting the value, the generation of moire on the boundary surface of the voxel image is suppressed. The voxel image creation range in step S40 is the range of the voxelized region based on the opaque circumscribed region calculated in step S30. Therefore, the number of voxels to be created is greatly reduced compared to the case where all the pixels of the tomographic image group are created as voxels, and the amount of data can be reduced. As a result, the number of voxels to be processed is reduced in the rendering process, which will be described later, so that the processing load can be reduced.

Details of the processing of steps S510 and S530 described above will be described. FIG. 11 is a flow chart showing the process of creating an elemental image in steps S510 and S530. First, an elemental image corresponding to the range of [Xs, Xe] and [Ys, Ye] of the target z-coordinate tomographic image is created (step S511). Specifically, using a color map as shown in FIG. 5, for each voxel of the voxel image corresponding to each pixel in the voxelized area of the tomographic image, R, An element image is created by giving each color component of G and B and opacity α as a voxel value. For example, when z=zu, an elemental image with voxel values V(x, y, zu, c) (c=0 (R), 1, (G), 2 (B), 3 (α)) is created. be done.

Next, the α value, which is opacity, is corrected (step S512). Specifically, an opacity correction table Sα(x, y, z) (0≦Sα(x, y, z)≦1; 0≦x≦Sx-1, 0≦ y≤Sy-1, 0≤z≤Sz-1), the opacity α is corrected by performing the processing according to [Equation 5] below. The opacity correction table Sα(x, y, z) is a multidimensional transfer function in which a correction magnification for opacity is defined in advance in a three-dimensional space (x, y, z). The same function sets the correction magnification to the maximum in the central area where the internal organs exist, decreases the correction magnification toward the periphery by a predetermined function, Outside the body: bed, headrest/fixing jig) is defined so that the correction magnification is 0. By executing the processing according to [Equation 5] using this, it is possible to obtain a volume rendering image in which the outside of the bone region becomes transparent and only the visceral region is exposed.

[Formula 5]
V' (x, y, zu, 3) = V (x, y, zu, 3) Sα (x + Xs, y + Ys, zu + Zs)

In [Equation 5], V(x, y, zu, 3) indicates the opacity α before correction since it is the case of c=3. V'(x, y, zu, 3) indicates the opacity α after correction by the opacity correction table Sα(x, y, zu). V(x, y, zu, 3) and V'(x, y, zu, 3) are narrowly defined in the ranges Xs≤x≤Xe, Ys≤y≤Ye, Zs≤zu≤Xe. On the other hand, Sα(x, y, z) is defined in the same range as the original tomographic image, 0 ≤ x ≤ Sx-1, 0 ≤ y ≤ Sy-1, 0 ≤ z ≤ Sz-1 Therefore, V(0,0,0,3) corresponds to Sα(Xs, Ys, Zs).

The voxel data created in this embodiment corresponds to each pixel of a plurality of two-dimensional tomographic images, and is a voxel structure composed of each voxel determined with reference to a predefined color map. . This voxel data has the total number of tomographic images as M, the total number of elemental images as M'=Ze-Zs+1, and the predetermined number of divisions as N (both M' and N are integers, N>1 and N<M'≤ M is satisfied), M′ element images created corresponding to M′ tomographic images corresponding to the range [Zs, Ze] are divided into N element image groups, The N elemental image groups are processed in parallel and created by voxel image creating means for creating a voxel image in which the N elemental image groups are integrated.

Details of the shadow adding process in step S540 described above will be described. In step S540, the shading means 62 performs shading processing. This is a process of adding a shadow corresponding to a predetermined lighting environment to a voxel image. Specifically, the voxel opacity α , the color components of voxels in the range of Xs<x<Xe, Ys<y<Ye and Zs<z<Ze are updated. For this purpose, a certain parallel light source and ambient light components are set, and shadows to be added are calculated. Specifically, first, a light source vector (Lx, Ly, Lz) and an ambient light component Ab are set as unit vectors indicating the directions of the parallel light sources. The numerical values of the light source vector (Lx, Ly, Lz) and the ambient light component Ab can be appropriately set according to the situation of the image to be expressed. .57735, 0.57735) and set Ab=0.2. Lx, Ly, and Lz are real numbers with absolute values less than 1 (including negative values), and Ab is a real number of 0 or more and 1 or less. Then, the gradient vector (Gx, Gy, Gz) at the voxel (x, y, zu-1) of the voxel image is calculated by performing the processing according to [Equation 6] below.

[Formula 6]
Gx=(V(x+1,y,zu-1,3)-V(x-1,y,zu-1,3)).(Rxy/Rz)
Gy = (V (x, y+1, zu-1, 3) - V (x, y-1, zu-1, 3)) (Rxy/Rz)
Gz = (V (x, y, zu, 3) - Vα (x, y, zu - 2, 3))
G=(Gx ² +Gy ² +Gz ² ) ^1/2

In [Formula 6], V (x, y, zu-1, 3) is the opacity α in the voxel (x, y, zu-1) of the voxel image, and Rxy is the resolution of the tomographic image (reciprocal of the pixel interval is the number of pixels per mm, and the x-axis direction and the y-axis direction are the same), Rz is the slice resolution of the tomographic image (reciprocal of the slice interval of the tomographic image, the number of slices of the tomographic image per mm), G is the gradient vector ( Gx, Gy, Gz). Rxy/Rz is defined as the z-axis direction scaling factor. You may multiply the reciprocal of the magnification. As shown in [Equation 6], the shading means 62 calculates the difference value of the opacity of two neighboring voxels in each of the x-axis direction, the y-axis direction, and the z-axis direction of the voxel of the voxel image as the x-axis of the gradient vector. A unit vector calculated as a direction component, a y-axis direction component, and a z-axis direction component, and corrected to the x-axis direction/y-axis direction component or the z-axis direction component based on a predefined z-axis direction scaling factor, It is calculated as a gradient vector at the voxel. Further, when G≧1, the processing according to [Equation 7] below is executed to calculate the luminance value S(x, y, zu−1) including only the diffuse reflection component. S(x, y, zu-1) is a real number of 0 or more and 1 or less.

[Formula 7]
S (x, y, zu-1) = (1-Ab) (|Gx · Lx + Gy · Ly + Gz · Lz |) / G + Ab

On the other hand, when the value of G calculated in [Equation 6] is G<1, it is a uniform region (in the background air, lung fields, etc.) where there is no opacity gradient, or it is the outermost region. Therefore, the brightness value S(x, y, zu-1) is set to 0 without adding shading. In [Equation 7], |Gx·Lx+Gy·Ly+Gz·Lz| represents an absolute value, and the light source vector in the opposite direction (-Lx, -Ly, -Lz) is also virtually irradiated at the same time, and the subject is in the opposite direction. Diffuse reflected light with the same brightness can be obtained even at an angle.

Then, using the calculated brightness value S (x, y, zu-1), the processing according to the following [Equation 8] is performed, and each voxel value V (x, y, zu-1) of the voxel image , c) to obtain corrected voxel values V′(x, y, zu−1, c).

[Formula 8]
V'(x,y,zu-1,c)=S(x,y,zu-1)V(x,y,zu-1,c)

In [Equation 8], c=0, 1, 2, corresponding to the R, G, and B color components, respectively. That is, in [Formula 8], only the color component is corrected. In this way, an elemental image V'(x, y, zu-1, c) that has undergone shading processing for the tomogram of z=zu-1 is obtained. As a result, the shading means 62 performs the processing according to [Equation 6] to [Equation 8], so that for each pixel of each element image (zu-1) created, an element adjacent to the element image By referring to the opacity α of the pixels located near the xyz axis direction including the image (zu, zu-2), the color component value of the pixel of the element image (zu-1) is corrected, and shading is added. Is going. In this embodiment, as a preferred example, adjacent pixels are used as pixels for which the opacity α is referred to, but the pixels are not limited to adjacent pixels, and other nearby pixels may be referred to.

At step S540 shown in FIG. 10, the above-described shading process is performed for each z-coordinate of (Z1+1)≤z≤(Z2-2). Also, in step S560 shown in FIG. 10, the above-described shading process is performed on the three elemental images of z=Z1, z=Z2−1, and z=Z2.

A voxel image is obtained as described above. After completing the voxel image creation process in step S40, the voxel image is registered in the 3D texture image (step S50). Specifically, the rendering means 70 registers a voxel image as a 3D texture image used when executing the 3D texture mapping method. As mentioned above, voxel images are created only in voxelized regions. Therefore, by registering a voxel image as an OpenGL 3D texture, the voxel creation load and the transfer load to the texture memory are reduced.

FIG. 12 is a diagram showing correspondence in 3D texture mapping. FIG. 12(a) shows a 3D texture map in the texture coordinate system, and FIG. 12(b) shows the stacked rectangles defined in the world coordinate system. The 4 vertices of each quadrangle forming the stacked quadrangle are associated with the 4 voxels of the 3D texture image, and each quadrangle is the basis of the texture map when it is not rotated as shown in FIG. The range of [Xs, Xe] [Ys, Ye] [Zs, Ze] for the tomographic image corresponds one-to-one with the elemental image cut out, and the four vertices of each quadrangle correspond to the four corners of the elemental image. Pixels are associated with voxels. A 3D texture image registered in OpenGL physically has (Xe-Xs+1)×(Ye-Ys+1)×(Ze-Zs+1) voxels, and the sizes in the xyz-axis directions are different. is defined uniformly in the range of [0, 1] in the xyz-axis direction, and the size is normalized to 1. Therefore, in the scaling described later, in addition to the z-axis direction scaling processing due to the difference in resolution between the xy-axis direction and the z-axis direction of the tomographic image, 1/(Xe−Xs+1) and 1/(Ye− Ys+1) and 1/(Ze-Zs+1).

Specifically, as shown in FIG. 12(b), 3D texture mapping is performed on the defined stacked quadrilaterals. In this case, the 3D world coordinates of the four vertices of each quadrangle in the world coordinate system (corresponding to the viewpoint coordinate system) are associated with the 3D texture coordinates of the 3D texture image in the texture coordinate system (corresponding to the voxel coordinate system). As shown in FIG. 12, a quadrangle composed of four vertices (-1, -1, z), (-1, 1, z), (1, -1, z), (1, 1, z) , the world coordinates (-1, -1, z) are associated with the texture coordinates (0, 0, r).

As will be described later, when there is an angle change instruction, a converted 3D texture image is generated in which the texture map in FIG. 12A is three-dimensionally rotated in the texture coordinate system. Correspondence breaks down. Voxel images are registered as 3D texture images via OpenGL (graphic API).

Subsequently, rendering processing is performed (step S60). FIG. 13 is a diagram showing the configuration of the rendering means 70. As shown in FIG. As shown in FIG. 13, the rendering means 70 has a 3D texture registration means 71, a coordinate conversion means 72, a stacked rectangle setting means 73, and a pixel value calculation means 74. As described above, in the case of the 3D texture mapping method, the rendering means 70 is realized mainly by executing a program in the GPU 7 with the assistance of the CPU 1. However, the coordinate transformation means 72 and the pixel The value calculation means 74 is preferably executed using the GPU and frame memory installed in the video card of the general purpose computer.

In the processing of the 3D texture mapping method, first, a projection screen is set. FIG. 14 is an explanatory diagram showing an example of projection screen setting by the rendering means 70 of this embodiment. The rendering means 70 sets the screen size (number of vertical and horizontal pixels, vertical and horizontal aspect ratio) of the rendered image. The rendering means 70 sets either parallel projection (normal external rendering) or perspective projection (endoscopic mode). A "rendered image" is synonymous with a "volume rendered image", and the "volume rendered image" may be referred to as a "rendered image". When the perspective projection is set, the rendering means 70 sets perspective projection parameters. The perspective projection parameters are composed of, for example, the viewing angle (focal length) of the camera and the viewpoint position (viewpoint and gaze point. The former is the eye position and the latter is the position on the object being viewed, both on the z-axis. In general, the gaze point is fixed at the origin of the world coordinate system), the clipping position (the distance on the z-axis from the viewpoint, two points near and far), and the like. Note that the clipping position may be limited to near objects. As shown in FIG. 14, in the case of parallel projection, all lines of sight from the viewpoint are assumed to be parallel to the z-axis, and the viewpoint is assumed to be infinitely far away in the leftward direction. All stacked rectangles defined along the z-axis are 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 viewpoint enters the inside of the stacked quadrilateral. far rectangle).

FIG. 15 is a flow chart showing details of rendering processing by the 3D texture mapping method. First, the 3D texture registration means 71 in FIG. 13 registers the voxel image as a 3D texture image in OpenGL. Next, the coordinate transformation means 72 performs predetermined coordinate transformation on the 3D texture image to generate a post-transformation 3D texture image (S61 to S63).

Specifically, the coordinate conversion means 72 includes a 4×4 matrix defining rotation in a predetermined three-dimensional coordinate system, a viewing angle, a viewpoint position, a clipping position, an offset value in each of the xyz axis directions, and an expansion in each of the xyz axis directions. Alternatively, predetermined coordinate transformation parameters including reduction ratio and z-axis direction scaling ratio are acquired, and the 3D texture image is subjected to predetermined coordinate transformation using the acquired parameters to generate a post-transformation 3D texture image.

Describing individual steps, the coordinate transformation means 72 performs scaling and z-axis direction scaling processing on the 3D texture image (step S61). Scaling is performed with the same magnification Scale in the xyz-axis direction. Enlargement or reduction is performed only at the specified magnification Scz (=Rxy/Rz). Furthermore, as described above, since the 3D texture images to be registered are normalized to the same size in the xyz axis directions, Add scaling processing with magnification. Then, the coordinate transformation means 72 performs rotation processing on the 3D texture image using a 4×4 matrix that defines rotation in a predetermined three-dimensional coordinate system (step S62).

Then, the coordinate conversion means 72 performs offset processing on the 3D texture image using the offset values (Xoff, Yoff, Zoff) in the directions of the xyz axes (step S63). A converted 3D texture image is generated by these processes (S61 to S63). Next, the stacked quadrangle setting unit 73 sets a stacked quadrangle composed of a plurality of quadrilaterals (step S64). Specifically, a stacked quadrangle is set by arranging quadrangles on the xy coordinate plane of the three-dimensional space in the z-axis direction. When setting 3D texture mapping, the stacked quadrangle setting means 73 sets a stacked quadrangle in which the same number of quadrangles as the number of xy plane images are arranged in the z-axis direction.

The number of quadrilaterals in the stacked quadrilateral can be set arbitrarily, but if the number of quadrilaterals and the number of element images (Ze-Zs+1) do not match, interpolation processing works in the z-axis direction, and coordinate values are changed by coordinate conversion. Rounding errors accumulate, and striped or grid-like moiré occurs. Then, as shown in FIG. 12, by defining four voxel coordinates in the corresponding post-transformation 3D texture image for the four vertices of each quadrangle, correspondence with the post-transformation 3D texture image is performed ( step S65). The voxel coordinates to be defined are based on the texture coordinate system of OpenGL (real values ranging from 0 to 1 with respect to the three axial directions of r, s, and t), not the coordinate system of the voxel image.

Specifically, as shown in FIG. 12(b), 3D texture mapping is performed on the defined stacked quadrilaterals. In this case, the three-dimensional world coordinates of the four vertices of each quadrangle in the world coordinate system are associated with the three-dimensional texture coordinates of the 3D texture image in the texture coordinate system. As shown in FIG. 12, a quadrangle composed of four vertices (-1, -1, z), (-1, 1, z), (1, -1, z), (1, 1, z) , the world coordinates (-1, -1, z) are associated with the texture coordinates (0, 0, r).

Next, the pixel value calculation means 74 performs scan conversion (filling in the inside of the rectangle) in order from the farthest rectangle on the line of sight parallel to the z-axis direction from a predetermined viewpoint to the nearest rectangle to the viewpoint (step S66). At this time, the inside of the rectangle is painted based on the color components (R, G, B) of the voxels of the 3D texture image after conversion attached to each rectangle, and when the same pixel is painted over, the voxels of the 3D texture image after conversion are painted. Also do alpha blending based on opacity. In this way, the color component (R, G, B) composed of (R, G, B) of the voxel of the converted 3D texture image corresponding to the xyz coordinates of the stacked rectangle is obtained from the viewpoint based on the opacity of the voxel. Color components composed of (R, G, B) obtained by alpha blending in the order of the farthest rectangle are given as pixel values of the volume rendering image. A volume rendering image is thus obtained.

More specifically, all pixel values (RGB values) of the rendering image are initialized in advance with the background color (for example, R=G=B=0). The pixel value calculation means 74 performs scan conversion on the quadrangle furthest from the viewpoint while referring to the attached post-conversion 3D texture image. Parallel projection or perspective projection is performed on the rendered image defined in FIG. Let (i, j) be the coordinates of a pixel to be filled in the rendered image, and u be the interval between pixels in the rendered image (where u is the interval in world coordinates corresponding to one pixel in the rendered image. For example, u= 0.002) for scan conversion. At each scan-converted world coordinate (-1+iu, -1+ju, -1) (i, j are coordinate values of the rendered image, for example, integers 0≤i, j≤511), the corresponding texture coordinates (2iu , 2ju, 0), refer to the converted 3D texture image based on the calculated texture coordinates, obtain the color component values (RGB values) and opacity, and calculate the corresponding rendered image coordinates (i, j ) and the acquired opacity, alpha blending is performed, and the pixel value at the coordinates (i, j) of the corresponding rendering image is updated.

In this way, when the alpha blending process is completed for all corresponding pixels of the rendered image on which the single rectangle furthest from the viewpoint is projected, the z coordinate is changed to v (the world coordinate where the rectangle is located in the z-axis direction is increased by v. For example, v=2/Sz), and the same process is repeated for the next closest quadrangle (-1+iu, -1+ju, -1+v) to the viewpoint direction to update the rendered image. Rendered image after scan conversion for all quadrilaterals and alpha blending processing for world coordinates (-1+iu, -1+ju, -1+kv) (k is the quadrilateral number and an integer of 0≤k≤Sz-1) is completed. The alpha blending process is performed as a process according to [Formula 9] below.

[Formula 9]
R′=R・α+(1−α)・Rb
G′=G・α+(1−α)・Gb
B′=B・α+(1−α)・Bb

In [Formula 9], R', G', and B' are the pixel values (RGB values) of the rendering image updated on the projection plane. R, G, B are the color components (RGB values) in the voxel (2iu, 2ju, 2kv) of the 3D texture image after conversion corresponding to the world coordinates (-1+iu, -1+ju, -1+kv) of the rectangle. corresponds to the value of α is also the α value in the transformed 3D texture image voxel (2iu, 2ju, 2kv) corresponding to the quadrangle world coordinates (−1+iu, −1+ju, −1+kv), and controls the overlapping ratio. Rb, Gb, and Bb are pixel values (RGB values) already recorded at the coordinates (i, j) of the rendered image corresponding to the world coordinates (-1+iu, -1+ju, -1+kv). The values of R', G', and B' calculated based on the above [Equation 9] correspond to the quadrangle located on the side opposite to the viewpoint and the quadrangle farthest from the viewpoint. , the background color (for example, Rb=Gb=Bb=0) is given. A volume rendering image is recognized by outputting the rendering image obtained as described above by the rendering image output means 80 .

The boundary surface of the voxelized area, that is, the value of the opacity α of the outermost voxel of the voxelized area is set to 0, and the RGB values of the color components are the natural values defined by the color map based on the signal value of the pixel. By giving RGB values, it is possible to suppress the moiré that occurs on the boundary surface while suppressing the processing load. In the present embodiment, the voxel image itself is also created in the same range as the voxelized area based on the opaque circumscribed area, so moire can be suppressed while further suppressing the processing load.

As described above, in the volume rendering apparatus 100 according to the present embodiment, the opaque circumscribed area calculation means 50 and the voxel image creation means 60 are implemented by the CPU 1, which is a multi-core CPU, executing programs. This program divides the read tomographic image into a plurality of (for example, N) tomographic image groups, assigns each tomographic image group to each CPU thread, and causes each CPU thread to calculate a part of the opaque circumscribed area. A process of creating an elemental image group corresponding to each tomographic image group is performed in parallel. A case in which the volume rendering apparatus 100 is realized by a personal computer, which is a general-purpose computer, will be described. FIG. 16 is a diagram showing parallel processing when the volume rendering device 100 is implemented by a computer. As shown in FIG. 16, the opaque circumscribed area calculation means 50 and the voxel image creation means 60 are implemented by an application program. Processing for the image group is divided into four processing threads #n1 to #n4. Then, the job scheduler of the multitasking operating system in which the opaque circumscribed area calculation means 50 and the voxel image creation means 60 (application program) operate allocates them to CPU threads #1 to CPU threads #8 of the plurality of CPU cores.

FIG. 17 is a diagram showing a software flow of parallel processing when the volume rendering apparatus 100 is implemented by a computer. As shown in FIG. 17, the application program that realizes the opaque circumscribed area calculation means 50 and the voxel image creation means 60 divides the processing for each of the four elemental image groups into four processing threads #n1 to #n4. After that, each invokes a parallel processing function. Then, the job scheduler of the multitasking operating system assigns four threads #n1 to #n4 to CPU thread #3, CPU thread #5, CPU thread #7, and CPU thread #2 each having a plurality of CPU cores. assign. In this way, the operating system allocates the processing threads to be processed to the CPU threads of the CPU cores that are free, thereby enabling parallel processing. In the example of FIG. 17, the application program waits for end messages from each CPU thread, and when end messages are obtained from all CPU threads, a voxel image is created by integrating the created element image groups. Then, after rendering processing is performed collectively, the rendered image is output to the display unit 6 and displayed.

As a result of actual measurement on a Windows 10/64bits personal computer equipped with a 4-core CPU (Intel CORE-i7), the configuration in Fig. 17 shows that when executing a single thread (512 x 512 x 370 full resolution image voxel by CPU A processing speed of about 1.5 times faster than the generation time: about 0.9 seconds was obtained (same generation time: about 0.6 seconds). However, some CPU threads were waiting for processing. Therefore, when the number of divisions of the elemental image group was increased to N = 8 and parallel processing was performed with 8 processing threads, the processing performance was almost double that of single thread execution (same generation time: about 0.45 sec ). Furthermore, as a result of experiments with N > 8, the processing performance hardly changed and reached a ceiling at 2 times. rice field. On the other hand, the rendering processing S60 after voxel image creation was executed by the GPU 7, and the generation time of the rendering image with full resolution of 512×512 was about 0.015 seconds.

<3. Decimation of 1/K>
In the volume rendering device, the voxel image creating means 60 creates a voxel image with a density equivalent to that of a plurality of tomographic images. good too. In this case, the voxel image generating means 60 arranges the voxel images three-dimensionally in association with every K pixels in each of the three axial directions of the two-dimensional xy-axis direction of the plurality of tomographic images and the z-axis direction orthogonal to the tomographic images. For each voxel obtained, the signal value of that voxel is replaced with the color component and opacity α to create a voxel image in which the color component and opacity are defined as voxel values. For example, when K=2, each of the xyz-axis directions can be thinned out by 1/2, so the number of voxels to be processed becomes 1/8, and high-speed processing can be performed. As a result, even when the user continuously switches the color map interactively via the instruction input I/F 4, the volume rendering image with coarse image quality is sequentially displayed, and the color map is updated. A volume rendering image can be displayed while tracking.

<4. Gradation Reduction of Tomographic Image>
In the above embodiment, a tomographic image in which each pixel records a 16-bit signal value is used as it is to create a voxel image. A voxel image may be created later. By converting the tomographic image into 8-bit gradation, the load in processing each pixel can be reduced, contributing to high-speed generation of volume rendering images.

In this case, tomographic image tone conversion means is further provided, and tone conversion is performed each time the tomographic image reading means 10 reads a single tomographic image. Specifically, a series of DICOM format tomographic images Do (x, y, z) (-32768 ≤ Do (x, y, z) ≤ 32767, 0 ≤ x ≤ Sx-1, 0 ≤ y ≤ Sy-1 , 0 ≤ z ≤ Sz-1; resolution: Rxy, Rz), the Sz/2th intermediate tomographic image Do (x, y, z/2) is first read, and all pixels in the same tomographic image Let Dmin be the minimum value of and Dmax be the maximum value of . For Dmin and Dmax, a method of determining the minimum and maximum values from all tomographic images may be used, but if a method of determining only from intermediate tomographic images is used, all slices of the large-capacity 16-bit original tomographic image may be used. Only the gradation-compressed tomographic image D8(x, y, z), which is a gradation-reduced image, can be constructed directly on the memory without reading it and storing it on the memory. Subsequently, a process according to [Equation 10] below is executed to calculate the lower limit value Lmin and the upper limit value Lmax.

[Formula 10]
Lower limit Lmin = (Dmax-Dmin) γ + Dmin
Upper limit value Lmax=(Dmax−Dmin)・(1−γ)+Dmin

In [Equation 10], γ is the contrast adjustment width of the gradation-compressed image, and the closer to 0, the higher the contrast but the lower the luminance. γ is a real value less than 0.3, but is usually set to γ=0.1. Since the luminance contrast of the rendered image can be adjusted by various settings such as a color map, γ may be fixed. Then, the range between the lower limit value Lmin and the upper limit value Lmax is compressed in 256 steps to obtain a gradation-compressed tomographic image D8(x, y, z). Specifically, the gradation compression tomographic image D8 (x, y, z) (0≤D8 (x, y, z)≤255, 0≤x≤Sx-1, 0≤y≤Sy-1, 0≤ When z≦Sz−1; resolution: Rxy, Rz), the process according to the following [Equation 11] is executed to calculate a gradation compression tomographic image D8(x, y, z).

[Formula 11]
D8 (x, y, z) = (Do (x, y, z) - Lmin) 255 / (Lmax - Lmin)
If D8(x,y,z)>255: D8(x,y,z)=255 If D8(x,y,z)<0: D8(x,y,z)=0

As shown in the second formula of [Equation 11], when the signal value exceeds the upper limit value Lmax, it is set to 255, and when it is below the lower limit value Lmin, it is set to 0. Convert. When creating a voxel image using a gradation-compressed tomographic image, a color map for 8-bit signal values is prepared as a color map for a gradation-reduced image instead of the color map in FIG. . Specifically, a color map in which R, G, B, and α are recorded in association with signal values 0 to 255 is used.

Although the preferred embodiments of the present disclosure have been described above, the present disclosure is not limited to the above embodiments, and various modifications are possible. For example, in the above embodiment, a multi-core CPU having a plurality of CPU cores and capable of parallel processing is used as the CPU 1, but a single-core CPU in which one CPU core only processes one processing thread is used. may be used.

In the above-described embodiment, in order to prevent moiré, the voxel image is created by setting the opacity of the outermost voxel of the voxelized area to 0. However, if the moiré is not noticeable, the outermost voxel There is no need to create a voxel image with 0 opacity.

1 CPU (Central Processing Unit)
2 RAM (random access memory)
3... Storage device 4... Instruction input I/F
5 Data input/output I/F
6... Display unit 7... GPU
8 frame memory 10 tomographic image reading means 20 color map reading means 30 ROI clipping setting means 40 mask setting means 50 opaque circumscribed area calculation means 60 Voxel image creation means 61 Opacity correction means 62 Shadow addition means 70 Rendering means 71 3D texture registration means 72 Coordinate conversion means 73 Laminated rectangle setting means 74. ..Pixel value calculation means 80..Rendered image output means 100..Volume rendering device

Claims (15)

Rendering with reference to a color map defined by associating color component values and opacity with signal values, based on multiple 2D tomographic images of a given object taken at given intervals. A volume rendering device for generating an image, comprising:
calculating an opaque circumscribed region, which is a region that circumscribes in each three-dimensional axial direction the region in which pixels having non-zero opacity obtained by referring to the color map exist among the pixels of each of the tomographic images; opaque circumscribed area calculation means;
Among the tomographic images, a region corresponding to the opaque circumscribed region is defined as a voxelized region, and the signal value of each pixel in the voxelized region is converted by referring to the color map, and color components and non-color components are obtained as pixel values. voxel image creation means for creating a voxel image in which two-dimensional element images with defined transparency are superimposed at predetermined intervals;
rendering means for generating a rendered image using the voxel image;
volume rendering device. The voxel image creating means is
2. The volume rendering apparatus according to claim 1, wherein the voxel image is created such that, when the opacity corresponding to the outermost pixel in the voxelized area is not 0, the opacity corresponding to the pixel is replaced with 0. Xso is the coordinate corresponding to the left end in the x-axis direction of each elemental image, Xeo is the coordinate corresponding to the right end in the x-axis direction, Yso is the coordinate corresponding to the upper end in the y-axis direction, and Yso is the coordinate corresponding to the lower end in the y-axis direction. is Yeo, Zso is the z-coordinate at which the first elemental image is arranged, and Zeo is the z-coordinate at which the M-th elemental image, which is the total number of tomographic images, is arranged. When six coordinates Xs, Xe, Ys, Ye, Zs, and Ze are defined that specify a voxel image area for which a rendering image is to be generated, satisfying Ys<Ye≦Yeo and Zso≦Zs<Ze≦Zeo ,
3. The volume according to claim 1, wherein said voxel image creation means creates a voxel image in a range that satisfies x≦Xs, x≧Xe, y≦Ys, y≧Ye, z≦Zs, or z≧Ze. rendering device. When mask data having a value of 0 (not displayed) or 1 (displayed) is defined corresponding to the tomographic image,
4. The volume rendering apparatus according to any one of claims 1 to 3, wherein said opaque circumscribed region calculating means refers to said mask data to specify pixels whose opacity is not zero. The opaque circumscribing region calculating means divides each of the tomographic images into N tomographic image groups with a predetermined number of divisions set to N (N>1), and divides each of the N tomographic image groups into the opaque circumscribing region. 5. The volume rendering apparatus according to any one of claims 1 to 4, wherein processes for creating a part of a region are executed in parallel. When mask data with a value of 0 (not displayed) or 1 (displayed) is defined corresponding to the element image,
The voxel image creation means refers to the mask data corresponding to the element image when creating the element image while giving opacity and color components to each pixel, and the value of the mask data corresponding to the pixel is 6. The volume rendering apparatus according to any one of claims 1 to 5, wherein the opacity of the pixel of the elemental image is set to 0 when the pixel is 0. 7. The voxel image creation means according to any one of claims 1 to 6, comprising opacity correction means for correcting opacity using an opacity correction table determined by three-dimensional coordinate values of the voxel image. Volume rendering device. The voxel image creating means divides each of the tomographic images corresponding to the voxelized region into N tomographic image groups with a predetermined number of divisions being N (N>1). 8. The volume rendering apparatus according to any one of claims 1 to 7, wherein the process of creating the element images is executed in parallel. The voxel image creating means refers to the opacity of each pixel of each created elemental image, which is positioned at least in the xyz-axis direction, including the pixel, to determine the color component of the pixel of the elemental image. 9. A volume rendering apparatus according to claim 8, comprising shading adding means for correcting values. In the tomographic image group, when one end arranged in order is the head and the other end is the end,
After the element images are created for all the tomographic images from the top tomographic image to the last tomographic image of the group of tomographic images, the shading means adds the element image corresponding to the top tomographic image and the last tomographic image. 10. The volume rendering apparatus according to claim 9, wherein the color component values of the pixels of the elemental image corresponding to the tomographic image are corrected. In the tomographic image group, when one end arranged in order is the head and the other end is the end,
The voxel image creating means, after creating the element image for the last tomographic image of the tomographic image group, creates each tomographic image from the first tomographic image to the last tomographic image of the tomographic image group. In contrast, creating the element image,
The shading means adds color component values to an element image corresponding to the head tomographic image, an element image corresponding to the tomographic image immediately before the tail, and an element image corresponding to the tail tomographic image. 10. A volume rendering apparatus according to claim 9, wherein the correction is performed. The rendering means is
3D texture registration means for generating a 3D texture image composed of the voxel image;
coordinate transformation means for performing a predetermined coordinate transformation on the 3D texture image to generate a post-transformation 3D texture image;
Lamination in which quadrilaterals on the xy coordinate plane of the three-dimensional space are arranged in the z-axis direction, and each of the three-dimensional coordinates of the four vertices of each of the quadrilaterals is associated with each of the predetermined four three-dimensional coordinates of the post-transformation 3D texture image. a laminated quadrangle setting means for setting a quadrangle;
A color component composed of (R, G, B) of a voxel of the converted 3D texture image corresponding to the three-dimensional coordinates on the stacked quadrangle on the line of sight parallel to the z-axis direction from a predetermined viewpoint is Pixel value calculation means for giving color components obtained by alpha blending in order of rectangles farther from the viewpoint based on opacity as pixel values composed of (R, G, B) of the rendered image. When,
12. A volume rendering apparatus according to any one of claims 1 to 11, comprising: The coordinate transformation means is
A 4×4 matrix defining rotation in a predetermined three-dimensional coordinate system, a predetermined Get the parameters of the coordinate transformation,
13. The volume rendering device according to claim 12, wherein the 3D texture image is subjected to the predetermined coordinate transformation using the acquired parameters to generate the post-transformation 3D texture image. 14. A volume rendering apparatus according to claim 13, wherein said coordinate transformation means and said pixel value calculation means are executed using a GPU and frame memory mounted on a video card of a general-purpose computer. A program for causing a computer to function as the volume rendering device according to any one of claims 1 to 14.

Images