JP2020142003A - Mask generation device, three-dimensional reconstruction image generation device, mask generation method and program - Google Patents

Abstract

To easily generate mask data for visualizing a desired region.SOLUTION: A mask generation device 1 comprises: image acquisition means for acquiring a plurality of tomographic images; subject region calculation means for, for every tomographic image, calculating a subject region on the basis of a signal value of each pixel; subject shape calculation means for, for every tomographic image, approximating an outline of the calculated subject region by a prescribed geometry for calculating a parameter of the geometry; and mask data generation means for, for every tomographic image, calculating a mask value which defines whether each pixel of the tomographic image is included in inside or outside of the geometry, and generating mask data holding a mask value with respect to each pixel.SELECTED DRAWING: Figure 1

Description

The present disclosure relates to a mask generation device, a three-dimensional reconstruction image generation device, a mask generation method, and a program.

In medical image diagnosis, there are scenes where it is desired to obtain a three-dimensional reconstructed image suitable for observing a specific human body tissue. For example, when it is desired to observe predetermined human tissues such as internal organs and brain in the chest and head, the internal organs and brain are surrounded by bones, and the bone region rather hinders diagnosis. In a three-dimensional reconstructed image generated from a CT image, bones are generally clearly depicted and internal organs and blood vessels may be hidden. Therefore, it may be necessary to take some measures such as cutting for visualization.

For example, a three-dimensional mask was created by extracting the bone region to be unobserved using the region expansion method (region glowing method), and the bone region to be unobserved was removed by mask processing while referring to the three-dimensional mask. A method for generating a volume rendered image is disclosed (see Patent Documents 2 and 3). The area expansion method is a method in which pixels in a non-observation target area are designated, and neighboring pixels are extracted one after another with the pixels as a start point (expansion start point).

JP-A-2009-240569 Japanese Patent No. 4087517 Japanese Patent No. 5257958

However, creating a 3D mask using the region expansion method requires the user to instruct multiple expansion start points, and the user must also instruct the expansion termination stage, which makes automation difficult and results for each user. There is a problem that variations occur. Further, not limited to the area expansion method, the creation of a three-dimensional mask requires the skill and experience of the user, and it takes time and effort to create the mask.

The embodiment of the present disclosure has been made in view of such a point, and an object of the present disclosure is to provide a mask generator or the like capable of easily creating mask data for visualizing a desired region. And.

According to one embodiment of the present disclosure, an image acquisition means for acquiring a plurality of tomographic images, a subject area calculation means for calculating a subject area based on a signal value of each pixel for each tomographic image, and a tomographic image for each. The subject shape calculation means that approximates the calculated outline of the subject area with a predetermined geometric shape and calculates the parameters of the geometric shape, and for each tomographic image, each pixel of the tomographic image is either inside or outside the geometric shape. A mask generation device including a mask data generation means for calculating a mask value defining whether or not the image is included in the image and generating mask data for holding the mask value for each pixel is provided.

According to one embodiment of the present disclosure, an image acquisition means for acquiring a plurality of tomographic images, a subject area calculation means for calculating a subject area based on a signal value of each pixel for each tomographic image, and a tomographic image for each. The subject shape calculation means that approximates the calculated outline of the subject area with a predetermined geometric shape and calculates the parameters of the geometric shape, and for each tomographic image, each pixel of the tomographic image is either inside or outside the geometric shape. A mask data generation means that calculates a mask value that defines whether it is included in the image and generates mask data that holds the mask value of each pixel, and a three-dimensional structure that is configured by arranging each pixel of the tomographic image in a voxel space. A three-dimensional reconstruction image generation device including a voxel and a three-dimensional reconstruction image generation means for generating a three-dimensional reconstruction image based on the mask data is provided.

According to one embodiment of the present disclosure, an image acquisition step in which a computer acquires a plurality of tomographic images, a subject area calculation step in which a subject area is calculated based on a signal value of each pixel for each tomographic image, and a tomographic image. Each time, the outline of the calculated subject area is approximated by a predetermined geometric shape, and the subject shape calculation step of calculating the parameters of the geometric shape, and for each tomographic image, whether each pixel of the tomographic image is inside the geometrical shape. A mask data generation step of calculating a mask value defining which one is included in the outside and generating mask data holding the mask value of each pixel, and a mask generation method for executing the above are provided.

According to one embodiment of the present disclosure, a computer is used as an image acquisition means for acquiring a plurality of tomographic images, a subject area calculation means for calculating a subject area based on a signal value of each pixel for each tomographic image, and each tomographic image. , Subject shape calculation means that approximates the calculated outline of the subject area with a predetermined geometric shape and calculates the parameters of the geometric shape. For each tomographic image, each pixel of the tomographic image is either inside or outside the geometric shape. A program is provided that calculates a mask value that defines whether or not it is included in the image, and functions as a mask data generation means that generates mask data that holds the mask value of each pixel.

According to the present disclosure, mask data for visualizing a desired region can be easily created.

The figure which shows the functional structure of the mask generation apparatus 1. The figure which shows the hardware composition of the mask generation apparatus 1. Flow chart showing the operation of the mask generator 1 Flowchart showing the flow of drawing weight table creation process The figure explaining the rectangular area 41, the rectangular area 42, and the subject candidate area 50. Flowchart showing the flow of subject area calculation processing The figure which shows the outline of the subject area calculation process The figure explaining the calculation of the 2nd center of gravity The figure explaining the calculation of the 3rd center of gravity The figure which shows the calculated rectangular area 42 (A) Rectangular area 42 when the signal value inside the area is relatively large, (b) Rectangular area 42 when the signal value inside the area is relatively small. The figure which shows the content which corrects the rectangular area 42 based on the area ratio with the subject candidate area 50. The figure which shows the case where the correction based on the area ratio is applied to the upper chest The figure which shows the content which corrects (reduces correction) the rectangular area 42 based on the aspect ratio of the rectangular area 42. Flowchart showing the flow of subject area correction processing The figure which shows the two-dimensional distribution pattern of the weight value Sv (x, y, z) Flowchart showing the flow of mask data generation processing Flowchart showing the flow of mask data generation processing The figure which shows the outline which executes the drawing weight table creation process and the mask data creation process in parallel. The figure which shows the mechanism which execute the drawing weight table creation process and the mask data creation process in parallel. The figure which shows the outline of the reconstruction image generation processing executed by 3D reconstruction image generation apparatus 2. The figure which shows the functional structure of the 3D reconstruction image generation apparatus 2. The figure which shows the hardware composition of the 3D reconstruction image generation apparatus 2. A flowchart showing the operation of the three-dimensional reconstruction image generator 2. Flowchart showing the overall flow of the process of generating a volume rendered image The figure which shows the effective voxel area Vr A flowchart showing a specific process flow for generating a volume rendered image Flowchart showing the flow of ray casting process The figure which shows the state of calculating the intersection of the effective voxel region Vr and the line-of-sight vector. The figure which shows the state of calculating the intersection of the effective voxel region Vr and the line-of-sight vector. Flowchart showing the flow of the starting coordinate search process Flowchart showing the flow of voxel α value acquisition processing (without interpolation) Flowchart showing the flow of voxel α value acquisition processing (with interpolation) A flowchart showing the overall flow of the process of generating a MIP image The figure which shows the effective voxel area Vr A flowchart showing a specific processing flow for generating a MIP image Flowchart showing the flow of ray casting process Flowchart showing the flow of the starting coordinate search process Flowchart showing the flow of voxel signal value acquisition processing (without interpolation) Flowchart showing the flow of voxel signal value acquisition processing (with interpolation) Flowchart showing the flow of processing to generate an MPR image The figure explaining the MPR image (body axis cross section, coronal section, sagittal section) Display example of volume rendered image of chest and abdomen Display example of MIP image of chest and abdomen Display example of MPR image of chest and abdomen (body axis cross section, coronal section, sagittal section) Comparison display example of volume rendered image of chest and abdomen

Hereinafter, embodiments of the present disclosure will be described in detail with reference to the drawings.

(1. Functional configuration of mask generator 1)
FIG. 1 is a diagram showing a functional configuration of the mask generation device 1. As shown in the figure, the mask generation device 1 includes an image acquisition unit 21, a drawing weight table creation unit 22, and a mask data generation unit 23.

The image acquisition unit 21 acquires the tomographic image group Do taken by the CT apparatus. The tomographic image group Do is a plurality of tomographic images (CT images) in which a subject (human body) is continuously photographed at predetermined intervals along the cranio-caudal axis. Each tomographic image is two-dimensional image data in DICOM format. The DICOM format is an image format generally used for medical images including a header part and an image data part in one file, and parameters and diagnostic information at the time of image taking can be saved.

One tomographic image is, for example, a 512 × 512 pixel image (CT image). A signal value v is assigned to each pixel of the tomographic image, and in the case of a CT image, the signal value v is a CT value. In the present embodiment, the signal value v (CT value) is 16 bits (-32768 ≦ v ≦ 32767) data (however, the number of bits of the signal value v is not particularly limited).

The CT value is an X-ray attenuation coefficient of tissue expressed with water as a reference, and the type of tissue or lesion can be determined from the CT value (unit is HU (Hounsfield Unit)). The CT value is standardized by the X-ray attenuation coefficient of water and air, the CT value of water is 0, and the CT value of air is -1000. The CT value of fat is about −120 to -100, the CT value of normal tissue is about 0 to 120, and the CT value of bone is about 1000.

The tomographic image group Do is a stack of tomographic images, which are two-dimensional data of XY, in the Z-axis direction, and can be expressed as voxel data (three-dimensional voxels) that can be arranged in the voxel space R defined by the coordinate axes of XYZ. Is. For example, the tomographic image group Do (three-dimensional voxel) is defined as follows. In the present disclosure, the X-axis corresponds to the left-right axis of the human body, the Y-axis corresponds to the dorsoventral axis (front-back axis) of the human body, and the Z-axis corresponds to the cranio-caudal axis (upper and lower axis) of the human body. The XY axis is set, for example, at the time of photographing with a CT apparatus. Further, in the present disclosure, the X-axis direction may be referred to as a horizontal direction, and the Y-axis direction may be referred to as a vertical direction.

(Equation 1)
-32768 ≤ Do (x, y, z) ≤ 32767
0 ≦ x ≦ Sx-1, 0 ≦ y ≦ Sy-1, 0 ≦ z ≦ Sz-1
Resolution: Rxy, resolution Rz

In (Equation 1), Sx represents the number of pixels in the horizontal direction (number of voxels) of the tomographic image, Sy represents the number of pixels in the vertical direction (number of voxels) of the tomographic image, and Sz represents the number of slices (number of voxels) in the scanning direction. Arranged at predetermined intervals. Rxy is the resolution of horizontal and vertical tomographic images, and the horizontal and vertical resolutions are usually the same (different resolutions can be specified according to the DICOM standard, but there is no existing modality to scan at different resolutions). , The inverse number of the pixel spacing (Wxy), that is, the number of pixels per unit distance. Rz is the resolution of the slice and represents the reciprocal of the slice interval (Wz) (for example, 0.5 mm or 1 mm) of the tomographic image, that is, the number of slices per unit distance.

The signal value of the tomographic image group Do may be gradation-compressed to 8 bits as needed. The tomographic image group Do (x, y, z) whose gradation is compressed to 8 bits is defined as follows.

(Equation 2)
0 ≦ Do (x, y, z) ≦ 255
0 ≦ x ≦ Sx-1, 0 ≦ y ≦ Sy-1, 0 ≦ z ≦ Sz-1

The drawing weight table creation unit 22 creates a drawing weight table Sv that stores weight values for each pixel (x, y, z) of the tomographic image Do (z) for each tomographic image Do (z). .. The drawing weight table Sv is defined as follows.

(Equation 3)
0 ≦ Sv (x, y, z) ≦ 1
0 ≦ x ≦ Sx-1, 0 ≦ y ≦ Sy-1, 0 ≦ z ≦ Sz-1

As shown in FIG. 1, the drawing weight table creation unit 22 includes a subject area calculation unit 31, a subject area correction unit 32, a subject shape calculation unit 33, and a weight value calculation unit 34.

The subject area calculation unit 31 calculates the subject area for each tomographic image Do (z) based on the signal value of each pixel of the tomographic image Do (z).

The subject area correction unit 32 corrects the subject area calculated by the subject area calculation unit 31 in the case of a chest CT image.

The subject shape calculation unit 33 has a predetermined geometric shape (subject shape) of the subject area calculated by the subject area calculation unit 31 or the outline of the subject area corrected by the subject area correction unit 32 for each tomographic image Do (z). The parameter (geometric parameter P (z)) that defines the geometric shape is calculated by approximating with.

The weight value calculation unit 34 determines each pixel (x, y, z) of the tomographic image Do (z) based on the geometric parameter P (z) calculated by the subject shape calculation unit 33 for each tomographic image Do (z). ) Is calculated as the weight value Sv (x, y, z) for drawing, and a three-dimensional drawing weight table Sv storing the weight value Sv (x, y, z) is created.

Since the weight value Sv (x, y, z) is a three-dimensional data array and can be treated as it is as a data table (drawing weight table), the weight value and the drawing weight table are substantially the same. Therefore, the weight value and the drawing weight table are represented by the same reference numeral “Sv”.

The mask data generation unit 23 binarizes the weight values of each pixel (x, y, z) of the tomographic image group Do with a predetermined threshold value based on the drawing weight table Sv created by the drawing weight table creation unit 22. The mask value is calculated, and three-dimensional mask data Mask holding the mask value corresponding to each pixel (x, y, z) is generated. The mask data Mask is defined as follows.

(Equation 4)
Mask (x, y, z) = 0 (non-drawing) or 1 (drawing)
0 ≦ x ≦ Sx-1, 0 ≦ y ≦ Sy-1, 0 ≦ z ≦ Sz-1

(2. Hardware configuration of mask generator 1)
FIG. 2 is a block diagram showing a hardware configuration of the mask generation device 1 according to the present embodiment. As shown in the figure, the mask generation device 1 includes a control unit 11, a storage unit 12, a media input / output unit 13, a communication control unit 14, an input unit 15, a display unit 16, a peripheral device I / F unit 17, and the like. It is realized by a general-purpose computer connected via 18. However, the present invention is not limited to this, and various configurations can be adopted depending on the application and purpose.

The control unit 11 is a single processor or a processor composed of a multi-core CPU or a multi-processor composed of a plurality of processors (hereinafter, both are simply referred to as "CPU", and a core possessed by the CPU is referred to as a CPU core. ), ROM (Read Only Memory), RAM (Random Access Memory), frame memory (Frame Memory), etc. The CPU calls and executes a rendering program stored in the storage unit 12, ROM, recording medium, etc. in the work memory area on the RAM and frame memory, drives and controls each device connected via the bus 18, and masks. The processing described later performed by the generation device 1 is realized. In the present disclosure, processing using a plurality of CPU cores will be described later, but the plurality of CPU cores include not only physically a plurality of CPU cores (in the case of distinction, they are referred to as physical cores). , Logically, a plurality of CPU cores (referred to as logical cores when described separately) are also included.

The ROM is a non-volatile memory, and permanently holds a computer boot program, a program such as a BIOS, data, and the like. The RAM and frame memory are volatile memories, and are work areas used by the control unit 11 to perform various processes while temporarily holding programs, data, etc. loaded from the storage unit 12, ROM, recording medium, and the like. To be equipped.

The storage unit 12 is an HDD (Hard Disk Drive) or the like, and stores a program executed by the control unit 11, data necessary for program execution, an OS (Operating System), and the like. As for the program, a control program corresponding to the OS and an application program for causing a computer to execute a process described later are stored. Each of these program codes is read by the control unit 11 as necessary, transferred to the RAM and frame memory, read by the CPU, and executed as various means.

The media input / output unit 13 (drive device) inputs / outputs data, and for example, media such as a CD drive (-ROM, -R, -RW, etc.) and a DVD drive (-ROM, -R, -RW, etc.). It has an input / output device. The communication control unit 14 has a communication control device, a communication port, and the like, and is a communication interface that mediates communication between a computer and a network, and controls communication between other computers via the network. The network can be wired or wireless.

The input unit 15 inputs data and has, for example, a pointing device such as a keyboard and a mouse, and an input device such as a numeric keypad. Operation instructions, operation instructions, data input, and the like can be given to the computer via the input unit 15.
The display unit 16 has a display device such as a liquid crystal panel, a logic circuit (video adapter or the like) for realizing a video function of a computer in cooperation with the display device. The input unit 15 and the display unit 16 may be integrated as in the touch panel display.

