JP6506842B2 - Correction of beam hardening artefacts of samples imaged in a container in microtomography - Google Patents

Description

The present invention relates to a method of correcting beam hardening effects in microtomography, and more particularly to a non-referenced multi-material beam hardening correction that enables automated correction.
This application claims the benefit of US Provisional Patent Application No. 62 / 038,263, entitled “Tomographic Reconstruction for Material Characterization,” filed Aug. 16, 2014.

  Tomography or computed tomography (CT) is the process of forming a three-dimensional model of an object by combining the 2D projections of the object of interest, the 2D projection of the object of interest usually being of any type Obtained by using transparent particles or transparent waves. Tomography is a rapidly advancing imaging technology and has a wide range of applications in various fields such as, but not limited to, medicine, dentistry, biology, environment, toxicology, mineralogy and electronics. Tomographic processes use various tools such as X-ray systems, transmission electron microscopy (TEM), scanning transmission electron microscopy (STEM) and / or atomic probe microscopy (APM), for example atomic structure and chemical analysis of samples Get various types of information such as Usually, reconstruct the volume by backprojecting a series of 2D images acquired at different angles across the sample, or in the case of APM, from the sequence of field evaporated atoms that hit the position sensitive detector To obtain a 3D tomographic data set.

  During the past 30 years, computer-coupled tomography (CT) has been applied as a nondestructive imaging technique for the study of geological samples such as fossils, reservoirs and soils. The number of possible applications is increasing with increasing resolution, and computer-assisted microtomography (μCT) imaging techniques can now reach resolutions of less than 1 micron. One technique that is increasingly being applied is multi-state imaging. In multi-state imaging, the same sample is imaged under different conditions, such as different pressures, different saturation states or different weathering effects. When performing such an analysis, a cylindrical vessel is required to maintain pressure or to hold saturated fluid. In high magnification, high cone angle imaging, such as that performed by the inventor of the present invention at the μCT facility of the Australian National University (ANU), the sample is placed very close to the x-ray source. In this case, a container is used to dissipate the radiant heat from the source while at the same time preventing movement of the sample. Thus, there is an increasing use of containers in μCT imaging. With this increase in use, containers present new challenges and opportunities, particularly with regard to the correction of what are called "beam hardening" effects.

  These detrimental effects are due to the interaction of the scanned material with the scanned material that can not be properly considered by the particular CT process used for a particular procedure. The standard linearization pre-processing step performed prior to tomographic reconstruction is based on Lambert-Beer's law, which assumes an exponential relationship between beam intensity and attenuation path length.

上式で、Ｉ_０は初期ビーム強度、ｘは材料中における深さ、μ（ｘ）は深さｘにおける減衰係数、Ｉ（ｔ）は厚さｔの材料を通過した後のビーム強度である。この式は、ビームに対する材料の唯一の効果が、ビームの強度を低減させることであると仮定している。すなわち、それぞれの体積要素内における減衰は、ビームが既にどれくらい材料を通過したのかということから独立していると仮定している。しかしながら、減衰係数μは、Ｘ線ビーム・エネルギーの関数でもあり、減衰は一般に、エネルギーの低下とともに低下する。そのため、（実験室ベースのマイクロフォーカスｘ資源によって生成されるものなどの）多色ビームの低エネルギー成分は優先的に吸収され、それにより、ビームの平均エネルギーは、ビームが試料を通過するにつれて増大する（よりエネルギーが高くなる）。この効果は、ビーム・ハードニングとして知られており、カッピング（ｃｕｐｐｉｎｇ）（試料の縁の近くのより高い見かけ密度）およびストリーキング（ｓｔｒｅａｋｉｎｇ）（比較的に高密度の領域間の暗影）などのアーチファクトの原因となりうる。このようなアーチファクトは画像忠実度を低減させる。このようなアーチファクトは、後に行われる定量分析に対して、特に、異なる材料であることを示すために断層像の領域にラベルを付けるセグメント化プロセスに対して、劇的な影響を与えうる。

Where I ₀ is the initial beam intensity, x is the depth in the material, μ (x) is the attenuation coefficient at depth x, I (t) is the beam intensity after passing through the material of thickness t . This equation assumes that the only effect of the material on the beam is to reduce the intensity of the beam. That is, it is assumed that the attenuation in each volume element is independent of how much the beam has already passed through the material. However, the attenuation coefficient μ is also a function of x-ray beam energy, and the attenuation generally decreases with decreasing energy. As such, low energy components of polychromatic beams (such as those produced by laboratory based microfocus x resources) are preferentially absorbed, so that the average energy of the beam increases as the beam passes through the sample Yes (more energy). This effect is known as beam hardening and is an artifact such as cupping (higher apparent density near the edge of the sample) and streaking (dark shadows between relatively dense areas) Can cause Such artifacts reduce image fidelity. Such artifacts can have a dramatic impact on later quantitative analysis, particularly on the segmentation process that labels the area of the tomogram to indicate that it is a different material.

  Generally, even when the beam spectrum and detector response are known exactly, there is insufficient information to perform tomographic reconstruction from the polychromatic attenuation data. Thus, beam hardening presents the problem of appearing to be unsolvable computationally without making assumptions about the sample. Although a number of methods have been proposed to target specific cases (such as bone and tissue imaging), there is no widely adopted method for μCT imaging. Thus, beam hardening is currently the main cause of imaging artefacts in many situations, as other artefacts of X-ray imaging have been dramatically reduced in recent years. It continues to be one of the main challenges. Provided herein is an automated beam hardening correction method that works in situations of widespread interest in μCT imaging.

  Various methods have been proposed to correct for beam hardening artifacts in μCT images. A common technique to reduce this artifact is the beam filtering technique, which places the sheet of material in front of or behind the sample and which is low energy x-rays that will be absorbed rapidly by the sample Remove This method has several disadvantages. First, this method reduces the severity of beam hardening artifacts but does not eliminate beam hardening. Furthermore, this method causes a reduction in the overall flux reaching the detector, reduces the signal to noise ratio and increases the time required to obtain the same result. When the container is present, in particular if the attenuation of the container is high as required for high pressure imaging, the flux is further reduced, which results in an additional loss of contrast. The dual energy approach of imaging the same sample at two different energies can also be employed, and the results can be combined under the assumption of linear decomposition to basis functions. However, this requires an additional μCT scan, which can increase the time and expense involved, especially if two scans are already required for dual state imaging.

  Other correction techniques are performed after data acquisition. Beam-hardening curve linearization corrects projection data by assuming that the sample consists of only one type of material. This assumption is a reasonable assumption for nearly homogeneous samples, but can not be used for heterogeneous materials and can not be used even when a container is present. A suitable group of correction methods is the so-called post-reconstruction method, which relies on reprojection of the data to estimate polychromatic or monochromatic projections. Reprojection uses forward projection of radiation to generate simulated projections from the reconstructed tomogram based on the knowledge of the x-ray spectrum and material attributes. The reprojection data and the original data are the same for monochromatic data that does not contain other artifacts. However, reconstructions that include beam hardening artifacts are generally not completely consistent with experimental data, so reprojected data is different from experimental projections. Post reconstruction uses the difference between experimental data and simulated reprojection as a correction that can be applied to that data. However, even if the x-ray spectrum and detector response are known, the material attenuation must be completely reconstructed. In addition, even for simple samples, this method does not work well for cylindrical samples, and beam-hardened projections are completely consistent including corresponding beam-hardening artifacts It is reconstructed into a tomogram. One type of post-reconstruction beam hardening correction technique is described by Krumm et al. (Krumm, M., Kasperl, S. and Franz, M., "Reducing non-linear artifacts of multi-material objects in industrial 3d computed tomography", It is proposed by NDT & E International 41 (4), 242-251 (2008)). This method is called a no reference post reconstruction correction method or an RPC method. The RPC method allows the correction to function without knowledge of the material or x-ray spectrum, but requires that the sample can be resolved into homogeneous regions.

Krumm, M .; Kasperl, S., et al. And Franz, M .; , "Reducing non-linear artifacts of multi-material objects in industrial 3d computed tomography", NDT & E International 41 (4), 242-251 (2008) "Dual-Energy Iterative Reconstruction for Material Characterization", Proc. SPIE 9212, Developments in X-Ray Tomography IX, 921213 (September 11, 2014) "Correction of Beam Hardening Artefacts in Microtomography for Samples Images in Container", Proc. SPIE 9212, Developments in X-Ray Tomography IX, 92120A (September 11, 2014) Alvarez, R .; E. And Macovski, A. et al. , "Energy-selective reconstructions in x-ray computerized tomography", Physics in medicine and biology 21 (5), 733 (1976). Siddiqui, S. And Khamees, A. et al. A. "Dual-energy ct-scanning applications in rock characterization", [SPE Annual Technical Conference and Exhibition (Society of Petroleum Engineers, 2004)]. Heismann, B .; Leppert, J. et al. And Stierstorfer, K .; , "Density and atomic number measurements with spectral x-ray attenuation method", Journal of applied physics 94 (3), 2073 to 2079 (2003). Optical Engineering and Applications. Myers, G., accepted in (2014). "Improving dynamic tomography quality through a maximum a posteriori model", SPIE Optics and Photonics

  It is an object of the present invention to provide an improved non-referenced multi-material beam hardening correction method. This method focuses on maintaining data quality for 3D imaging of geological materials, and looks to automation. This technique can be applied to other fields as well. A referenceless post-reconstruction (RPC) correction technique is provided that applies this correction in integral attenuation space. A container only pre-correction technique is also provided that allows for the automation of the segmentation process required for the beam hardening correction method. The method can use various segmentation techniques, including techniques that reference material libraries, at various stages of the process. The method is preferably implemented automatically by the tomography system under control of the controller after the operator has activated the tomography system by entering appropriate parameters. However, this is not intended to be limiting, and any particular steps may include operator intervention.

  Some embodiments of the present invention provide a method of correcting tomographic reconstruction for beam hardening. The method reconstructs an initial tomogram from the tomographic projection data and segments the tomogram into two or more phases. The two or more phases can include an agglomerate phase. Based on the segmented tomograms, the method calculates a referenceless post-reconstruction correction (RPC) correction in integrated attenuation space. The method then applies this RPC correction to tomographic projection data to correct for beam hardening effects. Using the corrected projection data, the method reconstructs the final tomogram from the corrected tomographic projection data.

  In some embodiments, applying the RPC includes reprojecting each phase of the segmented tomogram to produce a simulated monochromatic projection for each material, the simulated monochromatic projection having integrated attenuation Calculated in space. Some embodiments apply container-only pre-correction to the tomogram to correct for beam hardening effects due to the container holding the sample material, correction for beam hardening effects due to the sample itself do not do. The container-only precorrection is generated by selecting a set of radiation that passes through the container only and calculating an integral attenuation for the container based on the selected set of radiation. Some embodiments may further receive predetermined container shape data and generate container only pre-correction based on the predetermined container shape data. Container shape data used for pre-correction can also be obtained by segmenting the tomogram.

  In some embodiments, this segmentation includes performing a k-means clustering technique to distinguish phases using radial distance and grayscale intensity. In another embodiment, the segmentation comprises part of an iterative segmentation and reconstruction process for initial reconstruction, the iterative segmentation and reconstruction process comprising (a) one of a number of volume elements in a sample Determining one or more estimated material properties using a tomographic reconstruction algorithm, and (b) determining, for each volume element, the most likely material in the material dictionary Modifying the estimated material properties of the block based on the probabilistic classification by referring to the stored material property data, and (c) one or more additional iterations of the tomographic reconstruction algorithm Performing, each iteration updating one or more estimated material properties of the volume element, and (d) the sample DOO Label the reconstructed volume shows, and classifying the region of material having similar properties.

  In some embodiments, reconstructing the final tomogram can also include an iterative process, which (a) estimates one or more estimated materials of multiple volume elements in the sample. Determining the properties using a tomographic reconstruction algorithm, (b) storing one or more estimated material properties to determine the most likely material in the material dictionary for each volume element Modifying based on the probabilistic classification by referring to the selected material property data, (c) performing one or several additional iterations of the tomographic reconstruction algorithm, each iteration being a volume element Updating one or more of the estimated material properties of (d) and (d) relabeling the reconstructed volume that is indicative of a sample, the region of material having similar properties It includes classifying to.

  Some embodiments of the present invention provide a method of correcting tomographic reconstruction for beam hardening. This method reconstructs the initial tomogram from the tomographic projection data. The method then applies container-only pre-correction to the tomographic image to correct for beam hardening effects due to the container holding the sample material, correction for beam hardening effects due to the sample material itself do not do. The method then segments the tomogram into two or more phases. Based on the segmented tomograms, the method calculates a no reference post reconstruction correction (RPC) correction. The method then applies this RPC correction to tomographic projection data to correct for beam hardening effects. From the corrected projections, the method reconstructs the final tomogram.

  In some embodiments, applying the RPC further comprises reprojecting each phase of the segmented tomogram to produce a simulated monochromatic projection for each material, the simulated monochromatic projection being integral Calculated in attenuation space. In some embodiments, the method generates a container-only pre-correction by selecting a set of radiation that passes only the container and calculating an integral attenuation for the container based on the selected set of radiation. To do.

  In some embodiments, the method includes receiving predetermined container shape data and generating a container only pre-correction based on the predetermined container shape data. In another embodiment, segmenting the tomogram provides container shape data that is used for the purpose of generating a container only precorrection.

  In some embodiments, this segmentation includes performing a k-means clustering technique to distinguish phases using radial distance and grayscale intensity.

  In some embodiments, reconstructing the final tomogram determines (a) one or more estimated material properties of multiple volume elements in the sample using a tomography reconstruction algorithm And (b) determining, by volume element, the most likely material in the material dictionary, by referencing one or more estimated material properties by referring to the stored material property data. Modifying based on the dynamic classification and (c) performing one or more additional iterations of the tomography reconstruction algorithm, each iteration being estimated one or more of the volume elements (D) labeling the reconstructed volume indicating that it is a sample and classifying it into regions of material having similar properties; Further comprising.

  In some embodiments, segmenting the tomogram into two or more phases further comprises separating into at least one phase of the agglomeration material.

  In some embodiments, the method includes referring to the agglomerate material library to identify agglomerate material.

  It should be appreciated that embodiments of the present invention may be implemented or implemented by computer hardware, a combination of both hardware and software, or computer instructions stored in a computer readable non-transitory storage device. . The method may be implemented as a computer program using standard programming techniques and based on the methods and figures described herein, the computer program being configured to include the computer program A computer readable non-transitory storage medium, configured as such, causes the computer to operate in a predetermined and predetermined manner. Each program may be implemented in a high level procedural or object oriented programming language to communicate with a computer system. However, if desired, the programs can also be implemented in assembler or machine language. In any case, the language can be the language to be compiled or interpreted. Furthermore, the program can be run on a dedicated integrated circuit programmed to run the program.

  Further, the methodology communicates with, but is not limited to, separate from the charged particle tool or other imaging device, integral with the charged particle tool or other imaging device, or the charged particle tool or other imaging device Any type operatively connected to a suitable tomographic data scanning device, including personal computers, minicomputers, mainframes, workstations, networked or distributed computing environments, computer platforms, etc. Can be implemented on any computing platform. Aspects of the invention, whether removable or integral with a computing platform, hard disk, optical read and / or write storage medium, non-transitory storage medium such as RAM, ROM or the like or Machine readable code stored on a storage device, wherein the programmable computer configures the computer as it reads the storage medium or storage device to perform the procedures described herein. And for operation, it can be implemented as machine readable code stored so that the computer can read it. Additionally, the machine readable code or a portion of the machine readable code may be transmitted via a wired or wireless network. These types of non-transitory computer-readable storage media and various other types of non-transitory computer-readable storage media include instructions or programs that implement the above-described steps in conjunction with a microprocessor or other data processing apparatus. In the case, the invention described herein includes these non-transitory computer readable storage media. The invention further includes the computer itself when programmed according to the methods and techniques described herein.

  A computer program can be used on the input data to perform the functions described herein and thereby convert the input data to produce output data stored in the non-volatile storage device . This output information may also be output to one or more output devices, such as a display monitor. In a preferred embodiment of the invention, the transformed data represents a physical reality object, which includes generating on the display a specific visual depiction of the physical reality object.

