CN1301083C - 1/4 path deviation interpolation method for CT system - Google Patents

Abstract

A 1/4 channel offset interpolation method for CT systems, belonging to the field of image processing technology, the method includes four steps: 1. Interpolate the original channel DATA to obtain DATA1, and use the original channel for the data that needs to be interpolated. Obtained by simple interpolation, the interpolated data and the original data together form new data; 2. Use the opposite channel to obtain the corresponding data DATA2 of DATA1; 3. The low frequency part of DATA1 is combined with the high frequency part of DATA2 to obtain data DATA3, respectively for DATA1 Perform low-pass filtering and high-pass filtering with DATA2, and combine the obtained data to obtain DATA3; 4. Use the original data DATA and DATA3 to combine the data that does not overlap with the original data to output the result ResultDATA. Since the data of the opposite channel of the direct 1/4 offset interpolation is collected by different channels, and the time is inconsistent, it is difficult for the general system to ensure this consistency, so the obtained image will have some artifacts, the present invention proposes This improved 1/4 interpolation method can suppress such artifacts very well.

Description

A 1/4 channel offset interpolation method for CT system

technical field

The invention belongs to the technical field of image processing, in particular to an improved 1/4 channel offset interpolation method applied to a CT system.

Background technique

Spatial resolution is an important indicator of CT system. There are many ways to improve spatial resolution, such as adding physical detectors, flying focus, changing convolution kernel or image post-processing, etc. wait. The method of increasing the physical detector can increase the cut-off frequency of the system, but it will greatly increase the system cost, and also increase the complexity of the system, such as increasing the amount of data transmission and increasing the processing time of data; the method of flying focus can also Improve the cut-off frequency of the system, but this requires special tube technology, which is technically difficult; changing the convolution kernel or image post-processing method will enhance the effect of a special frequency band in a certain part of the image, but it is impossible to increase the cut-off frequency of the system.

Another method to improve the spatial resolution of the system is 1/4 channel offset interpolation reconstruction, which can increase the cutoff frequency of the system without changing the hardware. For simplicity of description, only equiangular fan beams are analyzed and described in the present invention, and the same is true for other types of scans.

Since the general fan beam is scanned at 360 degrees, each ray is sampled twice, that is, each ray has an opposite ray corresponding to it. As shown in Figure 1, the dotted line is emitted by Tube Position 1. On Tube Position 2, there is a ray completely coincident with it at the same position, and these two rays are opposite channels.

Assuming that the ray is (β, γ), where β is the sampling angle View, and γ is the fan angle, then the opposite channel ray is (β′, γ′) satisfying formula (1):

β'=π+β-2γ

(1)

γ′=-γ

Since the actual sampling is limited, in general, if you want to obtain the data of the opposite channel, you need to obtain it through two-dimensional interpolation. If you want to obtain twice the number of original channels through interpolation of the opposite channel, that is, to obtain the data between the original sampling two adjacent channels through interpolation, the discrete representation method is as follows:

Assuming that the View angle interval is Δβ, the channel angle interval is Δγ, and the subscript of the center channel is Po, then (iΔβ, (j+0.5-Po)Δγ) these data need to be interpolated through the opposite channel, where i, j are integers, respectively Represents a certain sample and a certain channel, and 0≤i<V, 0≤j<N, V, N are the total number of samples and the total number of channels respectively, according to formula (1), the opposite sides of these channels to be interpolated can be obtained Channel (i′, j′) can be expressed as formula (2):

i′＝Vie/2+i-(2j-1-2Po)·Δγ/Δβ (2)

j'=2Po-j-0.5

Through formula (2), you can use two-dimensional interpolation to obtain the desired data, such as bilinear interpolation, and perform normal convolution back projection after obtaining the interpolated data. However, the two-dimensional interpolation has too much influence on the resolution, and the interpolation in the channel direction has a particularly large impact on the resolution, so the image resolution obtained by simply performing interpolation to expand the channel cannot be significantly improved. But if Po satisfies formula (3)

Po=z+0.5±0.25                                                                    , there is only interpolation in the View direction, so when Po satisfies formula (3), using formula (2) to interpolate and expand the channel can improve the spatial resolution of the system, as shown in Figure 2, where iso-center is the detection corresponding to the rotation center The detector channel, that is, the channel corresponding to the Po mentioned above, the detector width indicates the width of the detector, and the interleaved sample is the opposite channel. In fact, it is impossible for Po to satisfy formula (3) exactly, so generally only an approximate center channel can be obtained. Facts have proved that as long as Po is within a certain range of formula (3), the spatial resolution of the system can be obtained very clearly improvement.