The peripheral device I / F (Interface) unit 17 is a port for connecting the peripheral device to the computer, and the computer transmits / receives data to / from the peripheral device via the peripheral device I / F unit 17. The peripheral device I / F unit 17 is composed of USB (Universal Serial Bus), IEEE1394, RS-232C, or the like, and usually has a plurality of peripheral device I / Fs. The connection form with peripheral devices may be wired or wireless. The bus 18 is a route that mediates the transfer of control signals, data signals, and the like between the devices.

(3. Operation of mask generator 1)
The operation of the mask generation device 1 will be described with reference to the flowchart of FIG.

First, the control unit 11 of the mask generation device 1 acquires the tomographic image group Do captured by the CT device (step S1).

Subsequently, the control unit 11 stores weight values Sv (x, y, z) for drawing each pixel (x, y, z) of the tomographic image Do (z) for each tomographic image Do (z). The drawn drawing weight table Sv is created (step S2). The process of creating the drawing weight table Sv will be described later (see “4. Creating the drawing weight table”).

Then, the control unit 11 refers to the drawing weight table Sv created in step S2, and determines the weight values Sv (x, y, z) of each pixel (x, y, z) of the tomographic image group Do. The mask value Mask (x, y, z) binarized by the threshold value is calculated, and the three-dimensional mask data Mask holding the mask value corresponding to each pixel (x, y, z) is generated (step S3). .. The process of generating the mask data Mask will be described later (see “5. Mask data generation”).

(4. Creation of drawing weight table)
The process of creating the drawing weight table Sv executed in step S2 of FIG. 3 will be described with reference to FIGS. 4 to 16.
FIG. 4 is a flowchart showing a flow of processing for creating the drawing weight table Sv.

<Subject area calculation process>
First, the control unit 11 sets the rectangular region 41 (first rectangular region) or the rectangular region 42 as the subject region for each tomographic image Do (z) based on the signal value of each pixel of the tomographic image Do (z). (Second rectangular area) is calculated (step S21).

As shown in FIG. 5, the rectangular area 41 (first rectangular area) is a rectangular area surrounding the outermost part of the subject candidate area 50. The subject candidate region 50 is a pixel region of the tomographic image Do (z) whose signal value is larger than a predetermined threshold value (for example, a signal value that is not air or more). Since the present inventor has already proposed a method for calculating the rectangular area 41 as the subject area in Japanese Patent Application No. 2018-103767, here, the rectangular area 42 (second rectangular area) is mainly used as the subject area. The method of calculating is described.

The rectangular area 42 (second rectangular area) is a rectangular area defined inside the rectangular area 41 as shown in FIG. Since the rectangular region 41 includes all pixels whose signal values are larger than a predetermined threshold value, as shown in FIG. 5, an unnecessary region such as a bed 50b or a headrest other than the human body 50a may be included. Therefore, in the present disclosure, by calculating the rectangular area 42 defined inside the rectangular area 41 in consideration of the signal value distribution of the subject candidate area 50, unnecessary areas such as a bed and a headrest are not included. Extract the subject area.

The process of calculating the subject area (rectangular area 42) executed in step 21 of FIG. 4 will be specifically described with reference to FIGS. 6 to 14.
FIG. 6 is a flowchart showing a flow of processing for calculating the subject area (rectangular area 42). When the signal value of each pixel of the tomographic image Do (z) is 16 bits, it is converted into an 8-bit signal value in advance as follows.

(Equation 5)
D8 (x, y, z) = {Do (x, y, z) +700} · 255/1724
However, when D8 (x, y, z) <0, D8 (x, y, z) = 0, and when D8 (x, y, z)> 255, D8 (x, y, z) = 255. To do.

First, the control unit 11 calculates the coordinates of the center of gravity weighted by the signal values of the pixels of the subject candidate region 50 as the center coordinates of the rectangular region 42 (step S31, FIG. 7A).

Specifically, the control unit 11 sets the center coordinates of the rectangular region 42 to be calculated as (C (x), Cy (z)) and the sum of the weights as Ws, and Cx (z) = Cy (z) = Ws. Initialized to = 0, D8 (x, y) for each pixel (x, y) (0 ≦ x ≦ Sx-1, 0 ≦ y ≦ Sy-1) of the tomographic image D8 (x, y, z). , Z)> Vmin (for example, Vmin = 10), it is assumed that the pixels are in the subject candidate region 50, and Cx (z), Cy (z), and Ws are updated as follows.

(Equation 6)
As w = {D8 (x, y, z) -Vmin} ² ,
Cx (z) ← Cx (z) + x · w
Cy (z) ← Cy (z) + y ・ w
Ws ← Ws + w

The control unit 11 executes the above processing on Sx × Sy pixels in the range of 0 ≦ x ≦ Sx-1 and 0 ≦ y ≦ Sy-1, and then when Ws> 1,

(Equation 7)
Cx (z) = Cx (z) / Ws
Cy (z) = Cy (z) / Ws

And.
As described above, the coordinates of the center of gravity weighted by the signal values of the pixels of the subject candidate region 50 are obtained, and the coordinates of the center of gravity are calculated as the center coordinates (Cx (z), Cy (z)) of the rectangular region 42.

Subsequently, the control unit 11 determines the signal value within the range of the rectangular region 41 (first rectangular region) in the four directions of up, down, left, and right starting from the calculated center coordinates (Cx (z), Cy (z)). The center of gravity (hereinafter, referred to as “second center of gravity”) weighted by the signal value of the pixels larger than the threshold value (pixels in the subject candidate region 50) is calculated (step S32, FIG. 7B).

That is, as shown in FIGS. 8A and 8B, the control unit 11 sets the subject candidate region 50 to X with the center coordinates (Cx (z), Cy (z)) (center of gravity of the subject candidate region 50) as the starting point. The X coordinates Gx1 and Gx2 (Gx1 <Gx2) of the center of gravity (second center of gravity) weighted by the signal values of the pixels in each of the divided subject candidate regions 50 are calculated by dividing into two in the axial direction.

Similarly, as shown in FIGS. 8C and 8D, the control unit 11 sets the subject candidate area 50 from the center coordinates (Cx (z), Cy (z)) (center of gravity of the subject candidate area 50). The Y coordinates Gy1 and Gy2 (Gy1 <Gy2) of the center of gravity (second center of gravity) weighted by the signal values of the pixels in each of the divided subject candidate regions 50 are calculated by dividing into two in the Y-axis direction.

The specific process for obtaining the second center of gravity is as follows.
First, the control unit 11 sets the sum of the weights for each second barycentric coordinate to Wx1, Wx2, Wy1, and Wy2, and initializes them to Gx1 = Gx2 = Gy1 = Gy2 = Wx1 = Wx2 = Wy1 = Wy2 = 0.

Then, the control unit 11 has D8 (x, y, z) for each pixel (x, y) (0 ≦ x ≦ Sx-1, 0 ≦ y ≦ Sy-1) of the tomographic image D8 (x, y, z). When y, z)> Vmin (for example, Vmin = 10) is satisfied, it is assumed that the pixels are in the subject candidate region 50, and Gx1, Gx2, Gy1, Gy2, Wx1, Wx2, Wy1, and Wy2 are updated as follows.

(Equation 8)
As w = {D8 (x, y, z) -Vmin} ² ,
When x <Cx (z), Gx1 ← Gx1 + x · w, Wx1 ← Wx1 + w
When x ≧ Cx (z), Gx2 ← Gx2 + x · w, Wx2 ← Wx2 + w
When y <Cy (z), Gy1 ← Gy1 + y · w, Wy1 ← Wy1 + w
When y ≧ Cy (z), Gy2 ← Gy2 + y · w, Wy2 ← Wy2 + w

The control unit 11 executes the above processing on Sx × Sy pixels in the range of 0 ≦ x ≦ Sx-1 and 0 ≦ y ≦ Sy-1, and then performs the above processing.

(Equation 9)
Gx1 = Gx1 / Wx1
Gx2 = Gx2 / Wx2
Gy1 = Gy1 / Wy1
Gy2 = Gy2 / Wy2

Then, the second barycentric coordinates Gx1, Gx2 (Gx1 <Gx2), Gy1, and Gy2 (Gy1 <Gy2) in each of the divided subject candidate regions 50 are obtained.

Subsequently, the control unit 11 is a pixel whose signal value is larger than a predetermined threshold value in the range of the rectangular region 41 (first rectangular region) in four directions of up, down, left, and right, starting from the second barycentric coordinate calculated in step S32. The center of gravity weighted by the signal value of (pixels in the subject candidate region 50) (hereinafter, referred to as “third center of gravity”) is calculated (step S33, FIG. 7 (c)).

That is, as shown in FIGS. 9A and 9B, the control unit 11 divides the subject candidate area 50 into three in the X-axis direction starting from the second centers of gravity Gx1 and Gx2, and two subject candidate areas at the left and right ends. The X coordinates Gtx1 and Gtx2 (Gtx1 <Gtx2) of the center of gravity (third center of gravity) weighted by the signal value of the pixel at 50 are calculated.

Similarly, as shown in FIGS. 9C and 9D, the control unit 11 divides the subject candidate area 50 into three in the Y-axis direction starting from the second centers of gravity Gy1 and Gy2, and two subject candidates at the upper and lower ends. The Y coordinates Gty1 and Gty2 (Gty1 <Gty2) of the center of gravity (third center of gravity) weighted by the signal values of the pixels in the region 50 are calculated.

The specific process for obtaining the third center of gravity is as follows.
First, the control unit 11 sets the sum of the weights for each third barycentric coordinate to Wx1, Wx2, Wy1, and Wy2, and initializes them to Gtx1 = Gtx2 = Gty1 = Gty2 = Wx1 = Wx2 = Wy1 = Wy2 = 0.

Then, the control unit 11 has D8 (x, y, z) for each pixel (x, y) (0 ≦ x ≦ Sx-1, 0 ≦ y ≦ Sy-1) of the tomographic image D8 (x, y, z). When y, z)> Vmin (for example, Vmin = 10) is satisfied, it is assumed that the pixels are in the subject candidate region 50, and Gtx1, Gtx2, Gty1, Gty2, Wx1, Wx2, Wy1, and Wy2 are updated as follows.

(Equation 10)
As w = {D8 (x, y, z) -Vmin} ² ,
When x <Gx1, Gtx1 ← Gtx1 + x · w, Wx1 ← Wx1 + w
When x> Gx2, Gtx2 ← Gtx2 + x · w, Wx2 ← Wx2 + w
When y <Gy1, Gty1 ← Gty1 + y · w, Wy1 ← Wy1 + w
When y> Gy2, Gty2 ← Gty2 + y · w, Wy2 ← Wy2 + w

The control unit 11 executes the above processing on Sx × Sy pixels in the range of 0 ≦ x ≦ Sx-1 and 0 ≦ y ≦ Sy-1, and then performs the above processing.

(Equation 11)
Gtx1 = Gtx1 / Wx1
Gtx2 = Gtx2 / Wx2
Gty1 = Gty1 / Wy1
Gty2 = Gty2 / Wy2

Then, the third barycentric coordinates Gtx1, Gtx2 (Gtx1 <Gtx2), Gty1 and Gty2 (Gty1 <Gty2) in each of the divided subject candidate regions 50 are obtained.

Then, the control unit 11 calculates the horizontal size W (z) and the vertical size H (z) of the rectangular region 42 based on the third barycentric coordinates calculated in step S33 as follows. (Step S34, FIG. 7 (d)).

(Equation 12)
Horizontal size: W (z) = Gtx2-Gtx1
Vertical size: H (z) = Gty2-Gty1

Further, the control unit 11 calculates two types of the horizontal size and the vertical size starting from the center coordinates (Cx (z), Cy (z)) as follows.

(Equation 13)
First lateral size: W1 (z) = Cx (z) -Gtx1
Second lateral size: W2 (z) = Gtx2-Cx (z)
First vertical size: H1 (z) = Cy (z) -Gty1
Second vertical size: H2 (z) = Gty2-Cy (z)

FIG. 10 is a diagram showing the calculated rectangular region 42. As shown in the figure, the lateral size W (z) of the rectangular region 42 corresponds to the distance from the third center of gravity Gtx2 in the right direction to the third center of gravity Gtx1 in the left direction, and is the vertical size of the rectangular region 42. H (z) corresponds to the distance from the downward third center of gravity Gty2 to the upward third center of gravity Gty1.

Further, the first lateral size W1 (z) corresponds to the distance from the third center of gravity Gtx1 in the left direction to the center coordinate of the rectangular region 42, and the second lateral size W2 (z) is the right. It corresponds to the distance from the third center of gravity Gtx2 in the direction to the center coordinates of the rectangular region 42. The first vertical size H1 (z) corresponds to the distance from the upper third center of gravity Gty1 to the center coordinates of the rectangular region 42, and the second vertical size H2 (z) is downward. It corresponds to the distance from the third center of gravity Gty2 to the center coordinates of the rectangular region 42.

By the above processing, the subject area has center coordinates (Cx (z), Cy (z)), horizontal size W (z) (first horizontal size W1 (z), second horizontal size W1 (z)). The rectangular area 42 defined by the size W2 (z)), the vertical size H (z) (the first vertical size H1 (z), the second vertical size H2 (z)) is calculated. To.

In step S34, the control unit 11 determines the lateral size W (z) of the rectangular region 42 (first lateral size W1 (z), second) based on the second center of gravity calculated in step S32. Horizontal size W2 (z)) and vertical size H (z) (first vertical size H1 (z), second vertical size H2 (z)) may be calculated. ..

Further, the control unit 11 may further calculate the fourth center of gravity in four directions of up, down, left, and right, starting from the third center of gravity calculated in step S33. In this case, the control unit 11 has a lateral size W (z) of the rectangular region 42 based on the fourth center of gravity (first lateral size W1 (z), second lateral size W2 (z)). ) And the vertical size H (z) (the first vertical size H1 (z), the second vertical size H2 (z)) are calculated.

<Subject area correction processing>
By the way, in the case of chest CT, if the vertical and horizontal sizes of the rectangular region 42 are calculated from the second center of gravity and the third center of gravity in the vertical and horizontal directions from the center point, the rectangular region depends on the signal value distribution of the organs in the bone region. The size of 42 may not be suitable for removing the bone region and visualizing the organs and the like in the bone region.

For example, when the signal value inside the region is relatively large (for example, in the case of the abdominal liver region) as shown in FIG. 11A, the calculated center of gravity (second center of gravity, third center of gravity) tends to move toward the center. , The rectangular area 42 is calculated to be smaller. That is, the rectangular region 42 may be defined inside the bone region. When mask data is created based on such a rectangular area 42, the built-in area to be drawn on the main body may be removed.

Further, when the signal value inside the region is relatively small as shown in FIG. 11B (for example, in the case of the thoracic lung region), the calculated center of gravity tends to move outward, so that the rectangular region 42 is calculated to be large. Rectangle. That is, the rectangular region 42 may be defined outside the bone region. When mask data is created based on such a rectangular region 42, the bone region that is a non-drawing target may not be removed accurately.

Therefore, in order to solve the above problem, as shown in FIG. 12, the calculated horizontal and vertical sizes of the rectangular area 42 are multiplied by the area ratio of the subject candidate area 50 and the rectangular area 42 as a magnification. By doing so, the size of the rectangular area 42 is corrected.

As a result, when the signal value inside the region is relatively large as shown in FIG. 11A (for example, in the case of the abdominal liver region), the rectangular region 42 is enlarged and corrected, and the inside of the region is corrected as shown in FIG. 11B. When the signal value of is relatively small (for example, in the case of the thoracic lung region), the rectangular region 42 is reduced and corrected, and the rectangular region 42 is a size suitable for removing the bone region and visualizing the organs in the bone region. Is corrected to.

However, when the above correction is applied, another problem arises in which the originally unnecessary rectangular region 42 (subject region) of the upper chest is enlarged and corrected.
FIG. 13 shows the case where the above correction is applied to the upper chest. As shown in the figure, the upper chest is not surrounded by the bone region, but the high-intensity bone region is distributed independently. Therefore, the area ratio of the subject candidate region 50 to the rectangular region 42 becomes as large as the abdomen, and the enlargement correction is performed in the same manner as the abdomen. However, since the upper chest is mostly occupied by the bone region that is not drawn, it is not desirable that the rectangular region 42 (subject region) is enlarged in the upper chest region.