7 is a flowchart of a tomographic reconstruction and segmentation process according to some embodiments. FIG. 7 is a flow chart of a tomographic reconstruction and segmentation process according to another embodiment. FIG. 1 is a schematic block diagram of an exemplary system according to an exemplary embodiment. FIG. 5 is a diagram of a composite object examined using the techniques herein, showing a 2D central slice of the Z (x) volume of the object and a 2D central slice of the ρ (x) volume, and a 3D multiplanar rendering is there. FIG. 7 shows the projection of a simulated object calculated assuming monochromatic x-rays with energies of 60 keV and 120 keV, respectively. FIG. 7 shows the results of several probabilistic matching functions compared to a reference image (in the central slice of the simulated object). FIG. 7 shows the results of several probabilistic matching functions compared to a reference image (in the central slice of the simulated object). FIG. 7 shows a similar comparison of the results when the material dictionary is limited to the chemicals that make up the acquired sample. FIG. 7 shows a similar comparison of the results when the material dictionary is limited to the chemicals that make up the acquired sample. The results of using the probability function [rho [Zn _[rho] In constrained repeated update was compared side by side the results shown in FIG. 4 (e) and 5 (e), a diagram showing the same comparison, each complete database FIG. 8 shows a central slice of the resulting data reconstructed / labeled using a restricted database. 5 is a flow diagram of a process according to an embodiment to correct for beam hardening effects. FIG. 6 is a flow diagram of another beam hardening correction process according to another embodiment that uses container pre-correction. FIG. 6 is a flow chart of a tomographic reconstruction and segmentation process that considers beam hardening over the spectrum of beams used. FIG. 7 shows a series of cross sections of an experimental sandstone sample used to show the effect of this correction method on different segmentations, (a) original reconstruction, (b) cupping artifacts due to beam hardening (C) convex segmentation used to represent being a rock material and (d) resulting correction, (e) more closely fitted segmentation And (f) the resulting correction. FIG. 9 (a) shows a cross section and intensity plot of a simulated silicon sample in a titanium container, and FIG. 9 (a) shows a reconstructed original with gray values readjusted to show cupping artifacts in the sample. Fig. 9 (b) shows an image and Fig. 9 (b) shows a trial of a failed two material correction where the inside of the container was treated as a single material. Fig. 9 (c) shows a gray scale intensity diagram showing the full range of values, Fig. 9 (d) shows the region of interest. It is a figure which shows the intensity | strength figure squeezed and expanded. FIG. 7 is a graph of a three-dimensional point cloud used in the correction process, the black data points showing the fitted regions for the container only correction technique. FIG. 11 shows a cross section and intensity plot of the correction applied to the sample of FIG. 9 (a), the contrast being readjusted to show the effect in the sample, and FIG. 11 (a) shows the beam in the container 11 shows the resulting reconstruction after only the hardening has been removed, and FIG. 11 (b) shows the corrected final reconstruction after the segmentation of the air phase has been performed. FIG. 11 (c) is a cross-sectional and intensity plot of the correction applied to the sample of FIG. 9 (a), with contrast re-adjusted to show the effect within the sample, and FIG. Figure 11 shows a grayscale intensity diagram of a pre-corrected image showing the inverse cupping effect, and Figure 11 (d) shows the final correction compared to the original image. FIG. 6 is a schematic diagram showing the relative path lengths of x-rays with zero fan angle (d1) and x-rays with non-zero fan angle (d2) through objects of various shapes. FIG. 6 is a schematic diagram showing the relative path lengths of x-rays with zero fan angle (d1) and x-rays with non-zero fan angle (d2) through objects of various shapes. FIG. 6 is a schematic diagram showing the relative path lengths of x-rays with zero fan angle (d1) and x-rays with non-zero fan angle (d2) through objects of various shapes. Fig. 13 (a) shows a cross section of a simulated silicon sample in a titanium container containing highly damped iron inclusions, and Fig. 13 (a) shows an image of the original that has been readjusted to show the effect within the sample. Figure 13 (b) with streaking artifacts almost completely removed is a diagram showing the final image corrected using three material correction.

  A tomographic reconstruction and segmentation process is provided, thus providing a system for producing segmentation based on physical attributes of materials using input from a material library. In a preferred embodiment, both reconstruction and segmentation update Z (atomic number) and rho (density) of the volume element in multiple iterations of the reconstruction algorithm. The preferred segmentation process provides voxels with updated labels using best matching techniques for known materials in the material library, and clustering based on voxel values and local intensity variations in tomograms. Do not provide. In some embodiments "labeling" corresponds to some degree to segmentation, but this segmentation is not only based on this data (voxels or regions are locally observed among them) It also depends on "voxels" matching using additional known data of the material dictionary (as separated from others by different intensity differences) (so that each voxel can be processed independently). Dual energy acquisition facilitates the reconstruction of physical attributes to provide this particular labeling / segmentation. In addition, at the next iteration, the segmentation results can be verified. This is because the forward model used in tomographic reconstruction is also based on these same physical properties.

  Various embodiments of the present invention are described in the article "Dual-Energy Iterative Reconstruction for Material Characterization", Proc. SPIE 9212, Developments in X-Ray Tomography IX, 921213 (September 11, 2014) and "Correction of Beam Hardening Artefacts in Microtomography for Samples Images in Container", Proc. SPIE 9212, Developments in X-Ray Tomography IX, 92120 A (September 11, 2014). These articles are incorporated herein by reference.

  Some embodiments provide a dual energy ordered subset convex method for transmission tomography based on material matching with the material dictionary. Such reconstructions include constrained updates that force the reconstructed atomic number (Z) and density (ρ) volume material properties to follow the distribution based on the provided material database. Probabilistic classification techniques are also provided to manage this material distribution. This process as a whole produces chemically segmented volume data and is superior to the sequential labeling calculated after tomographic reconstruction.

In the field of material characterization, the density ρ and atomic number Z (or effective atomic number Z _eff ) of the specimen by combining the data acquired at the two energies E1 and E2 using dual energy tomography ) Can be restored. Such characterization is dependent on both 減 衰 and Z and is based on the decay μ _E (x) model at each point x of the acquired sample, according to the acquisition energy E.

上式で、Ｋ_１およびＫ_２は２つの定数であり、式（１）の第１の部分は、

Where K ₁ and K ₂ are two constants, and the first part of equation (1) is

で示される光電効果を指し、これは、低エネルギーＥに対して現れ、高エネルギーでは無視できる。第２の部分は、Ｃｏｍｐｔｏｎ散乱

Refers to the photoelectric effect, which appears for low energy E and can be ignored at high energy. The second part is Compton scattering

を指し、下式によって定義されるＫｌｅｉｎ−Ｎｉｓｈｉｎａ式ｆ_ＫＮ（Ｅ）に依存する。

, And depends on the Klein-Nishina formula f _KN (E) defined by the following formula.

提案された最初のアルゴリズムは、［Ｚ^ｎρ］投影およびρ投影を両方の取得から計算して、［Ｚ^ｎρ］体積およびρ体積を、標準断層撮影アルゴリズムを使用して直接に再構成することからなる。（この手法は、Ａｌｖａｒｅｚ，Ｒ．Ｅ．およびＭａｃｏｖｓｋｉ，Ａ．、「Ｅｎｅｒｇｙ−ｓｅｌｅｃｔｉｖｅ ｒｅｃｏｎｓｔｒｕｃｔｉｏｎｓ ｉｎ ｘ−ｒａｙ ｃｏｍｐｕｔｅｒｉｓｅｄ ｔｏｍｏｇｒａｐｈｙ」、Ｐｈｙｓｉｃｓ ｉｎ ｍｅｄｉｃｉｎｅ ａｎｄ ｂｉｏｌｏｇｙ ２１（５）、７３３（１９７６）に出ている。）多色Ｘ線を使用する手法が２件の研究論文に出ている（Ｓｉｄｄｉｑｕｉ，Ｓ．およびＫｈａｍｅｅｓ，Ａ．Ａ．、「Ｄｕａｌ−ｅｎｅｒｇｙ ｃｔ−ｓｃａｎｎｉｎｇ ａｐｐｌｉｃａｔｉｏｎｓ ｉｎ ｒｏｃｋ ｃｈａｒａｃｔｅｒｉｚａｔｉｏｎ」、［ＳＰＥ Ａｎｎｕａｌ Ｔｅｃｈｎｉｃａｌ Ｃｏｎｆｅｒｅｎｃｅ ａｎｄ Ｅｘｈｉｂｉｔｉｏｎ（Ｓｏｃｉｅｔｙ ｏｆ Ｐｅｔｒｏｌｅｕｍ Ｅｎｇｉｎｅｅｒｓ，２００４）］、ならびにＨｅｉｓｍａｎｎ，Ｂ．、Ｌｅｐｐｅｒｔ，Ｊ．およびＳｔｉｅｒｓｔｏｒｆｅｒ，Ｋ．、「Ｄｅｎｓｉｔｙ ａｎｄ ａｔｏｍｉｃ ｎｕｍｂｅｒ ｍｅａｓｕｒｅｍｅｎｔｓ ｗｉｔｈ ｓｐｅｃｔｒａｌ ｘ−ｒａｙ ａｔｔｅｎｕａｔｉｏｎ ｍｅｔｈｏｄ」、Ｊｏｕｒｎａｌ ｏｆ ａｐｐｌｉｅｄ ｐｈｙｓｉｃｓ ９４（３）、２０７３〜２０７９（２００３））。この手法は、ｉ）減衰体積μ（ｘ，Ｅ_１）およびμ（ｘ，Ｅ_２）を再構成し、ビーム・ハードニング補正を適用し、ｉｉ）２つの単色測定が使用されたと仮定して式１に基づく式系を解くことによって減衰体積から材料特性体積を回復する２ステップ・プロセスによって、［Ｚ^ｎρ］（ｘ）体積およびρ（ｘ）体積を推定する。しばしばＺ（ｘ）が提供される。これは、Ｚ（ｘ）が、取得された試料を構成する化学元素をより正確に規定するためである。しかしながら、［Ｚ^ｎρ］（ｘ）からＺ（ｘ）を抽出することは、ρ（ｘ）の雑音の影響を受け、材料特性評価の誤りにつながる。

The first algorithm proposed calculates the [Z ⁿ ]] and ρ projections from both acquisitions and reconstructs the [Z ⁿ ]] and 体積 volumes directly using a standard tomography algorithm It consists of things. (This method is described in Alvarez, R. E. and Macovski, A., "Energy-selective reconstructions in x-ray computerized tomography", Physics in medicine and biology 21 (5), 733 (1976).) Methods using multicolor X-rays have been published in two research papers (Siddiqui, S. and Khamees, A. A., “Dual-energy ct-scanning applications in rock characterization”, [SPE Annual Technical Conference and Exhibition (Society of Petroleum Enginee s, 2004)] and Heismann, B., Leppert, J. and Stierstorfer, K., "Density and atomic number measurements with spectral x-ray attenuation method", Journal of applied physics 94 (3), 2073-2079 ( 2003)). This approach i) reconstructs the attenuation volumes μ (x, E ₁ ) and μ (x, E ₂ ), applies beam hardening corrections and ii) assumes that two monochromatic measurements are used The [Z ⁿ ρ] (x) and ρ (x) volumes are estimated by a two-step process that recovers the material property volume from the damping volume by solving the equation system based on Equation 1. Often Z (x) is provided. This is because Z (x) more accurately defines the chemical elements that make up the acquired sample. However, extracting Z (x) from [Z ⁿ ]] (x) is affected by the noise of ρ (x), which leads to an error in material characterization.

  One embodiment herein is a dual energy iterative reconstruction method based on the ordered subset convex (OSC) algorithm for transmission tomography that directly uses the damping function equation (1) in the forward model of the reconstruction algorithm Provide a solution to use. During reconstruction, this process estimates the physical properties of the specimen by matching each voxel to the material dictionary. This estimation is based on probabilistic classification which allows this process to assign the most probable material in the dictionary to each voxel. This classification can be verified during subsequent iterations of the algorithm to correct for false matches. The result of this algorithm is a reconstruction that has already been segmented according to the physical attributes of the scanned object.

  This iterative process avoids the independence of successive steps in the existing sequence of segmentation and labeling. In fact, the approximations during reconstruction may act dramatically on the remaining sequences. Conversely, the iterative process used herein can correct for labeling errors in subsequent iterations. In addition, the inventors of the present invention and their research group at the Australian National University (ANU) have developed several acquisition optimization and dynamic tomography methods for the high cone angle helical scanning microCT facility of ANU. . Some of these methods are described in Optical Engineering and Applications. Myers, G., accepted in (2014). "Improving dynamic tomography quality through a maximum a posteriori model", proposed in SPIE Optics and Photonics. This document provides an iterative method that can be adapted to include material characterization in accordance with the techniques herein, along with the processing goals of the methods themselves.

