JP7167699B2 - volume rendering device - Google Patents

Description

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

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

On the other hand, in the field of 3D computer graphics, it is practiced to set a virtual light source in a predetermined direction and add shadows for more three-dimensional expression (see Patent Documents 1 to 3). Patent Document 1 proposes a method of extracting the boundary surface of the organ region based on the voxel opacity gradient and adding shadows to the boundary surface based on the normal vector of the boundary surface in order to provide realistic shading. It is

Further, Japanese Patent Application Laid-Open No. 2002-200001 proposes a method of calculating a gradient from the gradient of the signal value of each voxel and directly applying a shadow to each voxel. Further, Japanese Patent Application Laid-Open No. 2002-200001 proposes a method of calculating a gradient from the gradient of the signal value of each voxel and directly applying a shadow to each voxel.

Japanese Patent No. 5065740 Japanese Patent No. 3542799 Japanese Patent No. 2627607

However, the technique described in Patent Document 1 has the advantage of providing an image that is easy to read because the shadow is limited to the boundary surface compared to the method of directly applying shadows to each voxel. There is a problem that the load is large and specular reflection cannot be supported. In addition, in the techniques described in Patent Documents 2 and 3, gradient calculation is not performed based on opacity but directly based on signal values, so there is a problem that shadows cannot be controlled with a color map and specular reflection cannot be supported. .

Therefore, an object of the present disclosure is to provide a volume rendering apparatus capable of adding shadows while reducing the computational load when generating a rendering image by applying a color map to a plurality of tomographic images. do.

The present disclosure includes multiple means for solving the above problems, but if one example is given,
A color map defined by associating color component values and opacity with signal values based on a voxel image composed of a plurality of two-dimensional tomographic images taken at predetermined intervals of a predetermined target. A volume rendering device for generating a rendered image with reference to
rendering means for generating a rendered image using the voxel image and the color map;
The rendering means refers to a plurality of voxels of the voxel image while performing a predetermined coordinate transformation on the voxel image for each pixel of the rendering image, and refers to the color map for each voxel to be referred to. a ray casting means for obtaining the color component and opacity of the voxel while obtaining the opacity of a plurality of neighboring voxels of the voxel;
calculating a gradient vector in the voxel based on the opacity of a plurality of voxels near the voxel obtained by the ray casting means, and calculating a regular reflection vector of a predetermined light source vector with respect to the calculated gradient vector; calculating the diffusely reflected light of the voxel based on the inner product value of the gradient vector and the light source vector, and calculating the specular reflected light of the voxel based on the inner product value of the regular reflection vector and a predetermined line-of-sight vector; and a shadow adding means for updating the color component of the voxel by multiplying the color component of the voxel obtained by the ray casting means by a shadow value calculated based on the calculated diffuse reflection light and the specular reflection light. .

According to the present disclosure, when applying a color map to a plurality of tomographic images to generate a rendering image, it is possible to add shadows while reducing the computational load.

FIG. 3 is a diagram showing a plurality of tomographic images and a color rendering image generated as a volume rendering image; 1 is a hardware configuration diagram of a volume rendering device 100 according to an embodiment of the present disclosure; FIG. 1 is a functional block diagram showing the configuration of a volume rendering device according to an embodiment; FIG. 4 is a flowchart showing processing operations of the volume rendering device according to the embodiment; It is a figure which shows the color map used by this embodiment. FIG. 5 is a flow chart showing a case where opaque circumscribed area calculation processing in step S30 shown in FIG. 4 is performed in parallel by a plurality of CPU cores. FIG. FIG. 7 is a flow chart showing calculation processing of an opaque circumscribed region for a group of tomographic images in step S330 of FIG. 6. FIG. FIG. 7 is a flow chart showing calculation processing of an opaque circumscribed region of the entire group of tomographic images in step S360 of FIG. 6. FIG. FIG. 11 compares ROIs and opaque circumscribing regions; FIG. 5 is a flow chart showing a process in which each CPU core performs rendering processing in step S40 shown in FIG. 4 in parallel; FIG. 4 is a conceptual diagram of coordinate transformation; FIG. FIG. 11 is a flow chart showing details of a divided image generation process in step S440 shown in FIG. 10; FIG. FIG. 13 is a flowchart showing details of voxel reference processing in step S640 shown in FIG. 12; FIG. FIG. 14 is a flow chart showing the details of search processing for starting point coordinates in step S642 shown in FIG. 13; FIG. FIG. 4 is a diagram showing parallel processing when the volume rendering device 100 is implemented by a computer; FIG. 4 is a diagram showing a software flow of parallel processing when the volume rendering device 100 is implemented by a computer.

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

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

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

FIG. 3 is a functional block diagram showing the configuration of the volume rendering device according to this embodiment. 3, 10 is tomographic image reading means, 20 is color map reading means, 30 is ROI clipping setting means, 40 is mask setting means, 50 is opaque circumscribed area calculation means, 60 is rendering means, 61 is ray casting means, Reference numeral 62 denotes starting point coordinate searching means, 63 denotes shading means, and 70 denotes rendering image output means.

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

The rendering means 60 is means for setting a rendering image generation area, determining pixel values obtained by projecting a voxel image composed of a plurality of tomographic images with reference to a color map, and generating a rendering image. The rendering means 60 further comprises ray casting means 61 and shading means 63 . The ray casting means 61 calculates, for each pixel of the rendered image, an intersection point between an effective area, which will be described later, and a vector in the line-of-sight direction closest to the viewpoint. This is means for irradiating a virtual ray in the line-of-sight direction and referring to a plurality of voxels of a voxel image while performing a predetermined coordinate transformation based on the ray casting method. The shading means 63 applies the color component (R, G, B) of each voxel converted by referring to the color map based on the signal value of the voxel image referred to, to the neighborhood of the voxel. Based on the signal value of the voxel of , based on the opacity α converted by similarly referring to the color map, calculate the gradient vector at that voxel, use the calculated gradient vector to calculate the shadow value, and color This is means for replacing each value of the components (R, G, B) with a value obtained by multiplying the shadow value. Furthermore, the ray casting means 61 has a starting point coordinate searching means 62 . The starting point coordinate searching means 62 is a means for searching for the coordinates of the starting point of projection in the viewpoint coordinate system when the ray casting means 61 calculates pixel values.

The tomographic image reading means 10 and the color map reading means 20 are realized mainly in the data input/output I/F 5 with the assistance of the CPU 1 . The ROI clipping setting means 30 and the mask setting means 40 are realized mainly in the instruction input I/F 4 with the assistance of the CPU 1 . The opaque circumscribed area calculation means 50 , ray casting means 61 , shadow addition means 63 , rendering means 60 and origin coordinate search means 62 are implemented by the CPU 1 executing programs stored in the storage device 3 . The rendering image output means 70 is realized mainly in the frame memory 8 and the display section 6 with the assistance of the CPU 1 .

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

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

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

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

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

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

Next, the opaque circumscribed area calculation means 50 performs opaque circumscribed area calculation processing (step S30). An opaque circumscribing area is an area that circumscribes a range in which pixels whose opacity is not 0 exists in each three-dimensional axial direction. The circumscribing area is defined as an area surrounded by a predetermined side (line segment) circumscribing the range in which pixels with non-zero opacity exist. This is because the range intersecting with can be calculated at high speed by calculation. Therefore, the opaque circumscribing area is limited to a rectangular parallelepiped area perpendicular to the three-dimensional coordinate axis that circumscribes the range in which pixels whose opacity is not 0 exist. (With an oblique rectangular parallelepiped or a curved surface, the calculation of the circumscribed surface requires a processing load, and the intersection with the line of sight cannot be analytically calculated using a linear formula, so the desired speed-up effect cannot be achieved.) It is preferable that the region is an internal region of a rectangular parallelepiped defined by the direction of the coordinate axis and the side along the direction in which the tomographic images are stacked. As a specific process, using the color map read from the color map reading means 20, the opacity corresponding to the signal value of each pixel included in the ROI of the tomographic image after clipping setting is applied. Calculate the opaque bounding area based on the opacity.

