JP2020038592A - Volume rendering apparatus - Google Patents

Abstract

To provide a volume rendering apparatus which can improve speed of generating a rendering image.SOLUTION: A volume rendering apparatus according to the present invention has: voxel image generating means 30 for converting, based on a plurality of two-dimensional tomographic images obtained by capturing a predetermined object at a predetermined interval, with reference to a color map defined in association with values of color components and opacity with respect to a signal value, signal values of respective pixels of respective tomographic images to generate a voxel image having overlapped two-dimensional element images defining the color components and opacity as pixel values; and rendering means 70 for generating a rendering image by using the voxel image. Defining that a total number of the topographic images is M and that a predetermined division number is set as N (N>1 and N<M), the voxel image generating means 30 divides M topographic images into N topographic image groups and executes parallelly processing for generating element images for the respective N topographic image groups.SELECTED DRAWING: Figure 3

Description

  The present disclosure relates to a volume rendering technique for three-dimensionally visualizing a plurality of two-dimensional tomographic images.

  2. Description of the Related Art Conventionally, a plurality of medical images (modalities), such as CT, MRI, and PET, which are continuously imaged at predetermined slice intervals and stored in a medical information system such as a PACS (Picture Archiving and Communication Systems) in a DICOM format or the like. The three-dimensional visualization of organs and the like has been performed based on the tomographic images.

In the field of 3D computer graphics, a rendering image is generated in limited colors such as 256 colors, and a color map (color look-up table) implemented in hardware is changed interactively to change the color of the rendering image in real time. A technique is used (see Patent Document 1).

Further, 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) independent of a color map has been proposed (see Patent Document 2).

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

  However, in the technique described in Patent Document 1, only the RGB values can be changed, and control of an α value expressing opacity, which is important in volume rendering, is not disclosed. Further, in the technique described in Patent Document 2, the time required for processing can be reduced by only about 10% as compared with a case where a signal value is applied to a color map and re-rendering is performed.

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

The present disclosure includes a plurality of means for solving the above-described problems.
Based on a plurality of two-dimensional tomographic images photographed at predetermined intervals for a predetermined object, rendering is performed by referring to a color map defined by associating color component values and opacity with signal values. A volume rendering device for generating an image,
Referring to the color map, a signal value of each pixel of each of the tomographic images is converted, and a voxel image in which a two-dimensional element image in which a color component and an opacity are defined as pixel values is superposed is created. Means for creating voxel images,
Rendering means for generating a rendering image using the voxel image,
The voxel image creating means divides each of the tomographic images into N tomographic image groups by setting a predetermined number of divisions to N (N> 1), and divides the elemental image into N tomographic image groups. Execute the process to be created in parallel.

Also, in the present disclosure,
Based on a plurality of two-dimensional tomographic images photographed at predetermined intervals for a predetermined object, rendering is performed by referring to a color map defined by associating color component values and opacity with signal values. A volume rendering device for generating an image,
A voxel image in which a signal value of each pixel of each tomographic image is converted with reference to the color map, and a two-dimensional element image in which color components and opacity are defined as pixel values is arranged in the z-axis direction. Means for creating a voxel image,
Rendering means for generating a rendering image using the voxel image,
The voxel image creating means divides each of the tomographic images into N tomographic image groups by setting a predetermined number of divisions to N (N> 1), and divides the elemental image into N tomographic image groups. Provided is a program for causing a computer to function as a volume rendering device that executes processing to be created in parallel.

Also, in the present disclosure,
A voxel structure corresponding to each pixel of a plurality of two-dimensional tomographic images taken at predetermined intervals of a predetermined target object and configured by voxels determined with reference to a predefined color map. Voxel data,
Assuming that a predetermined number of divisions is N (N> 1), an element image created corresponding to each of the tomographic images is divided into N element image groups,
The N element image groups provide voxel data that is processed in parallel by a voxel image creating unit that creates a voxel image integrating the N element image groups.

  According to the present disclosure, it is possible to perform a process of generating a rendering image at high speed by applying a color map to a plurality of tomographic images.

FIG. 4 is a diagram illustrating 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. 2 is a functional block diagram illustrating a configuration of the volume rendering device according to the embodiment. 5 is a flowchart illustrating a processing operation of the volume rendering device according to the embodiment. FIG. 3 is a diagram illustrating a color map used in the embodiment. 5 is a flowchart showing a case where voxel image creation processing and the like in step S20 shown in FIG. 4 are performed in parallel. 7 is a flowchart illustrating a process of creating an element image in steps S510 and S530 illustrated in FIG. FIG. 4 is a diagram illustrating association in 3D texture mapping. FIG. 3 is a diagram illustrating a configuration of a rendering unit when performing processing of a 3D texture mapping method. FIG. 4 is an explanatory diagram illustrating an example of a projection screen setting by a rendering unit 70 of the embodiment. It is a flowchart which shows the detail of the rendering process by a 3D texture mapping system. FIG. 7A is a diagram illustrating a display example of a rendering image in a state where clipping processing is performed on both the color value and opacity, and (b) mask processing is performed on both the color value and opacity. FIG. 7A is a diagram illustrating a display example of a rendered image in a state where (a) clipping setting and only opacity are clipped, and (b) mask setting and only opacity are masked. FIG. 11 is a diagram illustrating parallel processing when the volume rendering device 100 is implemented by a computer. FIG. 9 is a diagram illustrating a flow of software for parallel processing when the volume rendering apparatus 100 is implemented by a computer.