In the following, the study of several proposed material characterization methods based on numerical simulation is described. The objective is to perform noise efficient efficient voxel material labeling, most likely function with respect to Z (x) and ρ (x) or [Z ⁿ ]] (x) and ρ (x) volumes It is to define. We then introduce a dual energy iterative reconstruction algorithm that includes material characterization during the iteration.

Material Dictionary and Numerical Simulation The studies described below are based on numerical analysis of synthetic objects that simulate rock samples. This sample is modeled using the material database described in Table 1. This numerical sample consists of a cylindrical container randomly filled with different spherical mineral particles. The mineral composition chosen is 60% quartz, 20% kaolinite and 20% calcite. Analytes are rasterized into N ³ voxel volumes for comparison with the reconstructed data (ie, N ³ corresponds to the reconstructed volume size). FIGS. 1 (a)-(c) are diagrams of a composite object examined using the techniques herein, and of the Z (x) volume of that object shaded with the labels shown in Table 1. Figure 2 shows 2D central slice and 2D central slice of ρ (x) volume, and 3D multi-planar rendering. In these figures of the simulated object, (a) shows the central slice of the object Z (x) volume, (b) shows the central slice of the object ((x) volume, and (c) shows 512 ³ 3D shows a 3D visualization of an object imaged into a voxel volume of.

この合成物体は、合成データを使用するときに起こりうる逆クライム効果を防ぐために選択された画素の数の４倍の検出器を使用した模擬取得によって取得される。（源から検出器のそれぞれの画素までの）Ｘ線ラインが、スキャナ・ジオメトリに関してトレースされる。物体を構成するそれぞれの材料内においてそれぞれのＸ線ラインが横切る距離が測定され、

This composite object is obtained by simulated acquisition using a detector that is four times the number of pixels selected to prevent possible counter-crime effects when using composite data. An x-ray line (from the source to each pixel of the detector) is traced with respect to the scanner geometry. The distance traversed by each x-ray line in each material making up the object is measured,

によってサブ画素ｉｉごとに総吸収が計算される。この式で、Ｍ_ｓは、試料を構成する材料の集合であり、μＥ，ｍは、材料ｍのμＥ（式１）に対応する。この実施形態のこのプロセスは、４つのサブ画素値を平均し、Ｐｏｉｓｓｏｎ雑音モデルを仮定することによって、検出器上の最終的な投影値を得る。

The total absorption is calculated for each sub-pixel ii by In this formula, M _s is a collection of materials that make up the sample, and μ E, m corresponds to μ E (formula 1) of the material m. This process of this embodiment obtains the final projection value on the detector by averaging the four sub-pixel values and assuming a Poisson noise model.

上式で、γ_０は、放出された平均光子カウントであり（ブランク走査）、ｉｉ∈ｉは、検出器画素ｉに含まれるサブ画素のインデックスを示す。例えば、図３は、６０ｋｅＶおよび１２０ｋｅＶのＸ線エネルギーを仮定して模擬物体から計算された投影を示す。この研究では、再構成された体積・サイズがＮ^３＝２５６^３ボクセル、検出器サイズが２５６×１９２画素である。したがって、取得は、エネルギー・レベルごとに５１２×３８４の検出器サイズを使用して実行され、２５６^３個の体積を再構成するために２５６×１９２にダウンサンプリングされた。

Where γ ₀ is the average photon count emitted (blank scan) and ii ∈ i is the index of the sub-pixel contained in the detector pixel i. For example, FIG. 3 shows projections calculated from simulated objects assuming x-ray energies of 60 keV and 120 keV. In this study, the reconstructed volume / size is N ³ = 256 ³ voxels, and the detector size is 256 × 192 pixels. Thus, obtaining, using the detector size of 512 × 384 for each energy level is performed, downsampled to 256 × 192 to reconstruct the 256 ³ volume.

Ordered Subset Convex Reconstruction The established iterative expectation maximization (EM) algorithm refers to a broad class of tomographic reconstructions that provides a 3D structural volume of a sample from a set of radiographs. An expected value represented by p (R | μ) observing a measurement value R based on a given sample μ follows a Poisson distribution which can be defined by the following equation.

上式で、ｉは、投影ラインのインデックスであり、検出器位置ｉにおいてカウントされる光子の予想される数Ｒ＾は、

Where i is the index of the projection line and the expected number R ^ of photons counted at detector position i is

である。上式で、Ｗ_ｉｊは、検出器測定値ｉに対するボクセルｊの寄与を規定する重み係数である。Ｗ_ｉｊは普通、ラインｉがボクセルｊ内を移動する距離に比例する。尤度ｐ（Ｒ｜μ）および透過モデルＲ＾に基づいて、この方法は、対数尤度Ｌ＝ｌｏｇ ｐ（Ｒ｜μ）を最大にすることによって、すなわち、偏導関数∂μＬのゼロ根を見つけることによりμに関してＬを最大にする解μを見つけることによって機能する。これは、ｔで反復するＮｅｗｔｏｎ−Ｒａｐｈｓｏｎ法によって解くことができる。

It is. _Where W _ij is a weighting factor that defines the contribution of voxel j to detector measurement i. W _ij is usually proportional to the distance that line i travels within voxel j. Based on the likelihood p (R | μ) and the transmission model R ^, this method maximizes the log likelihood L = log p (R | μ), ie the zero root of the partial derivative ∂μL. It works by finding a solution μ that maximizes L with respect to μ. This can be solved by the Newton-Raphson method repeated at t.

解を安定させ、反復中の発散を回避するため、ｐ（Ｒ｜μ）に加えてペナライゼーションｐ（μ）が導入される。ｐ（μ）は、再構成されたデータのガウス分布に関して下式によって定義される。

In addition to p (R | μ), a penalization p (μ) is introduced to stabilize the solution and to avoid divergence during the iterations. p (μ) is defined by the following equation with respect to the Gaussian distribution of the reconstructed data.

上式で、βは緩和パラメータであり、Ｎ（ｊ）は画素ｊの近傍を示し、Φ（・）はポテンシャル関数、ｙ_ｊｋは重み因子（普通はボクセルｊとｋの間の距離に反比例する）である。この実施形態では、正則化のポテンシャルとして以下のＬｏｇＣｏｓｈ関数が使用される。

Where β is a relaxation parameter, N (j) indicates the neighborhood of pixel j, Φ (•) is a potential function, y _jk is a weighting factor (usually inversely proportional to the distance between voxel j and k ). In this embodiment, the following LogCosh function is used as the regularization potential.

上式で、ｇは感度パラメータである。この関数は、領域のより良好な均質性につながる小さな変動を減衰させ、所与のしきい値ｔを超えるエッジを保存する。式（６）において置換Ｌ＝ｌｏｇ ｐ（Ｒ｜μ）ｐ（μ）を行うことにより、ペナライゼーションを用いた全体の順序付きサブセット凸（ＯＳＣ）アルゴリズムは、体積μのそれぞれのボクセルＪを解が収斂するまで更新するために、ｔおよびサブセットＳ＋１で反復することからなる。

Where g is the sensitivity parameter. This function attenuates small variations leading to better homogeneity of the area and preserves the edge beyond a given threshold t. By performing the permutation L = log p (R | μ) p (μ) in equation (6), the entire ordered subset convex (OSC) algorithm using penalization solves each voxel J of the volume μ To repeat on t and subset S + 1 in order to update until.

上式で、Ｒ＾^ｔ（ｉ）は式（５）を使用することによってμ^ｔから計算された予想される光子カウント、α^ｔは緩和パラメータ、Ｓ（ｓ）はサブセットｓ中の放射線写真、Φ^Ｉ（・）およびΦ^ＩＩ（・）はそれぞれ、Φ（・）の第１および第２の導関数である。サブセットｓを使用したこの更新によって得られる体積は、次のサブセットｓ＋１の開始体積として使用される。主反復ｔは、全てのサブセットが処理されたときに完了となる。

Where R ^ ^t (i) is the expected photon count calculated from μ ^t by using equation (5), α ^t is the relaxation parameter, S (s) is the radiograph in subset s, ^{I I} (.) And Φ ^II (.) Are the first and second derivatives of ((.) Respectively. The volume obtained by this update using subset s is used as the starting volume of the next subset s + 1. The main iteration t is complete when all subsets have been processed.

Material characterization: Simulation studies μE ₁ (x) and μE ₂ (x), respectively, are tomograms reconstructed from the two acquisitions at energies E ₁ and E ₂ using the standard OSC algorithm. The volume pair (Z,)) or ([Z ⁿ ]],)) can be recovered by solving one of the following equation systems at each voxel x.