FIG. 5 is a diagram showing a color map used in this embodiment. In FIG. 5, the signal value is the signal value recorded in each pixel of the tomographic image, R, G, and B are the color components of red, green, and blue, respectively, and α is the opacity of the voxel. . FIG. 5 is a color map when signal values of a tomographic image are recorded in 16 bits and take values from -32768 to 32767. FIG. In the case of a CT image, it is often normalized to a range from -2048 to 2048 with the water component set to 0. As shown in FIG. 5, in the color map, the signal value of each pixel of the tomographic image is associated with each of the R, G, and B color components and the opacity α. This color map is also used in rendering processing by the rendering means 60 (described later: step S40). In the process of calculating the opaque circumscribed area in step S30, only the opacity α of the color map is used. By using the opacity α obtained by referring to the color map, it is possible to calculate, as an opaque circumscribed region, an area in which there are opaque pixels where α=0 in the tomographic image. The area corresponding to this opaque bounding area is defined as the valid area. As the "area corresponding to the opaque circumscribing area", the opaque circumscribing area can be used as it is as an effective area to be processed. is the effective area corresponding to the opaque circumscribed area.

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

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

In this embodiment, it is possible to change the color map in use or load a different color map by operating via the instruction input I/F 4 . The volume rendering device then refers to the changed color map or the newly loaded color map to generate a rendered image to be presented as a volume rendered image. The volume rendering apparatus of this embodiment does not refer to a color map from a tomographic image to generate a voxel image in the RGBα format. Since the conversion to color components and opacity is performed while referring to the color map for limited voxels, there is no need to reconstruct the voxel structure in RGBα format for changes in the color map. can update the volume rendering image.

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

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

After completing the calculation of the opaque circumscribing region for the tomographic image group in step S330, the opaque circumscribing region is calculated for the entire tomographic image consisting of the NA tomographic image groups (step S360). Details of the processing in step S360 will also be described later.

In this embodiment, a plurality of divided tomographic image groups are processed in parallel by a plurality of CPU cores. In this embodiment, z=Z1 is the head and z=Z2 is the tail, but z=Z1 may be the tail and z=Z2 may be the head. Ultimately, the head and tail should be set so that all read tomographic images can be processed from one to the other without changing the layered structure.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Next, the target tomographic image group is changed (step S363). Specifically, the value of tid is incremented. If tid is smaller than NA as a result of incrementing the value of tid, it means that processing has not been completed for all tomographic image groups. Perform processing to change variables for determining values and maximum values. If tid is greater than or equal to NA as a result of incrementing the value of tid, it means that the processing has been completed for all tomographic image groups, and boundary processing of the opaque circumscribed area is performed (step S364). In step S364, as a countermeasure against moire, a boundary plane (x=Xis or x=Xie or y=Yis or y=Yie or z=Zis or z=Zie) composed of transparent pixels with an opacity α of 0 is made opaque. One pixel width is added to the edge of the circumscribed rectangular parallelepiped that is the outermost circumference of the circumscribed area. Specifically, a process of expanding the calculated area by repeating step S362 by one pixel in each of the xyz axis directions is performed. More specifically, the value of the variable is changed by executing the processing according to [Formula 4] below.

[Formula 4]
If Xis=Xmin−1 and Xis<Xs, then Xis=Xs
When Xie=Xmax+1 and Xie>Xe, Xie=Xe
If Yis=Ymin−1 and Yis<Ys, then Yis=Ys
If Yie=Ymax+1 and Yie>Ye, then Yie=Ye
If Zis=Zmin−1 and Zis<Zs, then Zis=Zs
If Zie=Zmax+1 and Zie>Ze, then Zie=Ze

In step S364, the boundary pixels are moved outward by one pixel. The area expanded outward by one pixel in this manner is used as an effective area to be processed. Since the pixels on the outermost edge of the effective area should have an opacity α of "0", the pixels on the outermost periphery of the effective area, so-called pixels located on the boundary surface of the effective area, are all transparent pixels. Therefore, the opaque circumscribed area calculation means 50 calculates the effective area so that the opacity corresponding to the outermost pixels in each axial direction is zero. However, the color components (RGB values) of the pixels located on the boundary surface of the effective area are the signals of the tomographic image regardless of whether the mask data Mask (x, y, z) is 0 or outside the range of the ROI. Based on the value, the RGB value defined by the color map is given as it is. As a result, the opacity becomes discontinuous inside and outside the effective area, but the continuity of the RGB values is maintained. It is possible to suppress moiré that occurs on the surface.

A range defined by rectangular parallelepipeds of [Xis, Xie], [Yis, Yie], and [Zis, Zie] obtained as a result of the processing in step S364 is defined as an effective region and specified as a processing target. In step S364 of FIG. 8, by setting the area expanded by one pixel outside the opaque circumscribed area as the effective area, the opacity corresponding to the outermost voxel (pixel) in each axial direction of the effective area is is set to 0. However, it is not limited to this, and the area defined by the opaque circumscribed area is used as it is as an effective area. may be replaced with 0 so that the opacity corresponding to the outermost voxel becomes 0.

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

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

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

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

In this way, even if volume rendering is performed by defining an effective area reduced to the range corresponding to the opaque circumscribed area, the rendered image obtained is the same as the conventional volume rendering using the entire tomographic image (voxel image) as the effective area. There is no difference from the case of doing. However, it is possible to quickly obtain the starting point coordinates in the effective area on the virtual ray that is searched for each pixel of the rendering image described later.

After completing the opaque circumscribed area calculation processing in step S30, the rendering means 60 performs rendering processing (step S40). FIG. 10 is a flow chart showing the process of performing the rendering process in step S40 in parallel by each CPU core.

Rendering images are generated by dividing a generation region determined by a set size into a plurality of regions, and processing the obtained divided generation regions in parallel to generate divided images. The size of the rendered image is the same in both the x-axis direction and the y-axis direction, and is set to the number of pixels Size×Size. First, the rendering means 60 initializes divided generation areas to be divided for parallel processing (step S410). In this embodiment, the y-coordinate is not divided, but the x-coordinate is divided into NB pieces. Specifically, the start coordinate of the x coordinate of each divided generation area is set to X1, and the end coordinate is set to X2. Allocate the generation area. As initial values, X1=0 and X2=Size/NB-1 are set.

Next, coordinate transformation parameters to be generated are set (step S420). The rendering means 60 sets coordinate transformation parameters that are commonly used when performing coordinate transformation in subsequent steps. The setting of the coordinate transformation parameters is performed by setting the three-dimensional coordinates in the viewpoint coordinate system based on the viewpoint so that the sequential coordinate transformation processing can be performed smoothly for each voxel when irradiating the voxel with a virtual ray. is converted into three-dimensional coordinates in the voxel coordinate system in which the voxel image is defined based on , and the signal values of the corresponding voxels in the voxel coordinate system are acquired.

Here, a conceptual diagram of coordinate transformation and specific examples of coordinate transformation parameters will be described. FIG. 11 is a conceptual diagram of coordinate transformation. Since a rendered image obtained by rendering is an image viewed from a viewpoint, the viewpoint coordinate system that serves as a reference for each pixel (x, y) of the rendered image does not necessarily match the voxel coordinate system. Therefore, the rendering means 60 needs to convert the voxel image defined in the voxel coordinate system to the viewpoint coordinate system. That is, originally, as shown in FIG. 11(a), a virtual ray is emitted from each point on the projection surface to the voxel image in an arbitrary direction, and the specified range in the voxel image where the virtual ray intersects A plurality of signal values are extracted at , and the pixel values of the rendering image are 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 the pixel values in the specified range in the negative z-axis direction from each point on the projection plane are calculated as shown in FIG. 11(b). However, the coordinate-transformed voxel image as shown in FIG. 11(b) is not actually created separately, and the pixel value of the voxel image is acquired each time while performing the coordinate transformation in the opposite direction for each voxel. I'm trying

Various coordinate transformation parameters are used when performing the coordinate transformation. If the positional relationship between the viewpoint coordinate system and the voxel coordinate system is the same, it is considered that the coordinate transformation parameters are also the same. That is, in the rendering means 60, coordinate transformation parameters for carrying out coordinate transformation to the voxel coordinate system are the same for all voxels referred to by the viewpoint coordinate system. Therefore, the rendering means 60 calculates coordinate transformation parameters in advance in step S420, and uses the coordinate transformation parameters as common coordinate transformation parameters in subsequent steps. Examples of coordinate transformation parameters include the following.

・3×3 matrix defining rotation in a predetermined three-dimensional coordinate system (the inverse matrix of the 3×3 matrix specified and defined by the instruction input I/F 4 because projection transformation is performed in the opposite direction)
・Offset value in xyz axis direction: Xoff, Yoff, Zoff (unit: voxel)
・Area set by ROI clipping: x-axis direction: Xs to Xe, y-axis direction: Ys to Ye, z-axis direction: Zs to Ze
・Opaque circumscribed area: x-axis direction: Xis to Xie, y-axis direction: Yis to Yie, z-axis direction: Zis to Zie
・Enlargement or reduction ratio: Scale (use the same value for each xyz axis)
z-axis direction scaling factor: Scz=Rxy/Rz (ratio between xy-axis direction resolution Rxy and z-axis direction resolution Rz)
・Coordinate transformation sub-sampling ・Small offset value: x-axis direction dx, y-axis direction dy, z-axis direction dz
・ Virtual ray subsampling magnification: Sray (usually 1, the larger the value, the coarser the value, the less than 1, the finer)