By the way, since the bone region of the upper chest extends in the lateral direction, the rectangular region 42 has a characteristic of being flattened. Therefore, in the present disclosure, as shown in FIG. 14, the aspect ratio of the rectangular region 42 (vertical size H (z) / horizontal size W (z)) is calculated, and this value is a predetermined threshold value (for example,). If it is smaller than 0.6), it is identified as the upper chest, and instead of multiplying the horizontal and vertical sizes of the rectangular area 42 by the area ratio of the rectangular area 42 and the subject candidate area 50, the vertical and horizontal directions are used. Correct by multiplying the ratio (reduction correction). On the other hand, if the aspect ratio of the rectangular area 42 is equal to or greater than a predetermined threshold value (for example, 0.6), it is determined that the rectangular area 42 is the chest or abdomen, and the rectangular area 42 is relative to the horizontal and vertical sizes of the rectangular area 42. And the area ratio of the subject candidate area 50 are multiplied to correct.

As described above, in the present disclosure, the region (upper chest / chest or abdomen) of the subject region is determined based on the aspect ratio (vertical size / horizontal size) of the rectangular region 42, and the correction method according to the region is determined. Corrects the rectangular area 42 by.

FIG. 15 is a flowchart showing the flow of the correction process of the subject area in step S22.
First, the control unit 11 calculates the aspect ratio Rvh of the rectangular region 42 as follows (step S41).

(Equation 14)
Aspect ratio Rvh = H (z) / W (z)

When the threshold value for the aspect ratio is Srvh (for example, 0.62) and Rvh <Srvh (step S42; No), the control unit 11 determines that it is the upper chest and corrects the rectangular region 42 as follows ( Reduction correction) (correction in step S43, FIG. 14).

(Equation 15)
Horizontal size: W (z) = W (z) · Rvh
Vertical size: H (z) = H (z) · Rvh

(Equation 16)
First lateral size: W1 (z) = W1 (z) · Rvh
Second lateral size: W2 (z) = W2 (z) · Rvh
First vertical size: H1 (z) = H1 (z) · Rvh
Second vertical size: H2 (z) = H2 (z) · Rvh

On the other hand, when Rvh ≧ Srvh (step S42; Yes), the control unit 11
Judging that it is the chest or abdomen, first, for each pixel (x, y) (0 ≦ x ≦ Sx-1, 0 ≦ y ≦ Sy-1) of the tomographic image D8 (x, y, z), The number of pixels satisfying D8 (x, y, z)> Vmin (for example, Vmin = 10) is counted, and the count value Cobj is calculated as the area of the subject candidate region 50 (step S44).

Subsequently, the control unit 11 sets the area ratio Robj with the rectangular region 42.

(Equation 17)
Robj = Cobj / (H (z) · W (z))

In the horizontal direction and the vertical direction, the magnifications Rox and Roy are set to Rox = Roy = Robj (step S45).

The upper limit of the variable magnification is set to Mox and Moy, Rox = Mox when Rox> Mox, Roy = Moy when Roy> Moy, and the control unit 11 corrects the rectangular region 42 as follows (step). S46, correction of FIG. 12).

(Equation 18)
Horizontal size: W (z) = W (z) · Rox
Vertical size: H (z) = H (z) · Roy

(Equation 19)
First lateral size: W1 (z) = W1 (z) · Rox
Second horizontal size: W2 (z) = W2 (z) · Rox
First vertical size: H1 (z) = H1 (z) · Roy
Second vertical size: H2 (z) = H2 (z) · Roy

<Subject shape calculation process>
Subsequently, the control unit 11 determines the outline of the subject area (rectangular area 42) calculated in step S21 or the subject area (rectangular area 42) corrected in step S22 for each tomographic image Do (z). Approximate with a geometric shape (subject shape) and calculate a parameter (geometric parameter P (z)) that defines the geometric shape (step S23).

In the present disclosure, each of the following regions (4 quadrants) is approximated by a special ellipse (hereinafter, also referred to as "deformed ellipse") in which different diameters can be set, and each of the deformed ellipses is configured as a geometric parameter P (z). The parameters (ellipse parameters P1 (z) to P4 (z)) that define the ellipse of the region (4 quadrants) are calculated.
(Region 1) x <Cx (z) and y <Cy (z)
(Region 2) x ≧ Cx (z) and y <Cy (z)
(Region 3) x <Cx (z) and y ≧ Cy (z)
(Region 4) x ≧ Cx (z) and y ≧ Cy (z)

(Equation 20)
(Ellipse parameter P1 (z) in region 1)
Center coordinates of ellipse Cx (z), Cy (z)
Horizontal division size of ellipse rx1 = W1 (z) · a1 (horizontal diameter)
Vertically divided ellipse size ry1 = H1 (z) · b1 (vertical diameter)
(Ellipse parameter P2 (z) in region 2)
Center coordinates of ellipse Cx (z), Cy (z)
Horizontal division size of ellipse rx2 = W2 (z) ・ a2 (horizontal diameter)
Vertically divided ellipse size ry1 = H1 (z) · b1 (vertical diameter)
(Ellipse parameter P3 (z) in region 3)
Center coordinates of ellipse Cx (z), Cy (z)
Horizontal division size of ellipse rx1 = W1 (z) · a1 (horizontal diameter)
Ellipse vertical division size ry2 = H2 (z) · b2 (vertical diameter)
(Ellipse parameter P4 (z) in region 4)
Center coordinates of ellipse Cx (z), Cy (z)
Horizontal division size of ellipse rx2 = W2 (z) ・ a2 (horizontal diameter)
Ellipse vertical division size ry2 = H2 (z) · b2 (vertical diameter)

a1, a2, b1 and b2 are correction coefficients corresponding to each of the vertical and horizontal diameters W1 (z), W2 (z), H1 (z) and H2 (z) of the ellipse (1.0 when not corrected by real values). ), Left direction (a1, W1 (z)), right direction (a2, W2 (z)), front direction (b1, H1 (z)), back direction (b2, H2 (z)) Set separately.

For example, it is set as follows.
-For chest CT: a1 = a2 = 0.95, b1 = 1.05, b2 = 0.55
-For head CT: a1 = a2 = 1.15, b1 = 1.2, b2 = 1.25

<Weight value calculation process>
Then, the control unit 11 sets each pixel (for each tomographic image Do (z), based on the geometric parameters P (z) (ellipse parameters P1 (z) to P4 (z)) of the deformed ellipse calculated in step S23. A weight value Sv (x, y, z) for drawing x, y, z) is calculated, and a drawing weight table Sv is created (step S24).
For example, the weight value Sv (x, y, z) of each pixel (x, y, z) of the z-th tomographic image Do (z) is calculated as follows.

(Equation 21)
Sv (x, y, z) = Svxy (x, y, z) · Svz (z)
(0 ≦ Sv (x, y, z), Svxy (x, y, z), Svz (z) ≦ 1)
Svxy = 1-k · r (x, y, z)
(Region 1) When x <Cx (z) and y <Cy (z):
r (x, y, z) = {(x-Cx (z)) / rx1} ² + {(y-Cy (z)) / ry1} ²
(Region 2) When x ≧ Cx (z) and y <Cy (z):
r (x, y, z) = {(x-Cx (z)) / rx2} ² + {(y-Cy (z)) / ry1} ²
(Region 3) When x <Cx (z) and y ≧ Cy (z):
r (x, y, z) = {(x-Cx (z)) / rx1} ² + {(y-Cy (z)) / ry2} ²
(Region 4) When x ≧ Cx (z) and y ≧ Cy (z):
r (x, y, z) = {(x-Cx (z)) / rx2} ² + {(y-Cy (z)) / ry2} ²
Svz (z) = 1.0-kz · {(z-Sz / 2) / rz} ²
rz = Sz · c / 2

k is the attenuation coefficient in the XY plane and kz is the attenuation coefficient in the Z direction. In the case of chest CT, k = 1.0, kz = 0.0, and in the case of head CT, k = kz = 1.0.

c is a correction coefficient (1.0 when not corrected by a real value) corresponding to the Z direction width Sz / 2 set in the case of head CT, and for example, c = 0.9.

FIG. 16 is a diagram showing a two-dimensional distribution pattern of weight values Sv (x, y, z) calculated by Equation 21. As shown in the figure, the weight value Sv at the center of the deformed ellipse (Cx (z), Cy (z)) becomes 1 (maximum), and the weight value Sv inside the deformed ellipse becomes smaller as the distance from the center increases. Specifically, inside the deformed ellipse, the weight value Sv on each concentric ellipse having the same horizontal diameter, vertical diameter ratio, and center coordinate of each ellipse constituting the deformed ellipse is the horizontal diameter and vertical diameter of each concentric ellipse. The larger the diameter, the smaller it becomes (it becomes smaller in inverse proportion to the square of the diameter). The weight value Sv outside the deformed ellipse is uniformly set to 0.

Further, in the case of head CT, due to Svz (z), as the slice position of the tomographic image Do (z) is located from the center to the end (as the difference between Sz / 2 and the slice rank z increases), the weight is increased. The value Sv (x, y, z) is uniformly attenuated (Sv (x, y, z) is less than 1 even at the center of the deformed ellipse). Specifically, assuming that the slice order of the tomographic image is z, the weight value Sv (x, y, z) is uniformly attenuated in inverse proportion to the square of (z-Sz / 2) (Sz / 2).

The control unit 11 has weight values Sv (x, y, z) (0 ≦ x ≦ Sx-1, 0 ≦ y ≦ Sy-1, 0 ≦ z ≦ Sz-1) for all tomographic images Do (z). Is calculated, and a three-dimensional drawing weight table Sv that stores the calculated weight values Sv (x, y, z) is created.

The process of creating the drawing weight table Sv has been described above with reference to FIGS. 4 to 16.

(5. Generation of mask data)
The process of generating the mask data Mask executed in step S3 of FIG. 3 will be described with reference to the flowchart of FIG.
The process of FIG. 17 is executed for all pixels (x, y, z).

The control unit 11 acquires the weight value Sv (x, y, z) of each pixel (x, y, z) from the drawing weight table Sv, and the weight value Sv (x, y, z) · 255> threshold Th ( For example, for a pixel (x, y, z) satisfying Th = 10) (step S51; Yes), the mask value corresponding to the pixel (x, y, z) is set to 1.0 (Mask (x)). , Y, z) = 1.0, step S52), for a pixel (x, y, z) that does not satisfy the weight value Sv (x, y, z) · 255> 10 (step S51; No), the pixel. The mask value corresponding to (x, y, z) is set to 0 (Mask (x, y, z) = 0, step S53). As described above, the weight value Sv (x, y, z) of each pixel (x, y, z) is binarized with a predetermined threshold value, and the mask value Mask (x, y, z) of each pixel (x, y, z) is binarized. The mask data Mask holding y, z) is generated.

FIG. 18 is another example of the process of generating the mask data Mask.
The process of FIG. 18 is executed for all pixels (x, y, z).
First, the control unit 11 corrects the signal value by multiplying the 8-bit signal value D8 (x, y, z) of each pixel by the weight value Sv (x, y, z) obtained from the drawing weight table Sv. (D8'(x, y, z) = D8 (x, y, z) · Sv (x, y, z), step S61).

Then, the control unit 11 has a pixel (x, y, z) in which the corrected signal value D8'(x, y, z) satisfies D8'(x, y, z)> threshold Th (for example, Th = 10). ) (Step S62; Yes), the mask value corresponding to the pixel (x, y, z) is set to 1.0 (Mask (x, y, z) = 1.0, step S63), D8. '(X, y, z)> For a pixel (x, y, z) that does not satisfy the threshold Th, (step S62; No), the mask value corresponding to the pixel (x, y, z) is set to 0. (Mask (x, y, z) = 0, step S64).

(6. Parallel drawing weight table creation process and mask data generation process)
In the processes of steps S2 and S3 (drawing weight table creation process and mask data creation process) of FIG. 3 described above, as shown in FIG. 19, the processing area of the tomographic image group Do is equally divided in the Z-axis direction (in the figure). In the example, it is desirable to divide the data into four and then execute each processing area (Do1 to Do4) in parallel. As a result, mask data Mask can be generated at high speed.

FIG. 20 shows an example of a mechanism for executing the drawing weight table creation process and the mask data creation process in parallel. In the example of FIG. 20, the motherboard (lower part of the figure) mounted on the mask generator 1 (computer) has four CPUs (physical cores), and each CPU has two logical cores (multi-core CPU). And. That is, the mask generation device 1 includes substantially 8 logical cores (threads # 1 to # 8). First, when the mask generation program stored in the mask generation device 1 is started, it is assigned to thread # 1 of CPU core # 1 and executed by the job scheduler of the Windows multitasking operating system. Subsequently, the mask generation program includes drawing weight table creation processing and mask data creation processing (parallel processing thread # 1 in the figure) related to the processing area Do1, drawing weight table creation processing and mask data creation processing related to the processing area Do2 (FIG. Parallel processing thread # 2), drawing weight table creation processing and mask data creation processing related to the processing area Do3 (parallel processing thread # 3 in the figure), drawing weight table creation processing and mask data creation processing related to the processing area Do4 (FIG. Start each thread of the parallel processing thread # 4) and register it in the job scheduler of the Windows multitasking operating system.

In the Windows multitasking operating system, each registered parallel processing thread # 1 to # 4 is assigned to one of eight logical cores (threads # 1 to # 8) depending on the availability of cores, etc., and each parallel processing thread is assigned. Execute # 1 to # 4 in parallel. In the example of FIG. 20, the parallel processing thread # 1, the parallel processing thread # 2, and the parallel processing thread # 3 are the thread # 3 of the CPU core # 2 and the thread # of the CPU core # 3 whose physical cores are freed by the job scheduler. 5. Assigned to threads # 7 of CPU core # 4, respectively. The remaining concurrency thread # 4 is assigned to the logical thread # 2 of the CPU core # 1 although the physical core is not free. The rendering program body is in use in the logical thread # 1 of the CPU core # 1, but the mask generation program body starts four parallel processing threads and is in a waiting state for their termination. Therefore, the CPU core # 1. The parallel processing thread # 4 assigned to the logical thread # 2 of 1 is mainly executed.

The operation of the mask generation device 1 has been described above. The mask generation device 1 acquires a tomographic image taken by the CT device (step S1 in FIG. 3), and for each tomographic image Do (z), the subject area (rectangle) is based on the signal value of each pixel (x, y). Region 42) is calculated (see step S21 in FIG. 4, (Equation 7) (Equation 12) (Equation 13)), and the outline of the calculated subject region (rectangular region 42) is formed into a predetermined geometric shape (deformed ellipse). Calculate the geometric parameters P (z) (four-quadrant rectangle parameters P1 (z) to P4 (z)) that approximate and define the geometric shape (deformed ellipse) (see step S23 in FIG. 4, (Equation 20)). ). Subsequently, for each tomographic image Do (z), a weight value Sv (x, y, z) for drawing each pixel (x, y) of the tomographic image is calculated based on the geometric parameter P (z) ( See step S24 of FIG. 4, (Equation 21). Then, a mask value obtained by binarizing the weight values Sv (x, y, z) of each pixel (x, y, z) with a predetermined threshold value is calculated, and the mask value Mask of each pixel (x, y, z) is calculated. The mask data Mask holding (x, y, z) is generated (step S3 in FIG. 3).
This makes it possible to easily create mask data for visualizing a desired region based on a geometric shape that approximates the outline of the subject region.

Further, by calculating the subject area (rectangular area 42) in consideration of the signal value distribution of the subject candidate area 50 (see FIG. 6), mask data for making unnecessary areas such as a bed and a headrest other than the human body invisible is generated. can do.

Further, in the chest CT, a region (upper chest / chest or abdomen) is determined based on the aspect ratio of the subject region (rectangular region 42), and the subject region (rectangular region 42) is corrected by a correction method according to the region. As a result (see FIG. 15), the subject area can be calculated with high accuracy.
Furthermore, by approximating the outer shell of the subject area (rectangular area 42) with a deformed ellipse that can change the diameter of the ellipse on the sternum side and the vertebrae side (see FIG. 16), the bone area including the ribs and the sternum is accurate. It is possible to generate mask data in which the visceral regions such as lungs, heart, blood vessels, liver, and kidneys, which are important for diagnosis, are well visualized.

In the present disclosure, it is decided to create a drawing weight table Sv for storing weight values Sv (x, y, z) for each pixel (x, y, z) to be drawn, but the drawing weight table Sv is created. Is not required. That is, each time the weight value Sv (x, y, z) calculated in step S24 of FIG. 4 (see (Equation 21)) is processed, the mask value Mask (x, y, z) is processed by FIG. 17 or FIG. ) May be calculated. In this case, it is not necessary to secure the data area of the drawing weight table Sv in the memory.