上式で、

In the above formula,

、Ｂ_１＝Ｋ_２ｆ_ＫＮ（Ｅ_１）、

, B ₁ = K ₂ f _KN (E ₁ ),

、およびＢ_２＝Ｋ_２ｆ_ＫＮ（Ｅ_２）である。

, And B ₂ = K ₂ f _KN (E ₂ ).

In this section, we define several probabilistic material matching functions that assume that one reconstructed voxel corresponds to one material in the material dictionary and compare those functions. Since 関 数, Z or [Z ⁿ ]] can be determined using Equation 10, the matching function compares one or several reconstructed material properties with their eigenvalues in the material database based on.

Reconstructed data Z (x), ((x) and [Z ⁿ ]] (x) are each considered to follow a Gaussian distribution around valid material values. Then, the probability that voxel x corresponds to material m ∈ M is

またはこれらの式の結合

Or a combination of these expressions

によって定義することができる。上式で、ρ_ｍ、Ｚ_ｍおよび［Ｚ^ｎρ］_ｍは、ディクショナリＭ内の材料ｍの化学特性であり、σρ_ｍ、σＺ_ｍおよびσ［Ｚ^ｎρ］_ｍはそれぞれ、それらの標準偏差である。大きな標準偏差は不正確な推定に帰着するであろう。これは、ディクショナリ内の全ての材料が良好な確率（滑らかな推定）を達成するためである。反対に、あまりに小さな偏差は、過度に判別的な推定につながるであろう。ａ_ｍを、材料ｍに対する比較値の１つとする（すなわち、ａ_ｍ＝｛ρ_ｍ、Ｚ_ｍまたは［Ｚ^ｎρ］_ｍ｝）。ａ_ｎは、ｎ≠ｍおよびａ_ｎ＝ａｒｇｍｉｎ｛ａ_ｍ−ａ_ｍｉ，ｍ_ｉ∈Ｍ∧ｍ_ｉ≠ｍ｝であるような材料ｎ∈Ｍの対応する値を示す。材料特性ａに関するｍの標準偏差σ_ａを下式によって定義する。

It can be defined by Where _{m m} , Z _m and [Z ⁿ ]] _m are the chemical properties of material m in dictionary M, and σ _{m m} , σ Z _m and σ [Z ⁿ ]] _m are their standard deviations respectively It is. Large standard deviations will result in inaccurate estimates. This is because all the materials in the dictionary achieve a good probability (smooth estimation). Conversely, too small a deviation will lead to an overly discriminatory estimate. Let a _{m be} one of the comparison values for material m (ie, a _m = {ρ _m , Z _m or [Z ⁿ ]] _m }). a _n represents n ≠ m and _{a n = argmin {a m -a} mi, m i ∈M∧m i ≠ m} corresponding values of the material n∈M as is. The standard deviation sigma _a of m regarding material properties a defined by the following equation.

上式で、ｐ_ｔ∈｛０，１｝は、材料ディクショナリ内のｍとそのより近い近傍ｎとの間の識別可能性を定義する。材料ごとにこのような標準偏差を設定することは、データベースを構成する化学物質に対する確率的マッチング関数の難しさを局所的に調整することを可能にする。３つの単一ベストマッチング関数（式１１）の要約が表２に示されている。しかしながら、挙げられた関数は、限定を意図したものではなく、適当な他の関数を使用することもできる。

Where p _t ∈ {0, 1} defines the distinguishability between m and its closer neighbor n in the material dictionary. Setting such a standard deviation for each material makes it possible to locally adjust the difficulty of the probabilistic matching function for the chemical substances constituting the database. A summary of the three single best matching functions (Equation 11) is shown in Table 2. However, the listed functions are not intended to be limiting and other suitable functions can be used.

The inventor of the present invention describes in the article entitled “Dual-Energy Iterative Reconstruction for Material Characterization”, which is incorporated herein, five matching functions (Equation 11 and Equation 12) (in this example p _t = 0). According to .8), data are provided for examining the probability material distribution in the material database of Table 1. These results are particularly relevant when using a restricted material database (containing only the chemicals used in the numerical samples), as compared to using a larger or complete database for the same probability function. It shows an advantageous effect on this distribution.

それぞれのマッチング関数の効率を定量化するため、最初に、上で説明したシミュレーションを使用して２つのエネルギーで数値取得を計算する。次いで、上で詳細に説明したＯＳＣアルゴリズムを使用して、体積（μ_１，μ_２）を再構成する。ボクセルごとに式１０を解くことによって体積（Ｚ，ρ）および（［Ｚ^ｎρ］，ρ）を推定する。

In order to quantify the efficiency of each matching function, we first calculate the numerical acquisition at two energies using the simulations described above. The volume (μ ₁ , μ ₂ ) is then reconstructed using the OSC algorithm described in detail above. The volumes (Z,)) and ([Z ⁿ ]],)) are estimated by solving equation 10 for each voxel.

Finally, apply classification labeling using each matching function. This classification labels each voxel with a material m∈M such that m maximizes the probability function used. The shaded images in FIG. 4 (b-f) show the results of each method as compared to the reference image (center slice of the simulated object). Correctly labeled numbers (Corr.% In FIG. 4) indicate that each function has correctly recovered from the complete dictionary in at least 50% of the voxels making up the sample. Even though it is observed that the difference between the different results is not so great, the use of pZp gives the best results. FIG. 5 (b-f) shows the result of the most likely material when the material dictionary is limited to the chemicals that make up the acquired sample. This labeling process achieves better results regardless of the matching function. However, it can be seen that the results obtained with pZ and p [Z ⁿ ]] are higher than others. Thus, one can expect to perform an exact match from the (Z,)) and ([Z ⁿ ]],)) volumes.

Dual Energy Reconstruction and Material Matching It has been found herein that material labeling can be performed from (Z, ρ) and ([Z ⁿ ]],)) volumes. The latter has the advantage that it is not limited by noise in the 体積 volume. However, although the matching process detailed in the previous section is performed after tomographic reconstruction, this process can also be performed during the iterative reconstruction technique and corrects the incorrect matching during subsequent iterations of the algorithm can do. Therefore, the modified dual energy OSC algorithm will now be detailed. This algorithm is based on the reconstruction of the ([Z ⁿ ]],)) volume and includes a constrained update of voxels at each iteration. This update brings each voxel in the reconstruction closer to the most likely material in the database.

  The forward model of this algorithm is

の直前の推定に基づく。したがって、式５のμＥ_１（ｊ）およびμＥ_２（ｊ）は下式によって与えられる。

Based on the previous estimate of Therefore, μE ₁ (j) and μE ₂ (j) of Equation 5 are given by the following equation.

次いで、ＯＳＣ更新ステップ式９をＥ１およびＥ２に対して別々に適用することによって、更新値Δμ^ｔ（ｘ，Ｅ_１）およびΔμ^ｔ（ｘ，Ｅ_２）を得る。Δ［Ｚ^ｎρ］^ｔ（ｘ）およびΔρ^ｔ（ｘ）を、

The updated values Δμ ^t (x, E ₁ ) and Δμ ^t (x, E ₂ ) are then obtained by applying OSC update step equation 9 separately to E ₁ and E ₂ . Let Δ [Z ⁿ ]] ^t (x) and Δρ ^t (x),

の更新値とする。それらは、それぞれのｘについて下式を解くことによって回復される。

The updated value of They are recovered by solving the following equation for each x.

それぞれの超反復ｔの終わりに、ベストマッチング関数に基づいて、制約付き更新を追加する。この更新は、最も確からしい材料ｍがｐα＝ｐ（ｘ＝ｍ）＞ｐβであるような材料である場合に、［Ｚ^ｎρ］^ｔ（ｘ）およびρ^ｔ（ｘ）を修正することからなる。ｐβは、確率しきい値である。そのような場合、以下のように制約付き更新が実行される。

At the end of each super-iteration t, add a constrained update based on the best matching function. This update is from modifying [Z ⁿ ]] ^t (x) and ρ ^t (x), where the most likely material m is such that pα = p (x = m)> pβ Become. pβ is a probability threshold. In such a case, constrained update is performed as follows.

上式で、［Ｚ^ｎρ］_ｍおよびρ_ｍは、ｍの化学特性である。この制約付き更新は、ｍが最も確からしく、ｐ（ｘ＝ｍ）が、高い（しきい値ｐαよりも高い）確率を既に提供している場合に、再構成されたデータを、ディクショナリ内の材料ｍの理論値により近づける。ｐαよりも低い最良の推定を提供するマッチングは、この制約付き更新によって修正されない。この制約は、解を、ディクショナリ内の材料により近づけるという利点を有する。

Where [Z ⁿ ]] _m and _{m m} are the chemical properties of m. This constrained update is such that if m is the most probable and p (x = m) already provides a high (higher than threshold pα) probability, then the reconstructed data will be in the dictionary It will be closer to the theoretical value of material m. The matching that provides the best estimate below pα is not corrected by this constrained update. This constraint has the advantage of bringing the solution closer to the material in the dictionary.

