JP2020038589A - Rendering device - Google Patents

Abstract

To provide a method of creating a volume rendering image that can reduce the calculation load and shorten the computation time.SOLUTION: A rendering device includes two-dimensional projection image generation means 70 that divides a two-dimensional projection image generation region into N divided generation regions, irradiates a voxel image having each pixel associated with a plurality of two-dimensional tomographic images with a virtual beam in the line-of-sight direction from each pixel in the divided generation regions, extracts a plurality of voxels of the voxel image in which respective virtual beams cross each other, calculates a representative signal value on the basis of signal values of the extracted plurality of voxels, and gives the representative signal value as a pixel value of pixels of the divided generation regions irradiated with virtual beams, thereby generating a divided image for each of the divided generation regions. The two-dimensional projection image generation means 70 concurrently executes processing of generating divided images for the N divided generation regions to generate a two-dimensional projection image.SELECTED DRAWING: Figure 3

Description

  The present disclosure relates to a rendering technique for visualizing three-dimensionally using 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.

  One of the images visualized as a two-dimensional projection image obtained by projecting three-dimensional information in two dimensions is a MIP (Maximum Intensity Projection) image. The MIP image is a method of volume rendering called MIP (also referred to as 3D-MIP to distinguish it from the MIP display mode of MPR), and is generated by an algorithm similar to the ray casting method. Creation of a volume rendering image such as an MIP image has a problem that the calculation load is high and the calculation time is long.

  On the other hand, the MIP pixel is searched for while thinning out the generated MIP image at a predetermined interval in the xy-axis direction. For the thinned-out pixel, the position in the depth direction of the already searched MIP pixel in the vicinity is determined. There has been proposed a method for increasing the speed by searching for MIP pixels within a range limited to the center (see Patent Document 1). In addition, if the method of recalculating at an arbitrary angle is used in rotating and displaying the MIP image, real-time display cannot be performed. Therefore, the MIP image is limited to eight directions, and the image projected on the regular octahedron is calculated in advance, and the user There is also proposed a method of rotating in any one of eight directions in real time by designation and presenting an angle of a display image (see Patent Document 2).

JP 2005-46207 A JP 2001-229399 A

  In the techniques described in Patent Literature 1 and Patent Literature 2, a problem is hardly generated in the MIP display mode because a method of narrowing the search range of MIP pixels is used. However, they are arranged in the ray direction like RaySum. There is a problem that a display mode of calculating an average value of all voxels cannot be realized. That is, in Patent Literature 1, the effect of operation suppression is limited to the MIP mode, and the effect of operation suppression does not work in the RaySum mode or the MinIP mode. Further, the technique described in Patent Literature 2 has a problem in that it is not possible to fine-tune the angle to an arbitrary angle because recalculation is required for angles other than eight directions. That is, in Patent Literature 2, the effect of calculation suppression is limited to a specific angle, and when the angle is finely adjusted to an arbitrary angle, the effect of calculation suppression does not work.

  In short, the techniques described in Patent Literature 1 and Patent Literature 2 limit the applicable display methods.

The present disclosure includes a plurality of means for solving the above-described problems.
A representative signal calculated based on a signal value of each pixel passing through the line of sight based on a plurality of two-dimensional tomographic images in which an image of an object is captured at predetermined intervals and a signal value is assigned to each pixel. A rendering device for generating a two-dimensional projection image composed of values,
Dividing a generation area, which is an area for generating the two-dimensional projection image, into N divided generation areas, and irradiating virtual rays in the line-of-sight direction from each pixel forming the divided generation areas,
A plurality of voxels of the voxel image in which each pixel constituting the plurality of two-dimensional tomographic images corresponds to each voxel image where each virtual ray intersects, and a single voxel is extracted based on the signal values of the plurality of extracted voxels. Two-dimensionally generating a divided image corresponding to each divided generation region by calculating the representative signal value of the divided generation region as a pixel value of a pixel of the divided generation region irradiated with the virtual ray. Having projection image generation means,
The two-dimensional projection image generation means generates the two-dimensional projection image by executing a process of generating a divided image in parallel with the N divided generation regions.

In addition, the present disclosure,
A representative signal calculated based on a signal value of each pixel passing through the line of sight based on a plurality of two-dimensional tomographic images in which an image of an object is captured at predetermined intervals and a signal value is assigned to each pixel. To generate a monochrome image consisting of values,
Computer
The generation region of the monochrome image is divided into N divided generation regions, and a voxel image in which the plurality of two-dimensional tomographic images and each pixel are associated with each other is imagined in a line-of-sight direction from each pixel of the divided generation region. Irradiate light rays, extract a plurality of voxels of the voxel image where each virtual ray intersects, calculate a single representative signal value based on the signal values of the extracted plurality of voxels, and calculate the calculated representative signal value. By giving as a pixel value of the pixel of the divided generation region irradiated with the virtual ray, for each divided generation region, function as a monochrome image generating means for generating a divided image,
The monochrome image generation means executes a process of generating a divided image in parallel with the N divided generation regions, and integrates the obtained divided images to generate the monochrome image. provide.

  According to the present disclosure, it is possible to improve the speed of a process of generating a two-dimensional projection image for three-dimensional visualization using a plurality of tomographic images.

FIG. 3 is a diagram illustrating a plurality of tomographic images and a generated MIP image that is a monochrome two-dimensional projected image. 1 is a hardware configuration diagram of a rendering device 100 according to an embodiment of the present disclosure. FIG. 2 is a functional block diagram illustrating a configuration of a rendering device according to the embodiment. 5 is a flowchart illustrating a processing operation of the rendering device according to the embodiment. 9 is a flowchart illustrating details of MIP image generation. It is a conceptual diagram of coordinate conversion. 6 is a flowchart illustrating a divided image generation process in step S540 in FIG. 8 is a flowchart showing a process of calculating a representative signal value Vm in step S620 of FIG. 6 is a flowchart illustrating a process in which a starting point coordinate search unit 77 searches for starting point coordinates in each pixel. FIG. 7 is a diagram illustrating an example of a generated MIP image. FIG. 9 is a diagram illustrating a parallel processing when the rendering device 100 is implemented by a computer. FIG. 3 is a diagram illustrating a flow of software for parallel processing when the rendering apparatus 100 is implemented by a computer. 6 is a flowchart illustrating a processing procedure of progressive rendering. It is a figure showing the example of the screen displayed by progressive rendering. It is a figure showing various modes at the time of dividing a generation field into a plurality of division generation fields. It is a figure showing the flow by the software of the parallelization process in the case of performing a display process each time each processing thread ends.

Hereinafter, preferred embodiments of the present disclosure will be described in detail with reference to the drawings.
The present disclosure is to generate a MIP image that is a two-dimensional projection image having a predetermined value based on a plurality of tomographic images. FIG. 1 is a diagram showing a plurality of tomographic images and a generated MIP image which is a two-dimensional projection image. 1A shows a plurality of tomographic images, FIG. 1B shows one tomographic image, and FIG. 1C shows a two-dimensional projected 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 rendering apparatus according to an embodiment to be described later, a rendering process is performed on a plurality of tomographic images as illustrated in FIG. 1A, and a single two-dimensional projected image as illustrated in FIG. Generate The two-dimensional projected image shown in FIG. 1C is called a MIP image in a broad sense, and each pixel has a single signal value. The MIP image in a broad sense includes “MIP in a narrow sense” in which the maximum value is recorded as a single signal value, MinIP in which the minimum value is recorded as a single signal value, and RaySum in which the average value is recorded as a single signal value. In. Hereinafter, in this specification, a MIP image in a broad sense, which is a two-dimensional projection image in which some single signal value is recorded, is simply referred to as an MIP image. Although the two-dimensional projected image generated in the present embodiment is a monochrome MIP image, a two-dimensional projected image other than a monochrome image can be generated. Further, the rendering process in the rendering device according to the present disclosure is not limited to “MIP in a narrow sense”, “MinIP”, “RaySum”, and the like, but is also applied to other rendering processes such as ray casting.