After the coordinate transformation parameters are set, the coordinate transformation sub-sampling is offset (step S430). Specifically, dx=dy=dz=0 is set. Also, set l=0. Next, division image generation processing is performed for the division generation area (step S440). As described above, the rendering image generation area is divided into NB pieces only on the x-coordinate. X1+1). Therefore, in step S440, a divided image of Size.times.(X2-X1+1) pixels is generated. The processing of step S440 will be described later.

Subsequently, the coordinate transformation sub-samples are updated (step S450). Specifically, execute l←l+1, dx←dx+1/L, dy←dy+1/L, dz←dz+1/L. The process of step S440 is performed L times while updating the coordinate transformation sub-samples.

The process of step S440 is performed in parallel by a plurality of CPU cores. Therefore, after assigning the divided image generation processing for the first divided range to one CPU core and starting the execution processing, the next divided generation area is set and the generation processing is started without waiting for the end of the processing (step S460). ). Specifically, the values of X1 and X2 are each increased by Size/NB, and the leading coordinate X1 and the trailing coordinate X2 of the x coordinate in the next divided generation area are newly set.

After setting a new segmented region in step S460, if the new segmented region exceeds the size of the entire rendered image, that is, if X2>Size-1 is satisfied here, all segmented regions , the process of step S40 shown in FIG. 10 is terminated. If X2≦Size−1 is satisfied, the process has not been completed for all divided generation areas, so the process returns to step S420 to allocate another CPU core and repeat the process of activating divided image generation.

The divided image generation processing in step S440 of FIG. 10 will be described. FIG. 12 is a flow chart showing details of the divided image generation processing in step S440. Here, the divided images of the rendered image to be generated are defined. The rendered image is an image having 24-bit pixel values in total of 8 bits for each color of RGB at pixels specified by two-dimensional coordinates (x, y) (X1≤x≤X2, 0≤y≤Size-1). Then, the background color is first set for this rendered image (step S610). Specifically, the background color bg(c), c=0(R), 1(G), 2(B), which takes values from 0 to 255 in 8 bits, is defined. Let x=X1 and y=0 be the pixel of the target rendering image.

Next, the virtual ray intensity and the cumulative brightness value are initialized (step S620). Specifically, the virtual light intensity Trans=1.0 and the cumulative luminance value Energy(c)=0.0 are initialized. The cumulative luminance value Energy(c) is set for each c=0, 1, 2 corresponding to R, G, B, respectively.

Next, the z coordinate of the intersection point between the effective area and the virtual ray is calculated (step S630). Specifically, with respect to the target pixel (x, y) on the rendering image that is the target of the ray casting process, in the range from z = Zo (=Size/Sray-1) to z = 0 in the viewpoint coordinate system, An intersection point Zc with the effective area ([Xis, Xie], [Yis, Yie], [Zis, Zie]) closest to the viewpoint z=Zo is calculated.

In calculating the intersection point Zc, first, a virtual ray vector is calculated. Specifically, coordinate transformation is performed on (x, y, Zo) and (x, y, 0) of the viewpoint coordinate system using the above coordinate transformation parameters, and real-valued voxel coordinates are converted to (x1, y1, z1), (x2, y2, z2), and each element vx=x2-x1, vy=y2-y1, vz=z2-z1 of the virtual ray vector is calculated.

Next, the points of intersection between the virtual ray vector and the six surfaces of the effective area are calculated. Specifically, a process according to [Formula 6] below is executed.

[Formula 6]
Initialize to tx = ty = tz = 10 If |vx|≧1, calculate tx1=(Xis-x1)/vx, tx2=(Xie-x1)/vx, and set the smaller one of tx1 and tx2 to tx set to
If |vy|≧1, ty1=(Yis−y1)/vy and ty2=(Yie−y1)/vy are calculated, and the smaller one of ty1 and ty2 is set as ty.
If |vz|≧1, calculate tz1=(Zis−z1)/vz, tz2=(Zie−z1)/vz, and set the smaller one of tz1 and tz2 to tz.

tx, ty, and tz are obtained by processing according to [Equation 6]. Subsequently, the intersection point Zc closest to the viewpoint (z=Zo) is determined. Specifically, a process according to [Formula 7] below is executed.

[Formula 7]
Initialize to Zc = -1 and execute the following 1) to 3) 1) If tx ≤ 1, calculate Y = vy tx + y1 and Z = vz tx + z1, Yis ≤ Y ≤ Yie and Zis ≤ Z ≤ Zie Then, Z is the z-coordinate of the intersection point in the voxel coordinate system, so Zv = Zo (1-tx) is calculated with the z-coordinate of the intersection point in the viewpoint coordinate system as Zv. ) If ty ≤ 1, calculate X = vx ty + x1, Z = vz ty + z1, if Xis ≤ X ≤ Xie and Zis ≤ Z ≤ Zie, calculate Zv = Zo (1-ty), Zv > Zc Then replace with Zc=Zv 3) If tz≦1, X=vx tz+x1 and Y=vy tz+y1 are calculated. If Xis≦X≦Xie and Yis≦Y≦Yie, Zv=Zo (1−tz ) is calculated, and if Zv>Zc, replace with Zc=Zv

In [Equation 7], 1) to 3) may be executed more than once or may not be executed at all. The intersection point Zc is obtained by executing the processing according to [Equation 7]. However, when the obtained result is in front of the viewpoint (z=Zo) (Zc>Zo), the position of the viewpoint is corrected to Zc=Zo and the position of the viewpoint becomes the intersection point. In step S630, the ray casting means 61 calculates the intersection point Zc between the effective area and the line-of-sight direction vector closest to the viewpoint for each pixel of the rendering image.

After the z-coordinate Zc of the intersection point between the effective area and the virtual ray is calculated by the process of step S630, the voxel reference process is executed (step S640). In the voxel reference process, the signal value of the referenced voxel is converted to opacity while referring to the color map, and the opaque voxel is searched for. Convert to RGB values while referring to it, and calculate the virtual light intensity Trans and the cumulative luminance value Energy(c) of each color of RGB (c=0, 1, 2). The processing of step S640 will be described later.

After calculating the virtual light intensity Trans and the cumulative luminance value Energy(c) of each color of RGB (c=0, 1, 2), the pixel value of the pixel (x, y) of the rendering image is determined (step S650). Specifically, a process according to [Formula 8] below is executed.

[Formula 8]
Image (x, y, c) ← Image (x, y, c) + K · {Energy (c) + Trans · bg (c)} / L

In [Equation 8], K is an intensity magnification, and in this embodiment, K=1.0. After the pixel value of the pixel (x, y) is determined in this manner, the target pixel of the rendered image is moved in the x-axis direction by one (step S660), and then the processing of steps S620 to S650 is repeated. Execute and calculate the pixel value for the next pixel (x, y). When the upper limit X2 in the x-axis direction is reached, the target pixel of the rendering image is moved to the initial value in the x-axis direction (step S670). At this time, x=X1. Then, the processing of steps S620 to S650 is repeatedly executed to calculate the pixel value for the next pixel (x, y). After calculating the pixel values for all the pixels in the range of X1≦x≦X2 and 0≦y≦Size−1 in this way, the ray casting process in step S440 shown in FIG. 12 ends.

Here, the voxel reference processing in step S640 will be described. FIG. 13 is a flowchart showing the details of the voxel reference process in step S640. First, search start coordinates are set (step S641). Specifically, the search start z-coordinate Zst=Zc is set. Next, a search is performed using the z-coordinate of the first opaque voxel from the search start coordinate Zst as the starting coordinate (step S642).

The origin coordinate search processing in step S642 will be described using FIG. FIG. 14 is a flow chart showing the process by which the starting point coordinate searching means 62 searches for the starting point coordinates in each pixel of the rendered image. The origin coordinate searching means 62 executes the processing according to the flowchart shown in FIG. 14 for each target pixel (x, y).

The starting point coordinate searching means 62 sets the target pixel (x, y) for which the starting point coordinates are to be searched, and the three-dimensional coordinates (x, y, z) of the viewpoint coordinate system (step S901). The initial value of the z coordinate, which is the search start coordinate, is set to z=Zst (upper limit). The z-coordinate has a lower limit of 0 and an upper limit of Zst. Zst, which is the upper limit, is the z-coordinate of the viewpoint.