Next, apply this algorithm using the same set of acquisitions at E ₁ and E ₂ used in the previous section. In this example, pβ = pt = 0.8, and the constrained update used the probability function [[Zn _ρ ]. This is because this function averaged the best results in preliminary studies. The central slices of the resulting data reconstructed / labeled using the complete and restricted databases, respectively, are shown in FIG. 6 (b-c). This technique correctly labeled over 66% of the voxels containing the sample from the complete material database, and correctly labeled over 87% of the voxels containing the sample from the restricted material database. On the other hand, they were 54% and 77%, respectively, when they were segmented / labeled after reconstruction. Thus, using constrained iterative reconstruction, gains of more than 10% are observed for the same probability function. These results demonstrate the effectiveness of combining tomographic reconstruction with material labeling. However, those results still show material labeling errors due to noise and partial volume effects.

The foregoing has described several probabilistic techniques for labeling each voxel of a tomogram according to physical properties contained in a material dictionary. The matching techniques are based on local identifiability to their closest neighbors of each material in the database. Thus, material labeling from (Z,)) or ([Z ⁿ ]],)) volumes with equivalent quantitative results can be achieved. However, when labeling is performed after segmentation, it is difficult to correct material labeling according to the physical model. Furthermore, since the iterative reconstruction can incorporate physical interactions with x-ray radiation based on material properties, based on the reconstruction of the ([Z ⁿ ]],)) volume, at the end of each iteration Provided dual energy OSC algorithm including constrained update. This additional update brings the value of ([Z ⁿ ]],)) closer to the theoretical value of the material m providing the best probability. Such a process achieves better object reconstruction and more accurate labeling together as compared to the process of sequentially performing tomographic reconstruction and material classification. By applying a beam hardening correction (as provided by Siddiqui et al. And Heismann et al.) And then assuming two monochromatic energies, the above embodiment is applied to two polychromatic spectra Extended to These techniques extract the ([Z ⁿ ]],)) projections and reconstruct the ([Z ⁿ ]],)) volume, with the intermediates Δμ ^t (x, E ₁ ) and Δμ Avoid ^t (x, E ₂ ) and correct for polychromatic spectra at once. Various embodiments using dual polychromatic x-ray spectra centered at different energies can include a combination of these techniques and dual reconstruction / labeling as described herein.

  FIG. 1 (a) is a flowchart of a tomographic reconstruction and segmentation process in accordance with some embodiments. Process 100 begins at block 101 and provides at block 101 a material library or data storage that contains the characteristics of the material expected to be present in the sample. In this embodiment, this property includes at least the atomic number (Z) and density (ρ) of the expected material. The data storage area can include data regarding any material that can be examined using tomography, such as, for example, medical materials (biological or tissue samples), minerals or semiconductors. As discussed above, the data storage area is limited to the material that is expected to be present in the sample or is likely to be present in the sample, selected based on the type of sample to be analyzed Is preferred. The data storage area can also include other materials, such as materials that may be present in the sample but not as likely in the sample. Some processes can include selecting data storage to use based on the type of sample to be analyzed. Referring to FIG. 1 (c), which shows a block diagram of an exemplary system, a material library or data storage 19 comprising a plurality of material properties is shown, each set of material properties as discussed above Are associated with the corresponding respective materials.

  Next, at blocks 102 and 103, the process performs tomographic data acquisition by directing x-rays from multiple directions towards the sample and measuring the intensity of x-rays across the sample. Any suitable tomographic x-ray process may be performed, including a tomographic x-ray process using multiple x-ray energy spectra, such as the dual energy acquisition process described above. Other embodiments may use more than two energy levels, such as three or four, and some embodiments use only one energy level. As mentioned above, this process can use radiography or transmission tomography, depending on the type of sample to be examined, such as, for example, a medical sample (biological sample or tissue sample), mineral or semiconductor.

  As shown in the schematic block diagram of FIG. 1 (c), the scanning and detection of the tomographic X-rays are preferably performed under the control of the system controller 15. The system controller 15 commands the dual energy x-ray beam 11 to irradiate the sample 13 as desired. A rotating sample holder or sample stage 14 holds the sample and rotates the sample to achieve the desired number of tomographic projections. The X-ray detector 12 measures the intensity of X-rays that have passed through the sample. The tomography system 10 of FIG. 1 (c) is described in US Patent Provisional Application No. 61 / 981,330, entitled "X-Ray Microscopy and Tomography Imaging Methods and Stictions," a scanning electron microscope (SEM). It can be incorporated into another scanning device such as a device. This provisional US patent application has been changed to US patent application Ser. No. 14 / 645,689 entitled “High Aspect Ratio X-Ray Targets And Uses Of Same” filed on March 12, 2015. Both of these applications are incorporated herein by reference. Such an apparatus for realizing an embodiment of the present invention includes an electron source, an electron focusing column that focuses electrons from the electron source into a beam, and a target that emits an X-ray when the electron beam collides. The x-rays pass through the sample and are detected by one or more x-ray detectors. In other embodiments, the x-ray emitter 11 can direct radiation from the synchrotron from multiple directions towards the sample. The system controller 15 is typically implemented as a microprocessor motherboard or computer with a user interface 16 having a display and user input devices. The processor 17 is operatively connected to a storage device which stores a material library 19, instruction codes for tomographic reconstruction and segmentation, and drivers and controller software for the devices 11, 12 and 14. Processor 17 may include a dedicated image graphics processor or one or more separate graphics cards to accelerate analysis and display of tomographic data.

  Refer back to FIG. 1 (a). After obtaining the tomographic data, process 100 proceeds to block 104 and initiates a tomographic reconstruction and segmentation algorithm at block 104 that determines the atomic number and density of the material in the sample by iterative tomographic reconstruction. Block 104 determines an estimated atomic number and an estimated density of the number of volume elements in the sample using a tomographic reconstruction algorithm such as one of the suitable algorithms described above. Next, at block 106, the process calculates a probabilistic classification to determine the most likely material in the material dictionary for one or more volume elements. This calculation in this step has any suitable probability, including the classification function described above as having the best performance, the classification function that is evaluated to not work the same way, and other suitable classification functions. Dynamic classification functions can be included.

  Next, at block 107, the process modifies the estimated atomic number and estimated density by referencing the atomic number and density of the known material from the material library or data storage. As discussed above, this step preferably uses the results of the classification function from the immediately preceding step, and uses the most likely material properties identified. In some embodiments, this step corrects these numbers using a constrained update that includes thresholds such as those described above, and the corrected numbers are determined in the material dictionary Including getting closer to the numerical value of the most probable material.

  Next, at step 108, the process labels the reconstructed volume elements that are indicated to be samples and classifies them into regions of material having similar atomic numbers and densities. The label preferably indicates that the area of the material consists of the most likely determined material, or indicates that the area of the material comprises the most probable material. If step 108 was reached in this further process iteration, step 108 may update the label of the volume element, which updates the previously identified region in the reconstructed model. It can be changed. Note that this label can be a data element separate from the estimated attribute of the corresponding respective volume element used in the reconstruction algorithm. These labels may be data elements that, for each corresponding updated volume element, generate an association with one or more surrounding volume elements (voxels), or the boundaries of the identified area The voxels can be labeled, or the voxels can be grouped according to known techniques representing voxel data previously used in the field during the segmentation process following tomographic reconstruction as discussed in the background section. Can be combined or associated with

  Next, at block 109, the process determines if the iteration of the tomographic reconstruction and segmentation algorithm is complete. This block can include checking if the data has converged (ie, whether the data has changed a lot from the previous iteration) or if a predetermined number of iterations have been completed Can include checking. The number of iterations is selected based on the type of sample or the observed complexity of the three-dimensional structure of the sample. Other suitable methods of determining whether the iteration is complete can also be used. If the iteration is not complete at block 109, the process proceeds to block 110 and proceeds to the next iteration of the reconstruction algorithm at block 110. Returning to block 104 in this embodiment, the voxel representation of the sample uses current data. If the iteration is complete at block 109, the process ends at block 112. The output of this reconstruction algorithm is the current atomic number, density and label of the volume element indicating that it is the sample. Preferably, each iteration performing the indicated steps of the tomographic reconstruction algorithm updates the corrected atomic number and density of the volume element.

  FIG. 1 (b) is a flowchart of a tomographic reconstruction and segmentation process according to another embodiment. Process 200 begins at block 201 and at block 201 determines and provides a material library or data storage area that contains the characteristics of the material expected to be present in the sample. The remainder of this process proceeds in the same manner as the process described with respect to FIG. 1 (a). This flow chart is provided to illustrate that various embodiments can use other data characteristics. In various embodiments, the material properties estimated at block 204 and adjusted at block 207 are measured directly or indirectly by any type of tomographic scan with reference to data stored in data storage Or any suitable property that can be derived. For example, in embodiments using radiation tomography data scanning, the intensity level of the scan may be related to the concentration of radioactive element or other tracer introduced into the sample or object under study. For example, in positron emission tomography (PET), the system emits gamma rays (or annihilation photons) indirectly emitted by a radionuclide (tracer) that emits positrons introduced onto biologically active molecules in the body. Detect a pair. In such an embodiment a PET scanning device is used. As another example, the material dictionary or data storage 19 can store mass attenuation coefficient data of x-rays, gamma rays or bremsstrahlung, and the use of tomographic scanning can measure these transmissions . As another example, estimated material properties may be X-ray attenuation at multiple energy levels, eg at one, two, three or four energy levels, or a complete attenuation curve as a function of X-ray energy It can be done. The additional material properties used in the estimation of block 204 and the classification and labeling of blocks 206 and 207 may be different properties or different property sets, changing such properties between different iterations of the algorithm and iterations. can do. For example, one embodiment assigns one label to its volume element using atomic number and density, but then in forward estimation or simulation of iterative reconstruction, or in probabilistic classification, X-ray energy Use different additional information about the material, such as the complete decay curve as a function. As discussed above, the data storage area is limited to the material that is expected to be present in the sample or is likely to be present in the sample, selected based on the type of sample to be analyzed Is preferred. In some embodiments, block 201 includes selecting a material library to use based on the type of sample involved.

  Although steps 204, 206, 207 and 208 are described in this embodiment, other embodiments may not include all of these steps or, in a particular iteration, all of these steps It should be noted that fewer steps may be included. For example, some embodiments may not include labeling. As discussed above, the tomographic reconstruction algorithm implemented at block 204 may be any suitable algorithm such as, for example, algebraic reconstruction techniques (ART), simultaneous ARTs, simultaneous iterative reconstruction techniques (SIRT), ordered subsets SIRT, Multiplicative Algebra Reconstruction Technique, Maximum Likelihood Expectation-Maximization, Ordered Subset Expectation-Maximization, Ordered Subset Convex Algorithm, Iterative Coordinate Descent (ICD), Ordered Subset ICD or Model Based Iterative Replay It can include configurations. Additionally, the process can include correcting the data using a tomography algorithm that corrects for beam hardening.

  Preferably, determining the composition of the sample volume during reconstruction comprises segmenting the sample into regions of common composition, said segmentation being a voxel characteristic determined at the completion of the iterative reconstruction Instead of being run based on it, it is run during iterative reconstruction. In some embodiments, the method further comprises labeling the reconstituted volume indicating being a sample and classifying it into regions of material having similar properties.

  Preferred embodiments perform one or more additional iterations of the tomographic reconstruction algorithm, each iteration updating one or more estimated material properties of the volume element.

  In some embodiments, modifying estimated material properties by reference to material properties of a known material is a probabilistic classification to determine the most likely material in the sample by volume element Including determining those known materials based on Modifying the estimated atomic number and the estimated density for at least one of the iterations corrects those numbers using constrained updating, and the corrected numbers are Including being closer to the determined most probable material value in the dictionary.

  Although the preferred embodiment described above uses dual energy ordered subset convex methods for transmission tomography, the invention is applicable to any type of iterative tomography reconstruction. The techniques herein may be used with emission tomography techniques and with transmission tomography. Additionally, the preferred embodiments described above include determining atomic number and density, and comparing their properties to the atomic number and density of known materials in the material dictionary. The invention is applicable to determining any material property. Furthermore, the reference material properties can also be stored in other forms, such as in a damping table.

  Another aspect of the invention is the article "Correction of Beam Hardening Artefacts in Microtomography for Samples Images in Container", Proc. SPIE 9212, Developments in X-Ray Tomography IX, 92120 A (September 11, 2014), provides a method for image correction with respect to beam hardening, and this paper is incorporated by reference into this patent application Part of

Beam Hardening Corrections Along with the tomographic segmentation techniques discussed above, some embodiments of the present invention provide an improved non-referenced multi-material beam hardening correction method. This method focuses on maintaining data quality for 3D imaging of geological materials, and looks to automation. This technique can be applied to other fields as well. This beam hardening correction technique can be applied in combination with the reconstruction and material identification techniques discussed above, or alternatively can be applied independently using other reconstruction techniques. it can. A referenceless post-reconstruction (RPC) correction technique is provided that applies this correction in integral attenuation space. A container only pre-correction technique is also provided that allows for the automation of the segmentation process required for the beam hardening correction method. The method can use various segmentation techniques, including techniques that reference material libraries, at various stages of the process.

  When a container is introduced, the resulting beam hardening artifacts can make it difficult to segment the material inside the container. In order to solve this problem and to allow this technique to be applied to the sample in the container, herein reference-free post-reconstruction (RPC), which allows accurate automatic segmentation of the sample in the container A modification of the law is provided. The section entitled No Reference Post Reconfiguration Correction shows an improved Referenceless Post Reconfiguration Correction (RPC), which is the primary technique detailed herein. The lower section entitled Correction on the assumption of a single-material sample examines the validity of the RPC method when the sample is treated as a single material, and then the section entitled Correction of the sample in the container Now, this technique is improved for the sample in the container. Correction Space: The terms entitled Intensity or Attenuation provide a discussion of the implications of performing correction in intensity space or attenuation space.

No Reference Post Reconstruction Correction Overview of the RPC Technique The RPC method used here as the basis of the following technique is different from the Krumm et al. Method discussed above. This method does not apply correction in intensity space (I / I ₀ ) as the Krumm et al. Technique does, but applies correction in integral attenuation space (R). The inventors of the present invention have found that there are several distinct advantages to using R-space. Their advantages are shown below. As used herein, when referring to the RPC method herein, it generally refers to a method that has been modified to apply in the attenuation space. However, this is not intended to be limiting, and some of the techniques described herein may also benefit from use with prior RPC techniques that use intensity space. . Thus, the description herein should be construed to support the use of any of the various techniques individually combined with the various techniques described below.

Given normalized projection data I (v, x) / I ₀ , the integrated attenuation R can be derived as R = −log (I / I ₀ ) (in accordance with Equation 17), This is reconstructed using standard filtered backprojection to produce a tomogram containing beam hardening artifacts. The RPC method proceeds using segmentation of the reconstructed tomogram to distinguish each material present. Each phase (material) in the segmented image is then reprojected to produce a simulated monochromatic projection for each material. Because each reprojection assumes a monochromatic beam, the value of each pixel in the resulting projection is directly related to the path length through the corresponding material. For each material i, this preferred process calculates beam hardening corrections in the integral attenuation space.

上式で、

In the above formula,

は、点ｒにおけるｉ番目の材料セグメント化の２進数表現である。あるいは、Ｉ空間で算出することもできる。

Is the binary representation of the ith material segmentation at point r. Alternatively, it can be calculated in I space.

材料ごとに１つの投影があるため、Ｎ個の材料では、画素ごとに、Ｎ個のベクトル

Since there is one projection per material, N materials have N vectors per pixel

からなる投影データが得られる。次いで、Ｒ^ｎ＋１において散布図を構築することができ、この図では、それぞれの点が単一の画素に対応し、垂直軸に、その画素のオリジナルの投影データ値が示され、残りの軸は、その画素の再投影ベクトルを示す。したがって、それぞれの点は、単一の放射線に沿った実験に基づく減衰およびそれぞれの材料を通過するその放射線の経路長を表す。オリジナルのビームが単色である場合、それらの点は、減衰が経路長に正比例したＮ次元の超平面となるであろう。