The process described above is a direct 1/4 offset interpolation reconstruction method. According to the principle of sampling twice on each ray in general scanning, the method uses the sampling of the opposite channel to expand and interpolate the original data, which doubles the number of channels on a view, thereby increasing the cut-off frequency of the system . Because the data of the opposite channel is used, the consistency between the opposite channel and the original channel is very important. If the data is inconsistent, the built image will have many artifacts, and in severe cases, the image cannot be diagnosed; and due to the original channel The data and the data of the opposite channel are obtained at different times and at different locations. In addition, due to mechanical and tube focus shifts, etc., the possibility of data inconsistency is very high, so the 1/4 channel offset interpolation method is used directly to obtain Images generally have a lot of artifacts.

Contents of the invention

Aiming at the deficiencies of the prior art, the present invention provides a 1/4 channel offset interpolation method for CT systems, which can effectively reduce image artifacts, and at the same time, the cut-off frequency can be guaranteed to be the same as the direct 1/4 interpolation .

The inventive method comprises the following four steps:

1. Interpolate the original channel DATA to get DATA1

The central channel of the system requires formula (3). For the data that needs to be interpolated, the original channel is used for simple interpolation first, and the interpolated data and the original data form new data together.

2. Interpolate the channel opposite to DATA1 to obtain DATA2.

3. The low frequency part of DATA1 is combined with the high frequency part of DATA2 to obtain data DATA3

Use one-dimensional high-pass filtering to perform low-pass filtering and high-pass filtering on DATA1 and DATA2 respectively, and combine the obtained data to obtain DATA3.

4. Combine the data in the original data DATA and DATA3 that do not overlap with the original data to output the result ResultDATA

Since the data of the opposite channel of the direct 1/4 offset interpolation is collected by different channels, and the time is not consistent, it is difficult for the general system to ensure this consistency, so the obtained image will have some artifacts, the invention proposes This improved 1/4 interpolation method can suppress such artifacts very well.

Description of drawings

Fig. 1 is the opposite passage schematic diagram of the present invention;

Fig. 2 is a schematic diagram of 1/4 channel offset of the present invention;

Figure 3 is a comparison of CT-C3000 phantom results, where Figure 3a is the original 1/4 offset interpolation result, and Figure 3b is the improved result;

Figure 4 is a comparison of CT-C3000Dual phantom results, in which Figure 4a is the original 1/4 offset interpolation result, and Figure 4b is the improved result;

Figure 5 is a comparison chart of CT-C3000 phantom and MTF curve results, in which Figure 5a is the original 1/4 offset interpolation result, and Figure 5b is the improved result;

Figure 6 is a comparison chart of CT-C3000Dual phantom and MTF curve results, in which Figure 6a is the original 1/4 offset interpolation result, and Figure 6b is the improved result.

Detailed ways

Its specific implementation steps of the method provided by the present invention are as follows:

1. Interpolate the original channel DATA to get DATA1

The central channel of the system requires the following formula (3). For the data (iΔβ, (j+0.5-Po)Δγ) that needs to be interpolated, first use the original channel DATA to perform simple interpolation. You can use cubic spline or Lagrange interpolation, or Obtained by the method of frequency domain migration, the data composed of the interpolated data and the original data is DATA1(iΔβ, (j-2Po)Δγ/2), where 0≤i<V, 0≤j<2N.

2. Use the opposite channel to obtain the corresponding data DATA2 of DATA1

All the data of DATA1 has its opposite channel, and the interpolation of the opposite channel can be obtained according to formula (2), that is, the result obtained by the interpolation of the opposite channel is DATA2(iΔβ, (j-2Po)Δγ/2).

3. The low frequency part of DATA1 is combined with the high frequency part of DATA2 to obtain data DATA3

Let I be the identity operator, Γ is a one-dimensional high-pass filter, and for each sample i, the following formula (4) is performed

HighFre(iΔβ, j-2Po)Δγ/2)＝Γ(DA TA2(iΔβ, (j-2Po)Δγ/2))

(4)

LowFre(iΔβ, (j-2Po)Δγ/2)=(I-Γ)(DA TA1(iΔβ, (j-2Po)Δγ/2)) to get HighFre(iΔβ, (j-2Po)Δγ/2), LowFre(iΔβ, (j-2Po)Δγ/2), and then through formula (5) DATA3(iΔβ, (j-2Po)Δγ/2) = HighFre(iΔβ, (j-2Po)Δγ/2)