<1. Device Configuration>
FIG. 2 is a hardware configuration diagram of the rendering device 100 according to an embodiment of the present disclosure. The 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, an instruction input I / F (interface) 4 such as a keyboard and a mouse; 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 as a display device such as a liquid crystal display; and a GPU as an arithmetic processing unit specialized in graphics (Graphics Processing Unit) 7 and a frame memory 8 for holding an image to be displayed on the display unit 6. It is connected 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 logically having a plurality of CPU cores.

  FIG. 3 is a functional block diagram illustrating a configuration of the rendering device according to the present embodiment. In FIG. 3, 10 is a tomographic image reading unit, 15 is a coordinate conversion parameter input unit, 20 is an ROI clipping setting unit, 25 is a mask setting unit, 30 is a voxel image creating unit, 70 is a two-dimensional projection image generating unit, and 80 is 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 coordinate conversion parameter input means 15 is a means for inputting parameters for coordinate conversion. 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 creation means 30 arranges a plurality of two-dimensional tomographic images at predetermined intervals in the z-axis direction, and, for each voxel in which each pixel of the tomographic image is arranged three-dimensionally, This is a means for creating a voxel image having a signal value as a value.

  The two-dimensional projection image generation unit 70 is a unit that generates a MIP image that is a monochrome volume rendering image using voxel images. The two-dimensional projection image generation unit 70 can generate three types of MIP, MInIP, and RaySum in a narrow sense, as MIP images that are two-dimensional projection images. The image output unit 80 is a unit that outputs the generated MIP image in a predetermined mode. Normally, display output is performed on the display unit 6 or the like.

  The tomographic image reading means 10 is realized mainly in the data input / output I / F 5 with the assistance of the CPU 1. The coordinate transformation parameter input means 15 is realized mainly in the instruction input I / F 4 with the assistance of the CPU 1. The ROI clipping setting means 20 is realized mainly at the instruction input I / F 4 with the assistance of the CPU 1. The mask setting unit 25 is realized mainly in the instruction input I / F 4 and the display unit 6 with the assistance of the CPU 1. The two-dimensional projection image generation unit 70 is realized by the CPU 1 that is a multi-core CPU executing a program. This program divides the MIP image generation region into a plurality of division generation regions, allocates each division generation region to each core, and performs processing in which each core generates a MIP image division image in parallel. The 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. In the present embodiment, the CPU is preferably a multi-core CPU. Even if the CPU is not a multi-core CPU, the CPU must have a configuration that enables parallel processing that can execute programs in parallel, such as a multiprocessor configuration including a plurality of CPUs. 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 processing data, and is not limited to a general-purpose computer such as a personal computer, but also a tablet equipped with a multi-core CPU. Mobile terminals and computers embedded in various devices.