The starting point coordinate searching means 62 determines whether or not the corresponding voxel of the voxel image by performing coordinate transformation on the three-dimensional coordinates (x, y, z) is transparent (step S902). Coordinate transformation is an operation for transforming the viewpoint coordinate system into the voxel coordinate system, which is the reverse of the transformation operation on the GUI side. On the GUI side, it is assumed that clipping by ROI, scaling, z-axis direction scaling processing, rotation, and offset (simultaneously in the xyz-axis direction) are performed in this order, and the three-dimensional coordinates (x, y, z) of the given viewpoint coordinate system are assumed to be performed. Real number coordinates (xr, yr, zr) corresponding to voxels of a voxel image composed of a plurality of DICOM-format tomographic images Do (x, y, z) corresponding to (integer values) are calculated.

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 image of Sx×Sy×Sz. Since the size Sx×Sy of a DICOM-format tomographic image is usually Sx=Sy=512, Size=512 is often used.

Considering that the z-axis direction of the viewpoint coordinate system is enlarged by the virtual ray subsample (1/Sray) times (Sray≤1), first, the processing according to [Equation 9] below is executed. offset processing is performed simultaneously in the xyz-axis directions.

[Formula 9]
xx=x−Size/2+dx+Xoff
yy=y−Size/2+dy+Yoff
zz=z.Sray-Size/2+dz+Zoff

Subsequently, using the generated 3×3 rotation inverse matrix, the rotation processing is performed by executing the processing according to [Equation 10] below.

[Formula 10]
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 11] is executed to perform scaling and z-axis direction variable magnification processing at the same time.

[Formula 11]
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. The coordinates (xr, yr, zr) are real values. Next, the origin coordinate searching means 62 converts the calculated coordinates (xr, yr, zr) of the voxel coordinate system into integer values. Specifically, a process according to [Formula 12] below is executed to convert to an integer value.

[Formula 12]
xi=INT[xr+0.5]
yi=INT[yr+0.5]
zi=INT[zr+0.5]

In [Numerical formula 12], INT means rounding off the decimal point of the numerical value in [ ] to make it an integer. Therefore, by the processing according to [Equation 12], each value of the real-valued coordinates (xr, yr, zr) is rounded off to obtain the integer-valued coordinates (xi, yi, zi).

When the coordinates (xi, yi, zi) of the voxel coordinate system corresponding to the coordinates (x, y, z) of the viewpoint coordinate system are thus obtained, the effective area, the mask data Mask, and the opacity correction table Sα is used to execute the processing according to the following [Equation 13] to determine whether the voxel corresponding to the coordinates (xi, yi, zi) in the voxel coordinate system is transparent, opaque, or outside the effective area. . The opacity correction table Sα(x, y, z) is a multidimensional transfer function in which a correction magnification for opacity is defined in advance in a three-dimensional space. The same function sets the correction magnification to the maximum in the central area where the internal organs exist, decreases the correction magnification toward the periphery by a predetermined function, Outside the body: bed, headrest/fixing jig) is defined so that the correction magnification is 0. When volume rendering is performed using this, a volume rendering image in which the outside of the bone region becomes transparent and only the visceral region is exposed can be obtained. If the mask data Mask and the opacity correction table Sα are not defined, Mask(x, y, z)=Sα(x, y, z)=1.

[Formula 13]
If xi<Xis or xi>Xie or yi<Yis or yi>Yie or zi<Zis or zi>Zie (i.e. outside the valid area): Vα=−1 (invalid value)
Otherwise (i.e. valid area): Vα=Vo(xi, yi, zi)
However, it is defined as Vo(x, y, z)=Cmap(Do(x, y, z), 3)·Mask(x, y, z)·Sα(x, y, z). Cmap(d, c) is the color map, d is the signal value of the tomographic image, c is the color and opacity component, and c=0 (R), 1 (G), 2 (B), 3 (α). be.

In [Equation 13], if the signal value is outside the effective range, -1 is set to Vα. can be any value.

If Vα=0, the process proceeds to step S903; if Vα<0, the process proceeds to step S904; and if Vα>0, the process proceeds to step S905. That is, if the voxel is transparent, the process proceeds to step S903, if the voxel is outside the valid area, the process proceeds to step S904, and if the voxel is opaque, the process proceeds to step S905.

In step S903, the origin coordinate searching means 62 skips a predetermined number of voxels in the z-axis direction. Specifically, a predetermined number m2 (for example, m2=4) is subtracted from the current z-coordinate. Note that when compared with the predetermined number m1 used in steps S740 and S760, m2>m1. In this embodiment, m1=1. As a result of the subtraction, if the z-coordinate is greater than or equal to 0, the process returns to step S902 and repeats the same processing. If the z-coordinate is less than 0, the process proceeds to step S904.

In step S904, the starting point coordinate searching means 62 obtains a value indicating that the starting point coordinate was not found for the z coordinate corresponding to the current target pixel (x, y) (the original lower limit value of z is 0). If so, a negative value such as -1 is set, and the search processing for the coordinates of the starting point is terminated.

When proceeding to step S905, the origin coordinate search means 62 initializes the variable i (step S905), and then adds m1 (step S906). If i has reached m2, or if adding i to the current z-coordinate reaches the upper limit value Zst, which is the viewpoint z-coordinate, the process advances to step S908 to add m1 to the current z-coordinate. Set the value obtained by adding i before If i has not reached m2 and if adding i to the current z-coordinate does not reach the upper limit value Zst, which is the viewpoint z-coordinate, the process proceeds to step S907.

In step S907, the starting point coordinate searching means 62 performs coordinate transformation on the three-dimensional coordinates (x, y, z+i) and determines whether or not the corresponding voxels are transparent. The determination procedure is the same as step S902. However, it is not determined whether or not it is outside the effective area. If Vα>0, the process returns to step S906, m1 is added to i, and the same processing is repeated. If Vα=0, the process proceeds to step S908.

In step S908, the starting point coordinate search means 62 sets the z coordinate (z+i−m1) immediately before step S907 corresponding to the current target pixel (x, y), and ends the search for the starting point coordinates. .

14, the starting point coordinate searching means 62 changes the z coordinate of each pixel (x, y) of the rendering image from the search starting coordinate Zst to the lower limit value of 0 by a predetermined number of m2 in the viewpoint coordinate system. while performing a predetermined coordinate transformation for each of the three-dimensional coordinates (x, y, z) of the viewpoint coordinate system, sequentially obtaining the voxel values of the corresponding voxel image, and obtaining the viewpoint coordinates of the opaque voxel found first. The z coordinate in the system is searched as the origin coordinate. More specifically, the z-coordinate is subtracted by a predetermined number m2 (m2>m1) larger than the predetermined number m1 within a range from a predetermined upper limit value to a predetermined lower limit value, and coordinate transformation is performed for each three-dimensional coordinate. is performed to determine whether or not the voxel is opaque (S902), and if the signal value of the voxel at the z-coordinate after the subtraction of the predetermined number m2 is opaque (S905), the predetermined number m1 is subtracted from the z-coordinate after the subtraction. After increasing (S906), each three-dimensional coordinate (x, y, z) is subjected to coordinate transformation to determine whether or not the voxel is opaque (S907).

When the origin coordinate z is found by executing the process shown in FIG. 14, the process of step S642 shown in FIG. 13 is completed. If the obtained starting point coordinate z is z<0, it means that there is no opaque voxel, and the process of step S640 shown in FIG. 13 is terminated. If the obtained origin coordinate z satisfies z≧0, the three-dimensional coordinates (x, y, z) are coordinate-transformed to acquire the voxel opacity (step S643).

Specifically, as in steps S902 and S907, the processing according to [Formula 9] to [Formula 11] is executed to obtain the three-dimensional coordinates (x, y, z) (integer values) of the given viewpoint coordinate system. ), real number coordinates (xr, yr, zr) corresponding to voxels of a voxel image corresponding to each pixel of a plurality of DICOM-format tomographic images Do (x, y, z) corresponding to ) are calculated.

Next, the calculated coordinates (xr, yr, zr) of the voxel coordinate system are converted to integer values. Unlike the processing in steps S902 and S907, first, the coordinates obtained by truncating the decimal point are set to (xia, yia, zia), and the truncated fraction is set to (wx, wy, wz). That is, xr=xia+wx, yr=yia+wy, and zr=zia+wz.

Next, using the effective area and the mask data Mask, the process according to the following [Equation 14] is executed to obtain the voxel opacity Vα corresponding to the coordinates (xia, yia, zia) in the voxel coordinate system. . In the following [Formula 14], the processing is divided into three types. In [Equation 14], Vo(x, y, z)=Cmap(Do(x, y, z), 3)·Mask(x, y, z)·Sα(x, y, z).