Hereinafter, preferred embodiments of the present disclosure will be described in detail with reference to the drawings.
The present disclosure is to generate a rendering image by referring to a color map defined by associating a color component value and an opacity with a signal value based on a plurality of tomographic images. FIG. 1 is a diagram showing a plurality of tomographic images and a rendering image. 1A shows a plurality of tomographic images, FIG. 1B shows one tomographic image, and FIG. 1C shows a rendering image. The tomographic images shown in FIGS. 1A and 1B are of the DICOM format. FIG. 1A shows a group of 370 tomographic images, and FIG. 1B shows an enlarged one of them. In a volume rendering apparatus according to an embodiment to be described later, rendering processing is performed on a plurality of tomographic images as shown in FIG. 1A, and one rendering image as shown in FIG. Generate. The rendering image shown in FIG. 1C has 512 × 512 pixels, and one pixel is expressed 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 the present embodiment can be realized by a general-purpose computer. As shown in FIG. 2, a CPU (Central Processing Unit) 1 and a RAM (Random Access Memory) 2 which is a main memory of the computer. A large-capacity storage device 3 such as a hard disk, a solid state drive (SSD), or a flash memory for storing programs and data to be executed by the CPU 1, and an instruction input I / F (interface) 4 such as a keyboard and a mouse. , A data input / output I / F (interface) 5 for performing 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 specialized 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 For example, they are connected to each other via a bus. Since the calculation result by the GPU 7 is written into the frame memory 8, the GPU 7 and the frame memory 8 are often mounted on a video card provided with an interface to the display unit 6 and mounted on a general-purpose computer via a bus.

  In the present embodiment, the CPU 1 is a multi-core CPU, has a plurality of CPU cores, and can perform parallel processing. Although only one RAM 2 is shown in the example of FIG. 2, each CPU core of the CPU 1 is configured to access one RAM 2. Note that a plurality of CPUs 1 may be provided. Further, the multi-core CPU may be a CPU in which each CPU core logically receives a plurality of execution threads.

  FIG. 3 is a functional block diagram illustrating a configuration of the volume rendering device according to the present embodiment. 3, 10 is a tomographic image reading means, 15 is a color map reading means, 20 is an ROI clipping setting means, 25 is a mask setting means, 30 is a voxel image creating means, 40 is an opacity correcting means, and 50 is a shading adding means. , 70 are rendering means, and 80 is a rendered image output means.

  The tomographic image reading unit 10 is a unit that reads a plurality of tomographic images from a PACS collected and accumulated by a medical image diagnostic device such as CT, MRI, or PET via a storage medium or an input unit. The tomographic image is obtained by photographing the 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 15 is means for reading a color map created in advance by a color map creating means (not shown) from a storage medium or an input unit. The ROI clipping setting unit 20 is a unit that sets a voxel image creation target among a plurality of read tomographic images as an ROI that is a region of interest. The mask setting means 25 is a means for setting a mask of an area not to be displayed in the plurality of read tomographic images.

  The voxel image creating means 30 refers to a color map defined by associating the color component values and the opacity with the signal values, and associates each of the three-dimensionally arranged pixels with a plurality of tomographic images. This is a means for creating a voxel image by giving a color component value and opacity to the voxel. The voxel image creating means 30 further includes an opacity correcting means 40 and a shading adding means 50. The opacity correction means 40 is means for correcting the opacity α of each pixel of the voxel image using a predetermined table. Based on the opacity α of the voxel in the vicinity of the voxel for the color component composed of (R, G, B) of each voxel of the color component voxel image, the shading adding means 50 Is calculated by using the calculated gradient vector, and a shadow value is calculated and replaced with a value obtained by multiplying each value of the color components (R, G, B) by the shadow value.

  The tomographic image reading means 10 and the color map reading means 15 are realized mainly in the data input / output I / F 5 with the assistance of the CPU 1. The ROI clipping setting means 20 and the mask setting means 25 are realized mainly in the instruction input I / F 4 with the assistance of the CPU 1. The voxel image creating unit 30, the opacity correcting unit 40, and the shading adding unit 50 are realized by the CPU 1 executing a program stored in the storage device 3. The rendering means 70 is realized by executing a program mainly in the GPU 7 with the assistance of the CPU 1. The rendering image output unit 80 is realized mainly in the frame memory 8 and the display unit 6 with the assistance of the CPU 1.

  Each component shown in FIG. 3 is actually realized by installing a dedicated program in hardware such as a computer and its peripheral devices as shown in FIG. That is, the computer executes the contents of each means according to the dedicated program. Note that in this specification, a computer means a device having an arithmetic processing unit such as a CPU and a GPU and capable of data processing, and includes not only a general-purpose computer such as a personal computer but also a tablet or the like equipped with a GPU. It also includes mobile terminals and computers embedded in various devices.