<2. Processing Operation>
Next, the processing operation of the rendering apparatus shown in FIGS. 2 and 3 will be described. FIG. 4 is a flowchart illustrating the processing operation of the 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. 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 the specification of the present application, means a region of interest in a three-dimensional space (scientifically, 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, an area where a MIP image and a voxel image are to be created is shown. By setting the ROI, a MIP image in which a subject is cut three-dimensionally at an arbitrary position can be generated, and is used to depict an internal organ or an internal organ hidden by a body surface or a 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. Assuming that the coordinates are Yeo, the coordinates corresponding to the tomographic image of the first slice are Zso, and the coordinates corresponding to the tomographic image of the last slice are Zeo, Xso ≦ Xs <Xe ≦ Xeo, Yso ≦ Ys <Ye ≦ Yeo, Zso ≦ Zs A region defined by a rectangular parallelepiped defined by six coordinates [Xs, Xe], [Ys, Ye], and [Zs, Ze] that satisfies <Ze ≦ Zeo is defined as an ROI and specified as a processing target. 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 an MIP image exposing the organ to be observed is obtained, the number of pixels to be processed is reduced, the processing load is suppressed, and the responsiveness is improved. In many cases, the setting is performed because there is a double merit that the setting can be performed, but it is possible to omit the processing instead of the essential processing. 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) rendering image is generated in the MIP image generation 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 mask setting method is not particularly limited, but in the present embodiment, the tomographic image group read in step S10 is displayed on the screen as an MPR image (three axial coronal sagittal screens), and the method described in the present specification is used. A MIP image created based on the image or a volume rendering image created by a known method is displayed, and an image that is most easily instructed is selected from among the two-dimensionally projected images, and the selected image 2 While specifying the dimensional area, a range not to be displayed by masking and an area to be displayed without masking are specified.

  Further, in the case of the MIP image, since control such as making a specific signal region transparent by a color map as in a known volume rendering image cannot be performed, only a bone region is clearly displayed without any modification. The organ is blurred. Therefore, in the present embodiment, a function of automatically generating mask data in which a bone region is masked by applying a known bone removal method is also provided. 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). . In the case of the MIP image, if no mask data is used, only the bone region is clearly seen and other organs are blurred. Therefore, the process of creating such a three-dimensional mask is almost indispensable, but may be omitted. It is possible.

  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, the signal value of each pixel included in the effective area (ROI, specified by the mask) of the tomographic image after clipping setting and mask setting is given to each corresponding voxel, and the voxel image is create.

  In step S20, the voxel image creating means 30 outputs the signal at each pixel (x, y, z) corresponding to the region of Xs ≦ x ≦ Xe, Ys ≦ y ≦ Ye and Zs ≦ z ≦ Ze of the plurality of tomographic images. Get the value. By performing this process on all the pixels in the plurality of (Sz) tomographic images, a three-dimensional voxel image V (x, y, z) corresponding to the plurality of tomographic images is obtained. In the present embodiment, the signal value of the DICOM format tomographic image group Do (x, y, z) is given as a voxel value as it is, so that each pixel of the voxel image V (x, y, z) has one pixel. It is recorded in 16 bits and takes values from -32768 to 32767. This voxel image does not need to be stored in the non-volatile storage device 3 and may be created temporarily in the RAM 2.

  Next, the two-dimensional projection image generation means 70 generates a MIP image (step S70). The MIP image, which is a monochrome two-dimensional projected image, is composed of representative signal values calculated based on the signal values of each pixel passing through the line of sight. FIG. 5 is a flowchart showing details of MIP image generation. First, the size (the number of pixels) of the MIP image to be generated is set (step S510). Since the MIP image, which is a rendering image, is a two-dimensional projection image, the size is set to the size of the x coordinate and the size of the y coordinate, and is often set to be the same as the size of the tomographic image. Can be set. In the present embodiment, the size of the x coordinate and the size of the y coordinate are the same, and are set to Size.

  Next, the generated coordinate conversion parameters are set (step S520). The two-dimensional projection image generation unit 70 sets a coordinate conversion parameter that is a parameter commonly used when performing coordinate conversion in the subsequent steps. The setting of the coordinate conversion parameters is such that, when irradiating the voxel with virtual light, the three-dimensional coordinates in the viewpoint coordinate system based on the viewpoint are set so that the coordinate conversion processing can be smoothly performed sequentially for each voxel. Is a process for converting the voxel image into three-dimensional coordinates in a voxel coordinate system in which a voxel image is defined, and obtaining a signal value of a corresponding voxel in the voxel coordinate system.

  Here, a conceptual diagram of coordinate conversion and a specific example of coordinate conversion parameters will be described. FIG. 6 is a conceptual diagram of the coordinate conversion. Since the MIP image obtained by rendering is an image viewed from the viewpoint, the viewpoint coordinate system serving as a reference for each pixel (x, y) of the MIP image does not always match the voxel coordinate system. Therefore, the two-dimensional projection image generation means 70 needs to convert the voxel image defined in the voxel coordinate system into the viewpoint coordinate system. That is, as shown in FIG. 6 (a), a virtual ray (ray) is originally emitted from each point on the projection plane to the voxel image in an arbitrary direction, and a designated range in the voxel image where the virtual rays intersect is designated. Are extracted, and a single representative signal value (maximum value, minimum value, average value) is calculated based on the extracted signal values of the plurality of voxels. However, in order to simplify the calculation, the voxel image is subjected to coordinate transformation, and as shown in FIG. 6B, representative signal values (maximum value / minimum value・ Average value is calculated.

  When performing the coordinate conversion, various coordinate conversion parameters are used. If the positional relationship between the viewpoint coordinate system and the voxel coordinate system is the same, it is considered that the coordinate conversion parameters are also the same. That is, in the two-dimensional projection image generating means 70, the coordinate conversion parameters for performing the coordinate conversion to the voxel coordinate system for all the voxels referenced in the viewpoint coordinate system are the same. Therefore, the two-dimensional projection image generating means 70 previously calculates the coordinate conversion parameters in step S520, and uses the coordinate conversion parameters as common coordinate conversion parameters in the subsequent steps. The following are examples of the coordinate conversion parameters.

A 3 × 3 matrix defining rotation in a predetermined three-dimensional coordinate system (an inverse matrix of the 3 × 3 matrix specified and defined by the instruction input I / F 4 because projection transformation is performed in the reverse direction)
-Offset value in xyz axis direction: Xoff, Yoff, Zoff (unit: voxel)
Enlargement or reduction magnification: Scale (the same value is used for each axis of xyz)
• z-axis direction scaling factor: Scz = Rxy / Rz (ratio between resolution Rxy in xy-axis direction and resolution Rz in z-axis direction)
Sub-sampling magnification of the virtual ray: Sray (usually an integer of 1 or more, a larger value makes it coarser, and a real number less than 1 can also be specified)

  After the size of the MIP image and the coordinate conversion parameters are set, an MIP image is generated based on the set size and the coordinate conversion parameters. The generation of the MIP image is performed by dividing the generation area determined by the set size into a plurality of parts, and processing the obtained divided generation areas in parallel to generate a divided image. First, the two-dimensional projection image generation unit 70 initially sets a divided generation area to be divided for performing parallel processing (step S530). In the present embodiment, the x coordinate is not divided but the y coordinate is divided into N pieces. Specifically, the start coordinate of the y coordinate of each divided generation area is Y1, the end coordinate is Y2, and the start coordinate Y1 and the end coordinate Y2 are the divided generation areas of the entire MIP image specified by the set size. Assign. As initial values, Y1 = 0 and Y2 = Size / N−1 are set.

  After the division generation area is set, the generation processing of the division image of the MIP image is performed (step S540). As described above, since the MIP image is divided into N only at the y coordinate, the number of pixels at the x coordinate of the divided image is Size, and the number of pixels at the y coordinate of the divided image is (Y2−Y1 + 1). Become. Therefore, in step S540, a divided image of Size × (Y2-Y1 + 1) pixels is generated. The processing in step S540 will be described later.

  The process of step S540 is performed by a plurality of CPU cores in parallel. For this reason, if the generation processing of the divided image is allocated to one CPU core for the first divided range and the execution processing is started, the next divided generation area is set and the start of the generation processing is advanced without waiting for the end of the processing (step S550). ). Specifically, the values of Y1 and Y2 are each increased by Size / N, and the start coordinate Y1 and the end coordinate Y2 of the y coordinate in the next divided generation area are newly set.

  After the new divided generation area is set in step S550, if the new divided generation area exceeds the size of the entire MIP image, that is, if Y2> Size-1 is satisfied, all divided generation areas are set. Assuming that the setting and the start of the divided image generation process have been completed, the process of step S70 shown in FIG. 5 is terminated. However, at this stage, the processing of generating the individual divided images in step S540 has not been completed. After setting the new divided generation area in step S550, if the new divided generation area does not exceed the size of the MIP image, that is, if Y2 ≦ Size-1 is satisfied, all the divided generation areas are set. Since the process has not been completed, the process returns to step S540, and the process of allocating to another CPU core and starting generation of the divided image is repeated.

  The process in step S540 will be described. FIG. 7 is a flowchart showing the divided image generation processing in step S540 of FIG. As described above, in the divided image generation processing, the processing is performed in the entire range of the number of pixels Size in the x-axis direction, and in the range from the start coordinate Y1 to the end coordinate Y2 in the y-axis direction. First, all the initial values of each pixel of the generated divided image are set to 0 (Image (x, y) = 0). Also, a representative signal value calculation mode mode (0: MIP, 1: MinIP, 2: RaySum) is set. Then, the process according to the flowchart of FIG. 7 is performed on each of the two-dimensional coordinates (x, y) (0 ≦ x ≦ Size−1, Y1 ≦ y ≦ Y2).

  First, initial setting of a pixel to be processed is performed (step S610). Initially, the pixel to be processed is set to y = Y1, x = 0 where the values of the x coordinate and the y coordinate are both minimum. Next, a representative signal value Vm to be given to the pixel to be processed is calculated (step S620). The processing in step S620 will be described later.

  When the representative signal value Vm is calculated by the processing in step S620, the representative signal value Vm is given and recorded as the value of the pixel (x, y) to be processed. Specifically, a process for setting Image (x, y) = Vm is performed (step S630). After recording the representative signal value corresponding to the pixel to be processed, the pixel to be processed is changed. Specifically, the value of the x coordinate is incremented by 1 (step S640). As a result of changing the pixel to be processed in step S640, if the pixel to be processed is within the x-coordinate size, that is, if x <Size is satisfied here, the process returns to step S620 to calculate the representative signal value Vm. I do. As a result of changing the pixel to be processed in step S640, if the pixel to be processed exceeds the size of the x coordinate, that is, if x ≧ Size is satisfied here, the value of the y coordinate is incremented by 1 (step S650). ). As a result of changing the pixel to be processed in step S650, if the pixel to be processed is within the size of the y coordinate, that is, if y ≦ Y2 is satisfied, the process returns to step S620 to calculate the representative signal value Vm. I do. As a result of changing the pixel to be processed in step S650, if the pixel to be processed exceeds the size of the y coordinate, that is, if y> Y2 is satisfied, it is assumed that the processing has been completed for the assigned divided range. The process of step S540 shown in FIG. 7 ends.

  Next, the process of step S620 in FIG. 7 will be described in detail. FIG. 8 is a flowchart showing the process of calculating the representative signal value Vm in step S620 of FIG. First, variables are initialized (step S710). Specifically, the representative signal value Vm, counters Cnt, and Zs are initialized. The initial value of Zs is calculated as Zs = Size / Sray-1. The representative signal value Vm differs depending on the mode. Of the three types of modes, Vm = −32768 when mode = 0, Vm = 32767 when mode = 1, and Vm = 0 when mode = 2. The counter Cnt is used only when mode = 2, and is initialized as Cnt = 0.

  Next, a search for starting point coordinates is performed (step S720). The processing in step S720 will be described with reference to FIG. FIG. 9 is a flowchart showing a process in which the starting point coordinate searching means 77 searches for starting point coordinates in each pixel. The starting point coordinate search unit 77 executes a process according to the flowchart shown in FIG. 9 for each target pixel (x, y).

  The starting point coordinate search means 77 sets a target pixel (x, y) for searching for starting point coordinates and three-dimensional coordinates (x, y, z) of a viewpoint coordinate system (step S301). The initial value of the z coordinate which is the search start coordinate is z = Zs (upper limit). The lower limit of the z coordinate is 0, and the upper limit is Zs. The upper limit value Zs is the z coordinate of the viewpoint.

  The origin coordinate search means 77 performs coordinate transformation on the three-dimensional coordinates (x, y, z) to calculate a signal value Vs of a voxel of a corresponding voxel image (step S302). The coordinate conversion is an operation for converting the viewpoint coordinate system to the voxel coordinate system, and is the reverse of the conversion operation on the GUI side. On the GUI side, it is assumed that clipping by ROI, scaling, scaling processing in the z-axis direction, rotation, and offset (simultaneous in the xyz-axis directions) are performed in this order, and the three-dimensional coordinates (x, y, z) of the given viewpoint coordinate system The coordinates (xr, yr, zr) of the real number corresponding to the voxel of the voxel image corresponding to each pixel of the plurality of tomographic images Do (x, y, z) of the DICOM format corresponding to (integer value) on a one-to-one basis. calculate.

  Specifically, first, the three-dimensional coordinates (x, y, z) of the viewpoint coordinate system, which are integer values, are converted into real values (xx, yy, zz). The viewpoint coordinate system is a rectangular parallelepiped of Size × Size × Size / Sray, and is defined independently of the voxel structure of Sx × Sy × Sz. Since the size Sx × Sy of the tomographic image in the DICOM format is usually Sx = Sy = 512, it is often set to Size = 512.

  Considering that the z-axis direction of the viewpoint coordinate system is reduced to the virtual ray subsample (1 / Sray) times (Sray> 1), first, processing according to the following [Equation 1] is executed. Offset processing is performed simultaneously in the xyz-axis direction.