[Formula 14]
1) If xia<Xis or xia>Xie or yia<Yis or yia>Yie or zia<Zis or zia>Zie (outside the effective area): Vα=−1
2) If the above 1) is not satisfied and xia+1>Xie or yia+1>Yie or zia+1>Zie, Vα=Vo(xia, yia, zia) without performing interpolation processing
3) If the above 1) and 2) are not satisfied, the following interpolation processing is executed Vα=(1−wz)(1−wy)(1−wx)·Vo(xia, yia, zia)
+(1−wz)(1−wy)・wx・Vo(xia+1, yia, zia)
+(1−wz)・wy・(1−wx)・Vo(xia, yia+1, zia)
+(1-wz).wy.wx.Vo(xia+1, yia+1, zia)
+wz・(1−wy)(1−wx)・Vo(xia, yia, zia+1)
+wz.(1-wy).wx.Vo(xia+1, yia, zia+1)
+wz.wy.(1-wx).Vo(xia, yia+1, zia+1)
+wz.wy.wx.Vo(xia+1, yia+1, zia+1)

In 3) of [Equation 14], eight integer voxel coordinates (xia, yia, zia), etc. located in the vicinity of the real voxel coordinates (xr, yr, zr) are used for interpolation processing. Each value of the tomographic image, mask data, and opacity correction table corresponding to these eight coordinate values is obtained, and each signal value obtained from the tomographic image is converted into opacity while referring to the color map. is multiplied by predetermined weights wx, wy, wz, (1-wx), (1-wy), and (1-wz), and the values added are given as the voxel opacity Vα corresponding to the three-dimensional coordinates. . In [Equation 14], when it is outside the effective area, that is, when it is outside the effective area, -1 is set to Vα. can be any value. Further, in step S643, the opacity Vα is divided by 255, which is the maximum value that can be taken by 8 bits, to calculate Alpha (=Vα/255) with a maximum value of 1 or less.

If Alpha=0 and z>0 as a result of the processing in step S643, it means a transparent pixel within the effective area, so the process proceeds to step S647. If Alpha>0, it means that the pixel is opaque, so the process proceeds to step S644. If Alpha<0 or Alpha=0 and z=0, it means that it is outside the effective area or the lower limit of the search range has been reached, so the process of step S640 shown in FIG. 13 is terminated.

In step S647, the ray casting means 61 subtracts the coordinates of the viewpoint by a predetermined number m1. This predetermined number m1 is an integer. As described above, in this embodiment, m1=1 and m1<m2. Then, the process returns to step S642 and repeats the same processing.

As described above, if Alpha>0, it means that the pixel is opaque, so the ray casting means 61 acquires the RGB value V(c) (step S644).

Specifically, as in step S643, the processing according to [Formula 9] to [Formula 11] is executed to obtain the three-dimensional coordinates (x, y, z) (integer values) of the given viewpoint coordinate system. Calculate the real number coordinates (xr, yr, zr) corresponding to the voxel of the voxel image composed of the corresponding multiple tomographic images Do (x, y, z) in DICOM format, and round the decimal point to the integer. is (xia, yia, zia), and the truncated fraction below the decimal point is (wx, wy, wz). That is, xr=xia+wx, yr=yia+wy, and zr=zia+wz.

Next, using the effective area, the processing according to the following [Equation 15] is executed to obtain the voxel RGB value V(c) corresponding to the coordinates (xia, yia, zia) in the voxel coordinate system. In the following [Formula 15], the processing is divided into three types. In [Formula 15], defined as Vc (x, y, z, c) = Cmap (Do (x, y, z), c) (c = 0 (R), 1 (G), 2 (B)) do. Note that when calculating the RGB value V(c). Unlike [Equation 14], the mask data Mask and the opacity correction table Sα are not referred to.

[Formula 15]
1) If xia<Xis or xia>Xie or yia<Yis or yia>Yie or zia<Zis or zia>Zie (outside effective area): V(c)=0 (c=0, 1, 2)
2) If the above 1) is not satisfied and xia+1>Xie or yia+1>Yie or zia+1>Zie, V(c)=Vc(xia, yia, zia, c) without performing interpolation processing
3) If the above 1) and 2) are not satisfied, execute the following interpolation processing V(c)=(1−wz)(1−wy)(1−wx)·Vc(xia, yia, zia, c)+ (1-wz) (1-wy) wx Vc (xia+1, yia, zia, c)
+ (1-wz) wy (1-wx) Vc (xia, yia+1, zia, c)
+(1-wz) wy wx Vc (xia+1, yia+1, zia, c)
+wz (1-wy) (1-wx) Vc (xia, yia, zia+1, c)
+ wz (1-wy) wx Vc (xia+1, yia, zia+1, c)
+wz.wy.(1-wx).Vc(xia, yia+1, zia+1, c)
+wz wy wx Vc (xia+1, yia+1, zia+1, c)

In 3) of [Equation 15], eight integer voxel coordinates (xia, yia, zia), etc. located in the vicinity of real-valued voxel coordinates (xr, yr, zr) are used for interpolation processing. Each signal value of the tomographic image corresponding to these eight coordinate values is obtained, and each signal value obtained from the tomographic image is converted into an RGB value while referring to the color map, and given weights wx and wy , wz, (1-wx), (1-wy), and (1-wz) are multiplied and added as the RGB value V(c) of the voxel corresponding to the three-dimensional coordinates. In 1) of [Equation 15], V(c) is set to 0 if it is outside the effective region.

In the present embodiment, when the RGB values V(c) are acquired in step S644, shadow addition processing is also performed. In step S644, the shadow addition means 63 performs shadow addition processing. This is a process of adding a shadow corresponding to a predetermined lighting environment to a voxel image. Specifically, the voxel opacity α , the color components of voxels in the range of Xs<x<Xe, Ys<y<Ye and Zs<z<Ze are updated. For this purpose, a certain parallel light source and ambient light components are set, and shadows to be added are calculated. Specifically, the line-of-sight vector (Ex, Ey, Ez) as a unit vector indicating the direction of the line of sight, the light source vector (Lx, Ly, Lz) as a unit vector indicating the direction of the parallel light source, the ambient light component Ab, the diffuse reflection coefficient Df (real number of 0 or more), specular reflection coefficient Sp (real number of 0 or more), and specular glossiness Sh (integer of 1 or more) are set.

The line-of-sight vector (Ex, Ey, Ez) is usually a fixed value of (Ex, Ey, Ez) = (0, 0, -1) in the negative z-axis direction, the light source vector (Lx, Ly, Lz), the diffuse reflection coefficient The numerical values of Df, specular reflection coefficient Sp, specular glossiness Sh, and ambient light component Ab can be appropriately set according to the situation of the image to be expressed. , 0.57735, −0.57735), Df=1.2, Sp=1.7, Sh=8, Ab=0.2. Lx, Ly, Lz are real numbers with an absolute value of less than 1 (including negative values), Df is a real number of 0 or more, Sp is a real number of 0 or more, Sh is an integer of 1 or more, Ab is 0 or more of 1 The following are real numbers.

First, using the 3×3 rotation inverse matrix generated as the transformation parameter, the processing according to the following [Equation 16] is executed to obtain the light source vector (Lx, Ly, Lz), the sight line vector Ez) is rotated and converted into a vector in the voxel coordinate system.

[Formula 16]
Lx←R11・Lx+R12・Ly+R13・Lz
Ly←R21・Lx+R22・Ly+R23・Lz
Lz←R31・Lx+R32・Ly+R33・Lz
Ex←R11・Ex+R12・Ey+R13・Ez
Ey←R21・Ex+R22・Ey+R23・Ez
Ez ← R31 · Ex + R32 · Ey + R33 · Ez

A light source vector (Lx, Ly, Lz) and a line-of-sight vector (Ex, Ey, Ez) in a voxel coordinate system are obtained as coordinates after rotation processing. Next, for the voxel (xia, yia, zia) to be shaded, the resolution in the xy-axis direction of the voxel image is set to Rxy (the reciprocal of the pixel interval, the number of pixels per mm, the x-axis and y-axis The direction is the same), the resolution in the z-axis direction of the voxel image is Rz (reciprocal of the slice interval of the tomographic image, the number of slices of the tomographic image per mm), and Rz/Rxy<1, the following [Equation 17] to calculate real-valued z-coordinates zr1 and zr2. This is because when the resolution in the z-axis direction is lower than the resolution in the xy-axis direction, a step occurs in the calculated gradient vector and shadow value along the z-axis direction. This is for smoothing.