Further, the calculation itself of the weight value Sv (x, y, z) represented by (Equation 21) is not essential. In the case of the embodiment shown in FIG. 17, when Sv (x, y, z) · 255> Th is satisfied, the mask value Mask (x, y, z) = 1. In chest CT, when Svxy (x, y, z) = 1-k · r (x, y, z)> Th / 255 is satisfied, the mask value Mask (x, y, z) = 1. This is the mask value Mask (x, y, z) if the voxel coordinates (x, y, z) are contained within a single deformed ellipse of any of the concentric deformed ellipses shown in FIG. It can be interpreted that = 1. That is, in the chest CT, the mask value Mask (x, y, z) defines whether it is inside (mask value: 1) or outside (mask value: 0) of any single deformed ellipse shown in FIG. It is possible to calculate the mask value Mask (x, y, z) by the inside / outside determination processing for a single deformed ellipse of voxel coordinates (x, y, z). Also in the head CT, the mask value Mask (x, y) is determined by the inside / outside determination processing for a single deformed ellipsoid defined by Sv (x, y, z) = Svxy (x, y, z) · Svz (z). , Z) can be calculated. However, especially in the case of head CT, a drawing weight table Sv is created based on (Equation 21) from the method of calculating the mask value Mask (x, y, z) by the inside / outside determination processing, and binarized by Th. The above-mentioned method for calculating the mask value Mask (x, y, z) has a smaller calculation load.

Further, in the present disclosure, in step S23 of FIG. 4, the case where the outer shell of the subject region (rectangular region 42) is approximated by a deformed ellipse (see FIG. 16) has been described, but the geometric shape to be approximated is not limited to the deformed ellipse. For example, as the present inventor has already proposed in Japanese Patent Application No. 2018-103767, the outer shell of the subject area (rectangular area 42) may be approximated by a single ellipse, or already proposed in Japanese Patent Application No. 2018-124125. As shown above, the outer shell of the subject area (rectangular area 42) may be approximated by a plurality of ellipses. Even when approximating with a plurality of ellipses, the mask value Mask (x, y, z) is obtained by performing internal / external judgment processing on the plurality of ellipses of voxel coordinates (x, y, z) and performing AND or OR calculation of the judgment result. Can be calculated.

Further, in step S21 of FIG. 4, the rectangular area 42 (second rectangular area, see FIG. 5) was calculated as the subject area, but the rectangular area 41 (first rectangular area, see FIG. 5) can also be calculated. Good. In this case, the subject area is not corrected in step S22 of FIG. 4, and in step S23 of FIG. 4, for example, the outer shell of the subject area is a single ellipse (see the previously proposed Japanese Patent Application No. 2018-103767) or a plurality of ellipses. (See the previously proposed Japanese Patent Application No. 2018-124125).

[Application example of mask data Mask]
Hereinafter, as an application example of the generated mask data Mask, a process of generating a three-dimensional reconstructed image using the mask data Mask will be described.

(Outline of 7.3D reconstruction image generation process)
FIG. 21 is a diagram showing an outline of a reconstruction image generation process executed by the three-dimensional reconstruction image generation device 2 that generates a three-dimensional reconstruction image. As shown in the figure, the three-dimensional reconstruction image generator 2 uses mask data Mask based on a plurality of tomographic images (the example of the figure is 370 chest CT images of 512 × 512 pixels). Volume-rendered image, MIP (Maximum Intensity Projection) image, MPR (Multi-Planar), which is generated and observed from a predetermined viewpoint.
Reconstruction) Generates a 3D reconstruction image such as an image.

(8. Functional configuration of the three-dimensional reconstruction image generator 2)
FIG. 22 is a diagram showing a functional configuration of the three-dimensional reconstruction image generation device 2. As shown in the figure, the three-dimensional reconstruction image generation device 2 includes an image acquisition unit 21, a drawing weight table creation unit 22, a mask data generation unit 23, a clipping area setting unit 24, an effective voxel area calculation unit 25, and a three-dimensional reconstruction. The configuration image generation unit 26 and the three-dimensional reconstruction image display unit 27 are provided. The image acquisition unit 21, the drawing weight table creation unit 22 (subject area calculation unit 31, subject area correction unit 32, subject shape calculation unit 33, weight value calculation unit 34), and the mask data generation unit 23 are the mask generation devices of FIG. Since it is the same as the functional configuration of No. 1, the description thereof will be omitted.

The clipping area setting unit 24 sets the clipping area ROI (interest area) of the subject. In the present disclosure, the clipping ROI is a rectangular parallelepiped, and is defined by the X-direction ROI (X-direction open section), the Y-direction ROI (Y-direction open section), and the Z-direction ROI (Z-direction open section) as follows. ..

(Equation 22)
X-axis direction ROI: Xs-Xe
Y-axis direction ROI: Ys-Ye
Z-axis direction ROI: Zs-Ze
When the clipping region is not set, Xs = 0, Xe = Sx-1, Ys = 0, Ye = Sy-1, Zs = 0, Ze = Sz-1.

The effective voxel area calculation unit 25 includes all voxels (hereinafter, may be referred to as “effective voxels”) that satisfy the predetermined conditions so that the voxel value defines the voxels to be drawn, and the effective voxels include. The tangent rectangular parallelepiped is calculated as the effective voxel region Vr. The effective voxel region Vr is defined by the X-axis direction ROI, the Y-axis direction ROI, and the Z-axis direction ROI as follows. Here, the "voxel value" is numerical information associated with each voxel, and is a signal value (CT value) measured by modality, an opacity or a color value obtained by passing the signal value through a color map Cmap, and the like. And so on.

(Equation 23)
X-axis direction ROI: Xis-Xie
Y-axis direction ROI: Yis-Yie
Z-axis direction ROI: Zis-Zie

By calculating the effective voxel area Vr, the process of generating a volume rendering image or MIP image (raycasting process) is executed by referring only to the voxels in the effective voxel area Vr, so that the calculation load is greatly reduced. Since the boundary is a rectangular parallelepiped, the starting point coordinates in the ray casting process can be quickly searched based on the intersection with the straight line (virtual ray).

The three-dimensional reconstruction image generation unit 26 generates a three-dimensional reconstruction image from the tomographic image group Do while referring to the mask data Mask. As shown in FIG. 22, the three-dimensional reconstructed image generation unit 26 includes a volume rendering image generation unit 51, a MIP image generation unit 52, and an MPR image generation unit 53.

The volume rendering image generation unit 51 sets a mask value (transparent / opaque) for the opacity of each voxel of the tomographic image group Do based on the mask data Mask in the ray casting process described later (FIG. 33 described later). In step S144), the color map Cmap set by the color map setting unit 510 is applied to each voxel of the tomographic image group Do (step S145 in FIG. 33 described later) to generate a volume rendered image ImageVR.

The color map Cmap defines the correspondence between the signal value v and the color value (specifically, the RGB value) and the opacity (α value), and the signal value v is the 24-bit color (RGB value) and It is expressed as a function (substantially a two-dimensional data table) that converts to 8-bit opacity (α value). For example, the color map Cmap applied to the 16-bit tomographic image group Do (see (Equation 1)) is defined as follows.

(Equation 24)
0 ≤ Cmap (v, n) ≤ 255
-32768 ≤ v ≤ 32767
n = 0 (R), 1 (G), 2 (B), 3 (α)

The color map Cmap (v, n) applied to the 8-bit tomographic image group Do (see (Equation 2)) is defined as follows.

(Equation 25)
0 ≤ Cmap (v, n) ≤ 255
0 ≦ v ≦ 255
n = 0 (R), 1 (G), 2 (B), 3 (α)

Of the color maps Cmap (v, n) (0 ≦ n ≦ 3) shown in (Equation 24) and (Equation 25), the color map Cmap (v, n) (0 ≦ n ≦ 2) has a signal value v. Corresponds to a function that converts to a color value (RGB value).

Further, among the color maps Cmap (v, n) (0 ≦ n ≦ 3) shown in (Equation 24) and (Equation 25), in particular, the color map Cmap (v, 3) converts the signal value v into opacity. Corresponds to the function that does, and is generally called the "opacity curve".

In the color map Cmap (v, 3) (opacity curve), for example, the opacity of the normal structure (16-bit signal value v = 0 to 120) is set to 255 (opacity = 255 is opaque (all light is reflected). ), And the opacity of the bone region (16-bit signal value v = around 1000) is set to an intermediate value between 1 and 254 (opacity = 1 to 254 is translucent (part of the incident light). Is reflected, and others are transmitted), and the opacity of other parts such as tissues is set to 0 (opacity = 0 indicates that all of the incident light is transmitted). Will be done.

The generated volume-rendered image ImageVR is a full-color (24-bit) three-dimensional reconstruction image (Size × Size image), and is defined as follows.

(Equation 26)
0 ≤ Image VR (x, y, n) ≤ 255
0 ≦ x ≦ Size-1, 0 ≦ y ≦ Size-1
n = 0 (R), 1 (G), 2 (B)

In the ray casting process described later, the MIP image generation unit 52 sets the mask value for each voxel signal value of the tomographic image group Do based on the mask data Mask (step S224 in FIG. 40 described later), and then the tomographic image. A signal value is acquired for each voxel of the group Do (step S225 in FIG. 40 described later), and a MIP image ImageMIP is generated.

The generated MIP image ImageMIP is a monochrome (16-bit or 8-bit) three-dimensional reconstruction image (Size × Size image), and is defined as follows. At the stage of generating the MIP image ImageMIP, it is common to calculate the gradation in monochrome (16 bits) with the signal value of the original DICOM image as it is, and the ImageMIP (x, y) is -32768 ≤ ImageMIP (x, y) It may have a value of ≦ 32767, but at the stage of being displayed, it is converted to monochrome (8 bits) as follows.

(Equation 27)
0 ≤ Image MIP (x, y) ≤ 255
0 ≦ x ≦ Size-1, 0 ≦ y ≦ Size-1

The MPR image generation unit 53 generates an MPR image ImageMPR having an arbitrary cross section based on the tomographic image group Do, and synthesizes a mask image generated based on the mask data Mask into the MPR image ImageMPR.

The three-dimensional reconstruction image display unit 27 displays a three-dimensional reconstruction image generated by the three-dimensional reconstruction image generation unit 26 (volume rendering image generation unit 51, MIP image generation unit 52, MPR image generation unit 53). It is displayed on (display unit 16).

(Hardware configuration of 9.3-dimensional reconstruction image generator 2)
FIG. 23 is a diagram showing a hardware configuration of the three-dimensional reconstruction image generation device 2. As shown in the figure, the three-dimensional reconstruction image generation device 2 includes a control unit 11, a storage unit 12, a media input / output unit 13, a communication control unit 14, an input unit 15, a display unit 16, and a peripheral device I / F unit 17. Etc. are realized by a general-purpose computer connected via the bus 18. Since each configuration is the same as in FIG. 2, the description is omitted.

(Operation of 10.3 dimensional reconstruction image generator 2)
The overall operation of the three-dimensional reconstruction image generation device 2 will be described with reference to the flowchart of FIG. 24.

First, the control unit 11 of the three-dimensional reconstruction image generation device 2 acquires a 16-bit tomographic image group Do (see (Equation 1)) captured by the CT device (step 71). At this time, the control unit 11 may perform gradation compression of the 16-bit tomographic image group Do to the 8-bit tomographic image group Do (see (Equation 2)), if necessary.

(1) The minimum value Dmin and the maximum value Dmax of all the pixels in the Sz / second intermediate slice Do (z / 2) of the tomographic image group Do are calculated. Note that Dmin and Dmax may be determined by obtaining the minimum and maximum values from all the slices, but by adopting the method of determining only the intermediate slices as in the present embodiment, a large-capacity 16-bit tomographic image group. It is no longer necessary to keep Do in memory.

(2) Set the lower limit value Lmin = (Dmax-Dmin) · γ + Dmin and the upper limit value Lmax = (Dmax-Dmin) · (1-γ) + Dmin. Here, γ is the contrast adjustment width of the gradation compressed image, and the closer it is to 0, the higher the contrast (however, the brightness becomes smaller). Normally, it is set to γ = 0.1.

(3) The tomographic image group Do is converted into an 8-bit tomographic image group Do (see (Equation 2)) by the following formula (compressed in 256 steps in the range of the lower limit value Lmin to the upper limit value Lmax).

(Equation 28)
Do (x, y, z)
= (Do (x, y, z) -Lmin) · 255 / (Lmax-Lmin)
However, when Do (x, y, z)> 255, Do (x, y, z) = 255, and when Do (x, y, z) <0, Do (x, y, z) = 0. Saturate.

By performing gradation compression in this way, the memory capacity for holding the tomographic image group can be suppressed to half. Even if the gradation of the signal value is 16 bits, the gradation of the converted color value (RGB value) and opacity (α value) is limited to 8 bits according to the gradation of the display by the color map Cmap. Therefore, there is almost no deterioration in image quality due to gradation compression.

Subsequently, the control unit 11 stores weight values Sv (x, y, z) for drawing each pixel (x, y, z) of the tomographic image Do (z) for each tomographic image Do (z). The drawn drawing weight table Sv is created (step S72), and the weight value Sv (x, y, z) of each pixel (x, y, z) of the tomographic image group Do is referred to with reference to the created drawing weight table Sv. Is calculated as a mask value Mask (x, y, z) obtained by binarizing the above with a predetermined threshold value, and the mask value Mask (x, y, z) corresponding to each pixel (x, y, z) is held in three dimensions. Mask data Mask is generated (step S73). The processing of step S72 and step S73 is the same as that of step S2 and step S3 of FIG.

Then, the control unit 11 generates a three-dimensional reconstruction image from the tomographic image group Do (step S74) while referring to the mask data Mask generated in step S73, and displays the generated three-dimensional reconstruction image (step S74). Step S75).

Hereinafter, in step S74, an example of generating a volume rendering image, a MIP image, and an MPR image will be described in order.

(11. Generation of volume rendered image)
The process of generating the volume rendered image will be described with reference to FIGS. 25 to 33.

FIG. 25 is a flowchart showing the entire flow of the process of generating the volume rendered image.
The control unit 11 receives the color map Cmap setting from the user (step S81). Here, when the signal value of the tomographic image group Do is 16 bits (see (Equation 1)), a 16-bit compatible color map Cmap (see (Equation 24)) is set, and the signal value of the tomographic image group Do is set. When it is compressed to 8 bits (see (Equation 2)), a color map Cmap corresponding to 8 bits (see (Equation 25)) is set.

The method for setting the color map Cmap is not particularly limited. For example, the control unit 11 allows the user to select a color map Cmap to be applied to the tomographic image group Do from a plurality of color map Cmaps prepared in advance. Can be set. Further, the control unit 11 may read the color map Cmap used at the time of the previous rendering from the storage unit 12 and set it as the color map Cmap to be applied to the tomographic image group Do.

Further, the control unit 11 may set a color map Cmap created by the user on the color map setting screen (not shown). In this case, the control unit 11 displays a typical signal value v (for example, characteristic signal values at about 10 locations), a color value (RGB value), and an opacity (α) with respect to the signal value v on the color map setting screen. A color map Cmap that causes the user to set the value) and converts the signal value v in a predetermined range into the color value and the opacity based on the correspondence between the typical signal value v set by the user and the color value and the opacity. Automatically generated.

Subsequently, the control unit 11 receives a setting of the clipping region ROI (see (Equation 22)) from the user, if necessary (step S82).

Subsequently, the control unit 11 calculates the effective voxel region Vr (step S83). Specifically, as shown in FIG. 26, the control unit 11 uses voxels (opaque voxels) having a non-opaque value (α> 0) as effective voxels, and includes all of these effective voxels and is effective. The rectangular parallelepiped circumscribing the voxel is calculated as the effective voxel region Vr. When the clipping region ROI and the mask data Mask are set, the control unit 11 has a voxel (α> 0) having a non-zero opacity value (α> 0) within the drawing range defined by the clipping ROI and the mask data Mask. An opaque voxel) is used as an effective voxel, and a square body including all of these effective voxels and circumscribing the effective voxel is calculated as an effective voxel region Vr.

