WO2017169232A1 - Image reconstruction method for interior ct - Google Patents

Abstract

[Problem] To provide a thorough image reconstruction method for interior CT which uses a priori knowledge which is more practical. [Solution] In this image reconstruction for interior CT, projection data is obtained from particle beams passing through an ROI within a subject being imaged, approximate reconstruction (first step) is performed by a CT image reconstruction method using the obtained projection data, a region where an image value representing a physical quantity is at least piecewise constant or represented at least piecewise by a polynomial equation is specified in the ROI on the basis of the reconstructed CT image, and reconstruction (second step) which is more precise than the first reconstruction step is performed using the position of the specified region where the image value representing a physical quantity is at least piecewise constant or represented at least piecewise by a polynomial equation, and the feature of the physical quantity therein which is at least piecewise constant or represented at least piecewise by a polynomial equation.

Description

Image reconstruction method for interior CT

The present invention relates to an image reconstruction method for measuring a line integral value of a physical quantity distribution inside an object and generating an image of the physical quantity distribution by data processing, and more particularly to an image reconstruction method for interior CT.

First, a CT method called interior CT (computer tomography) will be described. In general, in many situations of CT imaging, an image of only a small region of interest (ROI: Region of interest) in an object (sample) may be desired. For example, in the diagnosis of heart disease and breast cancer using medical CT, it is sufficient to have only a small ROI image including the heart and breast. Even if the configuration method and the data acquisition method of the current CT apparatus are sufficient in the case where such an ROI-only image is sufficient, an X-ray beam that completely covers the cross-section including the ROI is irradiated. Instead, the projection data on all straight lines passing through the cross section of the object are measured (see FIG. 1A). This is because, in the filtered back projection (FBP) method, which is a calculation procedure used for image reconstruction of the CT apparatus, projection data on a straight line that does not pass through the ROI is also required to generate the ROI image. It is. However, intuitively, since projection data on a straight line that does not pass through the ROI does not include any ROI information, it is expected to be unnecessary. Therefore, the interior CT is a method of CT imaging in which only the ROI is irradiated with X-rays and only the projection data on (all) straight lines passing through the ROI are measured to generate only the ROI image (FIG. 1). (See (b)).

This interior CT has various advantages over conventional CT that measures unnecessary projection data wastefully. For example, (1) a significant reduction in exposure outside the ROI (sample damage), (2) a reduction in detector size and X-ray beam width, and (3) a large object that does not fit in the field of view can be taken. (4) High-resolution CT imaging can be performed in which only a small field of view of an object is irradiated with X-rays for enlarged imaging.

On the other hand, in the interior CT, projection data on a straight line that does not pass through the ROI is not measured. Therefore, a method of performing image reconstruction from incomplete projection data partially missing is necessary. More precisely, the definition of the above-described image reconstruction problem in the interior CT will be described as follows. That is, consider the object f (x, y) and the ROI S to be imaged (see FIG. 2A). Then, it is assumed that only the projection data p (r, θ) (r is the moving radius and θ is the angle) in which the straight line passes through ROI S can be measured. However, for the sake of simplicity, it is assumed that projection data is collected by a parallel beam. In this case, since p (r, θ) where the straight line does not pass through ROI S is not measured, the projection data at each angle θ is truncated and lost. The problem of correctly reconstructing the image f (x, y) in ROI S from such truncated projection data is the so-called interior CT image reconstruction.

This problem has been studied for many years and is outlined below. First, in Non-Patent Document 1, Natterer mathematically proves that the image reconstruction of the interior CT does not determine the solution to be “unique” (here, unique means that the image reconstruction solution is derived from the projection data). Because this non-uniqueness is known, it has been studied many approximate image reconstruction methods.

For example, as a typical technique, (1) a method of reconstructing an image after extrapolating the left and right missing portions of each direction projection data with a smooth function, and (2) an image by successive approximation with incomplete projection data. Although methods for performing reconstruction have been studied, artifacts due to approximation errors have occurred and have not been put into practical use (Non-Patent Documents 2, 3, and 4). Examples of typical artifacts that occur in the interior CT are shading artifacts (see FIG. 3A) in which distortion occurs in the low frequency components of the image, and cupping effects in which values increase in the periphery of the ROI (FIG. 3). (See (b)) occurs, and it is known that the value of the image is not stably determined.

In contrast to these previous studies, exact image reconstruction methods for interior CT were discovered independently of each other (Non-Patent Documents 5 and 6). In these papers, “any arbitrary region B within ROI S (ie, If there is a priori knowledge that the value of the image f (x, y) is known in advance in the area of the circle inside ROI S in FIG. 2A, the image reconstruction of the interior CT We proved that the solution of is uniquely determined. Surprisingly, the region B for guaranteeing the uniqueness of the solution can be anywhere in the ROI S (but not at one point). This result is already known as US Pat.

On the other hand, Yu et al. Discovered a technique that enables strict image reconstruction using another a priori knowledge (Non-Patent Document 7). In this paper, a lack of measurement data called compressed sensing is disclosed. Based on the method of signal restoration with high accuracy from the above, it is proved that if the image f (x, y) is piecewise uniform throughout the ROI S, the image reconstruction solution of the interior CT is uniquely determined did. However, the piecewise uniform means that the image is composed of a finite number of regions having a completely constant value like a numerical phantom (see FIG. 2B). This result is already known as a US patent (Patent Document 2).

US Pat. No. 7,697,658 US Pat. No. 8,811,700

Natterer F: The Mathematics of Computerized Tomography. Wiley, 1986 Ogawa K, Nakajima M, Yuta S: A reconstruction algorithm from truncated projections. IEEE Transactions on Medical Imaging 3: 34-40, 1984 Ohnesorge B, Flohr T, Schwarz K, Heiken JP, Bae KT: Efficient correction for CT image artifacts caused by objects extending outside the scan field of view. Medical Physics 27: 39-46, 2000 Ohyama N, Shiraishi A, Honda T, Tsujiuchi J: Analysis and improvement in region-of-interest tomography. Applied Optics 23: 4105-4110, 1984 Ye Y, Yu H, Wei Y, Wang G: A general local reconstruction approach based on a truncated Hilbert transform. International Journal of Biomedical Imaging 2007: Article ID 63634, 2007 Kudo H, Courdurier M, Noo F, Defrise M: Tiny a priori knowledge solves the interior problem in computed tomography. Physics in Medicine and Biology 53: 2207-2231, 2008 Yu H, Wang G: Compressed sensing based interior tomography. Physics in Medicine and Biology 54: 2791-2805, 2009 Yang J, Yu H, Jiang M, Wang G: High order total variation minimization for interior tomography. Inverse Problems, 26: Article ID 35013, 2010 Courdurier M, Noo F, Defrise M, Kudo H: Solving the interior problem of computed tomography using a priori knowledge. Inverse Problems 24: Paper No. 065001, 2008

However, the conventional techniques described above still have the following problems.

First, in order to use the exact solution known from Non-Patent Document 5 and Non-Patent Document 6 described above, a priori knowledge about an object (an image of an image in an arbitrary small region B inside ROI S before imaging). Value) must be known. However, it is rare that the value of the image is known before shooting. In the exact solution of Non-Patent Document 7 described above, it is necessary to assume that the image is piecewise uniform over the entire ROI S, but this is not possible for a general CT image having a smooth density change. In this method, there is a risk that a smooth density change may be lost (an example is shown in FIG. 9 later).