Projection data is obtained. A scatter plot can then be constructed at R ^{n +1} , where each point corresponds to a single pixel, the vertical axis shows the original projection data values for that pixel, and the remaining axes are , Indicates the reprojection vector of the pixel. Thus, each point represents the experimental attenuation along a single ray and the path length of that ray through each material. If the original beam is monochromatic, those points will be N-dimensional hyperplanes whose attenuation is directly proportional to the path length.

多色ビームからのデータに関しては、それらの点がＮ次元の超曲面上にあり、減衰は、それぞれの材料の経路長の、ある関数である。

For data from polychromatic beams, those points lie on N-dimensional hypersurfaces, and the attenuation is a function of the path length of each material.

このＲＰＣ法は、最初に、この多色データを通る超平面を、最良の単色近似を表すように当てはめることにより進行する。これは、減衰係数μ_ｉの最良の推定を見つけることと等価である。このデータを通る超曲面を当てはめ、その超曲面と超平面との差が、それぞれの点に適用する補正の必要量を形成する。

The RPC method proceeds by first fitting a hyperplane through the polychromatic data to represent the best monochromatic approximation. This is equivalent to finding the best estimate of the damping coefficient μ _i . By fitting a hypersurface through this data, the difference between the hypersurface and the hyperplane forms the required amount of correction to apply to each point.

このＲＰＣ法は補正として機能すること、すなわち、このデータを変換せず、むしろこのデータを改良しようとすることは、強調されるべきである。したがって、オリジナル・データと投影されたデータの間の完全な一貫性を達成する必要はなく、このことは、完全な情報が入手できないときに有利であり、仮定における誤りまたはセグメント化における誤りに対してこの方法を強くする潜在性を有する。

It should be emphasized that this RPC method functions as a correction, ie does not convert this data, but rather tries to improve this data. Thus, it is not necessary to achieve complete consistency between the original data and the projected data, which is advantageous when complete information is not available, against errors in assumptions or errors in segmentation. It has the potential to make this method stronger.

Segmentation Considerations Krumm et al. Demonstrate that the Krumm et al. RPC method applied in intensity space works very well when all materials (phases) are correctly distinguished from one another for a range of sample types Showed that. However, often the goal of μCT analysis is the separation of the different phases, and no correction is necessary if accurate segmentation is possible on the original image. If incorrect segmentation is used for the correction, new artifacts may be introduced as a result. This occurs because if phase 1 requires a different amount of correction than phase 2, then the area of phase 1 incorrectly identified as phase 2 is corrected differently from the other areas of phase 1 This is because there is a possibility that the incorrectly identified phase 1 region is no longer similar to the other regions of phase 1. Therefore, since correction methods such as Krumm's RPC method may introduce artifacts indistinguishable from the actual features, care must be taken when using those correction methods.

  Although the Krumm method has a certain level of robustness for slightly incorrect segmentation, the auto-segmentation method is desirable for the techniques herein, and therefore, the Krumm method is used directly It will not provide quality control over segmentation to mitigate the detrimental effects of the introduced artifacts. Automation is clearly advantageous, especially when both data collection and reconstruction take very long time. By minimizing user input after the process has begun, the entire process of data collection, reconstruction and correction can be performed without assistance and can be continued overnight. Furthermore, unlike the more subjective user-driven segmentation, automation has the ability to give consistent results. The disadvantage is that it is difficult to ensure that the quality of segmentation is at the desired level of accuracy. In order to avoid potential loss of fidelity due to incorrect segmentation, some embodiments herein provide for grouping multiple materials that are not easily distinguished as a single phase. This will be further described later under the heading of correction assuming a single material sample.

  After grouping the materials which still have to segment certain phases, such as containers and samples, but which are not easily distinguished, auto-segmentation becomes relatively easy. The choice of an automatic segmentation method depends on the type and quality of the image to be processed. In some embodiments, dual energy iterative techniques using material matching as discussed above provide advantageous segmentation. In other embodiments, the use of watershed segmentation using gradients defined by Sobel filters may be sufficient to separate these phases. In the preferred embodiment described below, this method uses the k-means clustering technique, since the shape of the phase was approximately cylindrical. In this technique, radial distances and grayscale intensities are used to distinguish phases.

Correction Assuming a Single Material Sample The first sample to be considered is experimental data of a 5 mm diameter bentheim sandstone imaged at 80 kV using a 0.3 mm flat aluminum sheet as a filter. This filter is insufficient for the purpose of reducing the beam hardening of such samples. FIG. 8 shows a series of cross sections of experimental sandstone samples used to show the effect of this correction method on different segmentations. FIG. 8 (a) shows this reconstructed two-dimensional slice of sandstone, and FIG. 8 (b) shows the contrast adjusted to highlight the cupping artifact of this rock, which is an artifact of beam hardening. This same slice is shown. This embodiment treats the entire sandstone sample as one phase consisting of air and silicon, rather than performing beam hardening correction by performing internal segmentation to separate pore spaces. Perform beam hardening correction by. In this process, it is assumed that all radiation passing through a certain amount of air-silicon phase suffers a similar amount of beam hardening.

  FIG. 8 (c) shows the first possible segmentation, assuming that this rock is convex for modeling. The resulting corrected image shown in FIG. 8 (d) corrected for cupping artifacts except for the slightly bright narrow ring around the sample. We state that, first, there is a dispersion in the attenuation of X-rays passing through a segmented phase of the same length, and that the dispersion is related to the proportion of air and silicon along that path Explain by. For radiation passing through only a small segmented region, i.e. for near tangential radiation, the variance of the attenuation is greater than the value of the attenuation. The radiation that passes just the edge of the segmented region actually passes through 100% air or 100% silicon, but those radiations are handled the same, resulting in a slightly bright ring Will appear. The average solid volume fraction is much smaller at the edge of the sample, so the grazing radiation passes through a much smaller amount of solid material than the solid material through which the "internal" radiation passes with respect to the path length through the sample. This effect decreases rapidly as the damping increases, so it may be useful to ignore the edge of the sample.

  If necessary, as in FIG. 8 (e), this edge ring artifact can be improved by making the segmentation closer to the sample periphery. Although this segmentation correction shown in FIG. 8 (f) reduces this problem, there is still a very faintly dark ring at the edge of the sample. Such techniques effectively function to correct beam hardening in multi-material objects by segmenting the multi-material object as a single material and applying a referenceless post-reconstruction correction technique It has been shown. However, there is a limitation as to what materials can be grouped as the same phase, ie all radiation of the same path length through the segmented phase cause similar beam hardening. This assumption is satisfied if the damping of all the materials grouped into a group is similar, or if the materials are uniformly heterogeneous, as is the case in this example. is there.

Correction of the sample in the container Two materials: correction given the sample and the container Next, consider simulated data that represents the case of a quartz rock sample in a highly damped titanium container. Fig.9 (a)-(b) shows the cross section of the simulation silicon sample in a titanium container. FIG. 9 (a) is a simulation of this procedure generated by combining multiple reprojections of different energy into one using the attenuation coefficients of silicon and titanium obtained from the National Institute of Standards and Technology (NIST) Indicates As shown in FIG. 9 (a), this beam hardening produced a bright halo effect near the inner edge of the container. This halo presents challenges when segmenting the air around the rock, especially if the segmentation is attempted automatically. The container itself can be easily segmented, for example by using automatic k-means clustering on intensity values. The container is considered as a single phase to avoid segmentation inside the container. For clarity, this process segments the image into two regions. The container is the first phase and everything inside the container is grouped as one group and treated as the other phase. The RPC method is then applied. The subsequent correction is shown in FIG. 9 (b). The intensity diagram in FIG. 9 (c) shows that the cupping artifact has been removed in the container phase. This is shown by the fact that the gray value line in that area is flat.

  FIG. 9 (d) shows the scaled intensity diagram. This intensity diagram shows that the cupping of the rock sample itself is reduced towards the center, but is still present near the vessel. Furthermore, this correction process does not necessarily facilitate air segmentation. This is because bright halos are still present near the inner edge of the container. While segmentation of air can be attempted, and repeated application of this technique can improve segmentation of air, each segmentation reinforces itself in subsequent reconstructions, and so this iterative method Does not guarantee that the segmentation is correct. The next section details modifications of the RPC method similar to those performed by Van de Casteele et al. And in particular Paziresh et al. This requires a two step correction, the first correction step removing the inner halo artifact caused by the container.

  FIG. 7 (a) is a flow diagram of a process in accordance with one embodiment for correcting beam hardening effects using some of the techniques described above. This process block is preferably performed automatically by the tomography system 10 under the controller of the system controller 15 after activation of the tomography system 10. The activation of the tomography system 10 is achieved by the operator inputting the appropriate parameters remotely via the user interface 16 or via the network, or a host adapted to control the process via the network Achieved by the controller. The process begins at block 702 with obtaining or loading data from a series of tomographic data scans. This block can include controlling the tomographic scanning system 10 to actually perform a data scan, or simply include retrieving a previously acquired data scan from a data store it can. Next, at block 704, a series of tomographic data projections, usually acquired using a topography scanner, are reconstructed into an initial set of tomograms or slices to create a three-dimensional model of the scanned sample. This model contains several tomograms, or they are commonly referred to as slices. In this initial reconstruction, the tomogram may contain a considerable amount of the beam hardening effect discussed above and needs to be corrected.

  Next, at block 706, the process begins to perform a beam hardening correction procedure. In this embodiment, the beam hardening correction procedure is performed by repeating blocks 708-712 for each slice or tomogram in this reconstruction. At block 708, the process segments the image into different phases. The phases usually comprise different regions of different materials identified in the image. This segmentation can sometimes involve segmenting the entire sample material into a single phase and segmenting the sample container into another phase. As discussed further above, this block sometimes involves grouping the hard-to-distinguish areas as a collective phase with aggregate material attributes across the hard-to-distinguish areas. I will. Such a process, as well as the identification based on that property of a single material performed during the reconstruction / segmentation process described herein, serves to help identify agglomerated materials by their properties. Reference may be made to the agglomeration material library 21 (FIG. 1 (c)). In addition, in some embodiments, agglomeration material can be selected from this library based on operator input. The sample container can be segmented into different phases in various ways, such as a segmentation process that identifies the container material as a phase separate from the sample material and the energy-mediated void. Container material properties can also be stored and referenced in the material library 19 to aid in container identification and segmentation. Although described below with respect to techniques for correcting the container-only beam hardening effect, another suitable process provides container size and shape. Furthermore, in some embodiments, the segmentation block can further include segmenting air between the container and the sample. This will be further described later. Some embodiments, such as the embodiment of FIG. 7 (b), apply container only pre-correction before segmenting the image. However, this is not intended to be limiting, and in some embodiments, the container is considered as part of the same beam hardening correction to account for beam hardening due to the sample. Further, it should be noted that the segmentation performed at block 708 can be performed in conjunction with an initial reconstruction based on any suitable variation of the iterative segmentation and reconstruction process described above. For example, one preferred embodiment comprises: (a) determining one or more estimated material properties of multiple volume elements in the sample using a tomographic reconstruction algorithm; (b) the material Modifying one or more estimated material properties based on the probabilistic classification by reference to stored material property data to determine the most likely material in the dictionary for each volume element And (c) performing one or more additional iterations of the tomographic reconstruction algorithm, each iteration updating one or more estimated material properties of the volume element , (D) labeling the reconstructed volume indicating that it is a sample and classifying it into regions of material having similar properties to generate segments Including the.

  Next, at block 710, the process simulates an RPC referenceless post reconstruction correction, preferably for each identified segment of material, through the simulated path from the meter to the container of each material and to the detector Calculated by integrating radiation attenuation. This calculation can be performed using the various techniques described above for RPC, along with their formulas and equivalents or their appropriate permutations. Next, at block 712, the process applies the calculated post-reconstruction corrections to the tomographic projection data, usually by subtracting a correction amount from the intensity level at the end of each simulated radiation of those projections .

  Next, at block 714, the process reconstructs a series of final tomograms again from the corrected tomographic projection data. In some embodiments, this process can be iterative, where multiple correction steps are calculated and multiple reconstructions are performed. Furthermore, the reconstruction performed after this correction may also be iterative, applying other techniques to generate more accurate reconstruction and segmentation of the data. For example, using any suitable variant of the reconstruction and segmentation process described above with respect to FIG. 1 (a) or 1 (b), in this step with reference to the material library, the fault of block 714 The final segmentation of the imaging reconstruction can be improved. At block 716, the final reconstruction is also segmented, but, as discussed, this block can be combined with the reconstruction or separated in an iterative process. This series of steps is preferably performed, slice by slice, over the entire volume that one wishes to reconstruct to model the sample. In order to properly model radiation attenuation and beam hardening effects, the orientation of the slice is usually determined by the direction of radiation travel and the direction of the original tomographic scan.