Then, the control unit 11 is set in the tomographic image group Do (three-dimensional voxel) acquired in step S71 (FIG. 24), the mask data Mask generated in step S73 (FIG. 24), and step S81 (FIG. 25). A volume rendered image ImageVR is generated with reference to the color map Cmap, the clipping region ROI set in step S82 (FIG. 25), and the effective voxel region Vr calculated in step S83 (FIG. 25) (step S83).

In particular, in the present disclosure, in the ray casting process, the control unit 11 sets the mask value for the opacity of each voxel of the tomographic image group Do based on the mask data Mask (step S144 in FIG. 33 described later), and then the fault. A color map Cmap (opacity curve) is applied to each voxel of the image group Do (step S145 in FIG. 33 described later) to generate a volume rendered image ImageVR. As a result, a volume rendered image reflecting the drawing / non-drawing information of the mask data Mask generated based on the subject shape (geometric shape) is generated.

FIG. 27 shows a specific flow of processing for generating the volume rendered image ImageVR, which is executed in step S84.
First, the control unit 11 sets the coordinate conversion parameters (step S91). In the present disclosure, in the ray casting process (step S93, FIG. 28) described later, since the voxels are referred to for each three-dimensional coordinate of the viewpoint coordinate system, a method of sequentially performing coordinate conversion is adopted. Therefore, the control unit 11 sets and calculates the coordinate conversion parameters common to each coordinate conversion process in advance as follows.

<Coordinate conversion parameters>
Rotation parameter matrix R:
R = [R11 R12 R13;
R21 R22 R23;
R31 R32 R33]
(Inverse matrix of 3 × 3 matrix for performing coordinate conversion from voxel coordinate system to viewpoint coordinate system, GUI side instructs coordinate conversion from voxel coordinate system to viewpoint coordinate system, but rendering side is viewpoint coordinate system Performs coordinate conversion from to the voxel coordinate system.)
XYZ axis offset:
Xoff, Yoff, Zoff (normally Zoff = 0 because it is specified on the 2D screen)
Clipping area:
X-axis direction ROI: Xs-Xe
Y-axis direction ROI: Ys-Ye
Z-axis direction ROI: Zs-Ze
Effective voxel area (circumscribed rectangular parallelepiped area):
X-axis direction ROI: Xis-Xie
Y-axis direction ROI: Yis-Yie
Z-axis direction ROI: Zis-Zie
Scaling Scale (same in XYZ axis direction)
X-axis direction variable magnification: Scx, Y-axis direction variable magnification: Scy, Z-axis direction variable magnification Scz
Coordinate conversion subsample offset: X-axis direction dx, Y-axis direction dy, Z-axis direction dz
Size of generated volume rendered image: Size (same in XY axis direction)
Virtual ray subsampling magnification: Slay (usually 1, if the value is greater than 1, it becomes coarser, and if it is less than 1, it becomes high definition. Is set to Sray = 2, etc.)

The control unit 11 sets a unit matrix as an initial value for the rotation parameter matrix R described above. That is, R11 = 1, R12 = 0, R13 = 0, R21 = 0, R22 = 1, R23 = 0, R31 = 0, R32 = 0, and R33 = 1.
Then, according to the instructions of the GUI, one of X-axis center rotation Rx, Y-axis center rotation Ry, and Z-axis center rotation Rz (angle unit: radian) is sequentially specified, and a rotation matrix A is generated for each as follows. The rotation parameter matrix R is multiplied from the right to update the rotation parameter matrix R. As a result, the inverse matrix of the rotation matrix from the voxel coordinate system to the viewpoint coordinate system generated by the instruction of the GUI is calculated.

Rotation matrix A
A = [A11 A12 A13;
A21 A22 A23;
A31 A32 A33]
Then
Each element of the rotation matrix A in the case of the X-axis center rotation Rx is
A11 = 1, A12 = 0, A13 = 0
A21 = 0, A22 = cosRx, A23 = sinRx
A31 = 0, A32 = sinRx, A33 = cosRx
Each element of the rotation matrix A in the case of Y-axis center rotation Ry is
A11 = cosRy, A12 = 0, A13 = sinRy
A21 = 0, A22 = 1, A23 = 0
A31 = -sinRy, A32 = 0, A33 = cosRy
Each element of the rotation matrix A in the case of Z-axis center rotation Rz is
A11 = cosRz, A12 = sinRz, A13 = 0
A21 = -sinRz, A22 = cosRz, A23 = 0
A31 = 0, A32 = 0, A33 = 1
Will be.
The rotation parameter matrix R is updated as R ← R × A.

The coordinate conversion process based on the coordinate conversion parameters defined above is sequentially executed in each process of the ray casting process (step S93, FIG. 28).

Subsequently, the control unit 11 initializes the number of processes l = 0 and the subsample offset value to dx = dy = dz = 0 (step S92).
Then, the control unit 11 executes a ray casting process (a process of calculating a color value for each coordinate (x, y)) (step S93).

<Ray casting process>
The ray casting process will be described with reference to the flowchart of FIG. 28.
First, the control unit 11 sets all the initial values of the 24-bit (RGB) volume-rendered image ImageVR (x, y, n) to be generated to 0 (ImageVR (x, y, n) = 0, n = 0). (R), 1 (G), 2 (B)). Then, as the number of subsamples L, the following processing is executed for each of the two-dimensional coordinates (x, y) (0 ≦ x ≦ Size-1, 0 ≦ y ≦ Size-1).

The control unit 11 sets the background color pg (n) (n = 0 (R), 1 (G), 2 (B), 0 ≦ pg (n) ≦ 255), and sets x = 0 and y = 0. (Step S101).

Subsequently, the control unit 11 initializes the virtual light intensity Trans = 1.0 and the cumulative luminance value Energy (n) = 0.0 (0 ≦ n ≦ 2) (step S102).

Subsequently, the control unit 11 calculates the Z coordinate (Zc) of the intersection of the effective voxel region Vr and the virtual ray (step S103). The process of calculating the Z coordinate (Zc) of the intersection of the effective voxel region Vr and the virtual ray will be described later (see "11-1. Intersection calculation process").

Subsequently, the control unit 11 sets Zst to Zst = Zc (step S104), and searches for the Z coordinate z (starting point coordinate at which the irradiation of the virtual ray is started) of the first effective voxel (opaque voxel) from Zst (step S104). Step S105). That is, the starting point coordinates for starting the ray casting calculation from the intersection calculated in step S103 toward the line of sight are searched for. The process of searching for the starting point coordinate z will be described later (see “11-2. Starting point coordinate search process”).

If z <0, the process proceeds to step S111.

On the other hand, when z ≧ 0, the control unit 11 transforms the three-dimensional coordinates (x, y, z) to obtain the voxel α value (Vα), and obtains Integer = Vα / 255 (step S106). The process of acquiring the voxel α value (Vα) will be described later (see “11-4. Acquisition of voxel α value and RGB value (with interpolation)”).

On the other hand, when the Alpha calculated in step S106 is Alpha <0 (passing through the effective voxel region Vr), or Alpha = 0 and z = 0, the process proceeds to step S111 and the RGB pixel value in the pixel (x, y) is determined. To do. As a result, when the effective voxel region Vr is deviated (passed) (Alpha <0), the voxel to be drawn does not exist ahead of the effective voxel region Vr, so that the raycasting process is terminated early and redundant processing is omitted.

On the other hand, when the Integer calculated in step S106 is Integer> 0, the control unit 11 transforms the three-dimensional coordinates (x, y, z) to acquire the voxel RGB value Vc (n) (step S108). The process of acquiring the voxel RGB value (Vc (n)) will be described later (see “11-4. Acquisition of voxel α value and RGB value (with interpolation)”).

Subsequently, the control unit 11 sets the α value based on the virtual ray subsample.
Alpha'= 1- (1-Alpha) ^{1 / Sray}
And correct
Cumulative brightness,
Energy (n) = Energy (n) + Trans / Alpha · Vc (n) / 255
Transmitted light intensity,
Trans = Trans · (1.0-Alpha)
Is updated (step S109).

When Trans <0.001, the process proceeds to step S111, and the RGB pixel value in the pixel (x, y) is determined. When Trans ≧ 0.001, the control unit 11 updates to z = z-1 (step S110), returns to step S106 when z ≧ 0, shifts to step S111 when z <0, and pixels. The RGB pixel value at (x, y) is determined.

In step S111, the control unit 11 determines the RGB pixel values in the pixels (x, y) as follows.

(Equation 29)
ImageVR (x, y, n) = ImageVR (x, y, n) + k · Energy (n) · Light (n) / L
Here, k is the intensity magnification, and the initial value is set to k = 1.0.

Subsequently, the control unit 11 updates to x ← x + 1 (step S112), and if x <Size, returns to step S102 and executes a process of determining the RGB value of the next pixel x (steps S102 to S111). .. If x ≧ Size, update to y ← y + 1 (step S113), if y ≦ Size, set x = 0, return to step S102, and change the RGB values of the pixels (x, y) on the y-th row. The process of determining (steps S102 to S111) is repeated. When y> Size, the control unit 11 ends the process.
That is, the processes of steps S102 to S111 are repeated until the RGB values of all the pixels are calculated (step S112; x ≧ Size and step S113; y> Size).

Return to the flowchart of FIG. 27.
The control unit 11 updates the number of processes to l ← l + 1, and updates the subsample offset values to dx ← dx + 1 / L, dy ← dy + 1 / L, and dz ← dz + 1 / L (step S94). The control unit 11 repeats the ray casting process of step S93 until the number of processes l satisfies l> L-1 (step S94; l> L-1). When the ray casting process is completed, the control unit 11 outputs the generated volume rendered image ImageVR (step S95).

(11-1. Intersection calculation process)
The process of calculating the Z coordinate (Zc) of the intersection of the effective voxel region Vr and the virtual ray, which is executed in step S103 of FIG. 28, will be described.

The control unit 11 performs the following processing procedure with respect to the calculation target pixel (x, y) in the ray casting process in the range of z = Zoo (= Size / Sray-1) to z = 0 in the viewpoint coordinate system. , The intersection Zc with the effective voxel region Vr ([Xis, Xie], [Yis, Yie], [Zis, Zie]) closest to the viewpoint z = Zo is calculated.

(1) Calculation of line-of-sight vector (virtual ray vector) As shown in FIG. 29, the control unit 11 refers to the viewpoint coordinates (x, y, Zo) and the lower limit coordinates (x, y, 0) in the viewpoint coordinate system. Perform coordinate transformation for each, obtain the viewpoint coordinates (x1, y1, z1) and lower limit coordinates (x2, y2, z2) in the boxel coordinate system, and calculate each element of the line-of-sight vector (vx, vy, vz) as follows. To do.

vx = x2-x1
by = y2-y1
vz = z2-z1

(2) Calculation of the intersection of the line-of-sight vector and the effective voxel region Vr The control unit 11 initializes tx = ty = tz = 10.
1) When | vx | ≧ 1, tx1 = (Xis-x1) / vx and tx2 = (Xie-x1) / vx are calculated, and the smaller one is set to tx.
2) When | by | ≧ 1, ty1 = (Yis-y1) / by and ty2 = (Yie-y1) / vy are calculated, and the smaller one is set to ty.
3) When | vz | ≧ 1, tz1 = (Zis-z1) / vz and tz2 = (Zie-z1) / vz are calculated, and the smaller one is set to tz.

As a result, as shown in FIG. 29 (right figure), the coordinates of the intersection in the voxel coordinate system can be obtained. Specifically, the intersection coordinates are given by (vx · t + x1, vy · t + y1, vz · t + z1) (t is any value of tx, ty, tz calculated in (2)), and the example of FIG. 29. Then, there are intersections on the two surfaces of z = Zis and z = Zie, and the coordinates of the intersections at z = Zie can be calculated as (vx · tz + x1, vy · tz + y1, vz · tz + z1). The intersection coordinates are the intersections of the line-of-sight vector and the six planes constituting the effective voxel region Vr (rectangular parallelepiped) (usually, there are two or more intersections), which is the closest to the viewpoint.

(3) Determination of intersection Zc in the viewpoint coordinate system The control unit 11 initializes the intersection coordinate Zc in the viewpoint coordinate system to Zc = -1.
1) When tx ≦ 1 Y = vy · tx + y1 and Z = vz · tx + z1 are calculated.
In the case of Yis ≤ Y ≤ Yie and Zis ≤ Z ≤ Zie (when the intersection in the voxel coordinate system exists in the effective voxel region Vr), the intersection Z2 = Zo · (1-tx) in the corresponding viewpoint coordinate system is calculated. , If Z2> Zc, replace with Zc = Z2.
2) When ty ≦ 1 X = vx · ty + x1 and Z = vz · ty + z1 are calculated.
In the case of Xis ≤ X ≤ Xie and Zis ≤ Z ≤ Zie (when the intersection in the voxel coordinate system exists in the effective voxel region Vr), the intersection Z2 = Zo · (1-ty) in the corresponding viewpoint coordinate system is calculated. , If Z2> Zc, replace with Zc = Z2.
3) When tz ≦ 1 X = vx · tz + x1 and Y = vy · tz + y1 are calculated.
In the case of Xis ≤ X ≤ Xie and Yis ≤ Y ≤ Yie (when the intersection in the voxel coordinate system exists in the effective voxel region Vr), the intersection Z2 = Zo · (1-tz) in the corresponding viewpoint coordinate system is calculated. , If Z2> Zc, replace with Zc = Z2.

As a result, as shown in FIG. 29 (left figure), the intersection coordinates (x, y, Zc) closest to the viewpoint in the viewpoint coordinate system can be obtained.

(4) Correction of intersections calculated before the viewpoint tx, ty, and tz may take negative values, in which case the intersections are before the viewpoint. When Zc> Zo, the control unit 11 corrects Zc = Zo. That is, as shown in FIG. 30, when the intersection coordinates in the viewpoint coordinate system are calculated in front of the viewpoint (invisible region), the intersection coordinates are corrected to the viewpoint position.

(11-2. Origin coordinate search process)
The process of searching for the starting point coordinate z, which is executed in step S105 of FIG. 28, will be described with reference to the flowchart of FIG. 31.

First, the control unit 11 initializes z to z = Zst (= intersection coordinates Zc) and inputs the search target pixels (x, y) (step S121). Subsequently, the control unit 11 performs coordinate conversion on the three-dimensional coordinates (x, y, z) and acquires the α value of the voxel (step S122). The process of acquiring the voxel α value will be described later (see “11-3. Obtaining the voxel α value (without interpolation)”).

When the α value of the voxel acquired in step S122 is α = 0 (when the voxel is transparent), the control unit 11 updates to z ← z−m (for example, m = 4) (step S123). m is the skip width in the Z-axis direction at the time of search, and the larger m is, the faster the search can be performed. However, if the skip width is too large, the search accuracy of the starting coordinate may deteriorate.

When the updated z in step S123 is z <0 (outside the voxel space), the control unit 11 returns z = -1 (a value indicating that the starting point coordinates do not exist) (step S124), and ends the process. To do. When z ≧ 0, the control unit 11 returns to step S122.

When the α value of the voxel acquired in step S122 is α <0 (when passing through the effective voxel region Vr), the control unit 11 returns z = -1 (a value indicating that the starting coordinate does not exist) (a value indicating that the starting coordinate does not exist). Step S124), the process is terminated.

When the α value of the voxel acquired in step S122 is α> 0 (when the voxel is opaque), the process returns to the + direction of the Z axis (opposite direction of the line of sight) and searches for accurate starting point coordinates as follows.

First, the control unit 11 initializes i to i = 0 (step S125) and updates i ← i + 1 (step S126).

When i ≧ m or z + i ≧ Zst, the control unit 11 returns z = z + i-1 (step S128) and ends the process.

In the case of i <m and z + i <Zst, coordinate conversion is performed on the three-dimensional coordinates (x, y, z + i) to acquire the α value of the voxel (step S127). The process of acquiring the voxel α value will be described later (see “11-3. Obtaining the voxel α value (without interpolation)”).

When the α value of the voxel acquired in step S127 is α = 0 (when the voxel is transparent), the control unit 11 returns z = z + i-1 (step S128), and ends the process.

When the α value of the voxel acquired in step S127 is α> 0 (when the voxel is opaque), the process returns to step S126 and the processes of steps S126 to S127 are repeated.

As described above, the starting point coordinates can be specified at high speed by searching the Z coordinate z (starting point coordinate for starting irradiation of virtual light rays) of the leading opaque voxel with a predetermined step width (steps S122 and S123). In addition, when an opaque voxel is searched, the opaque voxel that may exist between the position of the voxel and the position returned by the step width in the opposite direction of the line of sight (+ direction of the Z axis) can be searched again. (Steps S126 and S127), the starting point coordinates can be accurately specified.