Accordingly, the present invention has been developed from a strict interior CT image reconstruction method using the above-mentioned a priori knowledge, and in the points described below, a more practical interior CT image reconstruction method. The purpose is to provide.

In order to achieve the above object, according to the present invention, first, there is provided an image reconstruction method for interior CT, in which projection data is acquired by a quantum beam passing through an ROI inside an imaging target, and the projection obtained above is obtained. First, approximate reconstruction is performed by CT image reconstruction using data, and an image numerical value representing a physical quantity in the ROI based on the CT image reconstructed above is at least piecewise uniform or A region that is represented by a piecewise polynomial is specified, and the physical quantity is at least piecewise uniform or a piecewise polynomial expression that indicates the position of the specified region and the physical quantity is piecewise uniform inside or An interior CT image reconstruction method is provided that performs second-stage reconstruction with higher accuracy than the first-stage reconstruction using the property expressed in a piecewise polynomial form.

In the present invention, in the interior CT image reconstruction method described above, the numerical value of the projection data representing the physical quantity may include absorption of the quantum beam by the imaging target, or the imaging CT The phase shift of the quantum beam may be included, the diffraction of the quantum beam by the imaging target may be included, or the scattering of the quantum beam by the imaging target may be included.

In the interior CT image reconstruction method, the numerical value of the projection data including the phase shift, the phase shift of the quantum beam, diffraction, or diffraction is added to the detector by adding an optical element or changing its position. It is also possible to reconstruct an image using a numerical value of the projection data extracted from a set of acquired intensity data of the plurality of quantum beams and including the phase shift, diffraction, or diffraction of the extracted quantum beams. is there.

In addition, in the present invention, in the image reconstruction method for interior CT described above, the region in which the physical quantity specified in at least the ROI is piecewise uniform or piecewise represented by a polynomial is The CT image value reconstructed by the approximate reconstruction in the first stage may be a region in which the values are piecewise uniform or piecewise represented by a polynomial, or the CT reconstructed in the first step. The at least piecewise uniform or piecewise polynomial region identified within the ROI of an image may be selectable within the ROI.

Further, according to the present invention, in the image reconstruction method for interior CT described above, the at least piecewise uniform or piecewise polynomial area specified in the ROI is an approximation of the first step. May be set manually by a human using a CT image obtained by simple reconstruction, or set by image processing using a CT image obtained by approximate reconstruction in the first stage. May be. Alternatively, it may be set in advance in the ROI by at least one of specifying from the previously acquired CT image of the imaging target, a model representing the structure of the imaging target, and specification from a priori information. Note that the at least piecewise uniform or piecewise polynomial area specified in the ROI using the CT image obtained by the approximate reconstruction in the first stage is the boundary of the ROI. A part may be formed.

In the interior CT image reconstruction method, the first-stage approximate reconstruction includes conventional CT images including a filtered back projection (FBP) method, a successive approximation method, and a statistical reconstruction method. The second stage reconstruction may be performed by a differential backprojection Hilbert transform method, a constrained successive approximation method, and a constrained statistical reconstruction method. You may perform by at least 1 of the image reconstruction method of CT containing.

According to the present invention described above, it is possible to provide a more practical interior CT image reconstruction method capable of strict image reconstruction with much less a priori knowledge than the prior art.

It is a figure explaining interior CT in which this invention relates by comparison with normal CT. It is a figure which shows the example in which the conventional exact reconstruction becomes possible, and the conventional exact reconstruction becomes possible. It is a figure which shows the example of the typical artifact which generate | occur | produces in interior CT. It is a figure explaining the piecewise uniform and piecewise polynomial in the present invention. FIG. 11 is a diagram showing an actual reconstruction example when a priori knowledge of [Result 1] to [Result 3] is used. It is a flowchart figure which shows the detail of the two step image reconstruction method of this invention. It is a figure which shows an example of the incomplete image and exact image in the two-step image reconstruction method of this invention. It is a figure which shows the specific example which sets a priori information area | region B to a part of ROI S in the a priori knowledge non-identification type image reconstruction method (Example 3). It is a figure which shows the specific example of a simulation experiment of the a priori knowledge identification type image reconstruction method (Example 1) and the a priori knowledge non-identification type image reconstruction method (Example 3). It is a figure explaining the method which combined DBP method and truncation Hilbert transform used by this invention. It is a figure which shows the external appearance structure of the general X-ray CT apparatus which is an example of the apparatus using the image reconstruction method of this invention. It is a figure which shows an example of the internal structure of the said X-ray CT apparatus.

Hereinafter, prior to describing embodiments of the present invention in detail, in the present invention, in order to solve the above-described problems and to construct a more practical interior CT image reconstruction method, New theories that enable strict image reconstruction with much less a priori knowledge than the prior art papers of Ye et al., Kudo et al. (Non-Patent Documents 5 and 6), and Wang et al. Built. Based on this, a method of identifying from the measured projection data (a priori knowledge identification-type image reconstruction method), and an identification so that it is not necessary to give the a priori knowledge in advance. 2 (a priori knowledge non-identification type image reconstruction method) is proposed, which is fixed to the periphery of the image (approx. Approximately any image).