[Formula 1]
xx = x-Size / 2 + Xoff
yy = y-Size / 2 + Yoff
zz = zSray-Size / 2 + Zoff

  Subsequently, using the generated 3 × 3 rotation inverse matrix, a process according to the following [Equation 2] is executed to perform a rotation process.

[Equation 2]
xx ′ = R11 · xx + R12 · yy + R13 · zz
yy '= R21 · xx + R22 · yy + R23 · zz
zz ′ = R31 · xx + R32 · yy + R33 · zz

  (Xx ', yy', zz ') are obtained as the coordinates after the rotation processing. Subsequently, processing according to the following [Equation 3] is executed, and scaling and scaling processing in the z-axis direction are simultaneously performed.

[Equation 3]
xr = xx / Scale + Sx / 2
yr = yy / Scale + Sy / 2
zr = zz / (Scale · Scz) + Sz / 2

  Thus, the coordinates (xr, yr, zr) of the voxel coordinate system are calculated. These coordinates (xr, yr, zr) are real values. Next, the starting point coordinate search means 77 converts the calculated coordinates (xr, yr, zr) of the voxel coordinate system into integer values. Specifically, a process according to the following [Equation 4] is executed to convert the value into an integer.

[Equation 4]
xi = INT [xr + 0.5]
yi = INT [yr + 0.5]
zi = INT [zr + 0.5]

  In [Equation 4], INT means that the value in [] is rounded down to the nearest integer to convert it into an integer. Therefore, by the processing according to [Equation 4], the values of the coordinates (xr, yr, zr) of the real value are rounded off to obtain the coordinates (xi, yi, zi) of the integer value.

  When the coordinates (xi, yi, zi) of the voxel coordinate system corresponding to the coordinates (x, y, z) of the viewpoint coordinate system are obtained in this way, the following [Equation [Math. 5] to obtain a voxel signal value Vs corresponding to the coordinates (xi, yi, zi) of the voxel coordinate system.

[Equation 5]
Mask (xi, yi, zi) = 0 (the place to be masked), or xi <Xs or xi> Xe or yi <Ys or yi> Ye or zi <Zs or zi> Ze (outside the ROI) ( That is, when the area is not an effective area): Vs = −99999 (invalid value)
In cases other than the above (that is, in the case of an effective area): Vs = V (xi, yi, zi)

  In [Equation 5], if it is outside the masked portion or ROI, -99999 is set to Vs, but it can be shown that the signal value is out of the valid range and invalid. Any value may be used as long as it is a value.

  Here, the effective area indicates an area that is not masked by the mask data (Mask (xi, yi, zi) = 1) and is within the ROI. That is, in both the ROI clipping and the mask setting, the area set as the area to be displayed is the effective area. The effective range is a range of signal values that the pixel can take, and in the case of 16 bits, the effective range is −32768 to 32767. In the above [Equation 5], if the converted coordinates are not in the valid area, values outside the valid range are set as invalid values.

  If Vs is smaller than −32768, the process proceeds to step S303. If Vs is −32768 or more, the process proceeds to step S305. This “−32768” is a lower limit value that a 16-bit voxel signal value can take. That is, if Vs is an invalid value, the process proceeds to step S303, and if Vs is a valid value, the process proceeds to step S305.

  In step S303, the starting point coordinate search means 77 skips a predetermined number of voxels in the z-axis direction. Specifically, a predetermined number m2 (for example, m2 = 8) is subtracted from the current z coordinate. Note that m2> m1 as compared with the predetermined number m1 used in steps S740 and S760. In the present embodiment, m1 = 1. As a result of the subtraction, if the z coordinate is 0 or more, the process returns to step S302 to repeat the same processing. If the z coordinate has become less than 0, the process proceeds to step S304.

  In step S304, the starting point coordinate searching unit 77 determines that the starting point coordinates have not been found for the z coordinate corresponding to the current target pixel (x, y) (the original lower limit value of z is 0). In this case, for example, a negative value such as -1 is set, and the search processing of the starting point coordinates is ended.

  When the process proceeds to step S305, the starting point coordinate search means 77 initializes the variable i (step S305), and then adds m1 (step S306). If i has reached m2, or if i has been added to the current z-coordinate to reach the upper limit Zs, which is the viewpoint z-coordinate, the flow advances to step S308 to add m1 to the current z-coordinate for the z-coordinate. Is set to the value obtained by adding i before the execution. If i has not reached m2, and even if i is added to the current z coordinate, it does not reach the upper limit Zs, which is the viewpoint z coordinate, the process proceeds to step S307.

  In step S307, the starting point coordinate search unit 77 performs coordinate conversion on the three-dimensional coordinates (x, y, z + i) to calculate the corresponding voxel signal value Vs. The calculation procedure is the same as step S302. If Vs is equal to or larger than -32768, the process returns to step S306, and m1 is added to i, and the same processing is repeated. If Vs is smaller than -32768, the process proceeds to step S308.

  In step S308, the starting point coordinate search unit 77 sets the z coordinate (z + i-m1) immediately before step S307 corresponding to the current target pixel (x, y), and ends the search for the starting point coordinates. .

  Eventually, the starting point coordinate search means 77 subtracts the z coordinate by a predetermined number m2 (m2> m1) larger than the predetermined number m1 in a range from a predetermined upper limit to a predetermined lower limit in accordance with the processing shown in FIG. For each three-dimensional coordinate, a coordinate conversion is performed to calculate a voxel signal value (S302). If the voxel signal value at the z coordinate after subtracting a predetermined number m2 is included in the effective range (S305), Then, a predetermined number m1 is increased from the z coordinate after the subtraction (S306), and coordinate conversion is performed on each three-dimensional coordinate (x, y, z) to calculate a voxel signal value (S307). This processing is not limited to the MIP mode, and can be used in the RaySum mode or the MinIP mode, and the effect of suppressing the calculation can be obtained.

  When the process shown in FIG. 9 is executed and the starting point coordinate z is searched, the process of step S720 shown in FIG. 8 has been completed. If the obtained starting point coordinate z is z <0, it means that there is no valid voxel, and the process proceeds to step S770. If the obtained starting point coordinate z satisfies z ≧ 0, the three-dimensional coordinates (x, y, z + i) are subjected to coordinate conversion to obtain a voxel signal value Vs (step S730).

  Specifically, the processing according to [Equation 1] to [Equation 3] in steps S302 and S308 is executed to convert the given three-dimensional coordinates (x, y, z) (integer value) of the viewpoint coordinate system. A real number coordinate (xr, yr, zr) corresponding to a voxel of a voxel image corresponding to each pixel of the corresponding plurality of tomographic images Do (x, y, z) in the DICOM format is calculated.

  Next, the calculated coordinates (xr, yr, zr) of the voxel coordinate system are converted into integer values. Different from the processing in steps S302 and S308, first, the coordinates obtained by rounding down the decimal point to an integer are set to (xia, yia, zia), and the fractions rounded down to the decimal point are set to (wx, wy, wz). That is, xr = xia + wx, yr = yia + wy, and zr = zia + wz.

  Next, using the ROI and the mask data Mask, a process according to the following [Equation 6] is executed, and a voxel signal value Vs corresponding to the coordinates (xia, yia, zia) of the voxel coordinate system is obtained. In the following [Equation 6], processing is performed in three ways.

[Equation 6]
1) Mask (xia, yia, zia) = 0 (the place to be masked) or xia <Xs or xia> Xe or yia <Ys or yia> Ye or zia <Zs or zia> Ze (outside the ROI range) Case: Vs = -99999 (invalid value)
2) If the above 1) is not satisfied and xia + 1> Xe or yia + 1> Ye or zia + 1> Ze, Vs = V (xia, yia, zia) without performing the interpolation process
3) When the above 1) and 2) are not satisfied, the following interpolation processing is executed. Vs = (1−wz) (1−wy) (1−wx) · V (xia, yia, zia)
+ (1-wz) (1-wy) · wx · V (xia + 1, yia, zia)
+ (1-wz) · wy · (1-wx) · V (xia, yia + 1, zia)
+ (1-wz) · wy · wx · V (xia + 1, yia + 1, zia)
+ Wz · (1-wy) (1-wx) · V (xia, yia, zia + 1)
+ Wz · (1-wy) · wx · V (xia + 1, yia, zia + 1)
+ Wz · wy · (1-wx) · V (xia, yia + 1, zia + 1)
+ Wz · wy · wx · V (xia + 1, yia + 1, zia + 1)

  In (3) of [Equation 6], when performing the interpolation process, eight integer-valued voxel coordinates (xia, yia, zia) and the like located near the real-valued voxel coordinates (xr, yr, zr) are used. Using this value, a value obtained by multiplying each signal value by a predetermined weight wx, wy, wz, (1-wx), (1-wy), (1-wz) is used as a three-dimensional coordinate. Is given as the signal value Vs of the voxel corresponding to. In [Equation 6], if it is outside the masked location or ROI, that is, if it is outside the effective area, Vs is set to -99999, but a value that can indicate that it is invalid is set to Vs. Any value may be used.

  If Vs is smaller than −32768, the process proceeds to step S740, and if Vs is −32768 or more, the process proceeds to step S750. This “−32768” is a lower limit value that a 16-bit voxel signal value can take. That is, if Vs is an invalid value, the process proceeds to step S740, and if Vs is a valid value, the process proceeds to step S750.

  In step S740, the starting point coordinate search means 77 reduces the coordinates of the viewpoint by a predetermined number m1. This predetermined number m1 is an integer. As described above, in the present embodiment, m1 = 1 and m1 <m2. Then, the process returns to step S720 to repeat the same processing.

  In step S750, the starting point coordinate search means 77 sets the representative signal value Vm. In step S750, processing differs depending on the set mode. Specifically, a process according to the following [Equation 7] is executed to set the representative signal value Vm.

