WO2023007846A1 - 補正装置、システム、方法およびプログラム - Google Patents

補正装置、システム、方法およびプログラム Download PDF

Info

Publication number
WO2023007846A1
WO2023007846A1 PCT/JP2022/014455 JP2022014455W WO2023007846A1 WO 2023007846 A1 WO2023007846 A1 WO 2023007846A1 JP 2022014455 W JP2022014455 W JP 2022014455W WO 2023007846 A1 WO2023007846 A1 WO 2023007846A1
Authority
WO
WIPO (PCT)
Prior art keywords
absorption coefficient
linear absorption
incident
image
ray
Prior art date
Application number
PCT/JP2022/014455
Other languages
English (en)
French (fr)
Inventor
卓見 太田
Original Assignee
株式会社リガク
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 株式会社リガク filed Critical 株式会社リガク
Priority to JP2023538266A priority Critical patent/JPWO2023007846A1/ja
Priority to DE112022003708.3T priority patent/DE112022003708T5/de
Priority to CN202280052823.1A priority patent/CN117716231A/zh
Publication of WO2023007846A1 publication Critical patent/WO2023007846A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N23/00—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
    • G01N23/02—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material
    • G01N23/04—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and forming images of the material
    • G01N23/046—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and forming images of the material using tomography, e.g. computed tomography [CT]
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2223/00—Investigating materials by wave or particle radiation
    • G01N2223/40—Imaging
    • G01N2223/419—Imaging computed tomograph