That is, the key elements of the present invention are (1) exact reconstruction theory using much less a priori knowledge, and (2) a priori knowledge identification type image reconstruction that identifies a priori knowledge from projection data. (3) An a priori knowledge non-identification type image reconstruction method that fixes a priori knowledge without identifying it. A newly proposed image reconstruction method for this purpose will be described in detail below.

As shown in FIG. 1 (b), the interior CT is compared with normal CT (see FIG. 1 (a)) that irradiates the entire object that is the subject with X-rays from an X-ray source, In this technique, a part of a large object (region of interest (ROI)) is enlarged and imaged by irradiating only the region of interest (ROI) of the inside of the object with X-rays. According to this, it is possible to reduce the exposure dose outside the region of interest (ROI) to reduce the damage caused by X-rays, and to reduce the beam width of the irradiated X-rays and simultaneously reduce the size of the X-ray detector. In particular, it is effective for X-ray CT for enlarging a small field of view, and further for micro CT, electron beam CT, and the like.

<Processing Flow of Image Reconstruction Method of the Present Invention>

Example 1 is an “a priori knowledge identification type image reconstruction method” in which a priori knowledge necessary for strict image reconstruction is automatically identified from projection data and used for strict image reconstruction. This will be described below with reference to FIGS. 6 and 7A and 7B.

FIG. 6 shows a processing flow of the image reconstruction method which is a method proposed in the present invention. This proposed method is a so-called two-stage image reconstruction method including a first step and a second step. In the first first step (S61), an incomplete image including artifacts is generated using a conventional FBP method, a successive approximation image reconstruction method, a statistical image reconstruction method, or the like without a priori knowledge. This incomplete image includes artifacts, but since the artifacts generated in the interior CT are low frequency components, information on the boundaries of structures such as organs and tissues is accurately reflected in most cases.

Therefore, from the above incomplete image shown in FIG. 7A, for example, the result of [Result 1] to [Result 4] will be described later by the user manually or by image analysis (processing) using software. A priori information area B which is an arbitrary small area in ROI S that can be used as experimental knowledge is identified (S62). Then, in the second step, the prior information area B obtained in the first step is used for a priori knowledge, and as shown in FIG. Image reconstruction is performed with higher accuracy than the one-stage reconstruction (S63). Note that which one of [Result 1] to [Result 4] is used is determined depending on what a priori information area B can be extracted in the first step. For example, [Result 1] is definitely a constant B, and [Result 3] is definitely a piecewise uniform B. Although it is not a constant, B is close to a change in the density of a polynomial. [Result 2] is selected, and [Result 4] is selected and used if B is close to the density change of the piecewise polynomial.

Here, the results obtained from the above incomplete images are summarized in the following [Result 1] to [Result 4].

<Result 1 (constant value a priori knowledge)>
As shown in FIG. 2 (a), if an arbitrary small a priori information area B exists inside ROIS and it is known that f (x, y) is a constant value C (constant) in B The image reconstruction solution for the interior CT is uniquely determined. However, the constant value C may be unknown in advance, and as a result, the a priori knowledge of the exact solution of Ye et al., Kudo et al. (Non-Patent Documents 5 and 6) is reduced.

<Result 2 (polynomial a priori knowledge)>
As shown in FIG. 2 (a), it is known that there is an arbitrarily small a priori information area B inside ROIS, and f (x, y) is an M-order polynomial in B. For example, the image reconstruction solution for the interior CT is uniquely determined. Here, the polynomial means that the density change f (x, y) of the image has the following form.

However, the degree M of the polynomial needs to be known, and the coefficient a _mn of the polynomial may be unknown. [Result 1] is a result of the restriction on the form of the function in [Result 2] with the order of the polynomial set to M = 0. That is, [Result 2] is a priori knowledge of [Result 1]. It has become less.

<Result 3 (Categorical uniform a priori knowledge)>
As shown in FIG. 2 (a), it is known that an arbitrary small a priori information area B exists inside ROI S, and f (x, y) is known piecewise constant in B. If so, the image reconstruction solution for the interior CT is uniquely determined. Here, the piecewise uniform means that, as shown in FIG. 4, B is composed of a finite number (L) of regions D ₁ , D ₂ ,..., D _L and each region has a constant value C ₁ , C _2. ,…, C _L. However, the number of regions L and the values of constant values C ₁ , C ₂ ,..., C _L may be unknown in advance, in other words, [Result 3] is obtained by reducing the a priori knowledge of [Result 1]. It has become.

<Result 4 (a piecewise polynomial a priori knowledge)>
As shown in FIG. 2 (a), an arbitrary small a priori information region B exists inside ROI S, and f (x, y) is an M-order piecewise polynomial in B. If known, the image reconstruction solution for the interior CT is uniquely determined. Here, the piecewise polynomial, region D ₁ of the B is a finite number (L number) as shown in FIG. 4, D _2, ..., density change in the image of the l-th region consists D _L f _l ( x, y) has the following form:

However, the degree M of the polynomial needs to be known, the number of regions L and the polynomial coefficient a _mn ^(l) may be unknown, and [Result 4] is a priori knowledge of [Result 2] and [Result 3]. Has become less.

In this way, a priori knowledge about an object necessary for performing strict image reconstruction is theoretically considered and described in the above prior arts (Non-Patent Documents 5 and 6, Patent Document 1). It can be seen that strict image reconstruction is possible with much less a priori knowledge than a priori knowledge.

Here, the precautions when using [Result 2] to [Result 4] above are described. In the simulation experiment of the inventors, the a priori knowledge of [Result 1] and [Result 3] operates relatively stably, but the a priori knowledge of [Result 2] and [Result 4] is the degree M of the polynomial. When is large, the numerical instability tends to increase and the error tends to increase, and it has been substantially successful until M ≦ 3. This is because, as the degree M of the polynomial increases, the restriction on the image f (x, y) decreases and the a priori knowledge decreases, and this is a phenomenon that is naturally expected.