(11-3. Obtain voxel α value (without interpolation))
With reference to the flowchart of FIG. 32, the process of acquiring the α value of the voxel, which is executed in step S122 and step S127 of the origin coordinate search process (FIG. 31), will be described.

First, the control unit 11 performs coordinate conversion on the three-dimensional coordinate values (x, y, z) (integer values) of the given viewpoint coordinate system, and performs coordinate values (xr, yr,) of the real numbers of the corresponding voxels. zr) is calculated (step S131). Here, the coordinate conversion is a process of converting the viewpoint coordinate system into a voxel coordinate system, which is the opposite of the conversion process on the GUI side. On the GUI side, it is assumed that clipping by the region of interest ROI, scaling, Z-axis direction scaling processing, offset (simultaneous in the XYZ-axis direction), rotation, and fluoroscopic conversion are performed in this order, and the control unit 11 controls the given viewpoint coordinate system. The real coordinate values (xr, yr, zr) of the voxels corresponding to the three-dimensional coordinate values (x, y, z) (integer values) are calculated as follows.

First, the control unit 11 converts the coordinate values (x, y, z) (integer values) of the viewpoint coordinate system into real values (xx, yy, zz) as follows. The viewpoint coordinate system is a rectangular parallelepiped of Size (size in the X-axis direction) x Size (size in the Y-axis direction) x Size / Sray (size in the Z-axis direction), and is independent of the tomographic image group Do of Sx x Sy x Sz. Is defined. Normally, if Sx = Sy = 512, Size = 512. Considering that the Z-axis direction of the viewpoint coordinate system is expanded and contracted by the virtual ray subsample 1 / Sray times, first, the offset processing (simultaneous in the XYZ axis directions) is performed as follows.

(Equation 30)
xx = x-Size / 2 + dx + Xoff
yy = y-Size / 2 + dy + Yoff
zz = z · Sray-Size / 2 + dz + Zoff

Subsequently, the control unit 11 performs rotation processing (coordinate conversion) using the rotation parameter matrix R as follows.

(Equation 31)
xx'= R11, xx + R12, yy + R13, zz
yy'= R21 ・ xx + R22 ・ yy + R23 ・ zz
zz'= R31, xx + R32, yy + R33, zz
Let (xx', yy', zz') after the rotation process be (xx, yy, zz).

Then, the control unit 11 performs scaling and scaling in the XYZ axis directions at the same time, and calculates the coordinate values (xr, yr, zr) (real values) of the voxels as follows.

(Equation 32)
xx = xx / Scale / Scx + Sx / 2
yr = yy / Scale / Scy + Sy / 2
zr = zz / Scale / Scz + Sz / 2

Subsequently, the control unit 11 adds 0.5 to each value with respect to the calculated coordinate values (xr, yr, zr) (real values) of the voxel, and rounds down the decimal point as shown below. It is converted into the converted coordinate values (xi, yi, zi) (rounded integer values) (step S132).

(Equation 33)
xi = INT [xr + 0.5]
yi = INT [yr + 0.5]
zi = INT [zr + 0.5]

Then, the control unit 11 acquires the voxel α value as follows in consideration of the effective voxel region Vr and the mask data Mask (step S133).

(Equation 34)
1) When any of xi <Xis, xi> Xie, y <Yis, y> Yie, zi <Zis, or zi> Zie is satisfied (when outside the effective voxel region Vr)
α = -1
2) When the above 1) is not satisfied α = Cmap (Do (xi, y, zi), 3) · Mask (xi, y, zi)

(11-4. Acquisition of voxel α value and RGB value (with interpolation))
The process of acquiring the voxel α value (Vα) and the voxel RGB value, which is executed in steps S106 and S108 of the ray casting process (FIG. 28), will be described.

(11-4-1. Acquisition of voxel α value)
FIG. 33 is a flowchart showing a process of acquiring a voxel α value.

First, the control unit 11 performs coordinate conversion on the three-dimensional coordinate values (x, y, z) (integer values) of the viewpoint coordinate system in the same manner as in step 131 of FIG. 32, and coordinates the real numbers of the voxels. The values (xr, yr, zr) are calculated (step S141).

Subsequently, the control unit 11 converts the calculated coordinate values (xr, yr, zr) (real values) of the voxels into coordinate values (xi, y, zi) (integer values) obtained by rounding down the decimal point and converting them into integers. (Step S142). Here, the truncated fraction after the decimal point is (wx, wy, wz) (0 ≦ wx, wy, wz <1) (that is, xr = xi + wx, yr = yi + wy, zr = zi + wz).

Then, the control unit 11 considers the effective voxel region Vr, and sets the signal value of the voxel to be interpolated corresponding to the three-dimensional coordinate values (x, y, z) (integer value) of the viewpoint coordinate system to the tomographic image group Do. Is extracted as follows (step S143).

(Equation 35)
1) When any of xi <Xis, xi> Xie, y <Yis, y> Yie, zi <Zis, or zi> Zie is satisfied (outside the effective voxel region Vr)
D000 = -32769 (invalid value)
2) When the condition of 1) above is not satisfied and either xi + 1> Xie, yi + 1> Yie, or zi + 1> Zie is satisfied (central voxel is extracted).
D000 = Do (xi, yi, zi)
3) When the conditions of 1) and 2) above are not satisfied (signal values of voxels in the vicinity of 8 are extracted).
D000 = Do (xi, yi, zi)
D100 = Do (xi + 1, yi, zi)
D010 = Do (xi, yi + 1, zi)
D110 = Do (xi + 1, yi + 1, zi)
D001 = Do (xi, yi, zi + 1)
D101 = Do (xi + 1, yi, zi + 1)
D011 = Do (xi, yi + 1, zi + 1)
D111 = Do (xi + 1, yi + 1, zi + 1)

Subsequently, the control unit 11 sets the mask values S000 to S111 (0 or 1) for the α value (opacity) of each voxel to be extracted based on the mask data Mask as follows (step S144). ). That is, it is set whether or not to use the α value (opacity) of each voxel in the calculation.

(Equation 36)
1) When any of xi <Xis, xi> Xie, y <Yis, y> Yie, zi <Zis, or zi> Zie is satisfied (outside the effective voxel region Vr)
S000 = 0
2) When the condition of 1) above is not satisfied and either xi + 1> Xie, yi + 1> Yie, or zi + 1> Zie is satisfied (the mask value of the central voxel is set).
S000 = Mask (xi, yi, zi)
3) When the conditions of 1) and 2) above are not satisfied (set the mask value of 8 voxels in the vicinity).
S000 = Mask (xi, yi, zi)
S100 = Mask (xi + 1, yi, zi)
S010 = Mask (xi, yi + 1, zi)
S110 = Mask (xi + 1, y + 1, zi)
S001 = Mask (xi, yi, zi + 1)
S101 = Mask (xi + 1, yi, zi + 1)
S011 = Mask (xi, yi + 1, zi + 1)
S111 = Mask (xi + 1, yi + 1, zi + 1)

Then, the control unit 11 sets the α value Vα of the voxel corresponding to the three-dimensional coordinate value (x, y, z) (integer value) of the viewpoint coordinate system as follows, based on the signal value of the extracted voxel to be interpolated. (Step S145).

(Equation 37)
1) When any of xi <Xis, xi> Xie, y <Yis, y> Yie, zi <Zis, or zi> Zie is satisfied (outside the effective voxel region Vr)
Vα = -1
2) When the condition of 1) above is not satisfied and any one of xi = Xis, xi = Xie, y = Yis, y = Yie, zi = Zis, or zi = Zie is satisfied (boundary surface of effective voxel region Vr).
Vα = 0
3) When the above conditions 1) and 2) are not satisfied and either xi + 1> Xie, yi + 1> Yie, or zi + 1> Zie is satisfied (no interpolation).
Vα = Cmap (D000, 3) · S000
4) When the above conditions 1) 2) 3) are not satisfied (interpolate)
Vα = (1-wz) (1-wy) (1-wx) ・ Cmap (D000, 3) ・ S000 + (1-wz) (1-wy) ・ wx ・ Cmap (D100, 3) ・ S100 + (1- wz), wy, (1-wx), Cmap (D010, 3), S010 + (1-wz), wy, wx, Cmap (D110, 3), S110 + wz, (1-wy) (1-wx), Cmap (D001, 3), S001 + wz, (1-wy), wx, Cmap (D101, 3), S101 + wz, wy, (1-wx), Cmap (D011, 3), S011 + wz, wy, wx, Cmap (D111, 3) ・ S111

(11-4-2. Acquisition of voxel RGB values)
The control unit 11 has a voxel RGB value (Vc (n)) corresponding to a three-dimensional coordinate value (x, y, z) (integer value) of the viewpoint coordinate system based on the signal value of the voxel to be interpolated extracted in step S143. ) (0 ≦ n ≦ 2) is acquired as follows. When acquiring the RGB value, the mask process for multiplying S000, S100, ..., S111 is not performed. As a result, the occurrence of moire on the mask boundary surface is suppressed.

(Equation 38)
1) When any of xi <Xis, xi> Xie, y <Yis, y> Yie, zi <Zis, or zi> Zie is satisfied (outside the effective voxel region Vr)
Vc (n) = 0 (0 ≦ n ≦ 2)
2) When the condition of 1) above is not satisfied and either xi + 1> Xie, yi + 1> Yie, or zi + 1> Zie is satisfied (no interpolation).
Vc (n) = Cmap (D000, n) (0 ≦ n ≦ 2)
3) When the above conditions 1) and 2) are not satisfied (interpolate)
Vc (n) = (1-wz) (1-wy) (1-wx) · Cmap (D000, n) + (1-wz) (1-wy) · wx · Cmap (D100, n) + (1 -Wz), wy, (1-wx), Cmap (D010, n) + (1-wz), wy, wx, Cmap (D110, n) + wz, (1-wy) (1-wx), Cmap ( D001, n) + wz · (1-wy) · wx · Cmap (D101, n) + wz · wy · (1-wx) · Cmap (D011, n) + wz · wy · wx · Cmap (D111, n) (0 ≦ n ≦ 2)

(11-5. Shading calculation)
The control unit 11 may perform the shadow calculation as follows, if necessary.
First, the light source vector (Lx, Ly, Lz) (unit vector) is set. For example, (Lx, Ly, Lz) = (0.57735, 0.57735, 0.57735) is set. Further, the ambient light component Ab (0 ≦ Ab ≦ 1, for example, Ab = 0.2) is set.

Subsequently, the control unit 11 acquires the opacity V100, V200, V010, V020, V001, and V002 of the 6-neighboring voxels at the three-dimensional coordinates (xi, yi, zi) calculated by the coordinate transformation as follows. .. However, when xi + 1> Xie, V100 = 0, when xi-1 <Xis, V200 = 0, when yi + 1> Yie, V010 = 0, when y-1 <Yis, V020 = 0, and when zi + 1 <Zis, V001 = 0. , Zi-1 <Zis, V002 = 0. That is, when the adjacent voxels are in the clipping range, the opacity is treated as 0.

(Equation 39)
V100 = Cmap (Do (xi + 1, y, zi), 3) · Mask (xi + 1, y, zi)
V200 = Cmap (Do (xi-1, yi, zi), 3), Mask (xi-1, yi, zi)
V010 = Cmap (Do (xi, y + 1, zi), 3) · Mask (xi, y + 1, zi)
V020 = Cmap (Do (xi, y-1, zi), 3) · Mask (xi, y-1, zi)
V001 = Cmap (Do (xi, y, zi + 1), 3) · Mask (xi, y, zi + 1)
V002 = Cmap (Do (xi, y, zi-1), 3) · Mask (xi, y, zi-1)

Subsequently, the control unit 11 calculates the gradient vector (Gx, Gy, Gz) in the three-dimensional coordinates (xi, yi, zi) calculated by the coordinate transformation by the following formula.

(Equation 40)
Gx = (V100-V200) · Scx
Gy = (V010-V020) · Scy
Gz = (V001-V002) · Scz
G = {Gx ² + Gy ² + Gz ² } ^1/2

Subsequently, when G ≧ 1, the control unit 11 calculates the diffuse reflection component as the luminance value (shadow value) S (0 ≦ S ≦ 1).

(Equation 41)
S = (1-Ab) | Gx · Lx + Gy · Ly + Gz · Lz | / G + Ab

give. When G <1 (when α does not change), the control unit 11 gives S = 0 as the luminance value (shadow value) S.

Then, the control unit 11

(Equation 42)
Vc'(n) = S · Vc (n)

Then, the component of the calculated RGB value Vc (n) (0 ≦ n ≦ 2) is modified.

The process of generating the volume rendered image has been described above with reference to FIGS. 25 to 33. After setting the mask value for the opacity of each voxel of the tomographic image group Do based on the mask data Mask, apply the color map Cmap (opacity curve) to each voxel of the tomographic image group Do to obtain the voxel α value. By (interpolation calculation), it is possible to generate a volume rendering image in which the drawing / non-drawing information of the mask data Mask is reflected.

(12. Generation of MIP image)
Next, the process of generating the MIP image will be described with reference to FIGS. 34 to 40.

FIG. 34 is a flowchart showing the entire flow of the process of generating the MIP image.

First, the control unit 11 receives a setting of the clipping region ROI (see (Equation 22)) from the user as needed (step S161).

Subsequently, the control unit 11 calculates the effective voxel region Vr (step S162). Specifically, as shown in FIG. 35, the control unit 11 sets voxels whose signal values are equal to or higher than a predetermined threshold value (for example, equal to or higher than non-air signal values) as effective voxels, includes all of these effective voxels, and The rectangular parallelepiped circumscribing the effective voxel is calculated as the effective voxel region Vr. When the clipping region ROI and the mask data Mask are set, the control unit 11 sets the signal value to be equal to or higher than a predetermined threshold value (for example, equal to or higher than the non-air signal value) within the drawing range defined by the clipping ROI or the mask data Mask. ) Is defined as an effective voxel, and a rectangular body including all of these effective voxels and circumscribing the effective voxel is calculated as an effective voxel region Vr.

Then, the control unit 11 is set in the tomographic image group Do (three-dimensional voxel) acquired in step S71 (FIG. 24), the mask data Mask generated in step S73 (FIG. 24), and step S161 (FIG. 34). The MIP image ImageMIP is generated with reference to the clipping region ROI and the effective voxel region Vr calculated in step S162 (FIG. 34) (step S163).

In particular, in the present disclosure, in the ray casting process, the control unit 11 sets the mask value for the signal value of each voxel of the tomographic image group Do based on the mask data Mask (step S224 in FIG. 40 described later), and then the fault. The signal value of each voxel of the image group Do is acquired (interpolation calculation) (step S225 in FIG. 40 described later), and a MIP image ImageMIP is generated. As a result, a MIP image that reflects the drawing / non-drawing information of the mask data Mask generated based on the subject shape (geometric shape) is generated.

FIG. 36 shows a specific flow of processing for generating a MIP image ImageMIP, which is executed in step S163.
First, the control unit 11 sets the coordinate conversion parameters (step S171). The setting of the coordinate conversion parameter is the same as in step S91 of FIG. 27 described above.

Subsequently, the control unit 11 executes a ray casting process (a process of calculating a representative signal value for each coordinate (x, y)) (step S172).

<Ray casting process>
The ray casting process will be described with reference to the flowchart of FIG. 37. The generated MIP image ImageMIP (x, y) has a monochrome 8-bit or 16-bit value based on the gradation of the tomographic image group Do, and the tomographic image group Do is subjected to gradation compression processing in step S71 of FIG. In many cases, it is monochrome 16-bit. Hereinafter, a case where a monochrome (16-bit) MIP image ImageMIP (x, y) is generated will be described.

First, the control unit 11 sets all the initial values of the monochrome (16-bit) MIP image ImageMIP (x, y) to be generated to 0 (ImageMIP (x, y) = 0). Further, the representative signal value calculation mode mode (0: MIP, 1: MinIP, 2: RaySum) is set, and the two-dimensional coordinates (x, y) (0 ≦ x ≦ Size-1, 0 ≦ y ≦ Size-1) are set. ), The following processing is executed.

The control unit 11 sets x = 0 and y = 0 (step S181), and calculates the Z coordinate (Zc) of the intersection of the effective voxel region Vr and the virtual ray (step S182). The process of calculating the Z coordinate of the intersection is the same as the process of step S103 of FIG. 28 described above (see “11-1. Intersection calculation process”).