[Formula 17]
zr1=(zia+1).Rz/Rxy
zr2=(zia−1)·Rz/Rxy

Next, the calculated real-value z-coordinates zr1 and zr2 are converted to integer values. First, let zi1 and zi2 be real-valued z-coordinates obtained by truncating the decimal point and converting them to integers, and wz1 and wz2 be the truncated fractional components. That is, zr1=zi1+wz1 and zr2=zi2+wz2.

Subsequently, the opacity at voxel (x, y, z) is Vα (x, y, z) = Cmap (Do (x, y, z), 3), and for voxel (xia, yia, zia) , the opacity of all 8 neighboring voxels Vα (xia+1, yia, zia), Vα (xia−1, yia, zia), Vα (xia, yia+1, zia), Vα (xia, yia-1, zia), Vα (xia, yia, zi1), Vα (xia, yia, zi1+1), Vα (xia, yia, zi2), Vα (xia, yia, zi2+1).

When Rz/Rxy≧1, [Formula 17] is not executed, so Vα(xia, yia, zi1), Vα(xia, yia, zi1+1), Vα(xia, yia, zi2), Vα(xia, yia, zi2+1), obtain the opacities Vα(xia, yia, zia+1) and Vα(xia, yia, zia−1) of voxels that are positive and negative 2 in the z-axis direction.

however,
Vα(xia+1, yia, zia)=0 when xia+1>Xie,
Vα(xia-1, yia, zia)=0 when xia-1<Xis,
Vα(xia, yia+1, zia)=0 when yia+1>Yie,
Vα(xia, yia-1, zia)=0 when yia-1<Yis,
If Rz/Rxy<1,
Vα(xia, yia, zi1)=0 when zi1>Zie,
Vα(xia, yia, zi1+1)=0 when zi1+1<Zie,
Vα(xia, yia, zi2)=0 when zi2<Zis,
Vα(xia, yia, zi2+1)=0 when zi2+1>Zis,
If Rz/Rxy≧1,
Vα(xia, yia+1, zia+1)=0 when zia+1>Zie,
Vα(xia, yia-1, zia-1)=0 when zia-1<Zis,
and

Subsequently, when Rz/Rxy<1, four neighboring voxels of Vα(xia, yia, zi1), Vα(xia, yia, zi1+1), Vα(xia, yia, zi2), Vα(xia, yia, zi2+1) is interpolated in the z-direction using the fractional components wz1 and wz2 to calculate the z-axis direction interpolation components V'a and V'b.

[Formula 18]
V′a=(1−wz1)·Vα(xia, yia, zi1)+wz1·Vα(xia, yia, zi1+1)
V′b=(1−wz2)·Vα(xia, yia, zi2)+wz2·Vα(xia, yia, zi2+1)

Then, the gradient vector (Gx, Gy, Gz) at the voxel (xia, yia, zia) of the voxel image is calculated by executing the process according to [Equation 19] below.

[Formula 19]
Gx = Vα (xia+1, yia, zia) - Vα (xia-1, yia, zia)
Gy = V α (xia, yia + 1, zia) - V α (xia, yia - 1, zia)
If Rz/Rxy<1, then Gz=V'a-V'b
If Rz/Rxy≧1,
Gz=Vα(xia, yia, zia+1)−Vα(xia, yia, zia−1)
G=(Gx ² +Gy ² +Gz ² ) ^1/2