Container Precorrection and Subsequent Sample Correction Container Only Correction In addition to techniques that treat the interior of the container as a single phase, or techniques that use an agglomerate phase as discussed above, if only the container is present, ie Considering the correction that would be applied if the same container was filled with air rather than sample is advantageous for the automation process. In order to do this, one must approximate the attenuation that would occur if the process imaged an empty container with the same polychromatic beam. Thus, this process should ignore x-rays passing through the center of the container. These x-rays may be high energy depending on the sample and do not represent an empty container. To that end, the process must have information on the size and shape of the container. Such information may be obtained by performing segmentation on the initial tomogram, or may be obtained by receiving predetermined container shape data describing the container used. Such data can also include attenuation data.

  12 (a)-(c) are schematic diagrams showing the relative path lengths of x-rays having a zero fan angle (d1) and x-rays having a non-zero fan angle (d2). FIG. 12 (a) shows a spherical filter in which all beams pass the same amount of filter. FIG. 12 (b) shows a flat sheet filter where the path length increases as the fan angle increases. FIG. 12 (c) shows a container in which the path length becomes longer as the fan angle increases. If the inner wall of the container is convex, as in the case of a cylinder, all possible path lengths through the container are captured by at least one radiation that does not pass through the interior of the container. Put another way, grazing radiation can probe all possible path lengths. Thus, from this container only data, a complete correction curve for the container material can be constructed.

上式で、Ｃは、容器に起因する減衰、μｃは、同じ放射線セットを通る直線を当てはめることによって近似された容器の減衰係数である。これらの２つの関数の差をとった場合、その投影中の残りの放射線に適用する補正を容器経路長の関数として得ることができる。

Where C is the attenuation due to the container and μ c is the attenuation coefficient of the container approximated by fitting a straight line through the same set of radiation. If the difference between these two functions is taken, the correction applied to the remaining radiation in the projection can be obtained as a function of container path length.

上式で、Ｒ_ｅは、容器でない材料に起因する減衰である。容器は補正されていないことに留意する。

Where R _e is the attenuation due to the non-container material. Note that the container is not calibrated.

To apply this principle, the process segments the container and the interior in the same manner as in the above-mentioned section "Two materials: correction given the sample and the container", which segmentation is shown in FIG. Generate the same point cloud shown. FIG. 10 is a graph of a three-dimensional point cloud used in the correction process. The black data points show the fitted area for the container only correction technique. This technique preferably ignores gray data points representing radiation other than radiation that has only passed through the container, fits a one-dimensional curve to find C _estimated , and fits a straight line through black data points to find μ c C _seg .

  An exemplary result of this correction is shown in FIG. 11 (a). FIGS. 11 (a)-(d) show cross sections and intensity diagrams of the correction applied to the sample of FIG. 9 (a). Again, the contrast is reset to show the effect within the sample. FIG. 11 (a) shows the resulting reconstruction after only beam hardening in the vessel has been removed. FIG. 11 (b) shows the corrected final reconstruction after the segmentation of the air phase has been performed. FIG. 11 (c) shows a gray scale intensity diagram of the original image and the pre-corrected image showing the reverse cupping effect. FIG. 11 (d) shows the final correction compared to the original image. It should be noted that (c) and (d) have been narrowed down to the sample area and thus have not shown a flattening effect on the gray values in the container. However, this effect in the container area is very similar to that shown in FIG. 9 (c).

  There are several points to note about the results of FIGS. 11 (a)-(d). Most important is the complete removal of the bright halo around the inside of the vessel, which allows a very simple segmentation of the air phase surrounding the rock. There are two interesting aspects of the intensity line profile that need to be described. They are the overall loss of strength seen in the scaled magnitude diagram of FIG. 11 (c) and the new back cupping effect in the rock. The overall strength loss is not so great, which is the result of a correction assuming no attenuation in the vessel. Back-capping can be described in terms of the implicit assumption that beam hardening only occurs in the container by ignoring the radiation passing inside the container. In fact, beam hardening also occurs in rock samples, which appear as two artifacts. First, the conventional cupping artefact that results in a darker area centered by changing the amount of passing rock results. Second, there is a reverse cupping artifact where the central region is brighter than the edge. This can be explained by considering the radiation pq and pr in FIG. 11 (a). The inner region of the sample is probed only by radiation like pq and the amount of radiation like pq that passes through the container due to the angle of incidence of those radiation is better than the radiation like pr probing the outer region There are also relatively few. Thus the outer area looks less dense. This is because the outer region is, on average, probed by the higher energy beam. The relative impact of these two beam hardening artifacts depends on the thickness of the container, the attenuation of the container relative to the sample and the fan angle of the incoming beam. In this case, the container has a much higher attenuation compared to the sample, so the anti-Copping effect predominates.

It is worthwhile to show the relationship between this technique and the use of filters during acquisition discussed in Section 2. The attenuation of X-rays is defined as R = ln (I / I ₀ ). Here, I is the intensity of the beam when the sample is present, and I ₀ is the intensity of the beam when the sample is not present, which is known as the clear field. When a filter is used, the direct effect of the filter on the beam intensity is eliminated by this clear field. This is because there is exactly the same filtering for I and I ₀ . However, there is still an indirect effect on the intensity by changing the beam spectrum, ie by increasing the energy of the beam.

  In fan-beam arrangements with flat filtering sheets, larger angles of radiation pass through more material, so there is also a back-capping artifact described for container correction to pre-filtering. The extent to which the cupping and anti-cupping effects affect the sample also depends on the attenuation, thickness and fan angle. The reverse cupping effect on flat sheets is not as dramatic as the reverse cupping effect on container correction. This is because, for flat filters, the radiation path length exhibits less variation. As such, the container correction technique described is similar to applying filtering and using the clear field to remove its effects. The difference is that beam hardening can be corrected if the filtering takes place in the field of view, for example when "filtering" in the container.

  FIG. 7 (b) performs the beam hardening effect correction using the container precorrection filter or precorrection step discussed above applied before the beam hardening effect in the sample is modeled Is a flow chart illustrating one process of doing. This process block is preferably performed automatically by the tomography system 10 under the controller of the system controller 15 after activation of the tomography system 10. The activation of the tomography system 10 is achieved by the operator inputting the appropriate parameters remotely via the user interface 16 or via the network, or a host adapted to control the process via the network Achieved by the controller. However, this is not intended to be limiting, and any particular steps may include operator intervention. This process begins in the same manner as the process of FIG. 7 (a), and at block 732 the process acquires or acquires a series of tomographic data projections and proceeds similar to FIG. 7 (a). Next, at step 734, the process performs an initial reconstruction from those projections. Next, at step 736, the process begins applying a beam hardening correction to each slice or tomogram in this reconstruction. In this embodiment, this correction includes blocks 737-744. Next, at step 737, the process identifies the size and shape of the container. This identification can be performed by a segmentation process on the image or by receiving known data describing the size and shape of the container. This known data may include the material of the container. In some cases, a pre-modeled container, for example with a filter, may already be calculated, the filter comprising beam effects originating from that container or one and the same container of the same material. Such prefilters and characteristics can be stored in the container library 21 (FIG. 1 (c)). The container library 21 has selected by the operator the container to be used, for example by image recognition of the container contour, or by automatically recognizing the bar code or identifier of the container, or when launching an automated process Can be automatically referenced to identify the container in response to After determining the size and shape of the container, the process takes out the correction only for the existing container or calculates the container only pre-correction. If calculated, the calculation is preferably performed by the above-described technique of identifying a series of radiation to determine radiation that passes only the container and not the sample, as described above. The beam hardening effect of the container is modeled based either on data extracted from the tomographic projection or on known data on container material attributes. This correction is calculated by any suitable technique discussed above, or other known methods of modeling such effects. This calculation is performed by using data on the radiation path length through the vessel and applying that data to the radiation passing through the vessel. At this time, it does not matter whether those radiations passing through the container pass through the sample or not. As discussed above, some embodiments fit this data to a curve to determine the attenuation coefficient of the container material. In other embodiments, this data can be applied directly to the projected data as a correction. Next, at step 738, the resulting beam hardening effect is applied to this data. This is preferably done by applying a filter to the tomographic projection and then performing another reconstruction as described further above.

  At block 740, the process segments the image into different phases. The phases usually comprise different regions of different materials identified in the image. This block proceeds similar to block 708 of FIG. 7 (a), except that the container is already taken into account by pre-compensation of the container so that the container can be removed or ignored in steps 740-744. Next, the rock sample is regarded as a single material using the method described in the section “correction assuming a single-material sample”, in which a pre-corrected tomographic image or image (image to which the pre-correction is applied) is used. Can be segmented in a manner that In the preferred embodiment with the results presented herein, this is done automatically by applying k-means clustering of intensities and then excluding all regions connected to previously segmented container regions. To be executed. While this auto-segmentation works well in this case, more sophisticated segmentation routines are easily adopted, if necessary, especially for sample types that are not suitable for k-means clustering or the single-material hypothesis discussed above. can do. In some embodiments, the segmentation process described above with respect to FIG. 1 (a) or (b) uses an appropriate material library based on expected materials in the sample, provided as input to the process. used. As discussed further above, this block is sometimes difficult to distinguish, and attributes of agglomerate material across regions that make it difficult to distinguish regions where segmentation may be very inaccurate into separate phases Including grouping as a collective phase with. Such a process, as well as the identification based on that property of a single material performed during the reconstruction / segmentation process described herein, serves to help identify agglomerated materials by their properties. Reference may be made to the agglomeration material library 21 (FIG. 1 (c)). In some embodiments, the segmentation block can further include segmenting air between the container and the sample. This will be further described later.

  Next, at block 742, the process performs RPC no reference post reconstruction correction to the newly segmented image, preferably from instrument to container of each material, and for each identified segment of material. Calculation by integrating the attenuation of the simulated radiation through the simulated path to the This calculation can be performed using the various techniques described above for RPC, together with their formulas and equivalents, and thus, appropriate substitutions for those techniques.

  Next, at block 744, the process applies the calculated post-reconstruction corrections to the tomographic projection data, usually by subtracting a correction amount from the intensity level at the end of each simulated radiation of those projections . As shown in the intensity line plots of FIG. 11 (b) and FIG. 11 (d), the resulting correction using the preferred equations herein almost completely eliminated beam hardening artefacts . Notice that there is a slight drop in intensity that is visible at the very edge of the intensity line plot. This is described in "correction assuming a single material sample".

  Next, at block 746, the process reconstructs a series of final tomograms again from the corrected tomographic projection data. As discussed above with respect to FIG. 7 (a), this final reconstruction process may be iterative, with multiple correction steps being calculated and multiple reconstructions being performed. Furthermore, the reconstruction performed after this correction may also be iterative, applying other techniques to generate more accurate reconstruction and segmentation of the data. For example, using any suitable variant of the reconstruction and segmentation process described above with respect to FIG. 1 (a) or 1 (b), in this step with reference to the material library, the fault of block 746 The final segmentation of the imaging reconstruction can be improved. At block 748, the final reconstruction is segmented and proceeds similar to block 716 of FIG. 7 (a).

  FIG. 7 (c) is a flow chart of a tomographic reconstruction and segmentation process that considers beam hardening across the spectrum of the x-ray beam used. The process of FIG. 7 (a) and the process of FIG. 7 (b) perform a beam hardening correction and then preferably perform a subsequent analysis assuming two monochromatic beams, which are shown here The process uses beam hardening as part of the physical order model of the reconstruction algorithm, preferably using a reconstruction technique similar to that described above with respect to FIGS. 1 (a) and 1 (b). Consider beam hardening by prediction. The illustrated process begins at block 750 and provides data storage of material properties at block 750. In this embodiment, the material properties include attenuation characteristics across the x-ray spectrum used in the tomographic scan. For example, if two polychromatic x-ray beams are used at different energies, the material property data storage area or library will contain attenuation data for the spectra of both beams for all materials in it. This attenuation data can be stored as an attenuation curve or as a discretized attenuation value. Next, at block 752, the process acquires tomographic scans using polychromatic beams, each at a desired energy level. In the illustrated embodiment, this energy level is a dual energy level. Next, at block 754, the process estimates material properties of multiple volume elements using a tomographic reconstruction algorithm across the beam spectrum to account for beam hardening. This estimation can use any of the suitable reconstruction techniques described above, which are applied across the spectrum of each beam. Thus, the iterative techniques discussed above take beam hardening into account as part of the reconstruction process. The remaining process blocks of FIG. 7 (c) proceed in the same manner as the process blocks of FIG. 1 (b), and the iterative algorithm preferably models each material data storage area containing spectral attenuation data in a forward model. Reference by repetition.