<2. Processing Operation>
Next, the processing operation of the volume rendering device shown in FIGS. 2 and 3 will be described. FIG. 4 is a flowchart illustrating the processing operation of the volume rendering device according to the present embodiment. First, the tomographic image reading means 10 reads a plurality of tomographic images. The plurality of tomographic images have an aspect as shown in FIG. In the present embodiment, the tomographic image group Do (x, y, z) in DICOM format (−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 laminating Sz tomographic images having the x-coordinate size (number of pixels) of Sx and the y-coordinate size (number of pixels) of Sy. In the following description, the number Sz of tomographic images is also expressed as M. The horizontal direction of one tomographic image in FIG. 1B corresponds to the x-axis, the vertical direction corresponds to the y-axis, and the stacking direction when a plurality of tomographic images illustrated in FIG. Corresponding.

  After reading a plurality of tomographic images in the DICOM format, the ROI clipping setting means 20 sets ROI clipping (step S10). ROI is an abbreviation of Region of Interest, and in this specification, means a region of interest in a three-dimensional space. (Academically, ROI indicates a two-dimensional region, and in the case of a three-dimensional region, VOI (Volume of Interest). However, in the medical imaging industry, ROI is frequently used even in a three-dimensional area). Here, the range in which the rendering image and the voxel image are created is shown. By setting the ROI, it is possible to generate a volume rendering image in which a subject is cut three-dimensionally at an arbitrary position, and is used to depict an internal organ or an organ hidden in the body surface or bone. In the processing after step S20, the ROI is set in step S10 so that the voxel image creation target is limited to a predetermined range.

  The coordinates corresponding to the left end in the x-axis direction of each read tomographic image correspond to Xso, the coordinates corresponding to the right end in the x-axis direction correspond to Xeo, the coordinates corresponding to the upper end in the y-axis direction correspond to Yso, and the lower end in the y-axis direction. If the coordinates are Yeo, the coordinates corresponding to the tomographic image of the first slice (head: first) are Zso, and the coordinates corresponding to the tomographic image of the last slice (end: Mth) are Zeo, Xso ≦ Xs <Xe A range defined by a rectangular parallelepiped of [Xs, Xe], [Ys, Ye], [Zs, Ze] that satisfies ≦ Xeo, Yso ≦ Ys <Ye ≦ Yeo, Zso ≦ Zs <Ze ≦ Zeo is defined and processed as ROI. Identify as target. That is, six coordinates Xs, Xe, Ys, Ye, Zs, and Ze that specify a voxel image area for which a rendering image is to be generated are defined. As a result of the ROI clipping setting in step S10, the number of pixels to be processed is changed from (Xeo-Xso + 1) × (Yeo-Yso + 1) × (Zeo-Zso + 1) pixels to (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 20, in addition to the effect that a volume rendering image exposing the organ to be observed is obtained, the number of pixels to be processed is reduced and the processing load is suppressed to reduce the responsiveness. There is a double advantage that it can be improved, so it is often set, but it is not essential processing and can be omitted. In the stage of setting the ROI clipping in step S10, a specific clipping process is executed in each process after step S20 only by defining a range to which the ROI clipping is applied. In each of the processes, the processing load is reduced by performing the calculation limited to the ROI range of the voxel image, and a clipped (cut) volume rendering image is generated in the rendering process in step S70.

  After the ROI clipping setting is performed by the ROI clipping setting means 20, the mask setting means 25 performs mask setting (step S15). The mask setting is a setting as to whether or not to display each pixel of the tomographic image. The method of setting the mask is not particularly limited. In the present embodiment, the group of tomographic images read in step S10 is displayed on the screen as an MPR image (three axial coronal sagittal screens), and the method described in this specification is used. Alternatively, a volume rendering image created by a known method is displayed, an image that is most easily indicated is selected from the two-dimensionally projected images, and a mask is displayed while indicating a two-dimensional area of the selected image. Then, an area not to be displayed and an area to be displayed without masking are specified. With this mask setting, mask data in a three-dimensional space corresponding to the tomographic image group is obtained. The mask data Mask (x, y, z) is one in which a value of “0” or “1” is set corresponding to each pixel of the tomographic image group Do (x, y, z). . The process of creating such a three-dimensional mask is close to an essential process, but can be omitted. As a method of delineating the organ to be observed, there is the ROI clipping described above. However, since the operation is simpler and the processing load is reduced as compared with the three-dimensional mask, the ROI clipping is frequently used together. A specific area to be displayed is set by the ROI clipping setting means 20 and the mask setting means 25. If both areas are set, the area set as the area to be displayed is the effective area. When either one is set, as a result, the area set as the area to be displayed becomes the effective area. If neither the ROI clipping setting means 20 nor the mask setting means 25 is set, the entire area of the read tomographic image group is the effective area.

  Next, the voxel image creating means 30 performs a voxel image creating process (step S20). The voxel image is a voxel structure in which each pixel of the tomographic image is arranged as a voxel in a three-dimensional space. As a specific process in step S20, a value corresponding to the signal value of each pixel included in the ROI range of the tomographic image after clipping is set using the color map read from the color map reading unit 15. , Create a voxel image.

  FIG. 5 is a diagram showing a color map used in the present embodiment. In FIG. 5, the signal value is a signal value recorded in each pixel of the tomographic image, R, G, and B are red, green, and blue color components, respectively, and α is the degree of display of the voxel. Opacity to be specified. FIG. 5 is a color map in a case where the signal value of the tomographic image is recorded in 16 bits and takes a value of -32768 to 32767. In the case of a CT image, the water component is often normalized to a range from -2048 to 2048, with the water component as 0. As shown in FIG. 5, in the color map, the R, G, and B color components and the opacity α are associated with the signal values of each pixel of the tomographic image. By using such a color map, for each voxel of the voxel image corresponding to each pixel of the tomographic image, the R, G, and B color components and the opacity α corresponding to the signal value of the pixel are represented by the voxel value. Thus, a voxel image for generating a full-color volume rendering image can be created.

  By using a color map as shown in FIG. 5, a portion having a predetermined signal value in a tomographic image can be represented by an arbitrary color component and an arbitrary opacity. In tomographic images, signal values differ depending on the organs, organs, and other parts. Therefore, by defining a color map in a volume rendering image, a specific organ is artificially color-coded based on the difference in signal values. Or make certain organs transparent, exposing the organs hidden behind them. If the opacity α of a specific part is set to the maximum value (α = 255), only that part is clearly displayed, and a part located behind the part is not drawn. Conversely, if the opacity α of a specific part is set to the minimum value (α = 0), the part becomes invisible (transparent), and the part located behind the part is exposed. If the opacity α of a specific part is set to a half tone (0 <α <255), the part is displayed translucent, and the part located in the back is displayed by watermark synthesis. Further, as the color components R, G, and B, arbitrary values can be set to make the target portion easily visible. Therefore, a color map can be set for each region to be observed, and different volume rendering images can be obtained for the same tomographic image group by using different color maps. That is, if a color map for the lung is used, a volume rendering image in which the lung portion is conspicuous is obtained, and if a color map for the 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 not possible to create a color map depicting only a specific organ. For example, in the color map for the heart, the ribs and the sternum are simultaneously drawn, and the heart is hidden by these bones. Therefore, it is necessary to set ROI clipping and virtually cut the ribs.

  In the present embodiment, it is possible to change the color map in use or read a different color map by operating via the instruction input I / F 4. Then, the volume rendering device refers to the changed color map or the newly read color map, and generates a rendering image to be presented as a volume rendering image. Since the volume rendering device of the present embodiment can generate a rendering image at high speed, the volume rendering image can be updated in real time in response to a color map change.

  In step S20, the voxel image creating unit 30 outputs the signal at each pixel (x, y, z) corresponding to the range of Xs <x <Xe, Ys <y <Ye, and Zs <z <Ze of the plurality of tomographic images. The values of the color components R, G, and B corresponding to the signal values and the value of the opacity α are obtained by referring to the color map using the values, and the values are set as the values of the corresponding voxels. For values of voxels (x, y, z) in the range of x ≦ Xs or Xe ≦ x or y ≦ Ys or Ye ≦ y or z ≦ Zs or Ze ≦ z, R = G = B = α = 0 is set. By performing this processing for all the 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, the value of each voxel of the voxel image becomes R, G, B, and α. Each of the color components R, G, B, and α is recorded with 8 bits per pixel, so that in a voxel image, it is recorded with 32 bits per pixel.

  In the present embodiment, the voxel image creation processing in step S20 is performed by a plurality of CPU cores in parallel. FIG. 6 is a flowchart illustrating a case where the voxel image creation processing and the like in step S20 are performed in parallel. The voxel image creation unit 30 divides the tomographic image group read in step S10 into N (N ≧ 2 integers), and performs parallel processing in each CPU core. Since the direction in which the tomographic images having the (x, y) coordinates are stacked is the z-axis direction, the tomographic images are divided into N according to the value of the z coordinate. Specifically, the coordinates of the leading tomographic image of each tomographic image group are Z1, the coordinates of the trailing tomographic image are 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 of the tomographic image groups, the tomographic images are arranged in order, and one end is set at the top and the other end is set at the end.

  When the tomographic image group is divided into N and each is assigned to each CPU core, each CPU core executes a process according to the flowchart shown in FIG. First, the outline configuration of the processing of the core part from step S520 to S550 in FIG. 6 will be described. In this process, for each tomographic image of the divided tomographic image group, an element image is basically created in the order from the top Z1 to the end Z2 (step S530), and the created element image is shaded. (Step S540) is performed to update the element image. The element image is an image obtained by converting each pixel value (signal value) of the tomographic image into a color value (R, G, B) and an opacity (α). The z-coordinate is specified in a three-dimensional voxel image. And is part of a three-dimensional image. That is, a voxel image is obtained by arranging the two-dimensional element images in an overlapping manner. Therefore, the value of each voxel of the element image is R, G, B, and α as described above. The element image at the coordinate z is created in one-to-one correspondence with each pixel of the tomographic image at the coordinate z. To execute the shading adding process (step S540) on the element image at the coordinate z, Since it is necessary to refer to the element image at the coordinate z-1 and the element image at the coordinate z + 1, there is a first problem that the processing cannot be executed as it is. At the stage of executing the shading adding process (step S540), since there is no element image of the coordinate z + 1, the element image of the coordinate z-1 created immediately before is compared with the coordinate z− created before the element image. The shadow adding process is performed with reference to the element image of the coordinate No. 2 and the element image of the coordinate z created this time, and the element image of the coordinate z-1 created immediately before is updated. That is, the shadow addition processing (step S540) is always executed with the z coordinate being delayed by one.

  By repeating the processing from step S520 to step S550 in this order from the top Z1 to the end Z2, all the element images from the top Z1 to the end Z2 are created (in the configuration of FIG. 6, one before the end Z2). Therefore, in this processing, the element image of the end Z2 is not created), and the shading addition processing is not performed on the element images of the head Z1 and the end Z2. In order to execute the shading adding process on the element images at the head 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. Other CPU cores. Therefore, after a series of processes from step S520 to S550 is completed, the element images at the head Z1 and the end Z2 are shaded with reference to the element images created by other CPU cores arranged in the shared memory. Execute additional processing. This additional processing corresponds to step S560 in FIG.

  Here, another problem arises. When performing the shading adding process on the element image at the end Z2, the element image at the coordinate Z2 + 1 created by another CPU core arranged in the shared memory is referred to. Since the element image corresponding to the element image at the head Z1 is created first in order, the element image at the coordinates Z2 + 1 is obtained when the CPU core executes the shading adding process on the element image at the end Z2. There is no problem because there is a high possibility that it has already been created. However, when the shading processing is performed on the element image at the head Z1, the element image at the coordinates Z1-1 created by another CPU core is referred to. Since the element image corresponds to the element image and becomes the last element image in order, the CPU core executes the shading adding process on the element image at the head Z1 and creates the element image at coordinates Z1-1. May not be in time (more than half of the experiments were not in time). Therefore, each CPU core gives the highest priority to the creation of the last Z2 element image (step S510). Thereafter, the processes from steps S520 to S550 are executed. However, since the creation of the element image at 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. I do. Then, in addition to the element images at the head Z1 and the end Z2, the shading addition processing is not executed for the element image at the coordinate Z2-1. Therefore, as shown in step 560, only the shadow addition processing is additionally performed on the three element images of the head Z1, the end Z2, and the coordinate Z2-1.

  Based on the above, the processing content will be described in detail according to the flowchart shown in FIG. First, an element image corresponding to the tomographic image serving as the end z coordinate is created (step S510). In step S510, a process of correcting the opacity of the created element image is also performed as a non-essential process. Details of the processing in step S510 will be described later. After the specification of the voxel value of the element image corresponding to the last tomographic image and the correction of the opacity are completed in step S510, the processing target is shifted to the tomographic image having the leading z coordinate (step S520). Specifically, the processing target is shifted to the first z = Z1. Z1 is an integer satisfying Zs ≦ Z1 <Z2 ≦ Ze.

  Next, an element image corresponding to the tomographic image of the z coordinate to be processed is created (step S530). Also in step S530, a process of correcting the opacity of the created element image is performed as a non-essential process. The specific processing in step S530 is exactly the same as step S510 except that the element image to be processed is different. Therefore, the details of the process in step S530 will be described later.

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

  After the shading addition processing in step S540 is completed, z is incremented (z ← z + 1). If z <Z2, the process returns to step S530, and the processing up to step S550 is repeated. On the other hand, if z ≧ Z2, a shadow is applied to the element image at the beginning (z = Z1), the element image at the end (z = Z2), and the element image immediately before (z = Z2-1) at the end. An additional process is performed (step S560). The specific processing in step S560 is exactly the same as step S540 except that the processing target element image is different. Therefore, the details of the process in step S560 will also be described later.

  After all, in the processing shown in the flowchart of FIG. 6, first, the last (z = Z2) element image is created (step S510), and then the last (z = Z1) element image is sequentially processed. The generation of the element image up to one before (z = Z2-1) and the shading adding process are repeatedly performed (steps S530 to S550). However, shading addition processing is not performed on the first element image and the last two element images. Then, finally, the shadow addition processing is performed on the three element images for which the 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 first tomographic image to the last tomographic image in the tomographic image group. Shading is added by correcting the color component values of the pixels of the element image corresponding to the image and the element image corresponding to the last two tomographic images. More specifically, in step S510 to step S550, the voxel image creation unit 30 creates an element image for the tomographic image at the end of the tomographic image group, and then immediately before the end of the tomographic image at the beginning of the tomographic image group. An element image is created for each tomographic image up to the tomographic image of step S560. In step S560, the shading adding unit 50 sets the element image corresponding to the leading tomographic image and the element image corresponding to the last tomographic image immediately before the end. In addition, the value of the color component is corrected for the element image corresponding to the last tomographic image, and shading is added.

  In the present embodiment, a plurality of divided tomographic image groups are processed in parallel by a plurality of CPU cores, and finally integrated to create a voxel image. Further, it is necessary to perform the shading adding process on a certain element image by using an adjacent element image. Therefore, the shading process for the first element image requires the last element image of another group of tomographic images, and the shading process for the last element image requires the first element image of the other group of tomographic images. Becomes However, for the first element image corresponding to z = Zs and the last element image corresponding to z = Ze, the luminance value is not calculated in the shading adding process, and 0 is given as the luminance value, and Correct all the color component values to zero. In the present embodiment, z = Z1 is described as the head and z = Z2 is described as the tail. However, z = Z1 may be the tail and z = Z2 may be the head. Eventually, the head and tail may be set so that processing can be performed from one to the other without changing the stack configuration of all the read tomographic images.

  The details of the processing in steps S510 and S530 will be described. FIG. 7 is a flowchart showing the element image creation processing in steps S510 and S530. First, an element image corresponding to the tomographic image of the target z coordinate is created (step S511). Specifically, 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 of R, G, B corresponding to the signal value of that pixel An element image is created by giving the components and the opacity α as voxel values. For example, when z = zu, an element image having voxel values V (x, y, zu, c) (c = 0 (R), 1, (G), 2 (B), 3 (α)) is created. Is done.

  Next, the opacity α value is corrected (step S512). Specifically, an opacity correction table Sα (x, y, z) (0 ≦ Sα (x, y, z) ≦ 1; 0 ≦ x ≦ Sx−1, 0 ≦) determined by the three-dimensional coordinate values of the voxel image By using y ≦ Sy−1, 0 ≦ z ≦ Sz−1), the opacity α is corrected by executing processing according to the following [Equation 1]. 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 (x, y, z) three-dimensional space. This function sets the correction magnification in the central region where the internal organs are present to the maximum, attenuates the correction magnification toward the periphery by a predetermined function, and outside the bone region (body surface: fat, muscle, skin layer, (Outside the body: bed, headrest / fixing jig) is defined so that the correction magnification becomes zero. When the processing according to [Equation 1] is executed using this, a volume rendering image in which the outside of the bone region is transparent and only the internal organ region is exposed can be obtained.

[Formula 1]
V ′ (x, y, zu, 3) = V (x, y, zu, 3) · Sα (x, y, zu)

  In [Equation 1], V (x, y, zu, 3) indicates the opacity α before correction because c = 3. V ′ (x, y, zu, 3) indicates the opacity α after correction by the opacity correction table Sα (x, y, zu).

  The voxel data created in the present embodiment is a voxel structure corresponding to each pixel of a plurality of two-dimensional tomographic images and composed of voxels defined with reference to a predefined color map. . This voxel data corresponds to M tomographic images, where M is the total number of tomographic images and N is a predetermined number of divisions (M and N are both integers and satisfy N> 1 and N <M). The created M element images are divided into N element image groups, and the N element image groups are parallelized by a voxel image creating unit that creates a voxel image integrating the N element image groups. Created by processing.

  The details of the shading adding process in step S540 will be described. In step S540, the shadow adding means 50 performs a shadow adding process. This is a process of adding a shadow corresponding to a predetermined lighting environment to a voxel image. Specifically, the opacity α of a voxel in the range of Xs ≦ x ≦ Xe, Ys ≦ y ≦ Ye, and Zs ≦ z ≦ Ze , 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 an ambient light component are set, and a shadow to be added is calculated. Specifically, first, a light source vector (Lx, Ly, Lz) is set as a unit vector indicating the direction of the parallel light source, and an ambient light component Ab is set. 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. For example, (Lx, Ly, Lz) = (0.577735, 0) .57735, 0.57775), and Ab = 0.2. Lx, Ly, and Lz are real values (including negative values) whose absolute values are less than 1, and Ab is a real value that is 0 or more and 1 or less. Then, the gradient vector (Gx, Gy, Gz) in the voxel (x, y, zu-1) of the voxel image is calculated by executing a process according to the following [Equation 2].