[Equation 7]
mode = 0 (MIPmode): If Vs> Vm, Vm = Vs
mode = 1 (MinIPmode): If Vs <Vm, Vm = Vs
mode = 2 (RaySummode): Vm ← Vm + Vs, cnt ← cnt + 1

  After the representative signal value Vm is set, z is reduced by a predetermined number m1 (step S760). In the present embodiment, m1 = 1 and m1 <m2. If z is equal to or greater than 0, the process returns to step S730 to repeat the same processing. If z is less than 0, there are no more voxels of interest, and the process proceeds to step S770. As shown in the third equation of [Equation 1], z is multiplied by the sub-sampling magnification Sray of the virtual ray at the time of the offset processing in step S302. Therefore, in the z-axis direction, the voxel is referenced at every predetermined reference interval (m1 · Sray), and the voxel signal value Vs is acquired. Therefore, if (m1 · Sray)> 1, in the z-axis direction, processing is performed not by referring to all voxels but by referring to thinned voxels. Therefore, the number of voxels to be processed in the z-axis direction is reduced, and efficiency is improved. In this embodiment, since m1 = 1, it is preferable that Sray> 1. When Sray ≦ 1, it is preferable to set m1 such that the value satisfies (m1 · Sray)> 1. The reference interval (m1 · Sray) may be (m1 · Sray) = 1, but more preferably (m1 · Sray) ≧ 2.

  In steps S770 and S780, processing according to the mode is performed. First, in step S770, it is determined whether or not mode = 2. If the mode is not mode = 2, that is, if the mode is 0 or 1, the process ends. If mode = 2, Vm is divided by cnt to obtain a new Vm (step S780). Then, the process ends.

  The process of obtaining the representative signal value Vm as described above is performed by the N CPUs executing the process of step S540 shown in FIG. 5 on each of the N divided images. By integrating the divided images, a MIP image, which is a set two-dimensional projected image of Size × Size, is generated. It is not always necessary to integrate the set two-dimensional projected images of Size × Size. When displaying the divided images sequentially, the divided images are not integrated.

  After the generation of the MIP image, a process of displaying the MIP image is performed. Display devices that display MIP images often have 8 bits per pixel. Therefore, when displaying a MIP image recorded with 16 bits per pixel, it is necessary to perform display processing with the gradation of the MIP image reduced to match the bit number of the display unit 6. Specifically, by executing the processing according to the following [Equation 8], the gradation of the MIP image is converted to 8-bit.