High Density Inclusions The discussion herein considers a procedure involving a fairly uniformly heterogeneous sample that can be treated as one material. However, at times, despite the drawback of changing the relative contrast between the phases, the materials that need to be segmented due to the obscuring effect they have on the rest of the image are samples It may be inside. This is particularly true for highly dampening materials that produce streaking artifacts when imaged, such as high density pyrite inclusions in sandstone. In the inventor's simulation of the present invention, this was simulated using the above technique by including five iron inclusions. This will cause significant streaking effects as shown in FIG. 13 (a) without correction.

  FIGS. 13 (a)-(b) show cross sections of simulated silicon samples in a titanium container containing highly damped iron inclusions. FIG. 13 (a) shows an image of the original scaled to show the effect in the sample. Note the streaking artifacts that are characterized by the dark spots between the bright spots. FIG. 13 (b) shows the final image corrected using three material correction. The streaking artifacts are almost completely removed. The high damping of these inclusions again simplifies their separation from the rock matrix using simple k-means segmentation. This sample is then treated as three materials: a container, a pre-corrected segmented rock material as discussed in the section "container-only correction", and a dense inclusion. By applying this correction method including these segmentations, streaking is almost completely eliminated, as shown in FIG. 13 (b).

Correction Space: Intensity or Decay As mentioned in the section "Referenceless Post Reconstruction Correction (RPC)", this RPC correction is applied in the intensity space (I-space) as implemented in Krumm et al. Not in the damping space (R space). If the sample is perfectly segmented into its component materials, the choice of correction space is not important and both corrections produce the same result. The main difference between these two methods is the effect on the correction that the slope of the hyperplane fit has. In either method, the gradient in each direction determines the average grayscale value for each phase. For R space correction, the contrast within the phase remains constant for any slope. That is, when the pore space is grouped into the same phase as the rock, the average absolute gray scale difference between the rock and the pore remains nearly constant. In the case of I-space correction, the contrast in the phase changes to match the contrast change between the phases. That is, the overall consistency of the contrast is maintained. However, this also means that there is a "correct" slope for the R _mono fit, and the result degrades, especially with regard to contrast changes within the phase, if any other plane is fitted. This means that R-space correction is more robust against inaccurate segmentation, but R-space correction has the disadvantage (or potential advantage) of intra- and inter-phase contrast Comes with. To ensure that there is no change in contrast within the sample phase, in the above several embodiments, where many materials are grouped as a single phase and the relative contrast between the container and the sample is not taken into account. R space correction was more appropriate. However, in other cases, such as when all the material is segmented, the I-space correction potentially has a higher fidelity to the original relative contrast value of the entire image.

For the pre-correction case of correcting only the container, it is preferable to apply the correction in R space. This is because this process ignores the R-values of radiation other than radiation that has only passed through the container when finding R _mono and R _poly , so this process compensates for this remaining radiation When applied, it is not guaranteed that a correction equation is defined. In I-space, the following correction equation 25 for each radiation is similar to equation 22 above.

ある種の放射線に対してこの対数が定義されないことが起こりうる。このことは特に、試料内に高密度包有物があり、そのため、それらの高密度領域を通過した放射線の

It may happen that this logarithm is not defined for certain types of radiation. This is particularly true of the presence of high density inclusions in the sample, so that the radiation passed through those high density areas.

が非常に小さいときに生じた。

Arose when it was very small.

  The above description has shown that the no reference post reconstruction correction technique can be effectively applied to reduce beam hardening artifacts of cylindrical samples in a container. These examples include substantially cylindrical samples, but there is no restriction on the sample shape, and this method should work equally well for non-cylindrical samples in containers. It has been shown that certain materials can be treated as a single phase and still give good artifact reduction. This correction requires accurate segmentation of all other phases, which is easily done by performing an initial pre-correction on the image by correcting only the beam hardening in the container. became. As such, the techniques described above can be used in many different workflows that allow beam correction artifacts in the sample to be automatically corrected.

  The foregoing has outlined rather broadly the features and technical advantages of the present invention in order that the detailed description of the invention that follows may be better understood. The following describes additional features and additional advantages of the present invention. Those of ordinary skill in the art should understand that the disclosed concepts and specific embodiments can be readily utilized as a basis for modifying or designing other structures to achieve the same purpose of the present invention. is there. Furthermore, it is to be understood by those skilled in the art that such equivalent constructions do not depart from the spirit and scope of the present invention as set forth in the appended claims.

  Having described the invention and the advantages of the invention in detail, various changes, substitutions and alterations are made herein without departing from the spirit and scope of the invention as defined by the appended claims. It should be understood that you can. The combinations of features described herein should not be construed as limiting, and the features herein are used in any effective or partial combination according to the present invention. can do. Accordingly, this description should not be construed as providing written support under US patent law and other related patent laws for any effective combination or partial combination of features of the present invention.

  Moreover, the scope of the present application is not intended to be limited to the particular embodiments of the process, machine, manufacture, composition of matter, means, methods and steps described in the specification. As those skilled in the art will readily appreciate from the disclosure of the present invention, existing or future developments performing substantially the same function as the corresponding embodiments described herein, or achieving substantially the same results. Any process, machine, manufacture, composition, means, method or steps may be utilized in accordance with the present invention. Accordingly, the appended claims are intended to include within their scope such processes, machines, manufacture, compositions of matter, means, methods or steps.

Claims (19)

A method of correcting tomographic reconstruction with respect to beam hardening, comprising
(A) reconstructing an initial tomographic image from tomographic projection data;
(B) segmenting the tomogram into two or more phases;
(C) N vectors per pixel based on the segmented tomogram

By determining, without reference post reconstruction compensation (RPC) correct, be calculated by the integrated attenuation space,
(D) applying the RPC correction to tomographic projection data to correct for beam hardening effects; and (e) reconstructing a final tomogram from the corrected tomographic projection data. Method.   Applying the RPC further comprises reprojecting each phase of the segmented tomogram to generate a simulated monochromatic projection for each material, the simulated monochromatic projection calculated in the integrated attenuation space The method of claim 1, wherein A method of correcting tomographic reconstruction with respect to beam hardening, comprising
(A) reconstructing an initial tomographic image from tomographic projection data;
(B) segmenting the tomogram into two or more phases;
(C) calculating a no reference post reconstruction correction (RPC) correction based on the segmented tomograms;
(D) applying the RPC correction to tomographic projection data to correct for beam hardening effects;
(E) reconstructing a final tomographic image from the corrected tomographic projection data
Including
It is further preferable to apply a container-only pre-correction to the tomogram in order to correct for the beam hardening effect due to the container holding the sample material and not for the beam hardening effect due to the sample itself How to include it .   The method further includes generating a pre-correction for only the container by selecting a set of radiation passing only the container and calculating an integral attenuation for the container based on the selected set of radiation. Item 3. The method according to Item 3.   4. The method of claim 3, further comprising receiving predetermined container shape data and generating only the container pre-correction based on the predetermined container shape data.   The method according to claim 3, wherein segmenting the tomogram provides container shape data used for the purpose of generating the container only precorrection. A method of correcting tomographic reconstruction with respect to beam hardening, comprising
(A) reconstructing an initial tomographic image from tomographic projection data;
(B) segmenting the tomogram into two or more phases;
(C) calculating a no reference post reconstruction correction (RPC) correction based on the segmented tomograms;
(D) applying the RPC correction to tomographic projection data to correct for beam hardening effects;
(E) reconstructing a final tomographic image from the corrected tomographic projection data
Including
The method comprising the segmentation executes distinguish k-means clustering technique the phase using a radial distance and gray scale intensity. A method of correcting tomographic reconstruction with respect to beam hardening, comprising
(A) reconstructing an initial tomographic image from tomographic projection data;
(B) segmenting the tomogram into two or more phases;
(C) calculating a no reference post reconstruction correction (RPC) correction based on the segmented tomograms;
(D) applying the RPC correction to tomographic projection data to correct for beam hardening effects;
(E) reconstructing a final tomographic image from the corrected tomographic projection data
Including
The segmentation forms part of an iterative segmentation and reconstruction process for the initial reconstruction, the iterative segmentation and reconstruction process comprising
( I ) determining one or more estimated material properties of multiple volume elements in the sample using a tomographic reconstruction algorithm;
( Ii ) probabilistic classification into the one or more estimated material properties by reference to stored material property data to determine the most likely material in the material dictionary for each volume element Correcting based on
( Iii ) performing one or more additional iterations of the tomographic reconstruction algorithm, each iteration updating the one or more estimated material properties of the volume element When,
(Iv) a method comprising the steps of the label a reconstructed volume indicates that it is a sample, classifies the area of materials with similar properties. A method of correcting tomographic reconstruction with respect to beam hardening, comprising
(A) reconstructing an initial tomographic image from tomographic projection data;
(B) segmenting the tomogram into two or more phases;
(C) calculating a no reference post reconstruction correction (RPC) correction based on the segmented tomograms;
(D) applying the RPC correction to tomographic projection data to correct for beam hardening effects;
(E) reconstructing a final tomographic image from the corrected tomographic projection data
Including
Reconstructing the final tomogram,
( I ) determining one or more estimated material properties of multiple volume elements in the sample using a tomographic reconstruction algorithm;
( Ii ) probabilistic classification into the one or more estimated material properties by reference to stored material property data to determine the most likely material in the material dictionary for each volume element To fix based on
( Iii ) performing one or more additional iterations of the tomographic reconstruction algorithm, each iteration updating the one or more estimated material properties of the volume element; iv ) A method further comprising labeling the reconstructed volume indicating that it is the sample and classifying it into regions of material having similar properties. A method of correcting tomographic reconstruction with respect to beam hardening, comprising
(A) reconstructing an initial tomographic image from tomographic projection data;
(B) A container-only pre-correction is applied to the tomogram in order to correct for beam hardening effects due to the container holding the sample material and not for beam hardening effects due to the sample material itself about,
(C) segmenting the tomogram into two or more phases;
(D) calculating a no reference post reconstruction correction (RPC) correction based on the segmented tomograms;
(E) applying the RPC correction to tomographic projection data to correct for beam hardening effects; and (f) reconstructing a final tomogram from the corrected tomographic projection data. Method.   Applying the RPC further comprises reprojecting each phase of the segmented tomogram to generate a simulated monochromatic projection for each material, the simulated monochromatic projection calculated in the integrated attenuation space The method according to claim 9, which is   The method further includes generating a pre-correction for only the container by selecting a set of radiation passing only the container and calculating an integral attenuation for the container based on the selected set of radiation. 10. A method according to item 9.   10. The method of claim 9, further comprising receiving predetermined container shape data and generating only the container pre-correction based on the predetermined container shape data.   10. The method of claim 9, wherein segmenting the tomogram provides container shape data used for the purpose of generating the container only pre-correction.   10. The method of claim 9, wherein the segmentation comprises performing a k-means clustering technique that distinguishes the phases using radial distance and grayscale intensity. Reconstructing the final tomogram,
(A) determining one or more estimated material properties of multiple volume elements in the sample using a tomographic reconstruction algorithm;
(B) In order to determine the most likely material in the material dictionary for each volume element, the one or more estimated material properties may be probabilistically classified by referring to the stored material property data To fix based on
(C) performing one or more additional iterations of the tomographic reconstruction algorithm, each iteration updating the one or more estimated material properties of the volume element; 10. The method of claim 9, further comprising: d) labeling the reconstructed volume indicative of the sample and classifying into regions of material having similar properties.   10. The method of claim 9, wherein segmenting the tomogram into two or more phases further comprises separating into at least one phase of agglomerating material.   17. The method of claim 16, further comprising identifying the mass material with reference to a mass material library. An x-ray beam source guided through the sample;
A rotating sample holder,
An x-ray detector for detecting energy from the x-ray beam source;
A computer controller for performing the method according to any one of claims 1 to 8, which is the X-ray beam source, the rotating sample holder and a system controller operatively connected to the X-ray detector. And a system controller including a program storage device.

Images