In the above Non-Patent Documents 5 and 6, the same a priori knowledge of an arbitrary small region B in ROI S as [Result 1] to [Result 4] is used, but the image f (x, y) Whereas the value itself is required, in the present invention, in [Result 1] to [Result 4], the value of f (x, y) is much smaller, such as (piecewise) uniform or (piecewise) polynomial. It differs greatly only in a priori knowledge. In Non-Patent Document 7, a priori knowledge of piecewise uniform type is used, but an unreasonable assumption is necessary that piecewise uniform throughout ROI S, not an arbitrary small region B in ROI S. It is different in point.

Also, FIGS. 5A to 5C show actual reconfiguration examples when the a priori knowledge of [Result 1] to [Result 3] is used. In any case, it can be seen that significant artifact reduction can be achieved with a little a priori knowledge.

Example 2 is an “a priori knowledge automatic estimation type image reconstruction method” for identifying a priori knowledge from projection data.

According to the theory described in the first embodiment, a priori knowledge necessary for strict image reconstruction of the interior CT can be much lessened. However, in actual CT imaging, it is rare if a priori knowledge about an object of interest is known before imaging, for example, a CT image previously captured of the same patient or an image captured with another modality. Limited to very few special cases, such as

Therefore, in the second embodiment, in identifying the a priori information area B that is an arbitrary small area in ROI S that can be used as a priori knowledge in step S62 in the process flow shown in FIG. It identifies automatically without going through. It will be apparent to those skilled in the art that the automatic identification of the a priori information area can be easily realized by using, for example, an image analysis (processing) technique.

The success or failure of the a priori knowledge automatic estimation type image reconstruction method described in the second embodiment depends on whether or not the a priori knowledge (a priori information area B) can be identified successfully in the first step. In addition, since the step of identifying the a priori information area is somewhat complicated, there may be a desire not to perform the identifying step.

Therefore, in the third embodiment, “a priori knowledge non-identification type image reconstruction method”, which is a practical technique for performing image reconstruction by fixing the prior information area B without identifying it, is proposed. Note that the method of the third embodiment is based on the following three findings regarding interior CT image reconstruction.

[Knowledge 1]
8A to 8C, even if some errors and artifacts are generated around the edge of ROI S, if the center portion of the image can be reconstructed well, there will be no significant problem in application.

[Knowledge 2]
In the image reconstruction of the interior CT, when the a priori information area B is set on both sides of the ROIS, the image reconstruction can be stably performed numerically, and artifacts and noise are reduced. For example, if the a priori information area B is arranged only at the center, an error increases around the ROI S or noise increases.

[Knowledge 3]
Assuming that the density change f (x, y) of the image around ROI S is piecewise uniform (a priori knowledge of [Result 3]) or piecewise polynomial (a priori knowledge of [Result 4]) In many cases, this is reasonable and generally applicable as the first approximation.

Based on the above knowledge, in the third embodiment, without identifying the position of the a priori information area B, as shown in FIGS. 8 (a) to 8 (c), the edge portion around ROI S A method of applying a strict image reconstruction method based on [Result 3] or [Result 4] is proposed. Of course, the assumption that the edge part around ROI S is piecewise uniform or piecewise polynomial is not strictly correct, so the method is only an approximate image reconstruction method, but many CT imaging In this situation, [Knowledge 3] holds, so that image reconstruction can be performed much better compared to other approximate image reconstruction methods that do not use a priori knowledge (studied before the discovery of the exact solution). Be expected.

Next, FIGS. 9A to 9F show specific examples of the a priori knowledge identification type image reconstruction method (Example 1) and the a priori knowledge non-identification type image reconstruction method (Example 3). An example of a simple simulation experiment is shown by comparison with the prior art. The chest CT image was used for the experiment, and the image was reconstructed with the heart located at the center as ROI S (see FIG. 9A).

In the first embodiment, the user looks at the reconstructed image by the FBP method in the first step, designates the a priori information area B in ROI S, and performs the strict image reconstruction in the second step (note that The a priori knowledge was piecewise uniform with B corresponding to [Result 3]). The result is shown in FIG. Further, in Example 3, image reconstruction was performed with the a priori information area B fixed to the periphery of the ROI （S (the frame type in FIG. 8C). 3] is piecewise uniform with B corresponding to). The result is shown in FIG.

As a comparative example, the result of the local FBP method in which the FBP method is applied by extrapolating the missing portion of the projection data with a smooth function is shown in FIG. 9B, and the image f (x, y in the same a priori information region B is shown. ) Is a piecewise uniform result for the entire ROI S, not just the small a priori information area B, in FIG. FIG. 9D shows the results obtained by the compression sensing method (Non-Patent Document 7) that applies a total variation (TV: Total Variation), which is a constraint of the above.

As is clear from these results, the local FBP method has a strong cupping effect and image degradation is significant, and the compression sensing method has lost considerable details and smooth density changes due to the influence of TV. On the other hand, it can be seen that in the methods of the first and third embodiments of the present invention, the image can be reconstructed quite well by reducing artifacts.

<Image reconstruction method>
Next, an image reconstruction method for generating an image from projection data based on the uniqueness of the solutions of [Result 1] to [Result 4] described above will be described. Of course, if the uniqueness of the solution is proved by using some a priori knowledge, any image reconstruction method that uses a priori knowledge as a constraint condition can generate an image precisely. There are countless options for image reconstruction methods. Specifically, an image reconstruction method can be constructed by the following procedure. First, vectors obtained by discretizing the image f (x, y) and the projection data p (r, θ) are represented by x and b, respectively, and a projection calculation matrix that associates the projection data with the image is represented by A. However, the image x includes not only the pixels in the ROIS but also all the pixels belonging to the object existence region in the cross section (caution is required), and the projection data vector b is created by arranging all the measured values in a line. Further, an evaluation function for evaluating whether a priori knowledge is satisfied in the a priori information area B is represented by F (x). At this time, the image reconstruction can be formulated as one of the following three optimization problems.

[Formulation 1]

[Formulation 2]

[Formulation 3]

However, x∈C represents a constraint condition that can be known in advance with respect to an image, and the following is often used.
(A) (Support Constraint) The image x becomes zero outside the support region Ω _OBJ that is known in advance.
(B) (Non-negative condition) The component of the image x does not take a negative value.
(C) (Projected data value on the Hilbert straight line) In the image reconstruction method using the Hilbert transform described later, a Hilbert straight line L (u described later) is used to compensate for the loss of information when rewriting Ax = b to Hx = c. The projection data values above are used.