In [Equation 19], the z-axis direction component Gz of the gradient vector is obtained by multiplying the opacity of two voxels near the coordinates shifted +Rz/Rxy in the z-axis direction with respect to the voxel and adding a predetermined weight. It is the difference between the opacity and the opacity obtained by multiplying the opacity of two voxels in the vicinity of the coordinates shifted by -Rz/Rxy in the z-axis direction with a predetermined weight and adding them. In [Formula 19], G is the magnitude of the gradient vector (Gx, Gy, Gz). As shown in [Equation 19], the shading means 63 shades the voxel of the voxel image in the x-axis direction, the y-axis direction, and the z-axis direction by 2 each (when Rz/Rxy<1, the z-axis direction is 4 (neighborhood) voxels are calculated as x-axis, y-axis, and z-axis components of the gradient vector. A unit vector obtained by correcting the directional component or the z-axis component is calculated as the gradient vector in the voxel. Further, if G≧1, divide the gradient vector (Gx, Gy, Gz) by G and replace the gradient vector (Gx, Gy, Gz) with a unit vector. When G<1, the shadow value cannot be calculated because there is no opacity gradient at the voxel (xia, yia, zia). Gives the light component Ab.

Once the gradient vectors have been calculated, the shading coefficient S is calculated using the well-known Phong shading model. Specifically, first, when the specular reflection coefficient Sp>0, using the gradient vector (Gx, Gy, Gz) and the light source vector (Lx, Ly, Lz), which are unit vectors, the following [Equation 20] is A specular reflection vector (Fx, Fy, Fz) is calculated by executing the processing according to the above. If the specular reflection coefficient Sp=0, then (Fx, Fy, Fz)=(0, 0, 0).

[Formula 20]
Fx=2.Gx.(Gx.Lx+Gy.Ly+Gz.Lz)-Lx
Fy = 2 · Gy · (Gx · Lx + Gy · Ly + Gz · Lz) - Ly
Fz=2.Gz.(Gx.Lx+Gy.Ly+Gz.Lz)-Lz

In [Equation 20], the angle formed by the specular reflection vector and the gradient vector is the same as the angle formed by the light source vector and the gradient vector, and the angle formed by the light source vector is equal to the angle formed by the light source vector and the gradient vector. Defined to double. Next, gradient vectors (Gx, Gy, Gz), light source vectors (Lx, Ly, Lz), specular reflection vectors (Fx, Fy, Fz), line-of-sight vectors (Ex, Ey, Ez) are all used as unit vectors, Using the reflection coefficient Df (real number of 0 or more), the specular reflection coefficient Sp (real number of 0 or more), the specular gloss Sh (integer of 1 or more), and the ambient light component Ab, the following [Formula 21] A shadow coefficient S(xia, yia, zia) at the voxel (xia, yia, zia) is calculated by executing the processing according to the above.

[Formula 21]
S (xia, yia, zia) = Df · (Gx · Lx + Gy · Ly + Gz · Lz) + Sp · (Fx · Ex + Fy · Ey + Fz · Ez) ^Sh + Ab

"Df.(Gx.Lx+Gy.Ly+Gz.Lz)" in [Equation 21] is diffusely reflected light, and "Sp.(Fx.Ex+Fy.Ey+Fz.Ez) ^Sh " is specularly reflected light. In [Equation 21], the diffuse reflection light of the voxel is calculated based on the inner product value of the calculated gradient vector (Gx, Gy, Gz) and the light source vector (Lx, Ly, Lz), and the specular reflection coefficient Sp is When a positive value is set, the specular reflection vector (Fx, Fy, Fz) of the predetermined light source vector (Lx, Ly, Lz) with respect to the gradient vector (Gx, Gy, Gz) is calculated, and the specular reflection vector Based on the inner product value of (Fx, Fy, Fz) and a predetermined line-of-sight vector (Ex, Ey, Ez), the specular reflection light of the voxel was calculated, and the calculated diffuse reflection light and specular reflection light were calculated. The color component of the voxel obtained by the ray casting means is multiplied by the shadow value to update the color component of the voxel as follows. As the specular reflection light, the inner product value of the regular reflection vector (Fx, Fy, Fz) and the line-of-sight vector (Ex, Ey, Ez) is multiplied by the specular reflection coefficient Sp to the power of the specular glossiness Sh. value. After the shadow coefficient S (xia, yia, zia) is calculated, the RGB value V ( c)=0 (c=0, 1, 2) is corrected by executing the process according to [Equation 22] below. That is, add shadows.

[Formula 22]
V(c)←S(xia, yia, zia) V(c)

In [Equation 22], c=0, 1, 2, corresponding to the R, G, and B color components, respectively. That is, in [Formula 22], only the color component is corrected. As a result, the shading means 63 executes the processing according to [Formula 16] to [Formula 22], and refers to the opacity α of the pixels located near the voxels in the xyz-axis direction. Values of color components of voxels are corrected and shading is added. In this embodiment, as a preferred example, adjacent pixels are used as pixels for which the opacity α is referred to, but the pixels are not limited to adjacent pixels, and other nearby pixels may be referred to.

In step S645, when the RGB value V(c) of the shaded voxel is obtained, the accumulated luminance value Energy(c) and the transmitted light intensity Trans are updated (step S646). Specifically, the process according to the following [Formula 23] is executed to correct the α value, and then the cumulative luminance value Energy(c) and the transmitted light intensity Trans are updated.

[Formula 23]
Alpha ← 1-(1-Alpha) ^1/Sray
Energy (c) <- Energy (c) + Trans Alpha V (c)/255
Trans←Trans・(1.0-Alpha)

After updating the cumulative luminance Energy(c) and the transmitted light intensity Trans, if Trans<0.001, the process of step S640 shown in FIG. 13 is terminated. If Trans≧0.001, z is reduced by a predetermined number m1 (step S646). In this embodiment, m1=1 and m1<m2. Then, if z is 0 or more, the process returns to step S643 to repeat the same processing. If z is less than 0, there is no voxel to be processed, so the process of step S640 shown in FIG. 13 ends.

As shown in the third formula of [Equation 9], z is multiplied by the sub-sampling magnification Sray of the virtual ray during the offset processing in step S902. Therefore, substantially, in the z-axis direction, voxels are referenced at predetermined reference intervals (m1·Sray) to obtain the voxel pixel value V(c). Therefore, if (m1·Sray)>1, processing is performed with reference to thinned voxels instead of all voxels in the z-axis direction. Therefore, the number of voxels to be processed in the z-axis direction is reduced and efficiency is improved, but in the case of volume rendering, image quality degradation becomes noticeable (in the case of 3D-MIP, which generates a rendering image based on signal values, the processing load is is larger than the volume rendering, so set (m1·Sray)>1, usually set m1=1, Sray=2). In this embodiment, it is preferable that m1=1 and Sray≦1, and normally Sray=1.

In steps S630 and S640 shown in FIG. 12, the ray-casting means 61 sets the upper limit of the z-axis, which is the line-of-sight direction, in the viewpoint coordinate system to , Zc corresponding to the intersection of the line-of-sight direction vector calculated for each pixel (x, y) of the rendered image and the cuboid, and the lower limit of the z-axis, which is the line-of-sight direction, is set to the z-axis direction of the viewpoint coordinate system. is set to 0, which is the lower limit of , and the starting point coordinate searching means 62 searches for the starting point coordinate Zst at which the irradiation of the virtual ray starts in the range from Zc to 0 on the z-axis of the viewpoint coordinate system, with the search start coordinate being Zc. ing.

Each of the NB CPU cores executes the process of step S440 shown in FIG. 10 for obtaining the color component of each pixel as described above for each of the NB divided images. Then, by integrating the obtained divided images, a rendering image, which is an RGB color two-dimensional projection image of the set Size×Size, is generated. Note that it is not necessary to integrate the set Size×Size two-dimensional projection image all at once, and all the divided images are finally integrated while sequentially displaying each divided image as soon as the processing is completed. good too.

The boundary surface of the effective area, that is, the value of the opacity α of the outermost voxel of the effective area is set to 0, and the RGB value of the color component is the natural RGB value defined by the color map based on the signal value of the pixel. By providing , it is possible to suppress moiré generated on the boundary surface while suppressing the processing load.

As described above, in the volume rendering device 100 according to the present embodiment, the opaque circumscribed area calculation means 50 and the rendering means 60 are implemented by the CPU 1, which is a multi-core CPU, executing programs. This program divides the loaded tomographic image into a plurality of (for example, NA or NB) tomographic image groups, assigns each tomographic image group to each CPU thread, and each CPU thread calculates a part of the opaque circumscribed area. Processing and processing for assigning pixel values to divided generation regions are performed in parallel. A case in which the volume rendering apparatus 100 is realized by a personal computer, which is a general-purpose computer, will be described. FIG. 15 is a diagram showing parallel processing when the volume rendering device 100 is implemented by a computer. FIG. 15 shows an example in which the opaque circumscribed region calculation means 50 is divided into NA=4 tomographic image groups and executed. will have a similar configuration. As shown in FIG. 15, the opaque circumscribed area calculation means 50 and the rendering means 60 are implemented by an application program. For example, when dividing into tomographic image groups of NA=4, the application program divides the processing for each of the four tomographic image groups into four processing threads #n1 to #n4. For example, when dividing into NB=4 divided generation areas, the application program divides the processing for each of the four divided generation areas into four processing threads #n1 to #n4. Then, the job scheduler of the multitasking operating system in which the opaque circumscribed area calculating means 50 and the rendering means 60 (application program) operate allocates them to the CPU threads #1 to #8 of the plurality of CPU cores.

FIG. 16 is a diagram showing a software flow of parallel processing when the volume rendering apparatus 100 is realized by a computer. FIG. 16 shows an example in which the opaque circumscribed region calculation means 50 is divided into NA=4 tomographic image groups and executed. will have a similar configuration. As shown in FIG. 16, the application program that realizes the opaque circumscribed region calculation means 50 and the rendering means 60 processes four tomographic image groups and four divided generation regions using four processing threads #n1. ∼ After dividing into processing threads #n4, each starts a parallel processing function. Then, the job scheduler of the multitasking operating system assigns four threads #n1 to #n4 to CPU thread #3, CPU thread #5, CPU thread #7, and CPU thread #2 each having a plurality of CPU cores. assign. In this way, the operating system allocates the processing threads to be processed to the CPU threads of the CPU cores that are free, thereby enabling parallel processing. In the example of FIG. 16, the application program waits for end messages from each CPU thread, and when the end messages are obtained from all CPU threads, the next process is performed.

As a result of actual measurement on a Windows10/64bits personal computer equipped with a 4-core configuration CPU (Intel CORE-i7), the rendering means 60 was set to NB=8 using a 512×512×370 voxel image with the configuration shown in FIG. When divided into divided generation areas and executed without shadow calculation, it takes about 3.5 times less than single thread execution (512 × 512 full resolution rendering image generation time by CPU: about 1.05 seconds) A processing speed was obtained (same generation time: about 0.3 seconds). However, some CPU threads were waiting for processing. Therefore, when shadows including specular reflection proposed in the present application are added and parallel processing is performed with 8 processing threads, the processing performance is almost 4 times that of single thread execution (same generation time: about 1.8 seconds). (Generation time: about 0.45 sec). Furthermore, as a result of experimenting with NA > 8, the processing performance hardly changed and reached a ceiling at 4 times. it was easy to get The above processing time includes the processing time by the opaque circumscribed area calculation means 50, which is much shorter than the processing time of the rendering means 60. The effect of treatment was difficult to measure.

<3. Gradation Reduction of Tomographic Image>
In the above embodiment, a tomographic image in which each pixel records a 16-bit signal value is used as it is to create a rendered image. A rendered image may be created later. By converting the tomographic image into 8-bit gradation, the load in processing each pixel can be reduced, contributing to high-speed generation of volume rendering images. As a result of actual measurement, a reduction effect of about 20% was obtained, and when rendering with shadows was executed with 8 threads on a 4-core CPU configuration, it was about 0.45 sec when using a tomographic image with a 16-bit signal value. On the other hand, when the 8-bit signal value gradation reduction tomographic image was used, the time was shortened to about 0.36 sec, and almost no difference was observed in the generated rendered image.

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

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

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

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

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

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

In the above embodiment, the opacity of the outermost voxels in the effective area is set to 0 in order to prevent moiré. need not be 0.

1 CPU (Central Processing Unit)
2 RAM (random access memory)
3... Storage device 4... Instruction input I/F
5 Data input/output I/F
6... Display unit 7... GPU
8 frame memory 10 tomographic image reading means 20 color map reading means 30 ROI clipping setting means 40 mask setting means 50 opaque circumscribed area calculation means 60 Rendering Means 61 Ray Casting Means 62 Starting Point Coordinate Searching Means 63 Shadow Adding Means 70 Rendered Image Output Means 100 Volume Rendering Apparatus

Claims (16)

A color map defined by associating color component values and opacity with signal values based on a voxel image composed of a plurality of two-dimensional tomographic images taken at predetermined intervals of a predetermined target. A volume rendering device for generating a rendered image with reference to
rendering means for generating a rendered image using the voxel image and the color map;
The rendering means refers to a plurality of voxels of the voxel image while performing a predetermined coordinate transformation on the voxel image for each pixel of the rendering image, and refers to the color map for each voxel to be referred to. a ray casting means for obtaining the color component and opacity of the voxel while obtaining the opacity of a plurality of neighboring voxels of the voxel;
calculating a gradient vector in the voxel based on the opacity of a plurality of voxels near the voxel obtained by the ray casting means, and calculating a regular reflection vector of a predetermined light source vector with respect to the calculated gradient vector; calculating the diffusely reflected light of the voxel based on the inner product value of the gradient vector and the light source vector, and calculating the specular reflected light of the voxel based on the inner product value of the regular reflection vector and a predetermined line-of-sight vector; and a shadow adding means for updating the color component of the voxel by multiplying the color component of the voxel obtained by the ray casting means by a shadow value calculated based on the calculated diffuse reflection light and the specular reflection light. Volume rendering device. The predetermined light source vector is obtained by subjecting a predefined light source vector to rotation processing in the predetermined coordinate transformation. 2. The volume rendering apparatus according to claim 1, wherein rotation processing is performed in coordinate transformation. The rotation process is
3. The volume rendering apparatus according to claim 2, wherein a rotation process is performed on the predefined light source vector or the predefined line-of-sight vector using a 3×3 matrix that defines rotation in a predetermined three-dimensional coordinate system. . Assuming that the resolution in the xy-axis direction of the voxel image is Rxy and the resolution in the z-axis direction is Rz,
The shading means includes:
The difference in opacity between voxels shifted +1 voxel and -1 voxel in the x-axis direction with respect to the voxel is defined as the x-axis component of the gradient vector,
The difference in opacity between voxels shifted +1 voxel and -1 voxel in the y-axis direction with respect to the voxel is defined as the y-axis direction component of the gradient vector,
Claims 1 to 3, wherein a value obtained by multiplying a difference in opacity of voxels shifted by +1 voxel and -1 voxel in the z-axis direction from the voxel concerned by Rz/Rxy is taken as the z-axis direction component of the gradient vector. A volume rendering device according to any one of Claims 1 to 3. Assuming that the resolution in the xy-axis direction of the voxel image is Rxy and the resolution in the z-axis direction is Rz,
The shading means includes:
The difference in opacity between voxels shifted +1 voxel and -1 voxel in the x-axis direction with respect to the voxel is defined as the x-axis component of the gradient vector,
The difference in opacity between voxels shifted +1 voxel and -1 voxel in the y-axis direction with respect to the voxel is defined as the y-axis direction component of the gradient vector,
For the voxel, the opacity obtained by multiplying and adding the opacity of two voxels in the vicinity of the coordinates shifted +Rz/Rxy in the z-axis direction by a predetermined weight, and the opacity of the coordinates shifted -Rz/Rxy in the z-axis direction. 4. The volume rendering according to any one of claims 1 to 3, wherein the difference between the opacities obtained by multiplying the opacities of two neighboring voxels by a predetermined weight and adding them is used as the z-axis direction component of the gradient vector. Device. The shading means includes:
Predetermined diffuse reflection coefficient, specular reflection coefficient, specular glossiness, and ambient light are defined in advance,
the specular reflection vector is defined so that the angle formed by the gradient vector is the same as the angle formed by the light source vector and the gradient vector;
When all of the gradient vector, specular reflection vector, light source vector, and line-of-sight vector are unit vectors,
The diffusely reflected light is a value obtained by multiplying the inner product value of the gradient vector and the light source vector by the diffuse reflection coefficient,
As the specular reflected light, a value obtained by multiplying the inner product value of the specular reflection vector and the line-of-sight vector by the specular reflection coefficient to a value raised to the power of the specular glossiness,
6. The volume rendering apparatus according to any one of claims 1 to 5, wherein the shadow value is a value obtained by adding the ambient light to the diffuse reflection light and the specular reflection light. dividing the generation area of the rendered image into NB (NB>1) divided generation areas in the x-axis direction;
3. The rendering means generates a rendering image composed of NB divided images by causing NB different processors to simultaneously execute NB divided generation regions that are different from each other. A volume rendering device according to any one of claims 1 to 6. calculating an opaque circumscribed region, which is a region that circumscribes in each three-dimensional axial direction the region in which pixels having non-zero opacity obtained by referring to the color map exist among the pixels of each of the tomographic images; and an opaque circumscribed area calculation means,
The rendering means defines an area corresponding to the opaque circumscribed area in each of the tomographic images as an effective area,
The ray casting means calculates, at each pixel of the rendered image, an intersection point between the effective area and a line-of-sight direction vector closest to a viewpoint, and calculates the voxel from the calculated intersection point at each pixel of the rendered image. 8. The voxel image according to any one of claims 1 to 7, wherein a virtual ray is applied to the image in a line-of-sight direction, and a plurality of voxels of the voxel image are referred to while performing a predetermined coordinate transformation based on a ray casting method. volume rendering device. When mask data having a value of 0 (not displayed) or 1 (displayed) is defined corresponding to the tomographic image,
9. The volume rendering apparatus according to claim 8, wherein said opaque circumscribing region calculating means further refers to said mask data to specify pixels whose value of said mask data is 1 and whose opacity is not 0. . The opaque circumscribing region calculating means divides each of the tomographic images into NA tomographic image groups with a predetermined division number NA (NA>1), and divides each of the NA tomographic image groups into the opaque circumscribing region. 10. A volume rendering apparatus according to claim 8 or 9, wherein processes for creating a part of the region are executed in parallel. Xso is the minimum coordinate in the first dimension of each of the two-dimensional tomographic images, Xeo is the maximum coordinate in the first dimension, Yso is the minimum coordinate in the second dimension, and Yso is the minimum coordinate in the second dimension. Let Yeo be the maximum coordinate, Zso be the z coordinate at which the first tomographic image is arranged, and Zeo be the z coordinate at which the Mth tomographic image is arranged. If six coordinates Xs, Xe, Ys, Ye, Zs, Ze are defined that satisfy Ye≤Yeo and Zso≤Zs<Ze≤Zeo,
11. The volume rendering apparatus according to any one of claims 8 to 10, wherein the opaque circumscribing area calculating means calculates the opaque circumscribing area in a range of Xs to Xe, Ys to Ye, and Zs to Ze. Assuming that a coordinate system obtained by projectively transforming the voxel image into the rendering image is a viewpoint coordinate system, the ray casting means sets the upper limit value of the z-axis, which is the viewing direction in the viewpoint coordinate system, to the pixel (x, y ) is set to Zc corresponding to the calculated intersection point,
12. The apparatus according to any one of claims 8 to 11, further comprising starting point coordinate searching means for searching for starting point coordinates Zo at which irradiation of the virtual ray starts within a predetermined lower limit range from Zc on the z-axis of the viewpoint coordinate system. The volume rendering device according to . The origin coordinate search means is
In the viewpoint coordinate system, the three-dimensional coordinates (x, y, The predetermined coordinate transformation is performed for each z) to sequentially acquire the voxel values of the corresponding voxel of the voxel image, and the z coordinate of the opaque voxel found first in the viewpoint coordinate system is searched as the starting point coordinate Zo. 13. A volume rendering apparatus according to claim 12, wherein: The raycasting means comprises:
when irradiating each pixel (x, y) of the rendered image with a virtual ray having a predetermined light intensity from the starting point coordinates searched by the starting point coordinate searching means toward the lower limit value along the z-axis; , performing the predetermined coordinate transformation for each of the three-dimensional coordinates (x, y, z) of the viewpoint coordinate system, and acquiring the opacity as the voxel value of the corresponding voxel of the voxel image while referring to the color map. , if the obtained opacity is not 0, further obtain the color components (RGB) referring to the color map as voxel values;
The shading means updates the color component by multiplying the obtained color component by a shading value calculated based on the opacity of a plurality of neighboring voxels of the voxel;
The rendering means is
performing a process of attenuating the light intensity based on the opacity of the voxel and calculating a cumulative luminance value based on the opacity of the voxel, the updated color component, and the attenuated light intensity; 13. The volume rendering apparatus according to claim 12, wherein a pixel value (RGB) corresponding to the pixel (x, y) of the rendered image is given based on the accumulated brightness value. 3×3 matrix defining rotation in a predetermined three-dimensional coordinate system, offset values and minute offset values in the directions of xyz axes, enlargement or reduction ratios in the directions of xyz axes, and the predetermined coordinates including scaling factors in the z-axis direction get the parameters of the transform,
calculating real-valued voxel coordinates of the voxel image by converting the integer-valued three-dimensional coordinates (x, y, z) of the viewpoint coordinate system into the coordinate system of the voxel image based on the parameters;
identifying a plurality of voxels of the voxel image corresponding to a plurality of integer-valued voxel coordinates near the calculated real-valued voxel coordinates;
If voxel coordinates of integer values of the identified voxels are included in the effective area, the signal values of the plurality of voxels are converted into color components or opacity with reference to the color map, and the converted color components 15. The volume rendering apparatus according to any one of claims 12 to 14, wherein a color component or opacity is calculated as a voxel value by adding a value obtained by multiplying each opacity with a predetermined weight. . A program for causing a computer to function as the volume rendering device according to any one of claims 1 to 15.

Images