[Equation 8]
Image ′ (x, y) = (Image (x, y) −Lmin) · 255 / (Lmax−Lmin)
If Image '(x, y)> 255: Image' (x, y) = 255, If Image '(x, y) <0: Image' (x, y) = 0

  As shown in the second expression of [Equation 8], 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 changed from 0 to 255. Convert.

  The lower limit Lmin and upper limit Lmax used in [Equation 8] can be arbitrarily changed and adjusted by the user through the GUI during display. In the medical field, there is a custom to use a window width (window width: WW) and a window level (window level: WL) of the DICOM standard instead of the lower limit value Lmin and the upper limit value Lmax. Therefore, the GUI is also designed based on this. Is done. Since WL = (Lmax + Lmin) / 2 and WW = Lmax−Lmin, when WL and WW are given, Lmax = WL + WW / 2 and Lmin = WL−WW / 2 are set.

  FIG. 10 is a diagram illustrating an example of the generated MIP image. Each shows an example in which the ROIs are x: 0-511, y: 170-250 (0-511), and z: 0-369, based on a gradation-compressed tomographic image compressed to 8 bits, which will be described later. It is a generated 8-bit gradation MIP image. 10A is a MIP image of MIP (mode = 0), FIG. 10B is a MIP image of MinIP (mode = 1), and FIG. 10C is a MIP image of RaySum (mode = 2) Is shown. In the MIP image of FIG. 10A, a range of 50 to 180 (8 min gradation: Lmin = −406, Lmax = 574, WL = 84, WW = 980) is displayed in 256 steps in 8 bit gradation. It is a mode to do. In the MIP image of FIG. 10B, the range of 0 to 120 in 8-bit gray scale (Lmin = −784, Lmax = 121, WL = −331, WW = 905 in 16-bit gray scale conversion) is set in 256 steps. This is a mode of displaying. In the MIP image of FIG. 10C, the range of 0 to 130 in the 8-bit gradation (Lmin = −784, Lmax = 196, WL = −294, WW = 980 in 16-bit gradation conversion) is set in 256 steps. This is a mode of displaying.

  As described above, in the rendering apparatus 100 according to the present embodiment, the two-dimensional projection image generation unit 70 is realized by the CPU 1 that is a multi-core CPU executing a program. This program divides the MIP image generation region into a plurality of regions, allocates each divided generation region to each core, and performs processing in which each core generates a partial image of the MIP image in parallel. A case where the rendering apparatus 100 is realized by a personal computer which is a general-purpose computer will be described. FIG. 11 is a diagram illustrating a parallelization process when the rendering device 100 is implemented by a computer. As shown in FIG. 11, when the two-dimensional projection image generation unit 70 is realized by an application program, and divides the MIP image into N = 4, for example, the application program performs processing for each of the four divided generation regions by four. It is divided into two processing threads # n1 to # n4. Then, the job scheduler of the multitasking operating system on which the two-dimensional projection image generating means 70 (application program) operates allocates the threads to the CPU threads # 1 to # 8 of the plurality of CPU cores.

  FIG. 12 is a diagram illustrating a flow of software for parallel processing when the rendering apparatus 100 is implemented by a computer. As shown in FIG. 12, the application program for realizing the two-dimensional projection image generation unit 70 divides the processing for each of the four divided generation areas into four processing threads # n1 to # n4, and then executes the parallel processing. Invoke the function. 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. 12, the application program waits for an end message from each CPU thread, and when the end messages are obtained from all the CPU threads, generates an integrated image of the generated divided images and collectively displays the display unit 6. And display it. 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. 12 shows a single thread execution (generation time of a 512 × 512 full-resolution image: about 1.6). (Second), a processing speed of less than about three times (same generation time: about 0.6 sec) was obtained, and some CPU threads had to wait for processing. Therefore, if the number of divisions of the MIP image generation processing is increased to N = 8 and parallel processing is performed by eight processing threads, processing performance almost four times that of the single thread execution can be obtained (same generation time: about 0.4 sec). , N> 8, the processing performance hardly changed and reached 4 times.

  The example shown in FIG. 12 is an example in which display processing is started by integrating generated divided images after waiting for the end of all processing threads. A flow in which display processing is performed without integration may be adopted. FIG. 16 is a diagram showing a flow of software for parallel processing when display processing is performed every time each processing thread ends. As in the example shown in FIG. 12, the application program that implements the two-dimensional projection image generation unit 70 divides the processing for each of the four divided generation regions into four processing threads # n1 to # n4, Invokes a parallel processing function. Then, the job scheduler of the multitasking operating system executes the CPU thread # 3, the CPU thread # 5, the CPU thread # 7, and the CPU thread # 2 having four processing threads # n1 to # n4 of the plurality of CPU cores. Respectively. In the example of FIG. 16, the application program waits for an end message from each CPU thread, and when an end message is obtained from each CPU thread, outputs the divided image to the display unit 6 for each of the generated divided images for display.

<3. Progressive rendering>
Next, a process of performing progressive rendering using the rendering device according to the present embodiment will be described. FIG. 13 is a flowchart illustrating the processing procedure of the progressive rendering. First, initial settings of mouse coordinates and xy rotation angles are performed based on a press operation of a mouse button as an instruction input device (step S810). Specifically, mouse coordinates are set at predetermined positions, and the xy rotation angle (the rotation angle about the x axis, the rotation angle about the y axis, and the rotation angle about the z axis cannot be directly specified, but the xy rotation operation is repeated. Is set to a predetermined angle, and a 3 × 3 matrix defining rotation in a three-dimensional coordinate system is generated. Then, coordinate conversion parameters are set based on the generated 3 × 3 matrix, a MIP image is generated by thinning out 1/4 in each of the x-axis direction and the y-axis direction, and displayed on the display unit 6 (step S820). ). That is, MIP images thinned out by 1/4 in the x-axis direction and the y-axis direction are displayed as initial images. In step S820, a process according to the flowchart of FIG. 5 is executed to generate a MIP image. However, the pixel value calculation processing by rendering is not performed for all the pixels, but is performed for 1/4 pixels in the x-axis direction and the y-axis direction, and 1/16 pixels as a whole. For the remaining 15/16 pixels, pixel values are obtained by interpolating the pixels obtained by rendering. When the processing according to the flowchart of FIG. 5 is executed, the xy rotation angle initially set in step S810 is acquired in step S520, and a 3 × 3 matrix defining the rotation in the three-dimensional coordinate system is generated. , And set the coordinate conversion parameters.

  When a mouse drag operation (movement while pressing a button) is performed, the xy rotation angle is updated based on the amount of movement of the mouse coordinates (step S830). Specifically, when the mouse movement is detected at a predetermined timing, the amount of movement of the acquired mouse coordinates from the previously acquired mouse coordinates is obtained, and the rotation about the y axis is determined based on the amount of movement of the x coordinates. The angle is updated, and the rotation angle about the x axis is updated based on the movement amount of the y coordinate. Then, the 3 × 3 matrix is updated based on the updated xy rotation angle, the coordinate conversion parameters are set based on the updated 3 × 3 matrix, and 1/4 is thinned out in each of the x-axis direction and the y-axis direction. Then, a MIP image is generated and displayed on the display unit 6 (step S840). At this time, since the generation of the MIP image is executed in about 1/16 of the generation time of the full resolution (about 0.4 seconds in the case of a 4-core CPU), the MIP image is thinned out and displayed following the drag operation of the mouse. Since the existing MIP image is updated in real time, from the user's point of view, an instruction to change the angle (rotation angle) is issued by performing a mouse drag operation. The process in step S840 is performed in the same manner as in step S820. The processes of steps S830 and S840 are repeatedly executed while the drag operation is continued with the mouse button pressed. Therefore, while the user continues the drag operation using the mouse, the MIP image thinned to 1/4 in the x-axis direction and the y-axis direction is displayed on the screen of the display unit 6 while changing the angle. Will continue.

  When the release operation of the mouse button is performed, the xy rotation angle is updated based on the movement amount of the mouse coordinates (step S850). Specifically, the amount of movement of the mouse coordinates obtained at the time of releasing the mouse button from the previously obtained mouse coordinates is obtained, and the rotation angle about the y-axis is updated based on the amount of movement of the x-coordinate. The rotation angle about the x-axis is updated based on the movement amount of the coordinates. Then, the 3 × 3 matrix is updated based on the updated xy rotation angle, the coordinate conversion parameters are set based on the updated 3 × 3 matrix, a MIP image is generated without thinning, and the display unit 6 Is displayed (step S860). That is, the process according to the flowchart of FIG. 5 is executed, the pixel values are calculated for all the pixels, and the MIP image is generated and displayed. For this reason, when the mouse drag operation is stopped and the mouse button is released, the screen of the display unit 6 displays a still-resolution full-resolution MIP image that is not thinned out. Therefore, from the user's point of view, by performing the release operation of the mouse button, the angle (rotation angle) is instructed.

  FIG. 14 is a diagram illustrating an example of a screen displayed by progressive rendering. FIG. 14A shows a screen at the time of display in steps S820 and S840, and FIG. 14B shows a screen at the time of display in step S860. At the time of mouse drag operation, since the thinning-out display is performed, the image becomes coarse as shown in FIG. 14A, but images having different rotation angles are sequentially displayed in real time following the mouse movement operation. be able to. Then, when the rotation angle is determined after the mouse drag operation is completed, a high-definition image can be displayed as shown in FIG.