If F (x) is a function that does not have a local optimal solution called local function (local minimum), there are many iterative or non-iterative solutions to solve the above problem in the mathematical optimization and image reconstruction fields. All of these techniques are available. For example, a statistical image reconstruction method with a constraint condition or a successive approximation method with a constraint condition can be used. As another class of image reconstruction methods, an image reconstruction method can be constructed based on a later-described framework called differential back projection (DBP). Although details of the DBP method will be described later, in the DBP method, the relational expression between the image x represented by Ax = b and the projection data b is not used as it is, but the image x and the Hilbert image (Hilbert image) are temporarily used by a technique called DBP. ) Is converted into a relational expression Hx = c of an incomplete image c, and is formulated as follows to perform image reconstruction. This class of methods is called by names such as differential backprojection + truncation Hilbert transform (Non-Patent Documents 5 and 6, Patent Document 1).

[Formulation 4]

[Formulation 5]

[Formulation 6]

Of course, the problem formulated as described above is solved using an iterative or non-iterative solution.

Next, the evaluation function F (x) for evaluating a priori knowledge in the a priori information area B, which is a key for performing strict or accurate image reconstruction in the present invention, will be described.

First, the idea for designing F (x) is as follows. First, unlike the methods of Non-Patent Document 7 and Non-Patent Document 8 in which a constraint condition based on a priori knowledge is applied to the entire ROI S, the image reconstruction method of the present invention is an a priori that is an arbitrary small region in S. Constraint conditions are imposed only on information area B. In the case of [Result 1], the norm of the first derivative of f (x, y) in B is minimized because the first derivative of f (x, y) is zero in B. Or the variation (dispersion) in density change in B may be minimized. In the case of [Result 2], since the M + 1 derivative of f (x, y) is a priori knowledge that becomes zero in B, the norm of the M + 1 derivative is minimized, or What is necessary is to minimize the error between f (x, y) and f (x, y) ((x, y) ∈B) and a function in which an M-th order polynomial is applied. In the case of [Result 3], f (x, y) is a priori knowledge that is piecewise uniform within B, and therefore, a total variation based on the L ₀ norm or L ₁ norm in B (TV: Total Variation) The norm should be minimized. In the case of [Result 4], since the M-th derivative of f (x, y) is a priori knowledge that is piecewise uniform in B, the TV defined based on the M + 1-th derivative It is sufficient to minimize the norm.

Table 1 shows typical F (x) examples, though not all can be given because there are various options. However, in Table 1, the parameter p is the norm order, and its value may be 0 ≦ p ≦ 2 in [Result 1] and [Result 2], but is uniform or M in [Result 3] and [Result 4]. It is necessary to use 0 ≦ p ≦ 1 in order to avoid that the influence of the boundary of a finite number (L) of partial areas having density changes of the following polynomial is evaluated too large. However, when 0 ≦ p <1, F (x) does not become a convex function, so it is necessary to devise an iterative solution or a non-iterative solution, and it can be said that p = 1 should be used. In addition, the F (x) shown in Table 1 with multiple candidates was tested by numerical experiments, but all of them worked well and there was no significant difference.

In the following, a new interior CT image reconstruction that enables strict ROI S image reconstruction with much less a priori knowledge than Non-Patent Document 5, Non-Patent Document 6, and Patent Document 1 that form the basis of the present invention. Explain the basics of construction theory. Specifically, far less a priori knowledge means that f (x, y) is uniform in any small region B inside ROI S described in [Result 1] to [Result 4] ([Result 1 ]), Polynomial ([result 2]), piecewise uniform ([result 3]), piecewise polynomial ([result 4]). Obviously, [Result 4] includes [Result 1] to [Result 3] as special cases, so in order to save space, the result here is strictly in the case of [Result 4] (piecewise polynomial a priori knowledge). We prove that accurate image reconstruction is possible.

(1) Basics First, the problem setting and terms / symbols used for the proof are defined, and the differential back projection (DBP) method and truncation Hilbert transform, which are the image reconstruction methods of previous research used for the proof. The combined method will be described (see Non-Patent Documents 5 and 6). As shown in FIG. 10A, an object f (x, y) and ROI S to be imaged are considered. Assume a situation of an interior CT in which only parallel beam projection data p (r, θ) (r: radius, θ: angle) that passes through a straight line passes through S. In this case, since data that does not pass through S is not measured, p (r, θ) is truncated, and reconstructing f (x, y) at S from such incomplete projection data is the image reconstruction of interior CT. It is a configuration. As mentioned above, Ye and Kudo have a priori knowledge that, independently of each other, “f (x, y) is known in an arbitrary small region B (which can be as small as possible) in S” The solution of interior CT image reconstruction has been shown to be unique (see Non-Patent Documents 5 and 6). In these papers, the DBP method is used as a mathematical framework for showing the uniqueness of the solution, and the proof of [Result 4] is also performed using the DBP method. In the DBP method, ROI S is first decomposed into a set of straight lines L (u); u∈U (u is a parameter representing a straight line) called a Hilbert line. Then, the DBP method is used to perform image reconstruction by reducing the image reconstruction to an inverse transformation of an integral transformation called a Hilbert transform for each Hilbert straight line L (u). However, the set of Hilbert straight lines L (u); u∈U is selected so as to satisfy the following two conditions.

(a) Each point (x, y) of S belongs to at least one L (u).
(b) All L (u) intersect with the a priori information area B in which the value of f (x, y) is known in advance.

10 (b) and 10 (c) show how to take a typical Hilbert straight line L (u); u∈U. After decomposing S into Hilbert straight lines in this way, image reconstruction is performed for each L (u) by the following processing procedure.

[First Step] (DBP)
First, for all points (x, y) where L (u) intersects with S, DBP of the following equation is calculated to obtain an intermediate image g _u (x, y) called a Hilbert image.

However, θ (u) represents the angle that L (u) makes with the x-axis. Note that in the region outside L (u) ∩S, angle truncation of projection data occurs due to truncation, and DBP cannot be calculated.

[Second step] (Hilbert inverse transform)
As shown in FIG. 10A, a one-dimensional coordinate t is defined on L (u) (the origin may be arbitrary), and an image f (x, y) is converted to a single variable function f (t ) And the Hilbert image g _u (x, y) is represented by a univariate function g (t). At this time, the relationship between f (t) and g (t) is expressed by the following one-dimensional Hilbert transform H (see Non-Patent Documents 5 and 6).