[Equation 2]
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 [Equation 2], 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). 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), and G is the gradient vector (the number of pixels per mm in the x-axis direction and the y-axis direction are the same). Gx, Gy, Gz). Rxy / Rz is defined as the scaling factor in the z-axis direction. In [Equation 2], instead of multiplying the formula for calculating Gx and Gy by Rxy / Rz, the formula for calculating Gz is Rz / Rxy, that is, the scaling factor in the z-axis direction. The reciprocal of the magnification may be multiplied. As shown in [Equation 2], the shading adding unit 50 calculates the difference value of the opacity of the voxel in the x-axis direction, the y-axis direction, and the z-axis direction of the voxel image adjacent to each other in the voxel image by 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 for 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 in the voxel. Further, when G ≧ 1, the processing according to the following [Equation 3] 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 value of 0 or more and 1 or less.

[Equation 3]
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 2] is G <1, it is determined that the area is a uniform area where there is no opacity gradient (such as the air in the background or the lung field) or the outermost area. Therefore, it is assumed that no shading is added and the luminance value S (x, y, zu-1) = 0. In [Equation 3], | Gx · Lx + Gy · Ly + Gz · Lz | indicates an absolute value, and a light source vector (−Lx, −Ly, −Lz) in the opposite direction is virtually simultaneously illuminated, so that the object is directed in the opposite direction. Even at an angle, diffuse reflection light of the same luminance can be obtained.

  Then, using the calculated luminance value S (x, y, zu-1), a process according to the following [Equation 4] is executed, and each voxel value V (x, y, zu-1) of the voxel image is obtained. , C) to obtain each corrected voxel value V ′ (x, y, zu−1, c).