<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 MIP 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, the processing according to the following [Equation 9] is executed to calculate the lower limit Lmin and the upper limit Lmax. The lower limit Lmin and the upper limit Lmax operate in the same manner as the lower limit Lmin and the upper limit Lmax used in [Equation 8], but are fixedly provided unlike in [Equation 8].

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

  In [Equation 9], γ is the 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 brightness contrast of the MIP image can be adjusted on the GUI side by a calculation formula similar to [Equation 10] described later, γ 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 10] is executed to calculate the gradation compressed tomographic image D8 (x, y, z).

[Equation 10]
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 10], when the signal value exceeds the upper limit value Lmax, 255, and when the signal value is lower than the lower limit value Lmin, the value is set to 0, and the range from the lower limit value Lmin to the upper limit value Lmax is linearly changed from 0 to 255. Convert. When the processing based on FIG. 5 is performed on the calculated gradation-compressed tomographic image D8 (x, y, z), a MIP image with 8-bit gradation is generated. Therefore, in principle, display processing can be performed without performing processing according to [Equation 8]. However, in the MinIP or RaySum mode, pixel values are often concentrated below 256 gradations, and it is difficult to display with appropriate contrast. Therefore, as in the case of displaying a 16-bit gradation MIP image, Window Width (window width: WW) and Window Level (window level: WL) are adjusted on the GUI side, and the lower limit Lmin is set in the range of 256 gradations. And the upper limit Lmax are set, and processing according to [Equation 8] is often performed. FIG. 10 described above is an example of generating an 8-bit MIP image for all three. FIG. 10A shows Lmin = 50 and Lmax = 180, FIG. 10B shows Lmin = 0 and Lmax = 120, and FIG. (C) is displayed by setting Lmin = 0 and Lmax = 130, respectively, and performing the processing according to [Equation 8].

  In the above embodiment, the MIP image generation region is divided into a plurality of regions only in the y-axis direction. However, as long as the region can be divided into a plurality of regions, the mode of division is not particularly limited. FIG. 15 is a diagram illustrating various modes when the generation area is divided into a plurality of division generation areas. In FIG. 15, a rectangle indicates an image generation area, and the horizontal direction is the x axis, and the vertical direction is the y axis. FIG. 15 shows an example in which N is set to N = 4 and divided into four divided generation regions R1 to R4. FIG. 15A shows an example in which the image is divided into a plurality of parts only in the y-axis direction. In the above-described embodiment, the description has been made assuming a division mode as shown in FIG. FIG. 15B shows an example in which the image is divided into a plurality of parts only in the x-axis direction. FIG. 15C shows an example in which the image is divided in both the x-axis direction and the y-axis direction. In the examples of FIGS. 15A to 15C, one divided generation region is one continuous region.

  In the examples shown in FIGS. 15A to 15C, as shown in FIG. 15B, it is preferable to divide the image into a plurality only in the x-axis direction. As shown in the above, it is preferable to divide into a plurality only in the y-axis direction. Normally, the CPU determines that the y-axis direction of the screen (display screen) corresponds to the z-axis direction of the voxel image at the front angle, and that each divided generation area accesses the voxel image at a discrete address in the z-axis direction. Become. At this time, if the image is divided only in the x-axis direction, the voxel image is accessed with almost the same z-coordinate for each divided generation area, and the access time is reduced. For this reason, division may be performed only in the x-axis direction where the speed can be increased at the most frequently used front angle. Further, the dividing direction may be switched according to the angle.

  In the example of FIG. 15D, one divided generation area is not a single continuous area, and one divided generation area is separated. Also in the example of FIG. 15D, N = 4 and the image is divided into four divided generation regions R1 to R4, and each divided generation region R1 to R4 is further divided into four blocks. As described above, when processing is performed by dividing one divided generation area into a plurality of blocks, each CPU thread processes a plurality of distant blocks instead of processing one continuous area.

  As described with reference to FIG. 12, when the display process is started after all the processing threads are finished, even if one divided generation area is one continuous area, one divided generation area It does not change even if is divided into multiple blocks. However, as described with reference to FIG. 16, when the display process is performed every time each processing thread ends, a difference occurs in a display mode. When one divided generation area is continuous as in the examples of FIGS. 15A to 15C, display output of an image is performed every time the CPU thread is processed. Will be displayed. For example, when only the processing of the divided generation region R1 is completed first, only the upper part, only the left part, and only the upper left part are displayed first. On the other hand, when one divided generation area is divided into a plurality of blocks as in the example of FIG. 15D, the display is performed in a well-balanced manner as a whole. For example, even if only the processing of the divided generation region R1 is completed first, four separated blocks are displayed.

  Further, it is also possible to adopt a configuration in which the divided sections in the example of FIG. 15D are made finer, and each of the divided generation areas of R1 to R4 is one pixel. That is, this is a method in which four adjacent pixels of 2 × 2 pixels are executed by different processing threads. This method has the advantage that, once the execution of any one of the processing threads is completed, every other thinned image is completed, and the completed image can be generally confirmed during the processing.

  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.

  In the above embodiment, the generation of each divided image is assigned to each CPU core. However, when one CPU core can execute a plurality of processes in parallel, a plurality of divided images are assigned to one CPU core. Generation may be assigned.

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 Coordinate transformation parameter input means 20 ROI clipping setting means 25 Mask setting means 30 Voxel image creating means 70 Two-dimensional projection image generation means 77 ... origin coordinate search means 80 ... image output means 100 ... rendering device

Claims (17)