(5)+LowFre(iΔβ, (j-2Po)Δγ/2) to get DATA3(iΔβ, (j-2Po)Δγ/2), where j is an odd number, and 0<j<2N. The one-dimensional high-pass filter Γ mentioned in the above process can use an ideal frequency-domain high-pass filter, and the expression is as follows:

$\mathrm{<span\; class="notranslate">\; F\; Re\; sult\; Data</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\left\{\begin{array}{cc}\mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; \&le;</span>\mathrm{<span\; class="notranslate">\; D.</span>}\\ \mathrm{<span\; class="notranslate">\; FSourceData</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; ></span>\mathrm{<span\; class="notranslate">\; D.</span>}\end{array}\right.<span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 6</span>}<span\; class="notranslate">\; )</span>$

Among them, FSourceData is the Fourier transform of the original data, FResultData is the Fourier transform of the result data, and D is the cutoff frequency.

The following high-pass filter can also be used, whose expression is:

$\mathrm{<span\; class="notranslate">\; ResultData</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 2</span>}\mathrm{<span\; class="notranslate">\; NSourceData</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; 2</span>}\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; -</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 2</span>}\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}\underset{\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; N</span>}}{\overset{\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; \&NotEqual;</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}\mathrm{<span\; class="notranslate">\; SourceData</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 7</span>}<span\; class="notranslate">\; )</span>$

Among them, 2N+1 is the filter width, N is a positive integer, SourceData is the original data, and ResultData is the filtered data. The performance of the above filters on the image is basically the same.

4. Combine the data in the original data DATA and DATA3 that do not overlap with the original data to output the result ResultDATA

The data obtained by formula (5) and the original data form new data ResultDATA to complete the entire interpolation process, and then use the normal convolution back-projection process to obtain an image.

CT images are processed according to the improved 1/4 channel offset interpolation method provided by the present invention, Fig. 3 is a comparison chart of CT-C3000 phantom results, Fig. 4 is a comparison chart of CT-C3000Dual phantom results, as can be seen from the figure The artifacts have been well eliminated; Figure 5 is a comparison chart of CT-C3000 phantom and MTF curve results, and Figure 6 is a comparison chart of CT-C3000Dual phantom and MTF curve results. It can be seen from the figure that MTF Curves have little effect.

Claims (3)

1. A 1/4 channel offset interpolation method for a CT system, comprising the following four steps: Step 1: Perform simple interpolation on the original channel DATA to obtain DATA1, Step 2, interpolating the channel opposite to DATA1 to obtain DATA2, Step 3, use one-dimensional high-pass filtering to obtain the low frequency part of DATA1 and the high frequency part of DATA2, and combine the obtained data to obtain data DATA3, Step 4: Combine the data in the original data DATA and DATA3 that do not overlap with the original data to output the result ResultDATA. 2. A kind of 1/4 channel offset interpolation method for CT system according to claim 1, it is characterized in that the simple interpolation described in step 1 adopts cubic spline or Lagrangian interpolation or frequency domain offset method of moving. 3. A kind of 1/4 channel offset interpolation method for CT system according to claim 1, it is characterized in that the one-dimensional high-pass filtering described in step 3 comprises ideal frequency-domain high-pass filter, and expression is : $\mathrm{<span\; class="notranslate">\; F\; Re\; sult\; Data</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\left\{\begin{array}{cc}\mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; \&le;</span>\mathrm{<span\; class="notranslate">\; D.</span>}\\ \mathrm{<span\; class="notranslate">\; FSourceData</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; ></span>\mathrm{<span\; class="notranslate">\; D.</span>}\end{array}\right.$ Among them, FSourceData is the Fourier transform of the original data, FResultData is the Fourier transform of the result data, and D is the cut-off frequency; A high-pass filter is also included, expressed as: $\mathrm{<span\; class="notranslate">\; ResultData</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 2</span>}\mathrm{<span\; class="notranslate">\; NSourceData</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; 2</span>}\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; -</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 2</span>}\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}\underset{\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; N</span>}}{\overset{\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}\mathrm{<span\; class="notranslate">\; SourceData</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span>$ Among them, 2N+1 is the filter width, N is a positive integer, SourceData is the original data, and ResultData is the filtered data.

Images