[Equation 4]
V ′ (x, y, zu−1, c) = S (x, y, zu−1) · V (x, y, zu−1, c)

In [Equation 4], c = 0, 1, and 2, which correspond to the R, G, and B color components, respectively. That is, in [Equation 4], only the color components are corrected. In this way, an element image V '(x, y, zu-1, c) obtained by performing the shading adding process on the tomographic slice of z = zu-1 is obtained. Eventually, the shading adding means 50 executes the processing according to [Equation 2] to [Equation 4] to obtain, for each pixel of each created element image (zu-1), an element adjacent to the element image. Referring to the opacity α of a pixel located in the vicinity of the xyz-axis direction including the image (zu, zu-2), the value of the color component of the pixel of the element image (zu-1) is corrected, and shading is added. Is going. In the present embodiment, as a preferred example, an adjacent pixel is used as a pixel that refers to the opacity α. However, the present invention is not limited to the adjacent pixel, and another neighboring pixel may be referred to.

  In step S540 shown in FIG. 6, the above-described shading addition processing is performed for each z coordinate of (Z1 + 1) ≦ z ≦ (Z2-2). Further, in step S560 shown in FIG. 6, the above-described shading processing is performed on three element images of z = Z1, z = Z2-1, and z = Z2.

  When the voxel image creation processing in step S20 is completed, the voxel image is registered as a 3D texture image (step S60). Specifically, the rendering unit 70 registers a voxel image as a 3D texture image used when executing the 3D texture mapping method.

  FIG. 8 is a diagram showing association in 3D texture mapping. FIG. 8A shows a 3D texture map in a texture coordinate system, and FIG. 8B shows a laminated rectangle defined in a world coordinate system. Four voxels of the 3D texture image are associated with the four vertices of each of the quadrangles forming the stacked quadrilateral, and when no rotation is applied as shown in FIG. 8, each of the quadrangles is the basis of the texture map. It corresponds to the tomographic image and the element image on a one-to-one basis, and four vertices of each rectangle are associated with the tomographic image, the four corner pixels of the element image, and the voxel.

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

  As described later, when there is an instruction to change the angle, the texture map of FIG. 8A is three-dimensionally rotated in the texture coordinate system, and a converted 3D texture image is generated. The correspondence with the image is broken. The voxel image is registered as a 3D texture image via OpenGL (graphic API).

  Subsequently, a rendering process is performed (step S70). FIG. 9 is a diagram showing the configuration of the rendering means 70. As shown in FIG. 9, the rendering unit 70 includes a 3D texture registration unit 71, a coordinate conversion unit 72, a stacked rectangle setting unit 73, and a pixel value calculation unit 74. As described above, in the case of the 3D texture mapping method, the rendering unit 70 is realized mainly by executing a program in the GPU 7 with the assistance of the CPU 1. It is preferable that the value calculating means 74 is executed using a GPU and a frame memory mounted on a video card of a general-purpose computer.

  In the processing of the 3D texture mapping method, first, a projection screen is set. FIG. 10 is an explanatory diagram illustrating an example of a projection screen setting by the rendering unit 70 of the present embodiment. The rendering means 70 sets the screen size (vertical and horizontal pixels, vertical and horizontal aspect ratio) of the rendered image. The rendering means 70 sets either parallel projection (normal appearance rendering) or perspective projection (endoscope mode). Note that “rendered image” is synonymous with “volume rendered image”, and “volume rendered image” may be referred to as “rendered image” in some cases. Then, when the perspective projection is set, the rendering unit 70 sets the perspective projection parameters. The perspective projection parameters include, for example, the viewing angle (focal length) of the camera, the viewpoint position (a viewpoint and a gazing point, the former is the position of the eye, the latter is the position on the object being viewed, and both are on the z-axis. In general, the gaze point includes a fixed point at the origin of the world coordinate system), a clipping position (distance on the z-axis from the viewpoint, and two points near and far). The clipping position may be only near. As shown in FIG. 10, in the case of parallel projection, all lines of sight from the viewpoint are parallel to the z-axis, and the viewpoint is assumed to be virtually distant to the left in infinity. All stacked rectangles defined in the z-axis direction are to be rendered. On the other hand, in the case of perspective projection, since the line of sight expands according to the viewing angle from the viewpoint and the viewpoint enters the inside of the stacked rectangle, the rendering target is far clipping rather than near clipping (usually the most viewed from the stacked square viewpoint). Distant rectangle).

  FIG. 11 is a flowchart showing details of the rendering processing by the 3D texture mapping method. First, the 3D texture registration unit 71 in FIG. 9 registers a voxel image as a 3D texture image in OpenGL. Next, the coordinate conversion unit 72 performs predetermined coordinate conversion on the 3D texture image to generate a converted 3D texture image (S71 to S73).

  Specifically, the coordinate conversion unit 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 xyz axis direction, and an enlargement in each xyz axis direction. Alternatively, a predetermined coordinate transformation parameter including a reduction magnification and a z-axis direction magnification is acquired, and a predetermined coordinate transformation is performed on the 3D texture image using the acquired parameters to generate a transformed 3D texture image.

  The individual steps will be described. The coordinate conversion unit 72 performs scaling and scaling processing in the z-axis direction on the 3D texture image (step S71). Scaling enlarges or reduces the image in the xyz-axis direction at the same magnification Scale, but the z-axis direction scaling process corrects the difference between the resolution Rxy in the xy-axis direction and the resolution Rz in the z-axis direction. Only the specified magnification Scz (= Rxy / Rz) is enlarged or reduced. Then, the coordinate conversion unit 72 performs a rotation process on the 3D texture image using a 4 × 4 matrix defining the rotation in a predetermined three-dimensional coordinate system (step S72).

  Then, the coordinate conversion unit 72 performs an offset process on the 3D texture image using the offset values (Xoff, Yoff, Zoff) in the x, y, and z-axis directions (step S73). Through these processes (S71 to S73), a converted 3D texture image is generated. Next, the stacked rectangle setting means 73 sets a stacked rectangle composed of a plurality of squares (step S74). Specifically, a laminated rectangle in which rectangles on an xy coordinate plane in a three-dimensional space are arranged in the z-axis direction is set. When setting the 3D texture mapping, the stacked rectangle setting unit 73 sets a stacked rectangle in which the same number of squares as the number of the plurality of xy plane images are arranged in the z-axis direction.

  The number of squares in the stacked squares can be set arbitrarily, but if the number of squares does not match the number of tomographic images, interpolation processing works in the z-axis direction, and rounding errors of coordinate values due to coordinate conversion are accumulated, Moire in stripes and grids occurs. Then, as shown in FIG. 8, by defining the four voxel coordinates in the corresponding converted 3D texture image with respect to the four vertices of each square, the correspondence with the converted 3D texture image is performed ( Step S75). The voxel coordinates to be defined are not based on the coordinate system of the color component voxel image with opacity but based on the texture coordinate system of OpenGL (real values in the range of 0 to 1 with respect to the three axes of r, s, and t). is there.

  Specifically, as shown in FIG. 8B, 3D texture mapping is performed on the defined stacked rectangle. In this case, the three-dimensional world coordinates of the four vertices of each rectangle in the world coordinate system are associated with the three-dimensional texture coordinates of the 3D texture image in the texture coordinate system. As shown in FIG. 8, a quadrangle composed of four vertices of (-1, -1, z), (-1,1, z), (1, -1, z), (1,1, z) On the other hand, 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 of the inside of the rectangle) in order from the farthest rectangle on the line of sight parallel to the z-axis direction from the predetermined viewpoint to the rectangle closest to the viewpoint (step S76). At this time, the inside of the rectangle is filled based on the color components (R, G, B) of the voxels of the converted 3D texture image attached to each square, and when the voxels of the converted 3D texture image are overlaid on the same pixel. Alpha blending is also performed based on the opacity. As described above, the color component (R, G, B) composed of the (R, G, B) of the voxel of the converted 3D texture image corresponding to the xyz coordinates of the stacked rectangle is determined from the viewpoint based on the opacity of the voxel. A color component composed of (R, G, B) obtained by alpha blending in the order of distant rectangles is given as a pixel value of the volume rendering image. As a result, a volume rendering image is obtained.

  More specifically, all pixel values (RGB values) of the rendering image are initialized in advance with a background color (for example, R = G = B = 0). The pixel value calculation unit 74 performs scan conversion on the rectangle farthest from the viewpoint with reference to the converted 3D texture image attached thereto. Parallel rendering or perspective projection is performed on the rendering image defined in FIG. 10, and the projected image is projected in the x-axis direction and the y-axis direction from the projected xyz coordinates (−1, −1, −1) in the world coordinate system. The coordinates of the pixel to be filled in the rendered image are (i, j), and the interval u corresponding to the pixel interval in the rendered image (the interval of world coordinates corresponding to one pixel of the rendered image is u. For example, u = Scan conversion is performed at 0.002). In the scan-converted world coordinates (-1 + iu, -1 + ju, -1) (i and j are coordinate values of a rendering image, for example, integers of 0 ≦ i, j ≦ 511), the corresponding texture coordinates (2iu) , 2ju, 0), obtains the color component values (RGB values) and opacity by referring to the 3D texture image after conversion based on the calculated texture coordinates, and obtains the coordinates (i, j) of the corresponding rendering image. ), The alpha blending process is performed based on the pixel values (RGB values) already recorded and the obtained opacity, and the pixel values of the coordinates (i, j) of the corresponding rendering image are updated.

  In this manner, when the alpha blending process is completed for all the corresponding pixels of the rendering image on which the single rectangle farthest from the viewpoint is projected, the z coordinate is changed to v (the world coordinate where the rectangle is arranged in the z-axis direction). Is increased by v. For example, v = 2 / Sz), and the same processing is repeated for the next square (-1 + iu, -1 + ju, -1 + v) in the viewpoint direction to update the rendering image. Scan conversion is performed for all the rectangles, and the rendering image is displayed when the alpha blending process is completed for the world coordinates (-1 + iu, -1 + ju, -1 + kv) (k is the number of the rectangle and 0 ≦ k ≦ Sz−1). Is completed. The alpha blending process is performed as a process according to the following [Equation 5].