However, pv represents the main value of Cauchy of integration, and points a, b, c, d, e, and f are defined as shown in FIG. That is, the image reconstruction on L (u) is reduced to the inverse transform of the Hilbert transform of equation (10). In the interior CT, the observation interval [b, e] of the Hilbert transform measurement data g (t) is completely included in the support (non-zero) interval [a, f] of the image f (t), and g (t) is The range of [a, b] ∪ [e, f] on both the left and right sides is truncated. The Hilbert transformation with such truncation is called truncation Hilbert transformation, and the image reconstruction of the interior CT is reduced to the inverse transformation of the truncation Hilbert transformation. In Ye et al. (Non-Patent Document 5) and Kudo et al. (Non-Patent Document 6), the inverse transformation is not uniquely determined only by Equation (10), but the interval L (u) ∩B = [c, d ], There is a priori knowledge that f (t) is known, and it proves that f (t) is uniquely determined by L (u) ∩ S by combining this a priori knowledge with Equation (10) . Hereinafter, the section [a, f] corresponding to the object existence area (support) is set to Ω _OBJ , the section [b, e] corresponding to ROI S is set to Ω _ROI , and the section [c, d corresponding to the a priori information area B is set. ] _Is represented as Ω _PRI (see FIG. 10A).

(2) Uniqueness of solutions represented by [Result 1] to [Result 4] Based on the above basic matters, the uniqueness of the solutions of [Result 1] to [Result 3] is included as a special case in the following. Prove the uniqueness of the solution represented by the strongest [Result 4].

First, the most important main results are summarized in the theorem as follows.
[Theorem] If the Hilbert transform g (t) is known in Ω _ROI and there is a priori knowledge that the image f (t) is an Mth order piecewise polynomial in Ω _PRI , f (t) is unique in Ω _OBJ Determined.

This theorem is based on the results of Non-Patent Document 9 by Courdurier et al. And Patent Document 1 by Wang et al., Because the value of image f (t) itself is known as necessary a priori knowledge. It is positioned as showing that it can be reduced by being an Mth order piecewise polynomial. The fact that [Result 4] is derived using this theorem can be easily explained as follows. If f (x, y) is an Mth order (bivariate) piecewise polynomial in ROI S, f (t) is also an Mth order (univariate) piecewise polynomial. Further, g (t) in Ω _ROI is obtained by measurement. Therefore, applying the theorem, f (t) is uniquely determined at every point on the Hilbert straight line L (u). This holds for all Hilbert lines L (u), and the set of Hilbert lines L (u); u∈U is taken to cover ROI S, so f (x, y) is unique in ROI S Determined. One point of caution indicates that f (x, y) is uniquely determined up to the point covered by the Hilbert straight line outside ROI S, but the projection data p (r , θ) has an angular defect (not measured by 180 °), so that the value of f (x, y) is substantially stable and the image can be correctly reconstructed only at a point inside ROIS.

[Result 1] to [Result 3] are positioned as special cases in which the polynomial degree M of the a priori knowledge used in [Result 4] and the piecewise uniform region number L are as follows.
[Result 1] M = 0, L = 1
[Result 2] M is arbitrary, L = 1
[Result 3] M = 0 and L are arbitrary, and the following system is also established.

In [System] [Theorem], when the a priori knowledge in Ω _PRI is changed to the a priori knowledge of [Result 1] to [Result 3], exactly the same solution uniqueness is established.

In the following, we prove the above theorem. First, if both sides of the equation (10) are differentiated M + 1 times with respect to t, the differentiation and the Hilbert transform H are exchangeable in order, so that the following equation is obtained.

Where g ^{(M + 1)} (t) and f ^{(M + 1)} (t) are the M + 1 derivatives of g (t) and f (t), respectively, and g ^{(M + 1)} (t ) (T∈Ω _ROI ) can be calculated from g (t) (t∈Ω _ROI ). Also, from the a priori knowledge of [Result 4], f (t) (t∈Ω _PRI ) is an M-order piecewise polynomial, so f ^{(M + 1)} (t) (t∈Ω _PRI ) is Dirac. It is expressed in the following form using the δ function.

Where δ ^(M) (t) is the M derivative of the delta function, and t _l ∈Ω _PRI (l = 1,2, ---, L-1) is the function of the piecewise polynomial of order M Breakpoint coordinates that change shape. A priori that f ^{(M + 1)} (t) (t∈Ω _PRI ) is expressed in the form of equation (12) and that f (t) (t∈Ω _PRI ) is an M-order piecewise polynomial Knowledge is equivalent. The proof of the theorem is that the measurement data information g ^{(M + 1)} (t) (t∈Ω _ROI ) and f ^{(M + 1)} (t) (t∈Ω _PRI ) are expressed in the form of equation (12). This is based on the fact that f ^{(M + 1)} (t) (−∞ <t <∞) is uniquely determined from a priori knowledge. First, the lemma that shows it is described.

[Lemma 1]
^{f (M + 1) (t} ) (t∈Ω PRI) is represented in the form of Equation ^{(12), f (M +} 1) Hilbert transform Hf ^{(M + 1)} of (t) (t) is the following formula If f is satisfied, f ^{(M + 1)} (t) = 0 (t∈Ω _PRI ).

(Proof)
From Equation (13), the Hilbert inverse transform is expressed as follows.

From Eq. (14), f ^{(M + 1)} (t) is an analytic function at t∈Ω _PRI , and it is obvious that f ^{(M + 1)} (t) = 0 ( Only if t∈Ω _PRI ).

[Lemma 2]
If it is known that f ^{(M + 1)} (t) (t∈Ω _PRI ) is expressed in the form of equation (12), then g ^{(M + 1)} (t) (t∈Ω _ROI ) to f ^{(M + 1)} (t) (−∞ <t <∞) is uniquely determined.