Definitions

  • the present invention relates to a correction device, system, method and program for correcting artifacts.
  • a CT device reconstructs a CT image from multiple projection images acquired while rotating the sample or gantry.
  • a CT apparatus uses continuous X-rays, and since the linear absorption coefficient at each energy differs for each material, artifacts due to beam hardening effect occur in the reconstructed CT image.
  • Patent Document 1 discloses a technique for reducing artifacts by segmenting a reconstructed image, identifying a substance, and repeating reconstruction using the linear absorption coefficient of the substance.
  • Patent Document 2 discloses an original CT image in which metal artifacts are generated, forward projection calculation considering the energy dependence of the X-ray absorption coefficient of the material, and back projection calculation using an image reconstruction algorithm for a single wavelength. There is disclosed a technique for reducing artifacts by repeatedly performing forward projection calculation and back projection calculation so as to reduce the difference from the CT image obtained by performing the above.
  • the present invention has been made in view of such circumstances, and provides a correction device, system, method, and program capable of reducing calculation costs for correcting artifacts due to beam hardening effects in CT image reconstruction.
  • the purpose is to
  • the correction apparatus of the present invention is a correction apparatus for correcting artifacts due to beam hardening effects in reconstructing CT images, wherein the incident X-ray distribution is obtained by an acquisition unit, a linear absorption coefficient model acquisition unit that acquires a linear absorption coefficient model in which the energy dependence of the linear absorption coefficient is represented by a scale factor including a parameter, a projection image acquisition unit that acquires a projection image, and the projection image; and a correction unit that performs correction using the incident X-ray distribution and the linear absorption coefficient model.
  • the linear absorption coefficient model is characterized by being represented by the product of the linear absorption coefficient at the reference energy and the scale factor.
  • the parameter of the linear absorption coefficient model is one parameter.
  • the scale factor is characterized by being represented by a power function having a power exponent as the parameter.
  • the parameters of the linear absorption coefficient model are determined based on the energy range of the incident X-ray distribution.
  • the parameters of the linear absorption coefficient model are determined from the linear absorption coefficients of representative element groups.
  • the incident X-ray distribution acquisition unit is characterized by acquiring the incident X-ray distribution based on the number, intensity and energy value of monochromatic X-rays.
  • the correction apparatus of the present invention includes a reconstruction unit that performs reconstruction based on the projection data corrected by the correction unit and generates a CT image, and a display that displays the CT image on a display device. and a part.
  • the system of the present invention includes an X-ray source that generates X-rays, a detector that detects X-rays, and a rotation control unit that controls rotation of the X-ray source and the detector or a sample. and a correction device according to any one of (1) to (8) above.
  • the method of the present invention is a method of correcting artifacts due to the beam hardening effect in CT image reconstruction, comprising the step of obtaining the incident X-ray distribution, and the energy dependence of the linear absorption coefficient as a parameter obtaining a linear absorption coefficient model expressed by a scale factor comprising: obtaining a projection image; and correcting the projection image using the incident X-ray distribution and the linear absorption coefficient model It is characterized by
  • the program of the present invention is a program for correcting artifacts due to the beam hardening effect in reconstructing CT images, and is a program for acquiring the incident X-ray distribution and the energy dependency of the linear absorption coefficient as a parameter. a process of acquiring a linear absorption coefficient model represented by a scale factor including the It is characterized by executing
  • FIG. 4 is a graph showing an example of incident X-ray distribution; 4 is a graph showing the mass absorption coefficient of titanium; 4 is a graph showing the mass absorption coefficient of iron; 1 is a schematic diagram showing an example of the configuration of an entire system; FIG. FIG. 4 is a schematic diagram showing a modification of the configuration of the entire system; It is a block diagram which shows an example of a structure of a processing apparatus and a correction apparatus.
  • FIG. 11 is a block diagram showing a modification of the configuration of the processing device and the correction device;
  • FIG. 11 is a block diagram showing a modification of the configuration of the processing device and the correction device;
  • 4 is a flow chart showing an example of the operation of the correction device;
  • 4 is a flow chart showing an example of system operation;
  • 7 is a flow chart showing a modification of the operation of the system;
  • (a) and (b) are a CT image of the sample 1 reconstructed using the uncorrected projection image and a CT image of the sample 1 reconstructed using the corrected projection image, respectively.
  • 13(a) and 13(b) are graphs showing line profiles on the straight line AB of the CT images of FIGS.
  • a CT apparatus irradiates a sample with cone-shaped or parallel beams of X-rays from all angles, and obtains a distribution of X-ray absorption coefficients, ie, a projection image, with a detector.
  • the CT apparatus rotates the sample stage with respect to the fixed X-ray source and detector, or rotates the gantry in which the X-ray source and detector are integrated. is configured as
  • the distribution of the linear absorption coefficient f of the sample can be estimated from the density of the projected image of the sample obtained by projecting from various angles.
  • Obtaining a three-dimensional linear absorption coefficient distribution from a two-dimensional projected image is called reconstruction.
  • Reconstruction basically back-projects the projected image.
  • a CT device uses continuous X-rays as incident X-rays and irradiates them onto the sample to measure the projected image.
  • the linear absorption coefficient of a substance depends on energy, and high-energy X-rays tend to be less absorbed. Therefore, when continuous X-rays pass through the sample as incident X-rays, the energy distribution after passage is no longer similar to the original distribution, and the center of gravity shifts to the higher energy side. This phenomenon is called beam hardening.
  • Reconstruction is usually performed assuming monochromatic X-rays.
  • continuous X-rays with a wide energy distribution are used.
  • low-energy X-rays are more likely to be absorbed and attenuated by the sample than high-energy X-rays. Therefore, the linearity between the transmission distance of X-rays and the attenuation of X-rays is lost, and non-linearity appears. Linearity is preserved when monochromatic X-rays pass through matter.
  • artifacts arise due to the discrepancy between the reconstruction assumptions and the actual phenomenon. Such artifacts are called beam hardening artifacts.
  • the hardware correction method basically narrows the width of the incident X-ray distribution to make the artifacts due to the beam hardening effect less noticeable.
  • a filter is used to extract only X-rays near a specific energy, and a mirror is used to extract monochromatic X-rays.
  • the Helgason-Ludwig condition method reduces artifacts by nonlinearly correcting the intensity of a projected image using physical conditions such as the law of conservation of linear absorption coefficients.
  • correction is made based on the assumption that the projection image data is consistent. This method cannot be applied when phenomena other than the beam hardening effect occur, such as when there is a shift in the optical system such as a center shift.
  • the calculation of the correction amount uses the integration of the projected image, the compression of the information may cause the magnitude relationship of the X-ray intensity to be reversed or the pixel value to become negative. , the correction may not go well.
  • the present invention corrects the projection image using a linear absorption coefficient model that expresses the energy dependence of the incident X-ray distribution and the linear absorption coefficient with a scale factor including parameters (auxiliary variables). It can be done for each detector pixel. This eliminates the need for segmentation and the use of rigorous sample information. In addition, since the processing is performed before reconstruction, the calculation cost for correcting artifacts due to the beam hardening effect in reconstruction of the CT apparatus can be reduced. Furthermore, it can also be applied when phenomena other than the beam hardening effect occur.
  • the linear absorption coefficient model is a function selected to approximate the energy dependence of the linear absorption coefficient distribution. This model is expressed, for example, by multiplying the scale factor by the distribution of linear absorption coefficients at some energy E0.
  • the scale factor is a non-negative function s(E) with respect to the energy E that has a value of 1 at a certain reference energy E0.
  • FIGS. 1(a) and 1(b) respectively show a conceptual diagram of X-rays transmitted through a sample, and the relationship between a portion of discrete monochromatic X-rays and the measurement space (voxels) of the sample measured by CT measurement. It is a conceptual diagram.
  • the X-ray intensity I detected by the detector is defined by the thickness of the object l, the ray absorption coefficient ⁇ , and the incident X-ray Assuming that the intensity is I0 , the formula (1) is given by (Lambert-Beer law).
  • the nonlinearity of X-ray transmission distance and X-ray attenuation is expressed by summing the attenuation of each monochromatic X-ray with the total energy.
  • the intensity I detected by each detector pixel is represented by the following equation (2) as the sum of the intensities of monochromatic X-rays attenuated by substances.
  • the energy distribution of continuous X-rays can be represented as a distribution as shown in FIG. 2 by the intensity Ik of each of the N monochromatic X-rays and the energy Ek of each monochromatic X-ray.
  • FIG. 2 is a graph showing an example of incident X-ray distribution.
  • the energy value and number may be set to 10 at equal intervals such as 10, 15, 20, . . . , 55 keV.
  • the energy value and the number of monochromatic X-rays may be determined such that areas where the change is large are taken densely and areas where the intensity is low are taken sparsely.
  • the energy value, the number, and the intensity of monochromatic X-rays corresponds to setting Ik in Equation (2).
  • the number of monochromatic X-rays to be set must be 3 or more, preferably 5 or more. Since the incident X-ray distribution can be known to some extent from the tube voltage and filter information, such knowledge may be used. Alternatively, the user may arbitrarily select or designate.
  • the energy dependence of f(E k ) on the right side of Equation (2) is replaced with a linear absorption coefficient model expressed by a scale factor including parameters.
  • the energy dependence of the linear absorption coefficient is represented by a scale factor including parameters.
  • the energy range (domain of s(E k )) for calculating the scale factor s(E k ) is set to include the energy range in which the incident X-ray distribution is set.
  • f(E k ) is, as shown in the following equation ( 3 ), a scale factor s ( E k ).
  • Equation (2) Replacing f(E k ) on the right side of Equation (1) with the equation on the right side of Equation (3), Equation (2) becomes Equation (4) below.
  • Equation (2) computes the line integral of the total energy projection image. Equation (2) therefore requires calculating the value of the line integral for each energy.
  • equation (4) computes the line integral of the projection image at some energy E0 . Therefore, equation (4) gives information on other energies from the absorption coefficient of the reference energy by the scale factor. The number of unknowns is reduced, making the calculations easier. Also, it is preferable to set the parameter of the scale factor to one parameter. This makes the calculation even easier. However, multiple parameters may be set for accurate calculation.
  • a projected image of the distribution of linear absorption coefficients at a certain energy E0 can be obtained.
  • the value of this line integral corresponds to the projection of the distribution of the linear absorption coefficient at energy E0.
  • This is the corrected projected image. That is, the continuous X-ray projection image is converted into a projection image with a certain energy E0 . Reconstruction using such a corrected projection image makes it possible to obtain a CT image with reduced artifacts due to the beam hardening effect.
  • the reference energy E 0 may be selected in any way within the domain of s(E k ). For example, it may be the lower limit of the energy range of the incident X-ray distribution, or the average value of the energy range of the incident X-ray distribution.
  • the scale factor of the linear absorption coefficient model is preferably represented by a power function with a power exponent as a parameter. This further facilitates calculation for correction. That the scale factor of the linear absorption coefficient model is represented by a power function means that s(E k ) in Equation (3) is represented by Equation (5) below, for example.
  • the scale factor included in the linear absorption coefficient model is not limited to this. Any functional form that can approximate the energy dependence of the linear absorption coefficient may be used. For example, exponential functions, logarithmic functions, etc. may be used.
  • the scale factors included in the linear absorption coefficient model are preferably selected in terms of both the degree of approximation of the energy dependence of the linear absorption coefficient and the subsequent computational cost.
  • Equation (6) the linear absorption coefficient model is expressed as in Equation (6) below. Also, using a linear absorption coefficient model, the intensity I detected by each detector pixel is represented by Equation (7).
  • the linear absorption coefficient model is preferably determined based on the energy range of the incident X-ray distribution.
  • an appropriate linear absorption coefficient model can be determined according to the linear absorption coefficients of representative element groups, and artifacts due to the beam hardening effect can be corrected more appropriately. That the linear absorption coefficient model is determined based on the energy range of the incident X-ray distribution means that s(E k ) and E 0 in Equation (4) are changed according to the energy range of the incident X-ray distribution. means that you can When using the linear absorption coefficient model, the linear absorption coefficient of each monochromatic X-ray can be calculated from the behavior of the linear absorption coefficient of representative elements in the assumed incident X-ray distribution.
  • the energy dependence of the linear absorption coefficient of elements (Ti, Fe, etc.) considered to cause metal artifacts is Knowing that, from the ratio of the scale factor at the reference energy E 0 and the energy E k of each monochromatic X-ray, the linear absorption coefficient f (E k ) can be calculated. At this time, only parameter values may be changed in the same function.
  • the linear absorption coefficient model is represented by a power function, it is preferable to determine and change the parameter ⁇ as follows.
  • the parameter ⁇ in equations (6) and (7) is preferably set so that the attenuation of the mass absorption coefficient of the element applied to the correction is similar to that of the element contained in the sample. Since the mass absorption coefficient has characteristics such as a different slope for each element, a different position of the absorption edge, and a different slope before and after the absorption edge, it is preferable to set ⁇ by the following method.
  • the representative element is preferably an element contained in the sample.
  • Representative elements to be selected are preferably selected from light element group, medium element group and heavy element group.
  • the light element group can be H to Ne
  • the middle element group can be Na to Ca
  • the heavy element group can be Sc or later.
  • the boundaries of the light element group, medium element group, and heavy element group can be changed depending on the field.
  • the parameter ⁇ is preferably calculated by fitting a power function, but it does not necessarily have to match the fitting result, and may be arbitrarily set using the fitting result as a guideline.
  • the average slope is calculated from the slope before the absorption edge ⁇ R and the slope after the absorption edge ⁇ L and used as ⁇ . is preferred.
  • the linear absorption coefficient model As described above, by expressing the linear absorption coefficient model as a power factor and appropriately setting the value of the parameter ⁇ within a predetermined range, the calculation for correction becomes easier, and the beam Artifacts due to hardening effects can be effectively corrected. Moreover, the calculation cost for that can be reduced significantly.
  • FIG. 5 is a schematic diagram showing the configuration of the entire system 100 including the CT device 200 and the processing device 300 connected thereto, the correction device 400, the input device 510 and the display device 520.
  • the CT apparatus 200 shown in FIG. 5 has a configuration in which the sample is rotated with respect to the X-ray source 260 and the detector 270, but is not limited to this, and the X-ray source and the detector are integrated. A configuration in which the gantry is rotated may also be used.
  • the CT apparatus 200 can use both a parallel beam apparatus and a cone beam apparatus.
  • the processing device 300 is connected to the CT device 200 and controls the CT device 200 and processes acquired data.
  • Correction device 400 corrects the projected image.
  • the processing device 300 and the correction device 400 may be PC terminals or servers on the cloud.
  • the input device 510 is, for example, a keyboard and a mouse, and performs input to the processing device 300 and the correction device 400 .
  • the display device 520 is, for example, a display, and displays a projected image or the like.
  • FIG. 5 is a schematic diagram showing a modification of the overall system configuration. Using such a system can reduce the computational cost of correcting artifacts due to beam hardening effects in CT image reconstruction.
  • the CT apparatus 200 includes a rotation control unit 210, a sample stage 250, an X-ray source 260, a detector 270 and a driving section 280.
  • X-ray CT imaging is performed by rotating a sample stage 250 placed between an X-ray source 260 and a detector 270 .
  • the X-ray source 260 and detector 270 may be installed on a gantry (not shown), and the gantry may be rotated with respect to the sample fixed on the sample stage 250 .
  • the CT device 200 drives the sample stage 250 at the timing instructed by the processing device 300 and acquires a projected image of the sample.
  • the measurement data are transmitted to the processing device 300 .
  • the CT apparatus 200 is suitable for use in precision industrial products such as semiconductor devices, it can be applied not only to industrial equipment but also to animal equipment.
  • the X-ray source 260 emits X-rays toward the detector 270 .
  • the detector 270 has a light-receiving surface that receives X-rays, and can measure the intensity distribution of X-rays that have passed through the sample using a large number of pixels.
  • the rotation control unit 210 rotates the sample table 250 at a speed set during CT imaging by the drive unit 280 .
  • FIG. 7 is a block diagram showing the configuration of the processing device 300 and the correction device 400.
  • the processing device 300 is configured by a computer having a CPU (Central Processing Unit), a ROM (Read Only Memory), a RAM (Random Access Memory), and a memory connected to a bus.
  • a processing device 300 is connected to the CT device 200 to receive information.
  • the processing device 300 includes a measurement data storage unit 310, a device information storage unit 320, a reconstruction unit 330, and a display unit 340. Each unit can transmit and receive information through the control bus L.
  • FIG. The input device 510 and display device 520 are connected to the CPU via appropriate interfaces.
  • the measurement data storage unit 310 stores measurement data acquired from the CT device 200 .
  • the measurement data includes rotation angle information and corresponding projected images.
  • the device information storage unit 320 stores device information acquired from the CT device 200 .
  • the device information includes device name, beam shape, geometry at the time of measurement, scanning method, and the like.
  • the reconstruction unit 330 reconstructs a CT image from the target projection image.
  • the display unit 340 causes the display device 520 to display the reconstructed CT image. This allows the user to check the CT image based on the corrected projection image. Also, the user can give instructions and designations to the processing device, correction device, etc. based on the CT image.
  • the correction device 400 is composed of a computer having a CPU, a ROM, a RAM, and a memory connected to a bus.
  • the correction device 400 may be directly connected to the CT device 200 or may be connected to the CT device 200 via the processing device 300 . Further, the correction device 400 may receive information from the CT device 200 or may receive information from the processing device 300 .
  • the correction device 400 may be configured as a part of functions included in the processing device 300, or as shown in FIG. may be Also, the correction device 400 may have a part of the functions of the processing device 300 .
  • the correction device 400 includes an incident X-ray distribution acquisition unit 410, a linear absorption coefficient model acquisition unit 420, a projection image acquisition unit 430, and a correction unit 440. Each unit can transmit and receive information through the control bus L.
  • FIG. If the correction device 400 and the processing device 300 have different configurations, the input device 510 and the display device 520 are also connected to the CPU of the correction device 400 via appropriate interfaces. In this case, input device 510 and display device 520 may be different than those connected to processing device 300 .
  • the incident X-ray distribution acquisition unit 410 acquires the incident X-ray distribution.
  • the incident X-ray distribution acquisition unit 410 preferably acquires the incident X-ray distribution based on the number, intensity, and energy value of monochromatic X-rays. This makes it possible to acquire a more appropriate incident X-ray distribution according to the situation.
  • the incident X-ray distribution acquisition unit 410 preferably acquires the incident X-ray distribution based on user's designation. This allows the user to select a more appropriate incident X-ray distribution according to the situation.
  • the incident X-ray distribution is specified by the user, for example, one of a plurality of incident X-ray distributions stored in advance is specified by mouse operation, the point on the incident X-ray distribution graph is moved, and the graph is displayed. It is preferable to use a UI function that can set the incident X-ray distribution, such as increasing or decreasing the number of upper points.
  • the incident X-ray distribution may be preset. Also, the incident X-ray distribution may be automatically selected and acquired according to device information or the like.
  • the linear absorption coefficient model acquisition unit 420 acquires a linear absorption coefficient model in which the energy dependence of the linear absorption coefficient is represented by a scale factor including parameters.
  • the linear absorption coefficient model acquisition unit 420 preferably acquires the linear absorption coefficient model based on user designation.
  • a storage unit (not shown) of the correction device 400 or the processing device 300 stores a function form representing the energy dependence of the linear absorption coefficient. For example, it is stored as an equation representing a model of the linear absorption coefficient f(E k ), or an equation representing the reference energy (E 0 ) and the scale factor s(E k ). Also, the group of elements that can be contained in the sample and the parameter ⁇ of the linear absorption coefficient model are correlated and stored.
  • the linear absorption coefficient model is designated by the user, for example, one is designated from a plurality of functional forms representing linear absorption coefficients stored in advance by mouse operation, and the parameter ⁇ is set by selecting an element group. .
  • the parameter ⁇ may be configured so that the user can directly input a value. In this way, it is preferable to use a UI function that can set the linear absorption coefficient model.
  • a linear absorption coefficient model may be set in advance. Also, the linear absorption coefficient model may be automatically selected and acquired according to sample information and the like.
  • the projection image acquisition unit 430 acquires projection images from the CT device 200 or the processing device 300 .
  • the corrector 440 corrects the projected image using the incident X-ray distribution and the linear absorption coefficient model. As a result, a CT image with reduced artifacts due to the beam hardening effect can be obtained.
  • a sample is placed in the CT apparatus 200, and a projection image is acquired while irradiating the sample with X-rays by repeating movement of the rotating shaft and projection of X-rays under predetermined conditions.
  • the CT apparatus 200 transmits apparatus information such as the scanning method and the acquired projection image to the processing apparatus 300 or the correction apparatus 400 as measurement data.
  • FIG. 10 is a flow chart showing an example of the operation of correction device 400 .
  • the correction device 400 acquires the incident X-ray distribution (step S1).
  • a linear absorption coefficient model is obtained (step S2).
  • a projection image is obtained (step S3).
  • the projected image is corrected using the acquired incident X-ray distribution and linear absorption coefficient model (step S4).
  • the projected image can be corrected using a linear absorption coefficient model in which the incident X-ray distribution and the energy dependence of the linear absorption coefficient are represented by scale factors including parameters. Acquisition of the incident X-ray distribution, acquisition of the linear absorption coefficient model, and acquisition of the projection image may be performed in any order.
  • FIG. 10 shows the operation of the correction device 400 only. This is sufficient for correcting the projection image, but in order to clarify the difference from the prior art, the operation including the measurement and reconstruction of the projection image will also be explained.
  • FIG. 11 is a flowchart showing an example of the operation of system 100.
  • the correction device 400 acquires the incident X-ray distribution (step T1).
  • a linear absorption coefficient model is obtained (step T2).
  • CT apparatus 200 measures the projection image (step T3). Measurement of the projection image includes moving the rotation axis and repeating X-ray projection.
  • correction device 400 acquires a projection image (step T4).
  • the projected image is corrected using the acquired incident X-ray distribution and linear absorption coefficient model (step T5).
  • Processing device 300 or correction device 400 then reconstructs a CT image using the corrected projection image (step T6). In this way, a reconstructed CT image can be obtained using the corrected projection images.
  • acquisition of the incident X-ray distribution, acquisition of the linear absorption coefficient model, and acquisition of the projection image may be performed in any order.
  • FIG. 12 is a flowchart showing a modification of the operation of system 100.
  • the correction device 400 acquires the incident X-ray distribution (step U1).
  • a linear absorption coefficient model is obtained (step U2).
  • the CT apparatus 200 moves the rotating shaft (step U3).
  • X-rays are projected (step U4).
  • the correction device 400 acquires the projection image projected earlier (step U5).
  • the projected image is corrected using the obtained incident X-ray distribution and linear absorption coefficient model (step U6).
  • step U7 determines whether the measurement is completed. If it is determined in step U7 that the measurement has not been completed (step U7, NO), the process returns to step U3. On the other hand, if processing device 300 or correction device 400 determines in step U7 that the measurement has been completed (step U7, YES), it reconstructs a CT image using the corrected projection image (step U8). . By doing so, the measurement of the projection image and the correction of the measured projection image can be processed in parallel, and the overall time can be shortened. Acquisition of the incident X-ray distribution and acquisition of the linear absorption coefficient model can be performed in the reverse order without any problem.
  • Example 1 Using the system 100 configured as described above, a cross section of a hydroxyapatite phantom (Micro-CT HA Phantom manufactured by QRM) (Sample 1) was observed. Rigaku's CTLabHX was used as the CT device 200, and the measurement was performed at a tube voltage of 60 kV.
  • FIGS. 13A and 13B are a CT image of the sample 1 reconstructed using the uncorrected projection image and a CT image of the sample 1 reconstructed using the corrected projection image, respectively.
  • . 14(a) and (b) are graphs showing line profiles on the straight line AB of the CT images of FIGS. 13(a) and (b), respectively. Reconstruction used the FDK method.
  • FIGS. 15A and 15B are a CT image of the sample 2 reconstructed using the uncorrected projection image and a CT image of the sample 2 reconstructed using the corrected projection image, respectively.
  • FIGS. 15(a) and (b) are graphs showing line profiles on the straight line AB of the CT images of FIGS. 15(a) and (b), respectively.
  • FIGS. 15(a) and (b) By comparing FIGS. 15(a) and (b), it was found that streak artifacts and dark bands were reduced by the correction. Also, by comparing FIGS. 16(a) and 16(b), it was confirmed that the cupping artifact at the end of the aluminum rod was reduced.
  • the correction device, system, method and program of the present invention can effectively correct artifacts due to the beam hardening effect in CT image reconstruction, and can reduce computational costs.