[Equation 5]
R ′ = R · α + (1−α) · Rb
G ′ = G · α + (1−α) · Gb
B ′ = B · α + (1−α) · Bb

  Here, R ', G', and B 'are pixel values (RGB values) of the rendering image updated on the projection plane. R, G, and B are color components (RGB values) in the voxels (2iu, 2ju, 2kv) of the converted 3D texture image corresponding to the quadrangle world coordinates (-1 + iu, -1 + ju, -1 + kv), and the colors to be overlapped Corresponding to the value of α is also an α value in the voxel (2iu, 2ju, 2kv) of the converted 3D texture image corresponding to the square 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 rendering image corresponding to the world coordinates (-1 + iu, -1 + ju, -1 + kv), and The values of R ', G', and B 'calculated based on the above [Equation 5] for the rectangle located immediately before the viewpoint on the opposite side to the viewpoint correspond to the rectangle farthest from the viewpoint. In the case of calculating by the calculation, a background color (for example, Rb = Gb = Bb = 0) is given. The volume rendering image is recognized by the rendering image output unit 80 outputting the rendering image obtained as described above.

  FIG. 12 and FIG. 13 are diagrams illustrating display examples of rendered images in which either the clipping setting or the mask setting has been performed. FIGS. 12A and 13A are images in which clipping is set and viewed from the same viewpoint position. FIGS. 12B and 13B show images obtained by setting masks and viewed from the same viewpoint position. The rendering image of FIG. 12A is a case where the color values and opacity α of all voxels in the clipping region including the boundary surface are set to 0, and the rendering image of FIG. 13A is a clipping region including the boundary surface. In some cases, only the opacity α of all voxels is set to 0, and the color values are left unchanged. A comparison between FIG. 12A and FIG. 13A shows that when clipping is set up to the color value, moire occurs as shown in FIG. If not performed, it can be seen that moiré becomes inconspicuous as shown in FIG. There was almost no change in the processing load when clipping was set for color values. The rendering image of FIG. 12B is a case where the color values and opacity α of all voxels in the mask region including the boundary surface are set to 0, and the rendering image of FIG. 13B is a mask region including the boundary surface. In some cases, only the opacity α of all voxels is set to 0, and the color values are left unchanged. Further, comparing FIG. 12B and FIG. 13B, when mask setting is performed up to the color value, moire occurs as shown in FIG. When the setting is not performed, it can be seen that the moiré becomes inconspicuous as shown in FIG. There was almost no change in the processing load between when the mask setting was performed on the color value and when the mask setting was not performed. By setting the value of the opacity α outside the boundary surface to 0 as the moiré processing of the clipping boundary surface or the mask boundary surface, it is possible to suppress the moire while suppressing the processing load. As the moiré processing of the clipping boundary surface or the mask boundary surface, a method of interpolating the color value of a transparent voxel having an opacity α = 0 or applying a process of smoothing the opacity by referring to a nearby voxel may be used. Exists. However, when these processes are used, a fixed amount of processing load is applied. Therefore, when the speed of the process is to be increased, the setting of the opacity α outside the boundary surface, as proposed in the present application, is difficult. Is preferred.

  As described above, in the volume rendering device 100 according to the present embodiment, the rendering unit 70 is realized by the CPU 1 that is a multi-core CPU executing a program. 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 generate an element image group corresponding to each tomographic image group. Perform the process to be created in parallel. A case where the volume rendering apparatus 100 is realized by a personal computer which is a general-purpose computer will be described. FIG. 14 is a diagram illustrating parallel processing when the volume rendering device 100 is implemented by a computer. As illustrated in FIG. 14, when the rendering unit 70 is realized by an application program and divides into, for example, N = 4 tomographic image groups, the application program performs processing on each of the four tomographic image groups by four processing threads. ♯n1 to processing thread ♯n4. Then, the job scheduler of the multitasking operating system in which the rendering means 70 (application program) operates allocates the CPU threads # 1 to # 8 of the plurality of CPU cores.

  FIG. 15 is a diagram illustrating a flow of software for parallel processing when the volume rendering apparatus 100 is implemented by a computer. As shown in FIG. 15, the application program that implements the rendering unit 70 divides the processing for each of the four element image groups into four processing threads # n1 to # n4, and then activates the parallel processing functions. . Then, the job scheduler of the multitask operating system assigns four threads # n1 to # n4 to the CPU threads # 3, # 5, # 7, and # 2 of the plurality of CPU cores, respectively. assign. In this manner, the processing thread to be processed is allocated to the CPU thread of the CPU core that is free by the operating system, so that parallel processing is possible. In the example of FIG. 15, the application program waits for an end message from each CPU thread, and when end messages are obtained from all the CPU threads, creates a voxel image integrating the created element image groups. Then, after the rendering processing is performed in a lump, the rendered image is output to the display unit 6 and displayed.

  As a result of actual measurement on a Windows 10 / 64-bit personal computer on which a 4-core CPU (Intel CORE-i7) is mounted, the configuration shown in FIG. 15 shows that when executing a single thread (voxel by the CPU of a 512 × 512 × 370 full resolution image) A processing speed about 1.5 times lower than the production time (about 0.9 sec) was obtained (same production time: about 0.6 sec). However, some CPU threads waited for processing. Then, when the number of divisions of the element image group was increased to N = 8 and parallel processing was performed by eight processing threads, processing performance almost twice that of the single thread execution was obtained (same generation time: about 0.45 sec) ). Furthermore, as a result of an experiment with N> 8, the processing performance hardly changed and reached a plateau at double, and the access to the shared memory occurred between threads, and the desired parallelization effect was not sufficiently obtained. Was. On the other hand, the rendering process S70 after the creation of the voxel image was executed by the GPU 7, and the generation time of the 512 × 512 full-resolution rendered image was about 0.015 sec.