(Proof)
Let f ^{(M + 1)} (t) and g ^{(M + 1)} (t) be the pair of the image M + 1 derivative and its Hilbert transform. From the a priori knowledge that g ^{(M + 1)} (t) (t∈Ω _ROI ) and f ^{(M + 1)} (t) (t∈Ω _PRI ) are expressed in the form of Eq. (12), f ^(M Assuming that ⁺¹⁾ (t) (−∞ <t <∞) is not uniquely determined, an image derivative u ^{(M + 1)} (t) ≠ 0 that simultaneously satisfies the following (a) and (b) Exists.
(A) In t∈Ω _PRI , u ^{(M + 1)} (t) is expressed in the form of equation (12) (b)

From [Lemma 1], such u ^{(M + 1)} (t) satisfies the following equation.

Equation (15) satisfies equation (16) at the same time u ^{(M + 1)} (t) is, u ^{(M + 1)} from Lemma 2.1 (p.5) of Courdurier et al (Non-Patent Document 9) (t ) = 0 (-∞ <t <∞). This is true because u ^{(M + 1)} (t) ≠ 0 contradicts.

Return to the proof of the theorem. The M + 1 derivative g ^{(M + 1)} (t) (t∈Ω _ROI ) can be calculated from the measured projection data g (t) (t∈Ω _ROI ), and g ^{(M + 1) From the} a priori knowledge that (t) (t∈Ω _ROI ) and f (t) (t∈Ω _PRI ) are M-order piecewise polynomials, f ^{(M + 1)} (t) (−∞ <t < ∞) is uniquely determined. If f ^{(M + 1)} (t) (−∞ <t <∞) is uniquely determined, f (t) (t∈Ω _OBJ ) is uniquely determined by the following equation.

(3) Comparison with other a priori knowledge that enables rigorous reconstruction Reconstruction of interior CT since Ye et al. (Non-Patent Document 5) and Kudo et al. There have been several previous studies on the proposal of a priori knowledge that enables accurate image reconstruction. A priori knowledge of the present invention and those used in previous studies are summarized and shown in Table 2.

As is clear from the above detailed description, the present invention can be applied to any CT imaging apparatus based on the principle of measuring the line integral value of a physical quantity distribution inside an object and generating an image of the physical quantity distribution by data processing. is there.

Here, CT generally refers to absorption X-ray CT that generates an image of an X-ray absorption coefficient distribution in many cases, and therefore, the image reconstruction method according to the embodiment of the present invention will be described below. As an example of an apparatus to which is applied, an outline of an X-ray CT apparatus that obtains a cross-sectional image of the inside of a subject using X-rays will be described with reference to the drawings.

First, FIG. 11 attached will illustrate an embodiment of the present invention in which the above-described image reconstruction method is used to measure a line integral value of a physical quantity distribution inside an object and generate an image of the physical quantity distribution by data processing. The whole external appearance structure of a general X-ray CT apparatus is shown. That is, the X-ray CT apparatus accommodates components such as an X-ray irradiation unit, which will be described below, and a gantry unit 1 having a substantially cylindrical hollow portion in which a subject is positioned at the center thereof, A base unit 2 having a top plate (cradle) 4 on which an object is placed, a computer (not shown here) that is a data processing device, and an image obtained are displayed. A console unit 3 including a display device 5 and a keyboard 6 for performing necessary input is provided.

Subsequently, components constituting the X-ray CT apparatus are provided in the housing of the gantry unit 1 and the console unit 3 as shown in FIG. As an example, as shown in the figure, inside the housing, an X-ray generator 10 that irradiates a sample with X-rays in a fan shape, and a circle that detects X-rays irradiated from the apparatus and transmitted through the subject. Although not shown in the figure, the arc-shaped X-ray detection apparatus 20 is attached on, for example, a ring-shaped frame.

On the other hand, a space between the X-ray generator 10 and the X-ray detector 20 is provided with a top board 30 (corresponding to reference numeral 4 in FIG. 11) on which the subject is placed (set). . A member to which the X-ray generation device 10 and the X-ray detection device 20 are attached is a rotation drive mechanism such as a motor provided inside the gantry unit 1 via the rotation drive unit 50. Rotate in a predetermined direction at a predetermined rotation speed (see arrows in the figure). On the other hand, the top plate 30 for placing the subject is disposed so as to oppose a cylindrical space at a substantially central portion of the rotation surfaces of the X-ray generation device 10 and the X-ray detection device 20 and moves the sample placement table. It is moved by the unit 60. Furthermore, the X-ray generator 10 and the X-ray detector 20 are controlled by the rotation control of the motor, an X-ray high voltage unit 40 for generating and supplying a high voltage to the X-ray generator 10. A rotation driving unit 50 for rotating the member to which the motor is attached is provided.

The detection signal from the X-ray detection apparatus 20 described above is input to the data collection unit 70 and collected as image data, and further reproduced by the image reproduction unit 75 as a cross-sectional image or a three-dimensional image inside the sample. Is done. Reference numeral 76 in the drawing denotes a storage device (image memory) used when the image reproducing unit 75 reproduces a cross-sectional image or a three-dimensional image inside the sample. Further, the cross-sectional image or the three-dimensional image inside the sample reproduced by the image reproducing unit 75 is displayed on an image display unit 80 (corresponding to reference numeral 5 in FIG. 11) configured by a liquid crystal display device or the like, for example. . If a so-called touch panel (not shown) is incorporated in the image display unit 80, the image display unit 80 can perform input necessary for operating the apparatus. However, the present invention is not limited to this, and the apparatus may include a keyboard (corresponding to reference numeral 6 in FIG. 11), a numeric keypad, a mouse, and the like instead of the touch panel.

In addition, reference numeral 90 in the drawing indicates a control unit (corresponding to reference numeral 3 in FIG. 11) for controlling the operation of each unit constituting the above-described X-ray CT apparatus. More specifically, for example, it is constituted by a central processing unit (CPU), a storage device (memory) such as a RAM or a ROM, and an external storage device such as an HDD, and the like. Necessary control is executed based on software and firmware for controlling the operation of each unit stored in the.

The image reconstruction method for interior CT according to the present invention described above is stored in a storage device (memory) such as RAM or ROM as software in the image reproduction unit 75 constituting the X-ray CT apparatus, It is executed by a central processing unit (CPU).

The technology to which the present invention can be applied is not limited to this, for example, a phase X-ray CT for generating an image of a phase shift distribution from line integral data of a phase shift distribution when X-rays are irradiated, a radiopharmaceutical distribution administered into the body PET (Positron Emission CT) and SPECT (Single Photon Emission CT), CT, electron beam CT, and projection data using ultrasonic waves, microwaves, sound waves, seismic waves, etc. The present invention can also be applied to MRI (magnetic resonance imaging) using image reconstruction from images. That is, in the present invention, “object” or “image” refers to a spatial distribution of physical quantities to be imaged, and “projection data” refers to measurement data representing a line integral value on the straight line.

The numerical value of the projection data including phase shift, quantum beam phase shift, diffraction, or diffraction is extracted from a set of intensity data of a plurality of quantum beams acquired by a detector by adding an optical element or changing its position, It is also possible to reconstruct an image using numerical values of the projection data including phase shift, diffraction, or diffraction of the extracted quantum beam.

The above has described in detail the image reconstruction method of the interior CT according to various embodiments of the present invention. However, it is obvious that the present invention is not limited to the above-described embodiments, and includes various modifications. For example, the above-described embodiments are described in detail for the entire system in order to explain the present invention in an easy-to-understand manner, and are not necessarily limited to those having all the configurations described. Further, a part of the configuration of one embodiment can be replaced with the configuration of another embodiment, and the configuration of another embodiment can be added to the configuration of one embodiment. In addition, it is possible to add, delete, and replace other configurations for a part of the configuration of each embodiment.

Industrial availability

The present invention provides an image reconstruction method, particularly an interior CT image reconstruction method, that measures a line integral value of a physical quantity distribution inside an object and generates an image of the physical quantity distribution by data processing.

DESCRIPTION OF SYMBOLS 1 ... Gantry part, 3 ... Console part, 4, 30 ... Top plate, 10 ... X-ray generator, 20 ... X-ray detection apparatus, 90 ... Control part

Claims (15)

An image reconstruction method for interior CT, in which projection data is acquired by a quantum beam passing through an ROI inside the object to be imaged, and first-stage approximation is performed by the CT image reconstruction method using the projection data obtained above. A region in which the numerical value of the physical quantity representing the physical quantity in the ROI is at least piecewise uniform or piecewise represented by a polynomial is determined based on the reconstructed CT image. The numerical value of at least piecewise uniform or piecewise represented by a polynomial, and the position of the identified area and the numerical value of the image inside the piecewise uniform or piecewise represented by the polynomial, An image reconstruction method for interior CT, wherein a second-stage reconstruction is performed with higher accuracy than a one-stage reconstruction. An image reconstruction method for interior CT, in which projection data is acquired by a quantum beam passing through an ROI inside an imaging target, and a numerical value of an image representing a physical quantity in the ROI is at least piecewise uniform or piecewise polynomial The area represented by is fixed in advance in a part of the ROI, and the position of the area and the numerical value of the image are expressed in a piecewise uniform or piecewise polynomial form, An image reconstruction method for interior CT, characterized by performing reconstruction with high accuracy. 3. The interior CT image reconstruction method according to claim 1, wherein the numerical value of the projection data representing the physical quantity is obtained by absorption of the quantum beam by the imaging target. Configuration method. 3. The interior CT image reconstruction method according to claim 1, wherein the numerical value of the projection data representing the physical quantity is obtained by a phase shift of the quantum beam by the imaging target. Reconfiguration method. 3. The interior CT image reconstruction method according to claim 1 or 2, wherein the numerical value of the projection data representing the physical quantity is obtained by diffraction of the quantum beam by the imaging target. Configuration method. 3. The interior CT image reconstruction method according to claim 1, wherein the numerical value of the projection data representing the physical quantity is obtained by scattering of the quantum beam by the imaging target. Configuration method. The interior CT image reconstruction method according to any one of claims 4 to 6, wherein the numerical value of the projection data including phase shift, diffraction, or diffraction of the quantum beam is an addition of an optical element or a position thereof. Extracted from the set of intensity data of the plurality of quantum beams acquired by the detector by modification, and reconstructs an image using the numerical values of the projection data including phase shift, diffraction, or diffraction of the extracted quantum beams An image reconstruction method for interior CT. 2. The interior CT image reconstruction method according to claim 1, wherein the at least piecewise uniform or piecewise polynomial area specified in the ROI is obtained by the approximate reconstruction in the first stage. A method for reconstructing an image of an interior CT, which is manually set by a human using the obtained CT image. 2. The interior CT image reconstruction method according to claim 1, wherein the at least piecewise uniform or piecewise polynomial area specified in the ROI is obtained by the approximate reconstruction in the first stage. An image reconstruction method for interior CT, which is set by image processing using the obtained CT image. 2. The interior CT image reconstruction method according to claim 1, wherein the at least piecewise uniform identified in the ROI using a CT image obtained by the approximate reconstruction in the first stage. Alternatively, the piecewise polynomial region is set in advance in the ROI by at least one of specifying from the previously acquired CT image of the imaging target, a model representing the structure of the imaging target, and specifying from a priori information. A method for reconstructing an image of an interior CT. 2. The interior CT image reconstruction method according to claim 1, wherein the approximate reconstruction in the first stage includes a filtered back projection (FBP) method, a successive approximation method, and a statistical reconstruction method. An interior CT image reconstruction method, which is executed by at least one of CT image reconstruction methods. 2. The interior CT image reconstruction method according to claim 1, wherein the second-stage reconstruction includes differential backprojection Hilbert transform method, constrained successive approximation method, and constrained statistical reconstruction method. An image reconstruction method for interior CT, which is executed by at least one of the image reconstruction methods for CT. 10. An interior CT image reconstruction method according to claim 9, wherein the at least piecewise uniform identified in the ROI using a CT image obtained by the approximate reconstruction of the first stage. Alternatively, the interior CT image reconstruction method, wherein the piecewise polynomial region is formed including a part of the boundary of the ROI. 3. The interior CT image reconstruction method according to claim 2, wherein an area set in advance in the ROI includes a part of a boundary of the ROI. Reconfiguration method. 3. The interior CT image reconstruction method according to claim 2, wherein the position of a region set in advance in the ROI and the numerical value of the image within the ROI are expressed piecewise uniformly or piecewise by a polynomial expression. Interior CT, wherein the reconstruction is performed by at least one of a CT image reconstruction method including a differential backprojection Hilbert transform method, a constrained successive approximation method, and a constrained statistical reconstruction method Image reconstruction method.

Images