Subsequently, the control unit 11 initializes the representative signal value Vm and the counter ct as follows according to the representative signal value calculation mode mode, and sets Zst = Zc (step S183).

(Equation 43)
mode = 0 (MIP): Vm = -32768 (set the minimum value)
mode = 1 (MinIP): Vm = 32767 (set the maximum value)
mode = 2 (RaySum): Vm = 0, ct = 0

Subsequently, the control unit 11 searches for the Z coordinate z (starting point coordinate) of the leading effective voxel (the voxel whose signal value is equal to or higher than a predetermined threshold value) following Zst (step S184). That is, the starting point coordinates for starting the ray casting calculation from the intersection calculated in step S182 toward the line of sight are searched for. The process of searching for the starting point coordinates will be described later (see “12-1. Starting point coordinate search processing”, FIG. 38).

If z <0, the process proceeds to step S189.

On the other hand, when z ≧ 0, the control unit 11 transforms the three-dimensional coordinates (x, y, z) to acquire the voxel signal value Vs (step S185). The process of acquiring the voxel signal value Vs will be described later (see "12-3. Acquiring the voxel signal value Vs (with interpolation)").

When the voxel signal value Vs acquired in step S185 is Vs = -327769 (when passing through the effective voxel region Vr), the process proceeds to step S189. As a result, when the voxel region Vr has passed (Vs = -327769), the voxel to be drawn does not exist after that, so the ray casting process is terminated early and redundant processing is omitted.

When the voxel signal value Vs acquired in step S185 is Vs <-32769 (in the case of an invalid value), Zst = z-1 is set (step S186), the process returns to step S184, and the starting point coordinates are searched again.

When the voxel signal value Vs acquired in step S185 is Vs ≧ -32768 (when it is an effective value), the control unit 11 calculates the representative signal value Vm as follows according to the representative signal value calculation mode mode (when it is a valid value). Step S187).

(Equation 44)
mode = 0 (MIP mode):
If Vs> Vm, then Vm = Vs
mode = 1 (MinIP mode):
If Vs <Vm, then Vm = Vs
mode = 2 (RaySum mode):
Vm ← Vm + Vs, ct ← cnt + 1

It is updated to z = z-1 (step S188), and if z ≧ 0, the process returns to step S185, and if z <0, the process proceeds to step S189.

In step 189, when mode = 0 and 1 (step S189; mode <2), the control unit 11 records the representative signal value Vm obtained in step S187 as a pixel value (ImageMIP (x, y). = Vm, step S191). On the other hand, in the case of mode = 2 (RaySum mode) (step S189; mode = 2), the control unit 11 calculates the average of the signal values as Vm ← Vm / ct (step S190), and then represents the representative. The signal value Vm is recorded as a pixel value (ImageMIP (x, y) = Vm, step S191).

Subsequently, the control unit 11 updates to x ← x + 1 (step S192), and if x <Size, returns to step S182 and executes a process (steps S182 to S191) for obtaining the representative signal value Vm of the next pixel x. To do. If x ≧ Size, update to y ← y + 1 (step S193), if y ≦ Size, set x = 0, return to step S182, and return to the representative signal value of the pixel (x, y) on the y-th row. The process of obtaining Vm (steps S182 to S191) is repeated. When y> Size, the control unit 11 ends the process.
That is, the processing of steps S182 to S191 is repeated until the representative signal values Vm of all the pixels in the pixel region are obtained (step S192; x ≧ Size and step S193; y> Size).

Return to the flowchart of FIG.
The control unit 11 outputs the MIP image ImageMIP generated in step S172 (step S173).

(12-1. Origin coordinate search process)
The process of searching for the starting point coordinates, which is executed in step S184 of FIG. 37, will be described with reference to the flowchart of FIG. 38.

First, the control unit 11 initializes z to z = Zst and inputs the search target pixels (x, y) (step S201). Subsequently, the control unit 11 performs coordinate conversion on the three-dimensional coordinates (x, y, z) and acquires the voxel signal value Vs (step S202). The process of acquiring the voxel signal value Vs will be described later (see "12-2. Acquiring the voxel signal value Vs (without interpolation)").

When the voxel signal value Vs acquired in step S202 is Vs = -32769 (when passing through the effective voxel region Vr), the control unit 11 returns z = -1 (a value indicating that the starting coordinate does not exist) (a value indicating that the starting coordinate does not exist). Step S204), the process is terminated.

On the other hand, when the voxel signal value Vs <-327769 acquired in step S202 (in the case of an invalid value), the control unit 11 updates to z ← z-1 (step S203).

When z <0 is updated in step S203, the control unit 11 returns z = -1 (a value indicating that the starting point coordinates do not exist) (step S204), and ends the process. When z ≧ 0, the control unit 11 returns to step S202.

On the other hand, when the voxel signal value Vs acquired in step S202 is Vs ≧ -32768 (in the case of an effective voxel), the control unit 11 outputs z as the starting coordinate.

(12-2. Acquire voxel signal value Vs (without interpolation))
The process of acquiring the voxel signal value Vs, which is executed in step S202 of the origin coordinate search process (FIG. 38), will be described with reference to the flowchart of FIG. 39.

First, the control unit 11 performs coordinate conversion on the three-dimensional coordinate values (x, y, z) (integer values) of the viewpoint coordinate system in the same manner as in step S131 of FIG. 32, and coordinates the real numbers of the boxel. The values (xr, yr, zr) are calculated (step S211), and the calculated real number coordinate values (xr, yr, zr) (real values) of the boxel are converted into integer coordinates in the same manner as in step S132 of FIG. Convert to a value (xi, yi, zi) (step S212).

Then, the control unit 11 acquires the voxel signal value Vs as follows in consideration of the effective voxel region Vr and the mask data Mask (step S213).

(Equation 45)
1) When any of xi <Xis, xi> Xie, y <Yis, y> Yie, zi <Zis, or zi> Zie is satisfied (when outside the effective voxel region Vr)
Vs = -32769
2) When Mask (x, y, z) = 0 without satisfying 1) above (invalid value)
Vs = -99999
3) When the above 1) and 2) are not satisfied (valid value)
Vs = Do (xi, yi, zi)

(12-3. Acquisition of voxel signal value Vs (with interpolation))
The process of acquiring voxel signal values Vs, which is executed in step S185 of the ray casting process (FIG. 37), will be described with reference to the flowchart of FIG. 40.

First, the control unit 11 performs coordinate conversion on the three-dimensional coordinate values (x, y, z) (integer values) of the viewpoint coordinate system in the same manner as in step S141 of FIG. 33, and coordinates the real numbers of the boxel. The values (xr, yr, zr) are calculated (step S221), and the calculated coordinate values (xr, yr, zr) (real values) of the real numbers of the boxel are calculated in the same manner as in step S142 of FIG. Is converted into a coordinate value (xi, yi, zi) converted into an integer by truncating (step S222). The truncated fraction after the decimal point is (wx, wy, wz) (0 ≦ wx, wy, wz <1) (that is, xr = xi + wx, yr = yi + wy, zr = zi + wz).

Subsequently, the control unit 11 sets the signal values D000, D100, D010, D001, of the voxel to be interpolated corresponding to the three-dimensional coordinate values (x, y, z) (integer values) of the viewpoint coordinate system as follows. D101, D011, and D111 are extracted (step S223).

(Equation 46)
1) When any of xi <Xis, xi> Xie, y <Yis, y> Yie, zi <Zis, or zi> Zie is satisfied (outside the effective voxel region Vr)
D000 = -32769 (invalid value)
2) When the condition of 1) above is not satisfied and either xi + 1> Xie, yi + 1> Yie, or zi + 1> Zie is satisfied (central voxel is extracted).
D000 = Do (xi, yi, zi)
3) When the conditions of 1) and 2) above are not satisfied (signal values of voxels in the vicinity of 8 are extracted).
D000 = Do (xi, yi, zi)
D100 = Do (xi + 1, yi, zi)
D010 = Do (xi, yi + 1, zi)
D110 = Do (xi + 1, yi + 1, zi)
D001 = Do (xi, yi, zi + 1)
D101 = Do (xi + 1, yi, zi + 1)
D011 = Do (xi, yi + 1, zi + 1)
D111 = Do (xi + 1, yi + 1, zi + 1)

Subsequently, the control unit 11 sets the mask values S000 to S111 (0 or 1) for the signal value of each voxel to be extracted based on the mask data Mask as follows (step S224). That is, it is set whether or not the signal value of each voxel is used for the calculation.

(Equation 47)
1) 1) When any of xi <Xis, xi> Xie, y <Yis, y> Yie, zi <Zis, or zi> Zie is satisfied (outside the effective voxel region Vr)
S000 = 0
2) When the condition of 1) above is not satisfied and either xi + 1> Xie, yi + 1> Yie, or zi + 1> Zie is satisfied (the mask value of the central voxel is set).
S000 = Mask (xi, yi, zi)
3) When the conditions of 1) and 2) above are not satisfied (set the mask value of 8 voxels in the vicinity).
S000 = Mask (xi, yi, zi)
S100 = Mask (xi + 1, yi, zi)
S010 = Mask (xi, yi + 1, zi)
S110 = Mask (xi + 1, y + 1, zi)
S001 = Mask (xi, yi, zi + 1)
S101 = Mask (xi + 1, yi, zi + 1)
S011 = Mask (xi, yi + 1, zi + 1)
S111 = Mask (xi + 1, yi + 1, zi + 1)

Then, the control unit 11 sets the voxel signal values Vs corresponding to the three-dimensional coordinate values (x, y, z) (integer values) of the viewpoint coordinate system based on the extracted signal values of the voxels to be interpolated as follows. Acquire (step S225).

(Equation 48)
1) When any of xi <Xis, xi> Xie, y <Yis, y> Yie, zi <Zis, or zi> Zie is satisfied (outside the effective voxel region Vr)
Vs = -32769
2) When the condition of 1) above is not satisfied and S000 = 0 (invalid value)
Vs = -99999
3) When the above 1) and 2) are not satisfied
If either xi + 1> Xie, y + 1> Yie, or zi + 1> Zie is satisfied, or if any of S000, S100, S010, S110, S001, S101, S011, and S111 is 0 (no interpolation).
Vs = D000
4) When the above conditions 1) 2) 3) are not satisfied (interpolate)
Vs = (1-wz) (1-wy) (1-wx) · D000 + (1-wz) (1-wy) · wx · D100 + (1-wz) · wy · (1-wx) · D010 + (1 -Wz), wy, wx, D110 + wz, (1-wy) (1-wx), D001 + wz, (1-wy), wx, D100 + wz, wy, (1-wx), D011 + wz, wy, wx, D111

The case of generating a monochrome (16-bit) MIP image ImageMIP (x, y) has been described above, but the control unit 11 converts the generated MIP image ImageMIP (x, y) into 8 bits and then converts the generated MIP image ImageMIP (x, y) into 8 bits. indicate. Specifically, the DICOM standard Window Width (window width: WW) and Window Level (window level: WL) are set, and the upper limit value Lmax = WL + WW / 2 and the lower limit value Lmin = WL-WW / 2 are converted. The MIP image ImageMIP (x, y) is converted into 8 bits as follows.

ImageMIP'(x, y) = (ImageMIP (x, y) -Lmin) · 255 / (Lmax-Lmin)
If ImageMIP'(x, y) <0, then ImageMIP'(x, y) = 0.
When ImageMIP'(x, y)> 255, ImageMIP'(x, y) = 255.

The process of generating the MIP image has been described above with reference to FIGS. 34 to 40. After setting the mask value for the signal value of each boxel of the tomographic image group Do based on the mask data Mask (step S224 in FIG. 40), the signal value of each boxel of the tomographic image group Do is acquired (interpolation calculation). (Step S225 in FIG. 40), a MIP image reflecting the drawing / non-drawing information of the mask data Mask can be generated.

(13. Generation of MPR image)
Next, the process of generating the MPR image will be described with reference to FIGS. 41 and 42.

FIG. 41 is a flowchart showing a flow of processing for generating an MPR image.

First, when the signal value of the tomographic image group Do is 16 bits (see (Equation 1)), the control unit 11 converts the tomographic image group Do into 8 bits (step S241). Specifically, the control unit 11 controls the tomographic image group Do (-32768 ≦ Do (x, y, z) ≦ 32767, 0 ≦ x ≦ Sx-1, 0 ≦ y ≦ Sy-1, 0 ≦ z ≦ Sz. -1; Set DICOM standard Window Width (window width: WW) and Window Level (window level: WL) for the variable magnifications (Scx, Scy, Scz) in three directions, and set the upper limit value Lmax = WL + WW /. 2. Convert to the lower limit value Lmin = WL-WW / 2, and convert the tomographic image group Do to 8 bits as follows.

(Equation 49)
Do (x, y, z) = (Do (x, y, z) -Lmin) · 255 / (Lmax-Lmin)
When Do (x, y, z) <0, Do (x, y, z) = 0.
When Do (x, y, z)> 255, Do (x, y, z) = 255.

Subsequently, the control unit 11 generates an MPR image ImageMPR of an arbitrary cross section based on the 8-bit tomographic image group Do (step S242). For example, as shown in FIG. 42, in the case of a body axis cross section (Axial cross section) parallel to the XY plane, the (Scx, Scy) scaling is applied to the voxels on the XY plane at the specified slice position z. In addition, each pixel (x, y) is acquired to reconstruct a two-dimensional image. In the case of a coronal section parallel to the XZ plane, each pixel (x, z) is multiplied by (Scx, Scz) scaling for the voxels on the XZ plane at the specified slice position y. Acquire and reconstruct the two-dimensional image. In the case of a sagittal section parallel to the ZY plane (Sagittal section), each pixel (z, y) is multiplied by (Scz, Scy) with respect to the voxels on the ZY plane at the specified slice position x. To reconstruct the two-dimensional image. In this way, MPR images in three directions are generated.

Further, when generating an MPR image of an oblique cross section (Oblique cross section) obtained by diagonally cutting the XYZ space, an instruction to rotate in the three-dimensional space is instructed while rotating one of the MPR images in three directions in two dimensions. Then, based on the generated rotation matrix, the 8-bit tomographic image group Do is rotated in a three-dimensional space, and then each two-dimensional image is reconstructed. Further, when the MPR image of the body axis cross section (axial cross section) is generated by specifying the slab thickness, the slice position is made minute by the specified thickness (number of slices) in the vicinity of the specified slice position z. By moving, each pixel (x, y) is acquired in multiple ways while applying (Scx, Scy) scaling to the voxel on the XY plane, and each pixel is composed of representative signal values of a plurality of signal values. The resulting two-dimensional image is reconstructed. As a method of calculating the representative signal value from a plurality of signal values, any one of the maximum value (MIP), the minimum value (MinIP), and the average value (RaySum) can be set as in the case of generating the MIP image described above. ..

Then, the control unit 11 synthesizes (superimposes) the mask image generated based on the mask data Mask on the MPR image ImageMPR generated in step S242 and outputs the mask image (step S243). For example, the control unit 11 has Mask (x, y, z) = 0 (non-drawing target) pixels or Mask (x, y, z) = 1 (drawing) as a mask image to be combined (superimposed) with the MPR image ImageMPR. An image in which the pixels of the target) are filled with a predetermined color, or the pixels of Mask (x, y, z) = 0 (non-drawing target) and the pixels of Mask (x, y, z) = 1 (drawing target) are different. A color-coded image is generated and combined with an MPR image.

The process of generating the MPR image has been described above with reference to FIGS. 41 and 42. In the above example, the mask image generated based on the mask data Mask is combined with the MPR image, but as in the case of the volume rendered image and the MIP image, the drawing target defined by the mask data Mask is used. An MPR image displaying only pixels may be generated.

As described above, as an application example of the mask data Mask, an example of generating a three-dimensional reconstructed image (volume rendering image / MIP image / MPR image) using the mask data Mask has been described.

(Display example of 14.3D reconstruction image)
Finally, a display example of the three-dimensional reconstruction image (volume rendering image, MIP image, MPR image) generated by the three-dimensional reconstruction image generation device 2 of the present disclosure is shown.

FIG. 43 is a diagram showing a display example of a volume-rendered image of the thoracoabdominal region generated by the three-dimensional reconstructed image generation device 2 of the present disclosure.
FIG. 43 (a) is a volume-rendered image observed from the front, and FIG. 43 (b) is a volume-rendered image observed from the left side.
It can be seen that all bone regions including ribs and sternum are sufficiently removed, and visceral regions such as lungs, heart, blood vessels, liver, and kidneys, which are important for diagnosis, are clearly visualized.