<3. Decimation of 1 / K>
In the above volume rendering device, the voxel image creating means 30 creates a voxel image with the same number of voxels as the total number of pixels of a plurality of tomographic images, but creates a voxel image with a reduced number of voxels by thinning out voxels. You may do so. In this case, the voxel image creating unit 30 arranges the plurality of tomographic images three-dimensionally in association with every K pixels in each of three two-dimensional xy-axis directions and a z-axis direction orthogonal to the tomographic images. For each of the voxels, the signal value of the voxel is replaced with the color component and the opacity α, and an opacity voxel image in which the color component and the opacity are defined as the voxel value is created. For example, when K = 2, the xyz-axis direction can be thinned out to 各 々 each, so that the number of voxels to be processed becomes ８, and high-speed processing can be performed. Thereby, even when the user interactively switches the color map through the instruction input I / F 4 continuously, the volume rendering images with coarse image quality are sequentially displayed, and the color map is updated. A volume rendering image can be displayed while following.

<4. Degradation of tomographic image>
In the above-described embodiment, a voxel image is created using a tomographic image in which each pixel records a 16-bit signal value as it is. However, a gradation-reduced image is created by converting the tomographic image to 8 bits. Later, a voxel image may be created. By converting the tomographic image into 8-bit gradation, the load in processing each pixel can be reduced, which contributes to high-speed generation of a volume rendering image.

  In this case, a tomographic image gradation conversion unit is further provided, and the gradation conversion is performed every time the tomographic image reading unit 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 / 2-th intermediate tomographic image Do (x, y, z / 2) is read first, and all pixels in the tomographic image are read. Is Dmin, and the maximum value is Dmax. For Dmin and Dmax, a method of obtaining the minimum value and the maximum value from all the tomographic images may be used. However, if a method of determining only the intermediate tomographic images is used, all slices of a large-capacity 16-bit original tomographic image are obtained. It is possible to directly construct only the gradation-compressed tomographic image D8 (x, y, z), which is a gradation-reduced image, on the memory without reading and holding it on the memory. Subsequently, a process according to the following [Equation 6] is executed to calculate a lower limit value Lmin and an upper limit value Lmax.