A rendering apparatus for generating a two-dimensional projection image representing a state in which a plurality of two-dimensional tomographic images captured at predetermined intervals of a predetermined target object are observed from a predetermined line-of-sight direction,
Dividing a generation area, which is an area for generating the two-dimensional projection image, into N divided generation areas;
From each pixel constituting the divided generation area, irradiate a virtual ray in the line of sight direction,
A plurality of voxels of the voxel image in which the virtual rays intersect with a voxel image corresponding to each pixel constituting the plurality of two-dimensional tomographic images are extracted, and a single voxel is extracted based on the extracted signal values of the plurality of voxels. By calculating the representative signal values of the above and giving the calculated representative signal values as the pixel values of the respective pixels constituting the divided generation area irradiated with the virtual ray, a divided image corresponding to each divided generation area is generated. A two-dimensional projection image generating means,
A rendering apparatus for generating the two-dimensional projection image, wherein the two-dimensional projection image generation means executes a process of generating a divided image corresponding to the N divided generation regions in parallel.   2. The rendering according to claim 1, wherein the two-dimensional projection image generation unit is configured to calculate any one of a maximum value, a minimum value, and an average value of the signal values of the extracted voxels as the representative signal value. 3. apparatus. The two-dimensional projection image generating means includes:
When each of the divided generation areas is defined in the viewpoint coordinate system as an area that intersects a predetermined line-of-sight direction,
In the viewpoint coordinate system, from each pixel of the divided generation area, obtain signal values of a plurality of voxels corresponding to the voxel image at predetermined reference intervals in the line of sight direction, and obtain signal values of the obtained voxels. The rendering device according to claim 1, wherein the single representative signal value is calculated based on the value. The two-dimensional projection image generating means includes:
In a predetermined pixel of the divided generation area, an upper limit value of the line of sight direction is set as a search start coordinate, and a starting point coordinate search unit is searched with a coordinate value of a maximum among the line of sight directions included in the effective range as a start point coordinate. ,
While reducing the coordinates in the line-of-sight direction by a predetermined number m1 toward the lower limit value from the start point coordinates calculated by the start point coordinate search means, predetermined coordinate conversion is performed for each three-dimensional coordinate of the viewpoint coordinate system, and the voxel image The rendering device according to claim 3, wherein a signal value of a corresponding voxel is calculated. The two-dimensional projection image generating means includes:
For a predetermined pixel of the divided generation area, the signal value of the voxel calculated at the coordinates of the line of sight is included in the effective range, and the signal value of the voxel is calculated at the coordinates of the line of sight after subtracting the predetermined number m1. Voxel signal value is not included in the effective range,
5. The rendering device according to claim 4, wherein the starting point coordinate search unit searches for starting point coordinates included in the effective range, using coordinates in the line of sight direction after subtracting the predetermined number m1 as search start coordinates. 6. . The starting point coordinate search means includes:
The coordinates in the direction of the line of sight are subtracted by a predetermined number m2 (m2> m1) greater than the predetermined number m1 in a range of a predetermined lower limit from a predetermined upper limit, and the coordinate conversion is performed on each three-dimensional coordinate. To calculate the voxel signal value,
If the signal value of the voxel at the coordinates in the line-of-sight direction after the subtraction of the predetermined number m2 is included in the effective range, the three-dimensional coordinates are increased by increasing the predetermined number m1 from the coordinates in the line-of-sight direction after the subtraction. The rendering apparatus according to claim 4, wherein the coordinate conversion is performed on the image data to calculate a voxel signal value.   When the signal value of the voxel corresponding to the three-dimensional coordinates in the coordinates of the gaze direction after the increase of the predetermined number m1 is not included in the effective range, m1 is calculated from the coordinates of the gaze direction after the increase of the predetermined number m1. 7. The rendering apparatus according to claim 6, wherein a value obtained by subtracting the first value is given as the starting point coordinates.   If none of the voxel signal values at the coordinates in the line-of-sight direction after the subtraction of the predetermined number m2 is included in the effective range, a value smaller than the lower limit value of the coordinates in the line-of-sight direction is given as the starting point coordinate. The rendering device according to claim 6 or 7, wherein: The starting point coordinate search means includes:
Calculate the voxel coordinates of the real value by performing the coordinate transformation, to obtain the voxel coordinates of the integer value by rounding each value of the calculated voxel coordinates of the real value,
When the voxel coordinates of the obtained integer value are included in a predetermined effective area, a pixel value of the voxel image corresponding to the voxel coordinates is given as a signal value of a voxel corresponding to the three-dimensional coordinates,
9. The rendering apparatus according to claim 8, wherein when the voxel coordinates are not included in the effective area, a value not included in the effective range is given as a voxel signal value corresponding to the three-dimensional coordinates. The two-dimensional projection image generating means includes:
The coordinate conversion is performed to calculate real-valued voxel coordinates, and a plurality of integer-valued voxel coordinates located near the calculated real-valued voxel coordinates are obtained,
Among the obtained voxel coordinates of the plurality of integer values, the plurality of voxels of the voxel image corresponding to the voxel coordinates of the integer values included in the predetermined effective area are specified, and the signal values of the specified plurality of voxels are specified. The rendering device according to claim 4, wherein a value added while multiplying by a predetermined weight is given as a voxel signal value corresponding to the three-dimensional coordinates. The coordinates corresponding to the left end in the x-axis direction of each tomographic image constituting the voxel image are Xso, the coordinates corresponding to the right end in the x-axis direction are Xeo, the coordinates corresponding to the top end in the y-axis direction are Yso, If the coordinates corresponding to the lower end are Yeo, the z coordinate at which the first tomographic image is arranged is Zso, and the z coordinate at which the final tomographic image is arranged is Zeo, Xso ≦ Xs <Xe ≦ Xeo, Yso ≦ Ys <Ye ≦ When six coordinates Xs, Xe, Ys, Ye, Zs, Ze satisfying Yeo, Zso ≦ Zs <Ze ≦ Zeo are defined,
The rendering device according to claim 9, wherein the effective area is an area surrounded by the six coordinates Xs, Xe, Ys, Ye, Zs, and Ze.   10. The effective area is Mask = 1 when mask data Mask having a value of 0 (not displayed) or 1 (displayed) corresponding to the voxel image is defined in advance. The rendering device according to any one of claims 11 to 13. Coordinate conversion parameter input means for inputting the parameters for the coordinate conversion, and image output means for displaying and outputting a two-dimensional projection image generated by the two-dimensional projection image generation means,
If a predetermined change instruction is issued from outside,
By the coordinate conversion parameter input means, input of a coordinate conversion parameter,
The pixels of the two-dimensional projection image are subjected to the two-dimensional projection image generation means while the pixels of the two-dimensional projection image are thinned out at predetermined intervals. The image output means displays a two-dimensional projection image generated by interpolating based on the
When the predetermined confirmation instruction is given,
By the coordinate conversion parameter input means, input of a coordinate conversion parameter,
13. The image output means displays the two-dimensional projection image generated by the two-dimensional projection image generation means being executed for all the pixels of the two-dimensional projection image. The rendering device according to claim 1. When each pixel value of the two-dimensional projection image generated by the two-dimensional projection image generation means has a gradation exceeding 8 bits, a predefined lower limit value Lmin and an upper limit value Lmax of the gradation are used. ,
When the pixel value Image (x, y) satisfies Lmin ≦ Image (x, y) ≦ Lmax, the pixel value is calculated as (Image (x, y) ′ = (Image (x, y) −Lmin) · 255 / ( Lmax (x, y) -Lmin),
When the pixel value Image (x, y) is Image (x, y)> Lmax, the pixel value is set to Image ′ (x, y) = 255,
When the pixel value Image (x, y) is Image (x, y) <Lmin, the pixel value is set to Image ′ (x, y) = 0,
14. The rendering device according to claim 13, wherein each of the pixel values is converted into 8-bit gradation and displayed. When the generation region is a two-dimensional region in the x-axis direction and the y-axis direction,
The N divided generation areas are divided only in the x-axis direction or only in the y-axis direction,
The rendering device according to any one of claims 1 to 14, wherein each of the divided generation areas is one continuous area.   The rendering according to any one of claims 1 to 14, wherein when the generation area is a two-dimensional area in the x-axis direction and the y-axis direction, each of the divided generation areas is a plurality of separated blocks. apparatus. In order to generate a two-dimensional projected image representing a state in which a plurality of two-dimensional tomographic images taken at predetermined intervals of a predetermined target object are observed from a predetermined line of sight,
Computer
The generation area, which is an area for generating the two-dimensional projection image, is divided into N divided generation areas, and a virtual ray is irradiated from each pixel of the divided generation area in the line of sight direction, and the plurality of two-dimensional tomographic images are obtained. A plurality of voxels of the voxel image and the voxel image where the virtual ray intersects each pixel constituting the pixel are extracted, and a single representative signal value is calculated based on the signal values of the extracted voxels, By providing the calculated representative signal value as a pixel value of each pixel constituting the divided generation region irradiated with the virtual ray, a two-dimensional projection image generating unit that generates a divided image corresponding to each divided generation region Let it work,
A program for generating the two-dimensional projection image by executing, in parallel, a process of generating a divided image corresponding to the N divided generation regions.