FIG. 44 is a diagram showing a display example of a MIP image of the thoracoabdominal region generated by the three-dimensional reconstruction image generation device 2 of the present disclosure.
FIG. 44 (a) is a MIP image observed from the front, and FIG. 45 (b) is a MIP image observed from the left side.
Similar to the volume rendered image of FIG. 43, all bone regions including ribs and sternum are sufficiently removed, and visceral regions important for diagnosis such as lungs, heart, blood vessels, liver, and kidneys are clearly visualized. Can be confirmed.

FIG. 45 is a diagram showing a display example of an MPR image of the thoracoabdominal region generated by the three-dimensional reconstruction image generator 2 of the present disclosure. As shown in the figure, an MPR image in which mask images are combined (superimposed) for each cross section of Axial (body axis cross section), Coronal (coronal section), and Sagittal (sagittal section) is displayed. The mask image in the figure is an image in which the non-drawing target is filled with "red" (the figure is displayed in "gray" because it is a grayscale image). In this way, since the mask image is compositely displayed on the MPR image, the user can easily grasp the drawn / non-drawn area in the tomographic image.

As a comparative example, FIG. 46 shows a display example of a volume-rendered image of the thoracoabdominal region generated by using the opacity control method in Japanese Patent Application No. 2018-124125 already proposed by the present inventor. FIG. 44 (a) is a volume-rendered image observed from the front, and FIG. 44 (b) is a volume-rendered image observed from the left side. Japanese Patent Application No. 2018-124125 approximates the outline of the subject area with a geometric shape, obtains the correction magnification of each pixel from the parameters defining the geometric shape, and obtains the correction magnification of each pixel, and in the process of volume rendering (raycasting processing), the signal value of each pixel. Is a method of controlling the opacity by multiplying the opacity of each pixel obtained by passing through the opacity curve by the correction magnification.

In the volume rendered image (FIG. 43) obtained by the method of the present disclosure, it can be confirmed that the bone region is sufficiently removed and the internal organ region and the like inside the bone region are clearly visualized.
Since the method of Japanese Patent Application No. 2018-124125 corrects the opacity curve set by the user, when the opacity curve is changed, the form of bone removal changes as in the case of changing the CT image, and it is necessary to redo the parameter adjustment. It was. In addition, when the capacity curve is changed, the opacity of the organ to be observed is multiplied by the correction magnification, and the organ may become translucent (blurred) unlike the user's intention. In that case, the ray casting method has a problem that it leads to an increase in processing load.

On the other hand, in the method of the present disclosure, the opacity is not controlled in the process of volume rendering, and the mask data is directly generated based on the weight value (value corresponding to the correction magnification of Japanese Patent Application No. 2018-124125). Since the rendering process is performed by general mask reference, there is no need to adjust parameters due to changes in the opacity curve, and the organs are not translucent (blurred) unlike the user's intention. Furthermore, since it is not necessary to multiply the correction magnification in the process of volume rendering, the processing load is reduced.

Further, since the method of Japanese Patent Application No. 2018-124125 corrects the color map (opacity curve), it cannot be applied to MIP images and MPR images that do not use the color map. On the other hand, in the method of the present disclosure, as described above, the method of controlling the opacity in the process of volume rendering is not adopted, and the mask data is directly based on the weight value (value corresponding to the correction magnification of Japanese Patent Application No. 2018-124125). Is generated, and a three-dimensional reconstructed image generation process is executed by a general mask reference. Therefore, it can be applied to a MIP image or MPR image generation process that does not use a color map.

Further, in the method of Japanese Patent Application No. 2018-124125, since the subject area (rectangular area 41 in FIG. 5) is calculated based on all the pixel areas where the signal value is larger than the predetermined threshold value, the bed, headrest, etc. are included in the subject area. It was sometimes included. On the other hand, in the method of the present disclosure, the subject area (rectangular area 42 in FIG. 5) is calculated in consideration of the signal value distribution, so that the area not including the bed, headrest, etc. is extracted with high accuracy as the subject area. Can be done.

Further, in the method of Japanese Patent Application No. 2018-124125, the outer contour of the subject area is approximated by three ellipses, but since the diameter of the central ellipse for removing the vertebra is set to be small, the nearby aorta and the like are simultaneously removed. There was a problem that it was easy. In addition, the number of setting parameters increases, it takes time and effort to adjust the parameters, and it may be difficult to adapt them to various CT images. On the other hand, in the method of the present disclosure, the number of setting parameters can be suppressed and the ribs can be suppressed by approximating the outer shell of the subject area with one deformed ellipse whose diameter can be changed between the sternum side and the vertebrae side. , All bone regions including the sternum can be automatically removed, and visceral regions important for diagnosis such as lungs, heart, blood vessels, liver, and kidneys can be automatically extracted with high accuracy.

Although the preferred embodiments according to the present disclosure have been described above with reference to the accompanying drawings, the present disclosure is not limited to such examples. It is clear that a person skilled in the art can come up with various modifications or modifications within the scope of the technical ideas disclosed in the present application, and these also naturally belong to the technical scope of the present disclosure. Understood.

It is also possible to configure a subject area calculation device that executes a process of calculating a subject area from a tomographic image. This subject area calculation device includes the image acquisition unit 21, the subject area calculation unit 31, and the subject area correction unit 32 of FIG. 1, and is realized by the same hardware configuration as that of FIG. Then, the control unit of the subject area calculation device calculates the subject area from the tomographic image by executing the processes of step S1 in FIG. 3, step S21 in FIG. 4, and step S22.

1 ... Mask creation device 2 ... 3D reconstruction image generation device 21 ... Image acquisition unit 22 ... Drawing weight table creation unit 23 ... Mask data generation unit 24: Clipping area setting unit 25: Effective rectangle area calculation unit 26: 3D reconstruction image generation unit 27: 3D reconstruction image display Part 31: Subject area calculation unit 32: Subject area correction unit 33: Subject shape calculation unit 34: Weight value calculation unit 41: Rectangular area ( 1st rectangular area)
42 ... Rectangular area (second rectangular area)
50 ... Subject candidate area 51 ... Volume rendering image generation unit 52 ... MIP image generation unit 53 ... MPR image generation unit 510 ... Color map setting unit

Claims (24)

Image acquisition means for acquiring multiple tomographic images,
A subject area calculation means that calculates a subject area based on the signal value of each pixel for each tomographic image,
For each tomographic image, a subject shape calculation means that approximates the calculated outline of the subject area with a predetermined geometric shape and calculates the parameters of the geometric shape,
A mask data generation means that calculates a mask value that defines whether each pixel of the tomographic image is included inside or outside the geometric shape for each tomographic image, and generates mask data that holds the mask value for each pixel. When,
A mask generator comprising. For each tomographic image, a weight value calculating means for calculating a weight value for drawing each pixel of the tomographic image based on the parameter of the geometric shape is further provided.
The mask data generation means calculates a mask value obtained by binarizing the weight value calculated by the weight value calculation means for each pixel with a predetermined threshold value for each tomographic image, and holds the mask value for each pixel. The mask generation device according to claim 1, wherein the mask data to be generated is generated. The subject area calculation means calculates, as the subject area, a first rectangular area which is a rectangular area surrounding the outermost region in which pixels having a signal value larger than a predetermined threshold value exist.
The mask generator according to claim 1 or 2. The subject area calculation means is a second rectangular area defined inside a first rectangular area which is a rectangular area surrounding the outermost side of a region in which pixels having a signal value larger than a predetermined threshold value exist as the subject area. To calculate,
The mask generator according to claim 1 or 2. The subject area calculation means said that the center coordinates of the second rectangular area are barycentric coordinates weighted by the signal values of pixels whose signal values in the first rectangular area are larger than the predetermined threshold value. Calculate the second rectangular area,
The mask generator according to claim 4. The subject area calculation means
The second center of gravity, which is the coordinates of the center of gravity weighted by the signal values of pixels whose signal values are larger than the predetermined threshold value in the range of the first rectangular area in four directions of up, down, left, and right starting from the center coordinates of the second rectangular area. Is calculated in 4 places,
The lateral size of the second rectangular area is calculated based on the calculated second center of gravity in the left direction and the second center of gravity in the right direction.
The vertical size of the second rectangular region is calculated based on the calculated second center of gravity in the upward direction and the second center of gravity in the downward direction.
The mask generator according to claim 5. The subject area calculation means
A third barycentric coordinate whose signal value is weighted by the signal value of a pixel whose signal value is larger than the predetermined threshold value in the range of the first rectangular region in four directions of up, down, left, and right starting from the calculated four second centers of gravity. Calculate the center of gravity at 4 points
The lateral size of the second rectangular area is calculated based on the calculated third center of gravity in the left direction and the third center of gravity in the right direction.
The vertical size of the second rectangular region is calculated based on the calculated upper third center of gravity and the lower third center of gravity.
The mask generator according to claim 6. The subject area calculation means
The first lateral size, which is the distance from the calculated second center of gravity or the third center of gravity in the left direction to the center coordinates of the second rectangular region, is calculated.
The second lateral size, which is the distance from the calculated second center of gravity or the third center of gravity in the right direction to the center coordinates of the second rectangular region, is calculated.
The first vertical size, which is the distance from the calculated upper second center of gravity or the third center of gravity to the center coordinates of the second rectangular region, is calculated.
The second vertical size, which is the distance from the calculated downward second center of gravity or the third center of gravity to the center coordinates of the second rectangular region, is calculated.
The mask generator according to claim 6 or 7. A geometry is a geometry defined by center coordinates, horizontal size, and vertical size.
The subject shape calculation means
As the center coordinates of the geometric shape, the center coordinates of the first or second rectangular region are calculated.
As the horizontal size of the geometric shape, the horizontal size of the first or second rectangular region is calculated.
As the vertical size of the geometric shape, the vertical size of the first or second rectangular region is calculated.
The weight value calculating means calculates the weight value based on the calculated center coordinates of the geometric shape, the size in the horizontal direction, and the size in the vertical direction.
The mask generator according to any one of claims 3 to 7, wherein claims 3 and 4 cite claim 2. The subject shape calculating means multiplies the calculated horizontal size and vertical size of the geometric shape by a magnification with respect to a predetermined horizontal size and a magnification with respect to the vertical size to obtain the horizontal size and the vertical size. Correct the size of the direction,
The mask generator according to claim 9. The geometric shape is the center coordinate, the first horizontal size which is the distance from the center coordinate to the left end, the second horizontal size which is the distance from the center coordinate to the right end, and the distance from the center coordinate to the upper end. It is a geometric shape defined by the first vertical size and the second vertical size, which is the distance from the center coordinate to the lower end.
The subject shape calculation means
As the center coordinates of the geometric shape, the center coordinates of the second rectangular region are calculated.
As the first lateral size of the geometric shape, the first lateral size of the second rectangular region is calculated.
As the second lateral size of the geometric shape, the second lateral size of the second rectangular region is calculated.
As the first vertical size of the geometric shape, the first vertical size of the second rectangular region is calculated.
As the second vertical size of the geometric shape, the second vertical size of the second rectangular region is calculated.
The weight value calculating means is based on the calculated center coordinates of the geometric shape, the first horizontal size, the second horizontal size, the first vertical size, and the second vertical size. Calculate the value,
The mask generator according to claim 8, wherein claim 4 cites claim 2. The subject shape calculation means
Multiplying the calculated first lateral size and second lateral size of the geometric shape by a magnification for different lateral sizes to obtain the first lateral size and the second lateral size. Correct the size and
The calculated first vertical size and second vertical size of the geometric shape are multiplied by a magnification for different vertical sizes to obtain the first vertical size and the second vertical size. Correct the size,
The mask generator according to claim 11. For each tomographic image, a subject area correction means for correcting the subject area by a correction method according to the ratio of the vertical size and the horizontal size of the subject area calculated by the subject area calculation means is further provided.
The mask generator according to any one of claims 1 to 12. The subject area correction means calculates an aspect ratio obtained by dividing the vertical size of the second rectangular area calculated by the subject area calculation means by the horizontal size.
When the calculated aspect ratio is smaller than a predetermined threshold value, the aspect ratio is corrected by multiplying each of the horizontal size and the vertical size of the second rectangular area calculated by the subject area calculation means by the aspect ratio. And
When the calculated aspect ratio is equal to or greater than the predetermined threshold value, the area ratio obtained by dividing the number of pixels whose signal value is larger than the predetermined threshold value inside the first rectangular region by the area of the second rectangular region. Is calculated, and the area ratio is multiplied by each of the horizontal size and the vertical size of the second rectangular area calculated by the subject area calculation means to correct the problem.
The mask generator according to claim 12, wherein claim 13 cites any of claims 4 to 7. The subject area correction means calculates an aspect ratio obtained by dividing the vertical size of the second rectangular area calculated by the subject area calculation means by the horizontal size.
When the calculated aspect ratio is smaller than a predetermined threshold value, the first horizontal size, the second horizontal size, and the first vertical of the second rectangular region calculated by the subject shape calculating means. The size in the direction and the size in the second vertical direction are each multiplied by the aspect ratio to be corrected.
When the calculated aspect ratio is equal to or greater than the predetermined threshold value, the area ratio obtained by dividing the number of pixels whose signal value is larger than the predetermined threshold value inside the first rectangular region by the area of the second rectangular region. The first horizontal size, the second horizontal size, the first vertical size, and the second vertical direction of the second rectangular area calculated by the subject area calculation means. The area ratio is multiplied by each of the sizes of
The mask generator according to claim 12, wherein claim 13 cites claim 8. The weight value calculating means calculates the weight value so that the center of the geometric shape is the maximum, the outside of the geometric shape is the minimum, and the inside of the geometric shape becomes smaller as the distance from the center increases.
The mask generator according to claim 2. The weight value calculating means corrects so that the weight value corresponding to the tomographic image becomes smaller as the slice position of the tomographic image is located from the center to the end.
The mask generator according to claim 2. The mask data generation means calculates a mask value obtained by binarizing the signal value of each pixel weighted by the weight value with a predetermined threshold value for each tomographic image, and obtains mask data holding the mask value of each pixel. Generate,
The mask generator according to claim 2. Image acquisition means for acquiring multiple tomographic images,
A subject area calculation means that calculates a subject area based on the signal value of each pixel for each tomographic image,
For each tomographic image, a subject shape calculation means that approximates the calculated outline of the subject area with a predetermined geometric shape and calculates the parameters of the geometric shape,
A mask data generation means that calculates a mask value that defines whether each pixel of the tomographic image is included inside or outside the geometric shape for each tomographic image, and generates mask data that holds the mask value of each pixel. When,
A three-dimensional voxel formed by arranging each pixel of the tomographic image in a voxel space, a three-dimensional reconstructed image generating means for generating a three-dimensional reconstructed image based on the mask data, and a three-dimensional reconstructed image generating means.
A three-dimensional reconstructed image generator including. The three-dimensional reconstructed image generation means sets whether or not to use the opacity of each voxel in the calculation based on the mask data.
The three-dimensional reconstructed image generator according to claim 19. The three-dimensional reconstruction image generation means sets whether or not to use the signal value of each voxel in the calculation based on the mask data.
The three-dimensional reconstructed image generator according to claim 19. The three-dimensional reconstructed image generating means generates an image obtained by synthesizing the three-dimensional reconstructed image generated based on the tomographic image with the mask image generated based on the mask data as the three-dimensional reconstructed image. To do,
The three-dimensional reconstructed image generator according to claim 19. The computer
Image acquisition steps to acquire multiple tomographic images and
For each tomographic image, a subject area calculation step that calculates the subject area based on the signal value of each pixel, and
For each tomographic image, the subject shape calculation step of approximating the calculated outline of the subject area with a predetermined geometric shape and calculating the parameters of the geometric shape,
For each tomographic image, a mask data generation step of calculating a mask value that defines whether each pixel of the tomographic image is included inside or outside the geometric shape and generating mask data that holds the mask value of each pixel. When,
How to generate a mask to perform. Computer,
Image acquisition means for acquiring multiple tomographic images,
Subject area calculation means that calculates the subject area based on the signal value of each pixel for each tomographic image,
A subject shape calculation means that approximates the calculated outline of the subject area with a predetermined geometric shape for each tomographic image and calculates the parameters of the geometric shape.
A mask data generation means that calculates a mask value that defines whether each pixel of the tomographic image is included inside or outside the geometric shape for each tomographic image, and generates mask data that holds the mask value of each pixel. ,
A program that functions as.