[Equation 6]
Lower limit value Lmin = (Dmax−Dmin) · γ + Dmin
Upper limit value Lmax = (Dmax−Dmin) · (1−γ) + Dmin

  In [Equation 6], γ is a contrast adjustment width of the gradation-compressed image. As the value approaches 0, the contrast increases but the luminance decreases. Although γ is a real value less than 0.3, γ 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), a process according to the following [Equation 7] is executed to calculate a gradation compressed tomographic image D8 (x, y, z).

[Equation 7]
D8 (x, y, z) = (Do (x, y, z) −Lmin) · 255 / (Lmax−Lmin)
When D8 (x, y, z)> 255: D8 (x, y, z) = 255, When D8 (x, y, z) <0: D8 (x, y, z) = 0

  As shown in the second equation of [Equation 7], when the signal value exceeds the upper limit value Lmax, 255, and when the signal value is lower than the lower limit value Lmin, it is set to 0, and the range from the lower limit value Lmin to the upper limit value Lmax is linearly from 0 to 255. Convert. When a voxel image is created using a gradation-compressed tomographic image, an 8-bit signal value color map 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 the signal values 0 to 255 is used.

  As described above, the preferred embodiments of the present disclosure have been described, but the present disclosure is not limited to the above embodiments, and various modifications are possible. Aspects changed without departing from the spirit of the present disclosure are included in the technical scope of the present disclosure.

1 ... CPU (Central Processing Unit)
2 ... RAM (Random Access Memory)
3. Storage device 4. Instruction input I / F
5 Data I / F
6 Display unit 7 GPU
8 ... Frame memory 10 ... Tomographic image reading means 15 ... Color map reading means 20 ... ROI clipping setting means 25 ... Mask setting means 30 ... Voxel image creating means 40 ... Not Transparency correction means 50 ... Shadow addition means 70 ... Rendering means 71 ... 3D texture registration means 72 ... Coordinate conversion means 73 ... Laminated square setting means 74 ... Pixel value calculation means 80 ...・ Rendered image output means 100 ・ ・ ・ Volume rendering device

Claims (11)

Based on a plurality of two-dimensional tomographic images photographed at predetermined intervals for a predetermined object, rendering is performed by referring to a color map defined by associating color component values and opacity with signal values. A volume rendering device for generating an image,
Referring to the color map, a signal value of each pixel of each of the tomographic images is converted, and a voxel image in which a two-dimensional element image in which a color component and an opacity are defined as pixel values is superposed is created. Means for creating voxel images,
Rendering means for generating a rendering image using the voxel image,
The voxel image creating means divides each of the tomographic images into N tomographic image groups by setting a predetermined number of divisions to N (N> 1), and divides the elemental image into N tomographic image groups. A volume rendering device that executes processing to be created in parallel. The coordinates corresponding to the left end in the x-axis direction of each of the element images are Xso, the coordinates corresponding to the right end in the x-axis direction are Xeo, the coordinates corresponding to the upper end in the y-axis direction are Yso, and the coordinates corresponding to the lower end in the y-axis direction. Let Xeo be Zeo, the z coordinate at which the first element image is arranged be Zso, and the z coordinate at which the Mth element image, which is the total number of tomographic images, be arranged be Zeo, Xso ≦ Xs <Xe ≦ Xeo, Yso ≦ When six coordinates Xs, Xe, Ys, Ye, Zs, and Ze that specify a voxel image region for which a rendering image is to be generated and satisfy Ys <Ye ≦ Yeo and Zso ≦ Zs <Ze ≦ Zeo are defined. ,
The voxel image creating means sets the opacity of each pixel of an element image satisfying x ≦ Xs or x ≧ Xe or y ≦ Ys or y ≧ Ye or z ≦ Zs or z ≧ Ze to 0. 2. The volume rendering device according to 1. When mask data having a value of 0 (not displayed) or 1 (displayed) is defined corresponding to the element image,
The voxel image creating 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 2. The volume rendering device according to claim 1, wherein when 0, the opacity of the pixel of the element image is set to 0.   The voxel image creating means refers to the opacity of at least pixels in the xyz-axis direction, including the pixels, of each of the created element images, and refers to the opacity of the pixels of the element image. The volume rendering device according to claim 1, further comprising a shade adding unit that corrects a value. In the tomographic image group, when one end arranged in order is a head and the other end is a tail,
The shading adding means, after the element image is created for all tomographic images from the first tomographic image to the last tomographic image of the tomographic image group, the element image corresponding to the first tomographic image and the last tomographic image 5. The volume rendering device according to claim 4, wherein a value of a color component of a pixel of the element image corresponding to the tomographic image is corrected. In the tomographic image group, when one end arranged in order is a head and the other end is a tail,
The voxel image creating unit creates the element image for the tomographic image at the end of the tomographic image group, and then converts the tomographic image from the leading tomographic image of the tomographic image group to the preceding tomographic image at the end. On the other hand, creating the element image,
The shading adding means calculates a color component value for an element image corresponding to the leading tomographic image, an element image corresponding to the last tomographic image immediately before the ending, and an element image corresponding to the ending tomographic image. The volume rendering device according to claim 4, wherein the correction is performed. The rendering means,
3D texture registration means for generating a 3D texture image composed of the color component voxel images;
Coordinate conversion means for performing predetermined coordinate conversion on the 3D texture image and generating a converted 3D texture image;
A stack in which quadrangles on an xy coordinate plane of a three-dimensional space are arranged in the z-axis direction, and three-dimensional coordinates of four vertices of each quadrilateral are associated with predetermined three-dimensional coordinates at four predetermined positions of the converted 3D texture image. Stacking rectangle setting means for setting a rectangle,
A color component composed of (R, G, B) of the voxel of the converted 3D texture image corresponding to the three-dimensional coordinates on the stacked rectangle on a line of sight parallel to the z-axis direction from a predetermined viewpoint is converted to the voxel of the voxel. A pixel value calculating means for providing a color component obtained by alpha blending in the order of rectangles farther from the viewpoint based on the opacity as a pixel value composed of (R, G, B) of the rendering image When,
The volume rendering device according to any one of claims 1 to 6, further comprising: The coordinate conversion means,
A predetermined 4 × 4 matrix defining rotation in a predetermined three-dimensional coordinate system, including a viewing angle, a viewpoint position, a clipping position, an offset value in the xyz axis direction, an enlargement or reduction magnification in the xyz axis direction, and a magnification in the z axis direction. Get the parameters of the coordinate transformation,
8. The volume rendering device according to claim 7, wherein the predetermined coordinate conversion is performed on the 3D texture image using the acquired parameters to generate the converted 3D texture image.   9. The volume rendering device according to claim 8, wherein said coordinate conversion means and said pixel value calculation means are executed using a GPU and a 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 9. A voxel structure corresponding to each pixel of a plurality of two-dimensional tomographic images taken at predetermined intervals of a predetermined target object and configured by voxels determined with reference to a predefined color map. Voxel data,
Assuming that a predetermined number of divisions is N (N> 1), an element image created corresponding to each of the tomographic images is divided into N element image groups,
The N element image groups are voxel data processed in parallel by a voxel image creating unit that creates a voxel image by integrating the N element image groups.