Landscapes

  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Pulmonology (AREA)
  • Radiology & Medical Imaging (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • General Health & Medical Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Immunology (AREA)
  • Pathology (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Analysing Materials By The Use Of Radiation (AREA)

Abstract

ＣＴ画像の再構成におけるビームハードニング効果によるアーチファクトを補正するための計算コストを低減できる補正装置、システム、方法およびプログラムを提供する。ＣＴ画像の再構成におけるビームハードニング効果によるアーチファクトを補正する補正装置４００であって、入射Ｘ線分布を取得する入射Ｘ線分布取得部４１０と、線吸収係数のエネルギー依存性をパラメータを含むスケール因子で表した線吸収係数モデルを取得する線吸収係数モデル取得部４２０と、投影像を取得する投影像取得部４３０と、前記投影像を前記入射Ｘ線分布および前記線吸収係数モデルを用いて補正する補正部４４０と、を備える。

Description

補正装置、システム、方法およびプログラム
　本発明は、アーチファクトを補正する補正装置、システム、方法およびプログラムに関する。
　ＣＴ装置は、試料またはガントリを回転させながら取得された複数の投影像からＣＴ画像を再構成する。ＣＴ装置では連続Ｘ線を使用しており、物質ごとに各エネルギーにおける線吸収係数が異なるため、再構成したＣＴ画像にはビームハードニング効果によるアーチファクトが生じる。
　このようなビームハードニング効果によるアーチファクトを低減するため、従来、ハードウェアやソフトウェアによる補正が行なわれてきた。ソフトウェアによる補正方法として、例えば、特許文献１は、再構成画像のセグメンテーションを行ない、物質を特定して、当該物質の線吸収係数を利用して再構成を繰り返すことでアーチファクトを低減する技術が開示されている。また、特許文献２は、メタルアーチファクトが発生しているオリジナルＣＴ画像と、物質のＸ線吸収係数のエネルギー依存性を考慮した順投影計算と単一波長用画像再構成アルゴリズムを用いて逆投影計算を行なって得たＣＴ画像との差分が小さくなるように、順投影計算と逆投影計算を繰り返し行なうことでアーチファクトを低減する技術が開示されている。
特開２０１０－０６８８３２号公報 特開２０１７－２２１３３９号公報
　しかしながら、特許文献１記載の技術は、事前に試料の質量密度や質量吸収係数といった物性値を知る必要がある。また、セグメンテーションと再構成を繰り返す必要があるため、計算コストがかかる。特許文献２記載の技術は、順投影計算と逆投影計算を繰り返す必要があるため、計算コストがかかる。
　本発明は、このような事情に鑑みてなされたものであり、ＣＴ画像の再構成におけるビームハードニング効果によるアーチファクトを補正するための計算コストを低減できる補正装置、システム、方法およびプログラムを提供することを目的とする。
　（１）上記の目的を達成するため、本発明の補正装置は、ＣＴ画像の再構成におけるビームハードニング効果によるアーチファクトを補正する補正装置であって、入射Ｘ線分布を取得する入射Ｘ線分布取得部と、線吸収係数のエネルギー依存性をパラメータを含むスケール因子で表した線吸収係数モデルを取得する線吸収係数モデル取得部と、投影像を取得する投影像取得部と、前記投影像を前記入射Ｘ線分布および前記線吸収係数モデルを用いて補正する補正部と、を備えることを特徴としている。
　（２）また、本発明の補正装置において、前記線吸収係数モデルは、基準エネルギーにおける線吸収係数と前記スケール因子の積で表されることを特徴としている。
　（３）また、本発明の補正装置において、前記線吸収係数モデルの前記パラメータは、１パラメータであることを特徴としている。
　（４）また、本発明の補正装置において、前記スケール因子は、冪指数を前記パラメータとする冪関数で表されることを特徴としている。
　（５）また、本発明の補正装置において、前記線吸収係数モデルの前記パラメータは、前記入射Ｘ線分布のエネルギー範囲に基づいて決定されることを特徴としている。
　（６）また、本発明の補正装置において、前記線吸収係数モデルの前記パラメータは、代表的な元素群の線吸収係数から決定されることを特徴としている。
　（７）また、本発明の補正装置において、前記入射Ｘ線分布取得部は、単色Ｘ線の本数と強度とエネルギー値に基づいて前記入射Ｘ線分布を取得することを特徴としている。
　（８）また、本発明の補正装置は、前記補正部により補正された前記投影データに基づいて再構成を行ない、ＣＴ画像を生成する再構成部と、前記ＣＴ画像を表示装置に表示させる表示部と、を備えることを特徴としている。
　（９）また、本発明のシステムは、Ｘ線を発生させるＸ線源と、Ｘ線を検出する検出器と、前記Ｘ線源および前記検出器、または試料の回転を制御する回転制御ユニットと、を備えるＣＴ装置と、上記（１）から（８）のいずれかに記載の補正装置と、を備えることを特徴としている。
　（１０）また、本発明の方法は、ＣＴ画像の再構成におけるビームハードニング効果によるアーチファクトを補正する方法であって、入射Ｘ線分布を取得するステップと、線吸収係数のエネルギー依存性をパラメータを含むスケール因子で表した線吸収係数モデルを取得するステップと、投影像を取得するステップと、前記投影像を前記入射Ｘ線分布および前記線吸収係数モデルを用いて補正するステップと、を含むことを特徴としている。
　（１１）また、本発明のプログラムは、ＣＴ画像の再構成におけるビームハードニング効果によるアーチファクトを補正するプログラムであって、入射Ｘ線分布を取得する処理と、線吸収係数のエネルギー依存性をパラメータを含むスケール因子で表した線吸収係数モデルを取得する処理と、投影像を取得する処理と、前記投影像を前記入射Ｘ線分布および前記線吸収係数モデルを用いて補正する処理と、をコンピュータに実行させることを特徴としている。
（ａ）は、試料に透過させるＸ線の概念図、（ｂ）は、離散化した単色Ｘ線の一部とＣＴ測定によって計測される試料の計測空間（ボクセル）の関係を示す概念図である。 入射Ｘ線分布の一例を示すグラフである。 チタンの質量吸収係数を示すグラフである。 鉄の質量吸収係数を示すグラフである。 全体のシステムの構成の一例を示す概略図である。 全体のシステムの構成の変形例を示す概略図である。 処理装置および補正装置の構成の一例を示すブロック図である。 処理装置および補正装置の構成の変形例を示すブロック図である。 処理装置および補正装置の構成の変形例を示すブロック図である。 補正装置の動作の一例を示すフローチャートである。 システムの動作の一例を示すフローチャートである。 システムの動作の変形例を示すフローチャートである。 （ａ）、（ｂ）、それぞれ補正を行なっていない投影像を用いて再構成した試料１のＣＴ画像および補正を行なった投影像を用いて再構成した試料１のＣＴ画像である。 （ａ）、（ｂ）、それぞれ図１３（ａ）および（ｂ）のＣＴ画像の直線ＡＢ上におけるラインプロファイルを示すグラフである。 （ａ）、（ｂ）、それぞれ補正を行なっていない投影像を用いて再構成した試料２のＣＴ画像および補正を行なった投影像を用いて再構成した試料２のＣＴ画像である。 （ａ）、（ｂ）、それぞれ図１５（ａ）および（ｂ）のＣＴ画像の直線ＡＢ上におけるラインプロファイルを示すグラフである。
　次に、本発明の実施の形態について、図面を参照しながら説明する。説明の理解を容易にするため、各図面において同一の構成要素に対しては同一の参照番号を付し、重複する説明は省略する。
　［原理］
　ＣＴ装置は、あらゆる角度からコーン状または平行ビームのＸ線を試料に照射し、検出器によりＸ線の吸収係数の分布、すなわち投影像を取得する。あらゆる角度からＸ線を照射するために、ＣＴ装置は、固定されたＸ線源および検出器に対して、試料台を回転させるか、Ｘ線源と検出器が一体となったガントリを回転させるように構成されている。
　このようにして様々な角度から投影を行ない得られた試料の投影像の濃淡で試料の線吸収係数ｆの分布を推測できる。そして、２次元的な投影像から３次元的な線吸収係数分布を求めることを再構成という。再構成は、基本的には投影像の逆投影を行なう。
　ＣＴ装置では入射Ｘ線として連続Ｘ線を使用して、これを試料に照射して投影像を測定している。物質の線吸収係数にはエネルギー依存性があり、高エネルギーのＸ線の方が吸収されにくい傾向にある。そのため、入射Ｘ線として連続Ｘ線が試料を通過すると、通過後のエネルギー分布は元々の分布と相似ではなくなり、高エネルギー側に重心がシフトする。この現象をビームハードニング（線質硬化）という。
　再構成は通常、単色Ｘ線を仮定して行なう。一方、実際の測定では、広いエネルギー分布を持つ連続Ｘ線が使用される。連続Ｘ線が物質を透過する際、低いエネルギーのＸ線は高いエネルギーに比べて、試料に吸収されて減弱しやすい。そのため、Ｘ線の透過距離とＸ線の減衰との線形性が失われ、非線形性が現れる。単色Ｘ線が物質を透過する際は、線形性が保たれる。このように、再構成の仮定と実際の現象が一致しないため、アーチファクトが生じる。このようなアーチファクトを、ビームハードニング効果によるアーチファクトという。
　従来、ビームハードニング効果によるアーチファクトを抑制するために、ハードウェアまたはソフトウェアによる補正が行なわれてきた。完全な単色Ｘ線ではビームハードニング効果によるアーチファクトは起きないため、ハードウェアによる補正方法では、基本的に入射Ｘ線分布の幅を狭めることで、ビームハードニング効果によるアーチファクトを目立ちにくくしている。例えば、フィルターを使って特定のエネルギー付近のＸ線だけを取り出したり、ミラーを使って単色Ｘ線を取り出したりすることが行なわれている。
　また、ソフトウェアによる補正では、画像処理を行なって補正する。例えば、Ｈｅｌｇａｓｏｎ－Ｌｕｄｗｉｇ条件の方法は、線吸収係数の保存則など物理的な条件を用いて投影像の強度を非線形に補正することで、アーチファクトを低減している。この方法では投影像のデータに整合性があるとの仮定に基づいて補正をしているため、例えば、測定時にＦＯＶから試料がはみ出している場合（インテリアＣＴ）、投影像に欠陥がある場合、センターシフトなどの光学系にズレがある場合等、ビームハードニング効果以外の現象が起きている場合には適用できない。また、補正量の計算は投影像の積分を用いているため、情報が圧縮されていることに起因して、Ｘ線強度の大小関係が逆転したり画素値が負になったりするなどして、補正がうまくいかないこともある。
　また、試料の既知条件をもとに再構成する方法もある。例えば、試料には一様な部分が多いはずであり、そしてアーチファクトにより非一様な部分が生じているはずであるという仮定の下、再構成像の全変動（ＴＶ）が小さくなるように逐次的に再構成を行なう方法などである。このような方法は、事前に試料に含まれる物質は何であるかを知る必要がある。含まれる物質が限られる医療や歯科用のＣＴではこのような方法がよく用いられ、効果を発揮しているが、何が含まれるか分からない試料に対しては、この方法を用いることができない。また、ＴＶ正則化を用いた逐次法では計算時間がかかってしまい、実用上の問題が生じる。
　また、入射Ｘ線分布や吸収係数のエネルギー依存性を考慮し、逐次法の中で連続Ｘ線を仮定した投影演算を行なって補正する方法もある。エネルギー情報を用いた逐次法では必然的に計算時間がかかってしまい、実用上の問題が生じる。
　本発明は、投影像を入射Ｘ線分布および線吸収係数のエネルギー依存性をパラメータ（補助変数）を含むスケール因子で表した線吸収係数モデルを用いて補正することで、補正を各投影角度における各検出器ピクセルに対して行なうことができる。これにより、セグメンテーションを必要とせず、試料の厳密な情報を用いることがなくなる。また、再構成前の処理で行なうために、ＣＴ装置の再構成におけるビームハードニング効果によるアーチファクトを補正するための計算コストを低減することができる。さらに、ビームハードニング効果以外の現象が起きている場合にも適用できる。なお、線吸収係数モデルは、線吸収係数の分布のエネルギー依存性を近似するように選択された関数である。このモデルは、例えば、スケール因子にあるエネルギーＥ０での線吸収係数の分布を掛けることで表される。ここで、スケール因子とは、ある基準エネルギーＥ０での値が１となる、エネルギーＥに関する非負関数ｓ（Ｅ）である。
　以下で、本発明の補正方法を詳細に説明する。本発明では、まず、連続Ｘ線を有限個の単色Ｘ線の集まりであると仮定する。図１（ａ）、（ｂ）は、それぞれ試料に透過させるＸ線の概念図、および離散化した単色Ｘ線の一部とＣＴ測定によって計測される試料の計測空間（ボクセル）の関係を示す概念図である。一般的に、Ｘ線の透過距離とＸ線の減衰との線形性に関しては、検出器で検出されるＸ線強度Ｉは、対象物の厚さをｌ、線吸収係数をμ、入射Ｘ線強度をＩ０とすると、式（１）で表される（Ｌａｍｂｅｒｔ-Ｂｅｅｒ
の法則）。
Figure JPOXMLDOC01-appb-M000001
　連続Ｘ線をＮ本の有限個の単色Ｘ線の集まりであると仮定すると、入射Ｘ線強度Ｉ０は、各単色Ｘ線の強度Ｉｋ（ｋ＝１，２，…，Ｎ）を積算した強度に置き換えられる。また、Ｘ線の透過距離とＸ線の減衰の非線形性は、各単色Ｘ線の減衰を全エネルギーで足し合わせることで表される。各検出器ピクセルで検出される強度Ｉは、物質によって減衰された各単色Ｘ線の強度の総和として以下の式（２）で表される。
Figure JPOXMLDOC01-appb-M000002
　連続Ｘ線のエネルギー分布は、Ｎ本の各単色Ｘ線の強度Ｉｋと各単色Ｘ線のエネルギーＥｋによって、図２のような分布として表すことができる。図２は、入射Ｘ線分布の一例を示すグラフである。例えば、図２のような入射Ｘ線分布が想定される場合、連続Ｘ線は、Ｎ＝１０の単色Ｘ線から成るエネルギー分布として表されるとして、入射Ｘ線分布のエネルギー範囲とＮ本の単色Ｘ線のエネルギー値Ｅｋ（ｋ＝１，２，…，Ｎ）を設定する。例えば、入射Ｘ線分布のエネルギー範囲を１０ｋｅＶ～５５ｋｅＶであるとして、エネルギー値と本数を、１０、１５、２０、…、５５ｋｅＶといったように等間隔に１０個設定してもよい。また、分布の形状から、変化が大きいところは密にとり、強度が低いところはスパースにとるといったように単色Ｘ線のエネルギー値と本数の決定をしてもよい。
　単色Ｘ線のエネルギー値と本数と強度を設定することは、式（２）のIｋを設定することに対応する。設定する単色Ｘ線の本数は、３以上である必要があり、５以上であることが好ましい。入射Ｘ線分布は、管電圧、フィルターの情報からある程度は分かるため、そのような知見を利用してもよい。また、ユーザが任意に選択または指定してもよい。
　本発明では、式（２）の右辺のｆ（Ｅｋ）のエネルギー依存性をパラメータを含むスケール因子で表した線吸収係数モデルに置き換える。線吸収係数のエネルギー依存性は、パラメータを含むスケール因子で表す。スケール因子ｓ（Ｅｋ）を算出するためのエネルギー範囲（ｓ（Ｅｋ）の定義域）は、入射Ｘ線分布を設定したエネルギー範囲を含むように設定される。例えば、ｆ（Ｅｋ）は以下の式（３）のように、スケール因子の定義域に含まれる基準となるあるエネルギーＥ０に対する線吸収係数ｆ（Ｅ０）とパラメータを含むスケール因子ｓ（Ｅｋ）の積として表すことができる。このような線吸収係数モデルを使用して補正をすることで、基準エネルギーにおける線吸収係数の分布の投影像を得ることができる。
Figure JPOXMLDOC01-appb-M000003
　式（１）の右辺のｆ（Ｅｋ）を式（３）の右辺の式で置き換えると、式（２）は、以下の式（４）となる。
Figure JPOXMLDOC01-appb-M000004
　式（２）は、全エネルギーの投影像の線積分を計算している。よって、式（２）は、それぞれのエネルギーに対して線積分の値を計算する必要がある。これに対し、式（４）は、あるエネルギーＥ０の投影像の線積分を計算している。よって、式（４）は、スケール因子によって、基準エネルギーの吸収係数から他のエネルギーの情報がわかる。未知数が減るため、計算が容易になる。また、スケール因子のパラメータは、１パラメータにすることが好ましい。これにより、計算がさらに簡単になる。ただし、正確に算出する場合は、複数パラメータを設定してもよい。
　この方程式をＮｅｗｔｏｎ法などで解くと、あるエネルギーＥ０の線吸収係数の分布の投影像を得ることができる。この線積分の値はエネルギーＥ０の線吸収係数の分布の投影に対応する。これを補正された投影像とする。つまり、連続Ｘ線の投影像を、あるエネルギーＥ０の投影像に変換していることになる。このような、補正された投影像を用いて再構成をすると、ビームハードニング効果によるアーチファクトが低減されたＣＴ画像を得ることができる。
　基準となるエネルギーＥ０は、ｓ（Ｅｋ）の定義域の中であれば、どのように選択してもよい。例えば、入射Ｘ線分布のエネルギー範囲の下限値としてもよいし、入射Ｘ線分布のエネルギー範囲の平均値としてもよい。
　次に、線吸収係数モデルを表す関数形を設定する。線吸収係数モデルのスケール因子は、冪指数をパラメータとする冪関数で表されることが好ましい。これにより、補正のための演算がさらに容易になる。線吸収係数モデルのスケール因子が冪関数で表されるとは、例えば、式（３）のｓ（Ｅｋ）が、以下の式（５）のように表されることをいう。ただし、線吸収係数モデルに含まれるスケール因子は、これに限られない。線吸収係数のエネルギー依存性を近似することができる関数形であれば、そのような関数を使用してもよい。例えば、指数関数、対数関数等を使用してもよい。線吸収係数モデルに含まれるスケール因子は、線吸収係数のエネルギー依存性の近似の度合いおよびその後の計算コストの両方の観点から選択されることが好ましい。
Figure JPOXMLDOC01-appb-M000005
　これにより、線吸収係数モデルは、以下の式（６）のように表される。また、線吸収係数モデルを用いて、各検出器ピクセルで検出される強度Ｉを表すと式（７）のように表される。
Figure JPOXMLDOC01-appb-M000006
Figure JPOXMLDOC01-appb-M000007
　また、線吸収係数モデルは、入射Ｘ線分布のエネルギー範囲に基づいて決定されることが好ましい。これにより、代表的な元素群の線吸収係数に応じた適切な線吸収係数モデルを決定でき、ビームハードニング効果によるアーチファクトをより適切に補正できる。線吸収係数モデルが、入射Ｘ線分布のエネルギー範囲に基づいて決定されるとは、入射Ｘ線分布のエネルギー範囲に応じて、式（４）のｓ（Ｅｋ）およびＥ０を変更することができることを意味する。線吸収係数モデルを使用する場合、仮定された入射Ｘ線分布における、代表的な元素の線吸収係数の振る舞いから、各単色Ｘ線の線吸収係数を算出することができる。例えば、入射Ｘ線分布が一般的なＣＴ装置で使用するエネルギー範囲（例えば、１０～１００ｋｅＶ）のとき、メタルアーチファクトの原因として考えられる元素（Ｔｉ、Ｆｅ等）の線吸収係数のエネルギー依存性がわかると、基準のエネルギーＥ０と各単色Ｘ線のエネルギーＥｋにおけるスケール因子の比率から、基準のエネルギーＥ０の線吸収係数ｆ(Ｅ０)を基に各単色Ｘ線の線吸収係数ｆ(Ｅｋ)を算出することができる。このとき、同一の関数でパラメータの値のみ変更してもよい。線吸収係数モデルが冪関数で表される場合、以下のようにパラメータαの決定および変更をすることが好ましい。
　式（６）および式（７）のパラメータαは、補正に適用する元素の質量吸収係数の減衰の仕方が、試料に含まれる元素と類似するように設定することが好ましい。質量吸収係数は、傾きが元素ごとに異なる、吸収端の位置が異なる、吸収端の前後で傾きが異なる、といった特徴があるため、以下のような方法でαを設定することが好ましい。
（i）代表的な元素を選択する。例えば、Ｔｉ（α=－２．６６７）、Ｆe（α＝－２．７０７）等の元素を代表的な元素として選択することができる。図３および図４は、それぞれ、ＴｉおよびＦｅの質量吸収係数を示すグラフである。代表的な元素は、試料に含まれる元素であることが好ましい。また、吸収端が使用するエネルギーの範囲に含まれない元素を選ぶことが好ましい。
（ii）選択する代表的な元素は、軽元素群、中元素群、重元素群の中から選ぶようにすることが好ましい。例えば、軽元素群はＨからＮｅ、中元素群はＮａからＣａ、重元素群はＳｃ以降とすることができる。軽元素群、中元素群、重元素群の境界は分野によって変更することができる。また、線吸収係数のエネルギー依存性が似ている元素を特定しておくことが好ましい。パラメータαは、冪関数のフィッティングにより算出することが好ましいが、必ずしもフィッティング結果と一致しなくてもよく、フィッティング結果を目安に任意で設定してもよい。
（iii）使用するエネルギー範囲に吸収端が含まれる元素を選んだ場合、吸収端の前の傾きαＲ、吸収端の後の傾きαＬから、平均の傾きを算出してαとして使用することが好ましい。
　上記のように、線吸収係数モデルを冪因子で表し、パラメータαの値を所定の範囲内で適切に設定することで、補正のための計算がさらに容易になり、ＣＴ装置の再構成におけるビームハードニング効果によるアーチファクトを効果的に補正することができる。また、そのための計算コストを大幅に低減することができる。
　［全体のシステム］
　図５は、ＣＴ装置２００とこれに接続された処理装置３００、補正装置４００、入力装置５１０および表示装置５２０を含む全体のシステム１００の構成を示す概略図である。ここで、図５に示すＣＴ装置２００は、Ｘ線源２６０および検出器２７０に対し試料を回転させる構成であるが、これに限定されることはなく、Ｘ線源と検出器が一体となったガントリを回転させる構成でもよい。また、ＣＴ装置２００は、平行ビームを使用する装置もコーンビームを使用する装置もいずれも使用できる。
　処理装置３００は、ＣＴ装置２００に接続され、ＣＴ装置２００の制御および取得されたデータの処理を行なう。補正装置４００は、投影像の補正を行なう。処理装置３００および補正装置４００は、ＰＣ端末であってもよいし、クラウド上のサーバであってもよい。入力装置５１０は、例えばキーボード、マウスであり、処理装置３００や補正装置４００への入力を行なう。表示装置５２０は、例えばディスプレイであり、投影像などを表示する。
　なお、図５では、補正装置４００の補正機能を強調するために、処理装置３００と補正装置４００を別の構成要素として表示しているが、図６のように、補正装置４００は処理装置３００に含まれる一部の機能として構成されてもよいし、補正装置４００と処理装置３００は一体的なものとして構成されてもよい。図６は、全体のシステムの構成の変形例を示す概略図である。このようなシステムを使用することにより、ＣＴ画像の再構成におけるビームハードニング効果によるアーチファクトを補正するための計算コストを低減することができる。
　［ＣＴ装置］
　図５に示すように、ＣＴ装置２００は、回転制御ユニット２１０、試料台２５０、Ｘ線源２６０、検出器２７０および駆動部２８０を備えている。Ｘ線源２６０と検出器２７０の間に設置された、試料台２５０を回転させてＸ線ＣＴ撮影を行なう。なお、Ｘ線源２６０および検出器２７０は、ガントリ（図示しない）に設置し、試料台２５０に固定された試料に対しガントリを回転させてもよい。
　ＣＴ装置２００は、処理装置３００により指示されたタイミングで試料台２５０を駆動し、試料の投影像を取得する。測定データは、処理装置３００に送信される。ＣＴ装置２００は、半導体デバイス等の精密な工業製品に用いることに適しているが、産業用装置のみならず動物用装置にも適用できる。
　Ｘ線源２６０は、Ｘ線を検出器２７０に向けて照射する。検出器２７０は、Ｘ線を受ける受光面を有し、多数のピクセルにより試料を透過したＸ線の強度分布を測定できる。回転制御ユニット２１０は、駆動部２８０によりＣＴ撮影時に設定された速度で試料台２５０を回転させる。
　［処理装置］
　図７は、処理装置３００および補正装置４００の構成を示すブロック図である。処理装置３００は、ＣＰＵ（Central Processing Unit／中央演算処理装置）、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）、メモリをバスに接続してなるコンピュータによって構成されている。処理装置３００は、ＣＴ装置２００に接続され情報を受け取る。
　処理装置３００は、測定データ記憶部３１０、装置情報記憶部３２０、再構成部３３０、および表示部３４０を備える。各部は、制御バスＬにより情報を送受できる。入力装置５１０および表示装置５２０は適宜のインターフェースを介してＣＰＵに接続されている。
　測定データ記憶部３１０は、ＣＴ装置２００から取得した測定データを記憶する。測定データには、回転角度情報とそれに対応する投影像が含まれる。装置情報記憶部３２０は、ＣＴ装置２００から取得した装置情報を記憶する。装置情報には、装置名、ビーム形状、測定時のジオメトリ、スキャン方式等が含まれる。再構成部３３０は、対象となる投影像からＣＴ画像を再構成する。表示部３４０は、再構成したＣＴ画像を表示装置５２０に表示させる。これにより、補正された投影像に基づいたＣＴ画像をユーザが確認することができる。また、ユーザがＣＴ画像に基づいて処理装置、補正装置等に指示、指定をすることができる。
　［補正装置］
　補正装置４００は、ＣＰＵ、ＲＯＭ、ＲＡＭ、メモリをバスに接続してなるコンピュータによって構成されている。補正装置４００は、ＣＴ装置２００に直接接続されてもよいし、処理装置３００を介してＣＴ装置２００に接続されてもよい。また、補正装置４００は、ＣＴ装置２００から情報を受け取ってもよいし、処理装置３００から情報を受け取ってもよい。なお、図８のように、補正装置４００が処理装置３００に含まれる一部の機能として構成されてもよいし、図９のように、補正装置４００と処理装置３００は一体的なものとして構成されてもよい。また、補正装置４００は、処理装置３００の機能の一部を備えていてもよい。
　補正装置４００は、入射Ｘ線分布取得部４１０、線吸収係数モデル取得部４２０、投影像取得部４３０、および補正部４４０を備える。各部は、制御バスＬにより情報を送受できる。補正装置４００と処理装置３００が別の構成である場合、入力装置５１０および表示装置５２０は適宜のインターフェースを介して補正装置４００のＣＰＵにも接続されている。この場合、入力装置５１０および表示装置５２０は、処理装置３００に接続されるものとは異なっていてもよい。
　入射Ｘ線分布取得部４１０は、入射Ｘ線分布を取得する。入射Ｘ線分布取得部４１０は、単色Ｘ線の本数と強度とエネルギー値に基づいて入射Ｘ線分布を取得することが好ましい。これにより、状況に応じたより適切な入射Ｘ線分布を取得することができる。また、入射Ｘ線分布取得部４１０は、ユーザ指定に基づいて入射Ｘ線分布を取得することが好ましい。これにより、状況に応じたより適切な入射Ｘ線分布をユーザが選択することができる。入射Ｘ線分布がユーザにより指定される場合、例えば、マウス操作で予め記憶されている複数の入射Ｘ線分布から１つを指定する、入射Ｘ線分布のグラフの上の点を移動させる、グラフ上の点の数を増減させる等、入射Ｘ線分布を設定できるＵＩ機能を用いることが好ましい。入射Ｘ線分布は、予め設定されていてもよい。また、入射Ｘ線分布は、装置情報等に応じて自動的に選択され、取得されてもよい。
　線吸収係数モデル取得部４２０は、線吸収係数のエネルギー依存性をパラメータを含むスケール因子で表した線吸収係数モデルを取得する。線吸収係数モデル取得部４２０は、ユーザ指定に基づいて線吸収係数モデルを取得することが好ましい。補正装置４００または処理装置３００の図示しない記憶部に線吸収係数のエネルギー依存性を表す関数形を記憶させておく。例えば、線吸収係数ｆ（Ｅｋ）のモデルを表す式、あるいは基準エネルギー（Ｅ０）とスケール因子ｓ（Ｅｋ）を表す式として記憶させておく。また、試料に含まれうる元素群と線吸収係数モデルのパラメータαを相関づけて記憶させておく。線吸収係数モデルがユーザにより指定される場合、例えば、マウス操作で予め記憶されている複数の線吸収係数を表す関数形から１つを指定し、元素群を選択することによってパラメータαを設定する。また、パラメータαは、ユーザによって直接値を入力できるようにしてもよい。このように、線吸収係数モデルを設定できるＵＩ機能を用いることが好ましい。線吸収係数モデルは、予め設定されていてもよい。また、線吸収係数モデルは、試料情報等に応じて自動的に選択され、取得されてもよい。
　投影像取得部４３０は、ＣＴ装置２００または処理装置３００から投影像を取得する。補正部４４０は、投影像を入射Ｘ線分布および線吸収係数モデルを用いて補正する。これにより、ビームハードニング効果によるアーチファクトを低減したＣＴ画像を得ることができる。
　［測定方法］
　ＣＴ装置２００に試料を設置し、所定の条件で回転軸の移動とＸ線の投影を繰り返すことで、Ｘ線を試料に照射しつつ投影像を取得する。ＣＴ装置２００は、スキャン方式のような装置情報および取得された投影像を測定データとして処理装置３００または補正装置４００に送信する。
　［補正方法］
　図１０は、補正装置４００の動作の一例を示すフローチャートである。まず、補正装置４００は、入射Ｘ線分布を取得する（ステップＳ１）。次に、線吸収係数モデルを取得する（ステップＳ２）。次に、投影像を取得する（ステップＳ３）。そして、取得した入射Ｘ線分布および線吸収係数モデルを用いて投影像を補正する（ステップＳ４）。このようにして、投影像を入射Ｘ線分布および線吸収係数のエネルギー依存性をパラメータを含むスケール因子で表した線吸収係数モデルを用いて補正することができる。なお、入射Ｘ線分布の取得、線吸収係数モデルの取得、および投影像の取得は、どのような順序で行なっても問題ない。
　［補正および再構成方法］
　図１０のフローチャートは、補正装置４００のみの動作を示している。投影像の補正としてはこれで十分であるが、従来技術との違いを明確にするため、投影像の測定および再構成を含む動作も説明しておく。
　図１１は、システム１００の動作の一例を示すフローチャートである。まず、補正装置４００は、入射Ｘ線分布を取得する（ステップＴ１）。次に、線吸収係数モデルを取得する（ステップＴ２）。次に、ＣＴ装置２００は、投影像を測定する（ステップＴ３）。投影像の測定は、回転軸の移動とＸ線の投影の繰り返しを含む。次に、補正装置４００は、投影像を取得する（ステップＴ４）。次に、取得した入射Ｘ線分布および線吸収係数モデルを用いて投影像を補正する（ステップＴ５）。そして、処理装置３００または補正装置４００は、補正された投影像を用いて、ＣＴ画像を再構成する（ステップＴ６）。このようにして、補正した投影像を用いて再構成したＣＴ画像を得ることができる。なお、上記と同様に、入射Ｘ線分布の取得、線吸収係数モデルの取得、および投影像の取得は、どのような順序で行なっても問題ない。
　従来技術では、再構成したＣＴ画像を得てから、アーチファクト低減のための画像処理を繰り返していた。これに対し、本発明では、再構成する前の投影像に対してアーチファクト低減のための補正をしている。このようにすることで、計算コストを低減することができる。
　図１２は、システム１００の動作の変形例を示すフローチャートである。まず、補正装置４００は、入射Ｘ線分布を取得する（ステップＵ１）。次に、線吸収係数モデルを取得する（ステップＵ２）。次に、ＣＴ装置２００は、回転軸の移動をする（ステップＵ３）。次に、Ｘ線の投影をする（ステップＵ４）。次に、次に、補正装置４００は、先程投影した投影像を取得する（ステップＵ５）。次に、取得した入射Ｘ線分布および線吸収係数モデルを用いて投影像を補正する（ステップＵ６）。
　次に、処理装置３００または補正装置４００は、測定が完了しているかを判断する（ステップＵ７）。ステップＵ７において、測定が完了していないと判断した場合（ステップＵ７、ＮＯ）は、ステップＵ３に戻る。一方、ステップＵ７において、処理装置３００または補正装置４００は、測定が完了したと判断した場合（ステップＵ７、ＹＥＳ）は、補正された投影像を用いて、ＣＴ画像を再構成する（ステップＵ８）。このようにすることで、投影像の測定と測定した投影像の補正を並列的に処理することができ、全体の時間を短縮することができる。なお、入射Ｘ線分布の取得、および線吸収係数モデルの取得は、逆の順序で行なっても問題ない。
　［実施例１］
　上記のように構成されたシステム１００を用いて、ハイドロオキシアパタイトファントム（ＱＲＭ社製：Micro-CT HA Phantom）（試料１）の断面を観察した。ＣＴ装置２００には、リガク製ＣＴＬａｂＨＸを用いて、６０ｋＶの管電圧で測定した。図１３（ａ）および（ｂ）は、それぞれ補正を行なっていない投影像を用いて再構成した試料１のＣＴ画像および補正を行なった投影像を用いて再構成した試料１のＣＴ画像である。また、図１４（ａ）および（ｂ）は、それぞれ図１３（ａ）および（ｂ）のＣＴ画像の直線ＡＢ上におけるラインプロファイルを示すグラフである。再構成は、ＦＤＫ法を用いた。
　図１３（ａ）および（ｂ）を比較することで、補正により、２個の高吸収体の間のストリークアーティファクトが低減していることが分かった。また、図１４（ａ）および（ｂ）を比較することで、高吸収体の端のカッピングアーティファクトが低減されていることが確認できた。
　［実施例２］
　次に、木材に径３ｍｍのアクリル棒を５本、径９ｍｍのアルミ棒を４本挿入し、径３ｍｍの中空の穴を２箇所、径１．５ｍｍの中空の穴を２箇所形成した試料２を作製した。上記システム１００を用いて、試料２の断面を観察した。管電圧は１００ｋＶとした。図１５（ａ）および（ｂ）は、それぞれ補正を行なっていない投影像を用いて再構成した試料２のＣＴ画像および補正を行なった投影像を用いて再構成した試料２のＣＴ画像である。また、図１６（ａ）および（ｂ）は、それぞれ図１５（ａ）および（ｂ）のＣＴ画像の直線ＡＢ上におけるラインプロファイルを示すグラフである。
　図１５（ａ）および（ｂ）を比較することで、補正により、ストリークアーティファクトおよびダークバンドが低減しているのが分かった。また、図１６（ａ）および（ｂ）を比較することで、アルミ棒の端のカッピングアーティファクトの低減が確認できた。
　以上の結果により、本発明の補正装置、システム、方法およびプログラムは、ＣＴ画像の再構成におけるビームハードニング効果によるアーチファクトを効果的に補正でき、計算コストを低減できることが確認された。
　なお、本出願は、２０２１年７月２７日に出願した日本国特許出願第２０２１－１２２５８５号に基づく優先権を主張するものであり、日本国特許出願第２０２１－１２２５８５号の全内容を本出願に援用する。
１００　システム
２００　ＣＴ装置
２１０　回転制御ユニット
２５０　試料台
２６０　Ｘ線源
２７０　検出器
２８０　駆動部
３００　処理装置
３１０　測定データ記憶部
３２０　装置情報記憶部
３３０　再構成部
３４０　表示部
４００　補正装置
４１０　入射Ｘ線分布取得部
４２０　線吸収係数モデル取得部
４３０　投影像取得部
４４０　補正部
５１０　入力装置
５２０　表示装置

Claims (11)

  1. 　ＣＴ画像の再構成におけるビームハードニング効果によるアーチファクトを補正する補正装置であって、
    　入射Ｘ線分布を取得する入射Ｘ線分布取得部と、
    　線吸収係数のエネルギー依存性をパラメータを含むスケール因子で表した線吸収係数モデルを取得する線吸収係数モデル取得部と、
    　投影像を取得する投影像取得部と、
    　前記投影像を前記入射Ｘ線分布および前記線吸収係数モデルを用いて補正する補正部と、を備えることを特徴とする補正装置。
  2. 　前記線吸収係数モデルは、基準エネルギーにおける線吸収係数と前記スケール因子の積で表されることを特徴とする請求項１記載の補正装置。
  3. 　前記線吸収係数モデルの前記パラメータは、１パラメータであることを特徴とする請求項１または請求項２記載の補正装置。
  4. 　前記スケール因子は、冪指数を前記パラメータとする冪関数で表されることを特徴とする請求項１から請求項３のいずれかに記載の補正装置。
  5. 　前記線吸収係数モデルの前記パラメータは、前記入射Ｘ線分布のエネルギー範囲に基づいて決定されることを特徴とする請求項１から請求項４のいずれかに記載の補正装置。
  6. 　前記線吸収係数モデルの前記パラメータは、代表的な元素群の線吸収係数から決定されることを特徴とする請求項１から請求項５のいずれかに記載の補正装置。
  7. 　前記入射Ｘ線分布取得部は、単色Ｘ線の本数と強度とエネルギー値に基づいて前記入射Ｘ線分布を取得することを特徴とする請求項１から請求項６のいずれかに記載の補正装置。
  8. 　前記補正部により補正された前記投影像に基づいて再構成を行ない、ＣＴ画像を生成する再構成部と、
    　前記ＣＴ画像を表示装置に表示させる表示部と、を備えることを特徴とする請求項１から請求項７のいずれかに記載の補正装置。
  9. 　Ｘ線を発生させるＸ線源と、Ｘ線を検出する検出器と、前記Ｘ線源および前記検出器、または試料の回転を制御する回転制御ユニットと、を備えるＣＴ装置と、
    　請求項１から請求項８のいずれかに記載の補正装置と、を備えることを特徴とするシステム。
  10. 　ＣＴ画像の再構成におけるビームハードニング効果によるアーチファクトを補正する方法であって、
    　入射Ｘ線分布を取得するステップと、
    　線吸収係数のエネルギー依存性をパラメータを含むスケール因子で表した線吸収係数モデルを取得するステップと、
    　投影像を取得するステップと、
    　前記投影像を前記入射Ｘ線分布および前記線吸収係数モデルを用いて補正するステップと、を含むことを特徴とする方法。
  11. 　ＣＴ画像の再構成におけるビームハードニング効果によるアーチファクトを補正するプログラムであって、
    　入射Ｘ線分布を取得する処理と、
    　線吸収係数のエネルギー依存性をパラメータを含むスケール因子で表した線吸収係数モデルを取得する処理と、
    　投影像を取得する処理と、
    　前記投影像を前記入射Ｘ線分布および前記線吸収係数モデルを用いて補正する処理と、をコンピュータに実行させることを特徴とするプログラム。
PCT/JP2022/014455 2021-07-27 2022-03-25 補正装置、システム、方法およびプログラム WO2023007846A1 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
JP2023538266A JPWO2023007846A1 (ja) 2021-07-27 2022-03-25
DE112022003708.3T DE112022003708T5 (de) 2021-07-27 2022-03-25 Korrekturvorrichtung, System, Verfahren und Programm
CN202280052823.1A CN117716231A (zh) 2021-07-27 2022-03-25 校正装置、系统、方法以及程序

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2021-122585 2021-07-27
JP2021122585 2021-07-27

Publications (1)

Publication Number Publication Date
WO2023007846A1 true WO2023007846A1 (ja) 2023-02-02

Family

ID=85087845

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2022/014455 WO2023007846A1 (ja) 2021-07-27 2022-03-25 補正装置、システム、方法およびプログラム

Country Status (4)

Country Link
JP (1) JPWO2023007846A1 (ja)
CN (1) CN117716231A (ja)
DE (1) DE112022003708T5 (ja)
WO (1) WO2023007846A1 (ja)

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2011203160A (ja) * 2010-03-26 2011-10-13 Tokyo Institute Of Technology Ｘ線ｃｔ画像再構成方法及びｘ線ｃｔ画像再構成プログラム
WO2014192889A1 (ja) * 2013-05-29 2014-12-04 地方独立行政法人東京都立産業技術研究センター Ｘ線エネルギー別画像再構成装置及び方法並びにｘ線三次元測定装置及び方法
US20150212015A1 (en) * 2014-01-24 2015-07-30 Institute For Basic Science Apparatus and method for computed tomography image processing
WO2016002034A1 (ja) * 2014-07-03 2016-01-07 株式会社ニコン Ｘ線装置、画像形成方法、構造物の製造方法、及び構造物製造システム
JP2017221339A (ja) * 2016-06-14 2017-12-21 国立大学法人信州大学 Ｘ線ｃｔ画像再構成方法およびコンピュータプログラム
JP2018515160A (ja) * 2015-03-18 2018-06-14 プリズマティック、センサーズ、アクチボラグＰｒｉｓｍａｔｉｃ Ｓｅｎｓｏｒｓ Ａｂ フォトンカウンティングマルチビン検出器からのエネルギー分解画像データに基づく画像再構成

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP5243160B2 (ja) 2008-09-16 2013-07-24 株式会社日立メディコ Ｘ線ｃｔ装置
JP7265171B2 (ja) 2020-02-07 2023-04-26 サミー株式会社 遊技機

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2011203160A (ja) * 2010-03-26 2011-10-13 Tokyo Institute Of Technology Ｘ線ｃｔ画像再構成方法及びｘ線ｃｔ画像再構成プログラム
WO2014192889A1 (ja) * 2013-05-29 2014-12-04 地方独立行政法人東京都立産業技術研究センター Ｘ線エネルギー別画像再構成装置及び方法並びにｘ線三次元測定装置及び方法
US20150212015A1 (en) * 2014-01-24 2015-07-30 Institute For Basic Science Apparatus and method for computed tomography image processing
WO2016002034A1 (ja) * 2014-07-03 2016-01-07 株式会社ニコン Ｘ線装置、画像形成方法、構造物の製造方法、及び構造物製造システム
JP2018515160A (ja) * 2015-03-18 2018-06-14 プリズマティック、センサーズ、アクチボラグＰｒｉｓｍａｔｉｃ Ｓｅｎｓｏｒｓ Ａｂ フォトンカウンティングマルチビン検出器からのエネルギー分解画像データに基づく画像再構成
JP2017221339A (ja) * 2016-06-14 2017-12-21 国立大学法人信州大学 Ｘ線ｃｔ画像再構成方法およびコンピュータプログラム

Also Published As

Publication number Publication date
JPWO2023007846A1 (ja) 2023-02-02
CN117716231A (zh) 2024-03-15
DE112022003708T5 (de) 2024-05-16

Similar Documents

Publication Publication Date Title
US10274439B2 (en) System and method for spectral x-ray imaging
Van Gompel et al. Iterative correction of beam hardening artifacts in CT
Schmidt et al. A spectral CT method to directly estimate basis material maps from experimental photon-counting data
JP5384521B2 (ja) 放射線撮像装置
Brabant et al. A novel beam hardening correction method requiring no prior knowledge, incorporated in an iterative reconstruction algorithm
Rinkel et al. A new method for x-ray scatter correction: first assessment on a cone-beam CT experimental setup
JP5264268B2 (ja) 単位面積質量画像の作成方法
Van de Casteele et al. A model-based correction method for beam hardening artefacts in X-ray microtomography
US20110168878A1 (en) Method and apparatus for empirical determination of a correction function for correcting beam hardening and stray beam effects in projection radiography and computed tomography
AU2019202047B2 (en) X-ray beam-hardening correction in tomographic reconstruction using the Alvarez-Macovski attenuation model
Schörner et al. Comparison between beam-stop and beam-hole array scatter correction techniques for industrial X-ray cone-beam CT
EP2652709B1 (en) Imaging system for imaging a region of interest
US9589373B2 (en) Monte carlo modeling of field angle-dependent spectra for radiographic imaging systems
Ingacheva et al. Polychromatic CT Data Improvement with One‐Parameter Power Correction
JP2017221339A (ja) Ｘ線ｃｔ画像再構成方法およびコンピュータプログラム
Grimmer et al. Empirical binary tomography calibration (EBTC) for the precorrection of beam hardening and scatter for flat panel CT
WO2023007846A1 (ja) 補正装置、システム、方法およびプログラム
JP2004357969A (ja) Ｘ線計測装置
EP2584532A1 (en) Empirical cupping correction for CT scanners with primary modulation
Mohapatra Development and quantitative assessment of a beam hardening correction model for preclinical micro-CT
Ingacheva et al. Methods of preprocessing tomographic images taking into account the thermal instability of the X-ray tube
Krumm et al. Referenceless beam hardening correction in 3d computed tomography images of multi-material objects
US20240202994A1 (en) Correction apparatus, method and program for correcting artifacts
Jakab et al. High quality cone-beam CT reconstruction on the GPU
Krumm et al. Beam hardening correction of multi-material objects

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 22848930

Country of ref document: EP

Kind code of ref document: A1

WWE Wipo information: entry into national phase

Ref document number: 2023538266

Country of ref document: JP

WWE Wipo information: entry into national phase

Ref document number: 18292823

Country of ref document: US

WWE Wipo information: entry into national phase

Ref document number: 202280052823.1

Country of ref document: CN