2.2.1

Pre-processing

From the image of a grid pattern, dots are segmented and grouped into horizontal lines and vertical lines as demonstrated in Fig 5.

Figure 5.

Dots are segmented and grouped into lines: (a) Horizontal lines; (b) Vertical lines.

2.2.2

Calculation of the center of distortion

Points on each group are fitted to parabolas in which the horizontal lines are represented by

and vertical lines by

where i, j are the index of the horizontal lines and vertical lines respectively. The rough estimate of the CoD, (x_c, y_c), is explained in the Fig. 6(a) where (x₀, y₀) is the average of the axis intercepts c of two parabolas between which the coefficient a changes sign. The slopes of the red and green line are the average of the b coefficients of these parabolas. For accurately determining the CoD coordinates, they are varied inside the bounds of these parabolas and metrics are calculated as the following steps: The point with minimum distance to the current CoD is located for each parabola. The horizontal and vertical parabolas yield two sets of points. Each set of points is fitted to a straight line. The sum of the intercepts of two fitted straight lines is the metric of the estimated CoD. The best CoD is the one with the minimum metric (Fig. 6(b)).

Figure 6.

Demonstration of estimating the center of distortion: (a) Coarse estimate; (b) Accurate estimate.

2.2.3

Calculation of the polynomial coefficients of the backward model

In the backward (BW) model, the relationship between an undistorted point (x_u, y_u) and a distorted point (x, y) is described as

where r_u and r are their distances from the CoD in the undistorted and distorted image, respectively; k_n are the coefficients. They are calculated by solving the following equation⁵

where and .

The undistorted intercepts, and , of undistorted lines can be determined without distortion correction. We use the assumption that the undistorted, uniform line spacing can be calculated and extrapolated from the area near the CoD which have negligible distortion. The undistorted intercepts are calculated as

where sgn() function returns the value of -1, 0, or 1 corresponding to its input of negative, zero, or positive value. i0 is the index of the line closest to the CoD. Δc is the average of the difference of c_i near the CoD. The same routine is used for the vertical lines as

In practice, using polynomial coefficients up to the fourth order is accurate enough, as there is no significant gain in accuracy by using higher orders.

Full stripe

A typical profile of a full stripe exhibits intensities that are offset at all angles compared with that of a neighboring good pixel (Fig. 13(a,b)). It gives rise to a half-ring artifact (in 180-degree tomography) in the reconstructed image (Fig. 13(c)).

Figure 13.

Demonstration of a full stripe artifact: (a) Stripe in the sinogram (arrowed); (b) Intensity profile along the stripe (red color) in comparison with an adjacent good pixel (blue color); (c) Ring artifact in the reconstructed image.

Partial stripe

Different to the full stripe, intensities of a partial stripe are offset only at certain ranges of angles as demonstrated in Fig. 14. As a result, it gives rise to a partial ring artifact in the reconstructed image (Fig. 14(c)).

Figure 14.

Demonstration of a partial stripe artifact: (a) Partial stripe in the sinogram (arrowed); (b) Intensity profile along the stripe (red color) in comparison with an adjacent good pixel (blue color); (c) Part-ring artifact in the reconstructed image.

Unresponsive stripe

Intensities of this type of stripe are independent of the angles. The pixel is not responsive to the variation of intensity versus the angle in the same way as its neighboring good pixel. This type of stripes may come from dead pixels of the sensor chip, light-blocking dusts or damaged scintillator (bright blobs in Fig. 9(a)) giving rise to stripes of constant brightness as clearly visible in Fig. 15(a). Missing information in these stripes strongly affects the reconstructed image; the constant intensity results in a prominent half-ring, and high-frequency edges of the stripes cause severe streak artifacts (Fig. 15(c)).

Figure 15.

Demonstration of an unresponsive stripe artifact: (a) Unresponsive stripe in the sinogram (red arrows; (b) Intensity profile along the stripe (red color) in comparison with an adjacent good pixel (blue color); (c) Ring artifact (red arrows) and streak artifacts (yellow arrows) in the reconstructed image.

Fluctuating stripe

Intensities fluctuate significantly between the angles. This type of stripes may come from defective pixels of the sensor chip rather than optical components. They are extremely small in number and their size is between 1 or 2-pixel as observed in our detector systems¹. Like unresponsive stripes, they give rise to both ring artifacts and streak artifacts as demonstrated in Fig. 16.

Figure 16.

Demonstration of a fluctuating stripe artifact: (a) Fluctuating stripe in the sinogram (zoomed in and arrowed); (b) Intensity profile along the stripe (red color) in comparison with an adjacent good pixel (blue color); (c) Ring artifact and streak artifacts in the reconstructed image.

Other types of stripe

There are other types of stripes which are the combination of the above defined stripes such as full stripes and partial stripes, or the extension of them such as large stripes or blurry stripes. Blurry stripes may come from the use of pre-processing methods such as phase retrieval or distortion correction. Large stripes may need a separate treatment to reduce the side effects of cleaning methods.

3.3.1

Equalization-based methods

Figure 18 shows three different ways of retrieving the underlying response curve of each pixel. In the first approach, each column of the sinogram (Fig. 18(a)) is fitted to a polynomial as shown in Fig. 18 (d). In the second approach, low-frequency components of each column are separated (Fig. 18(b, e)), and the grayscale values in each column are sorted in the last approach (Fig. 18(c,f). These underlying curves reveal the difference in response of adjacent pixels which is the cause of stripe artifacts. Correcting the difference in response of neighboring pixels is the way of removing stripe artifacts. This can be done by applying a smoothing filter along the horizontal direction of the retrieved images.

Figure 18.

Demonstration of different ways of retrieving the underlying response curves of each pixel. (a) and (d): Original sinogram and the result of fitting each column to a second order polynomial. (b) and (e) Original sinogram and the result of separating the low-frequency components of each column. (c) and (f) Original sinogram and its column-sorted intensities.

From the smoothed images, the stripe-cleaned images are obtained differently using above approaches. In the first approach, called the fitting-based approach, the corrected sinogram is the result of dividing the multiplication of the original sinogram with the smoothed image to the fitted sinogram. In the second approach, called the filtering-based approach, the corrected sinogram is the addition of the smoothed image and the high-frequency components of the original sinogram. In the last approach, called the sorting-based approach, simply re-sorting the smoothed image returns the cleaned sinogram. Detail of these implementations and python source code can be found in Ref. [8] and at https://github.com/nghia-vo/sarepy. The effect of applying these removal methods can clearly be seen in the reconstructed images before and after correction (Fig. 19).

Figure 19.

Reconstructed images before and after correction. (a) and (d) Without and with the fitting-based approach. (b) and (e) Without and with the filtering-based approach. (c) and (f) Without and with the sorting-based approach.

Each of these approaches has its own pros and cons as summarized in Table. 1. However, we can combine them in different ways to achieve superior approaches. This will be presented in the section 3.3.3.

Table 1.

Summary of pros and cons of each technique.

3.3.2

Detection of stripe artifacts

The above techniques are able to remove full and partial stripe artifacts which are of small or medium size without significantly affecting other areas. However, for large stripes, applying a stronger filter to the whole sinogram will degrade the final image. Furthermore, the large stripes are few in number. In this case, the correction could be selectively applied only on defective pixels. To do that, we propose a segmentation algorithm which works on a 1D array to locate stripes. Depending on the types of stripe artifacts, different pre-processing steps are used to generate the 1D array which then is input to the segmentation algorithm. Steps of this algorithm⁸, called SFTS (Sorting, Fitting, and Thresholding-based Segmentation) for short, is described as follows and demonstrated in Fig. 20.

R can be understood as a signal-to-noise ratio which controls the sensitivity of the algorithm. A smaller R is more sensitive to detect the stripes. A reasonable choice of R is around 3.0. This algorithm can be used as a binarization method. This is convenient for users because one does not need to know the absolute values of the array.

Figure 20.

Demonstration of the detection algorithm. (a) Normalized 1D array, i.e. non-uniform background is corrected. (b) Sorted array and fitted array using the middle part of the sorted array.

3.3.3

Combination of techniques for tackling all types of stripe artifacts

From the equalization-based techniques and the SFTS algorithm, we can derive various algorithms for removing different types of stripe artifacts.

Removal of large stripes

Large stripes (Fig. 21(a)) may come from partially defective regions (Fig. 10(b)) or the adjacent areas of the damaged scintillator which receive extra scattered light, i.e. the so-called halo effect (Fig. 9(a)). To detect them using the SFTS algorithm from the sinogram, the pre-processing steps are:

• 1 - Sort intensities in each column of the sinogram (Fig. 21(b)).

• 2 - Apply the strong median filter along each row to remove stripes (Fig. 21(c)).

• 3 - Average along the columns of the sorted sinogram where some percentage of pixels at the top and bottom are dropped. This simple technique helps to reduce the possibility of wrongly detecting stripes caused by high-frequency edges of the sinogram (Fig. 21(a)). It also can be used to improve other ring removal methods.

• 4 - Do the same for the smoothed sinogram.

• 5 - Divide the result of step 3 to the result of step 4 to get the normalized 1D array.

• 6 - Use the SFTS algorithm to get the locations of the stripes.

Then large stripes are removed by:

• 1 - Normalize each row of the sinogram using the result of the step 5 in the pre-processing algorithms (Fig. 22(b)). This step helps to correct the non-uniform background around the large stripes (Fig. 22(a)).

• 2 - Apply the sorting-based algorithm with the strong filter to remove large stripes.

• 3 - Replace intensities in the stripes of the normalized sinogram with the one of the result of step 2 (Fig. 22(c)).

Figure 21.

Demonstration of pre-processing steps used for detecting large stripes. (a) Original sinogram with large stripes and high-frequency edges (arrowed). (b) Averaging along columns of the sorted sinogram inside the blue box to avoid the false detection of stripes caused by the edges in (a). (c) Same as in (b) but on the smoothed sinogram.

Figure 22.

Results of removing large stripes. (a) Reconstructed image without ring removal; (b) Reconstructed image after the normalization step; (c) Final reconstructed image.

Removal of unresponsive stripes and fluctuating stripes

The unresponsive stripes and fluctuating stripes (Fig. 23(a)) give rise to both ring artifacts and streak artifacts in the reconstructed image. Their intensity profiles show opposite characteristics. The unresponsive stripe shows very little variation (Fig. 15(b)) while the fluctuating stripe shows excessive variation (Fig. 16(b)). Exploiting these features can help us detect them all together based on the SFTS algorithm by the following pre-processing steps:

• 1 - Apply the strong mean filter along each column of the sinogram (Fig. 23(b)).

• 2 - Take absolute values of the difference between the result of step 1 and the original sinogram (Fig. 23(c)).

• 3 - Average along each column of the result of step 2 resulting a 1D array.

• 4 - Apply the strong median filter to the result of step 3.

• 5 - Divide the result of step 3 to the result of step 4 to get the normalized 1D array.

• 6 - Use the SFTS algorithm to get the locations of the stripes.

Then the stripes are removed by:

• 1 - Interpolate values inside the stripes from neighboring pixels.

• 2 - Apply the algorithm of removing large stripes to the result of step 1. This step is needed because intensities around the unresponsive stripes are modulated by the scattered light resulting the large stripes.

Figure 23.

Demonstration of pre-processing steps used for detecting unresponsive and fluctuating stripes. (a) Original sinogram having both unresponsive stripe (right arrows) and fluctuating stripe (left arrows). (b) Smoothed sinogram using the strong mean filter along each column. (c) Absolute difference between (a) and (b).

Figure 24.

Results of removing the unresponsive and fluctuating stripes. (a) Reconstructed image without ring removal; (b) Reconstructed image after the interpolation step; (c) Final reconstructed image.

3.3.4

Applications of the sorting-based method

In our approaches, the sorting technique has proven its efficiency in tackling partial stripes, detecting stripes, and reducing extra artifacts. Indeed it can be combined with other ring removal methods to improve their efficiency. In the filtering-based method, the sorting-based method can be applied to the low-frequency components of the sinogram columns (Fig. 18(e)). This combination helps to remove void-center artifacts introduced by the filtering-based approach and reduce the streak artifacts introduced by the sorting-based approach. In the fitting-based approach, the sorting technique can be used before the fitting step. The sorted intensity profiles can easily be fitted to polynomial functions. This helps to extend the applications of the fitting-based method to more complex sinograms.

Well-known methods can be improved by using the average along the columns of the sorted sinogram where some percentage of pixels at the top and bottom are dropped (Fig. 21(b)). This helps to reduce extra stripe artifacts in the normalization-based and regularization-based methods. Applying stripe removal methods to the sorted sinogram help to remove void-center artifacts (Fig. 17(c) and (e)).

Here we present two more applications of the sorting-based technique which have not been introduced before.

Removing stripe artifacts using a 2D smoothing filter and the sorting-based technique

In the filtering-based method introduced in section 3.3.1, the filtering step which separates the low-frequency component from the high-frequency component is applied to each intensity profile of each sinogram column. Then stripes are removed by applying a smoothing filter perpendicular to the low-frequency components of the sinogram columns. Here we introduce another way (demonstrated in Fig. 25) where a 2D strong smoothing filter is used to separate two frequency components of the sinogram, then the stripes in the high-frequency component is removed by the sorting-based method. This approach uses the same assumption as the FFT-based methods in which stripes artifacts are corresponding to high-frequency components. The advantage is that it does not yield void-center artifacts.

Figure 25.

Demonstration of another method for removing stripe artifacts using a 2D smoothing filter and the sorting-based technique. (a) Low-frequency components of the sinogram in Fig. 18(c). (b) High-frequency component of the sinogram in Fig. 18(c) in which stripe artifacts are enhanced (arrowed). (c) Stripe-removed image of (b) using the sorting-based technique. (d) Combination of (a) and (c).

Removing ring artifacts in cone-beam tomography

In a cone-beam geometry, projections of a sample slice at different angles are not constrained in a single row of a detector. As a result, reconstruction is performed directly without the intermediate step of sinogram creation. Furthermore, partial rings are more popular in the reconstructed image because projections of a point in the reconstruction space do not always stay inside the field of view of defective areas. As there are no sinograms, methods of removing ring artifacts working on the reconstructed images are reasonable choices. However, here we demonstrate how the sorting-based method, which works in the projection domain, can be used to remove ring artifacts in the reconstructed image of a cone-beam system²⁵.

The idea is very simple. The steps of the sorting-based method are applied to the series of 2D projections: sorting intensities along the angular direction of each pixel; applying a 2D edge-preserving smoothing filter (e.g the median filter) to each sorted image; and resorting the smoothed images to the original positions to get the corrected images. This routine equalizes the responses of neighboring pixels in the projection domain which helps to remove ring artifacts in the reconstruction domain. The algorithm is quite memory-consuming because the indices of a 3D array need to be kept for the resorting step. Results of this approach in comparison to the wavelet-fft approach²³, which is applied to the polar coordinate transformation of the reconstructed image, are shown in Fig. 26 and Fig. 27. As can be seen in these images, our approach not only is able to remove all full and partial rings but also leave no artifacts along the vertical center of the reconstructed volume.

Figure 26.

Reconstructed image from a cone-beam system. (a) Without ring removal; (b) Using the wavelet-FFT approach but some partial rings still left (arrowed); (c) Using the sorting-based approach in 3D.

Figure 27:

Xz-slice of the 3D reconstructed volume. (a) Without ring removal. (b) Using the wavelet-FFT approach where artifacts are left at the middle. (c) Using the sorting-based approach in 3D.

Python implementations of presented methods can be found at https://github.com/nghia-vo/sarepy. These approaches are also available in other tomographic software packages such as Tomopy²⁶ and Savu²⁷.

Acquiring a calibration image

An ideal calibration image is the one can provide information about the intensity profiles inside the dark areas, which are not exposed to the X-ray, and sharp edges between the free spaces and the dark areas. Furthermore, any contribution from diffraction owning to the coherent property of the X-rays need to be removed or minimized. Figure 32 shows some examples of calibration images can be used to determine the MTF. The X-ray image in Fig. 32(a) is the projection of a long Tungsten rod held by a plastic holder. The rod is aligned to be parallel to the X-rays to block any refracted X-rays going inside the dark area. This type of calibration image allows to determine the 1D MTFs at different directions in the scintillator plane, then they can be combined to form a 2D MTF function. The calibration images in Fig. 32(b) and (c) can be used to determine the x and y-component of the MTF separately.

Figure 32.

Calibration image used for determining the MTF. (a) Projection of a long Tungsten rod. (b) Projection of two thick Tungsten plates. (c) Projection of a long Tungsten plate.

Determining the MTF function

Our approach aims to determine the MTF directly because it is most efficiently used for correcting a large number of projections of a tomographic dataset, instead of using the PSF (Point Spread Function) ²⁹ However, because the effect of the scintillation photon scattering depends on the quality of a scintillator such as the uniformity of thickness, material distribution, and scattering direction, it is challenging to design a generic MTF. For simplicity, here we apply two constraints: the MTF is space-invariant and separable.

The formula of the filter window is defined as

where is the unit function which can be changed to best describe an optical component; u is a variable in the range of [-1;1]; a_k are coefficients. To determine a_k we use a global optimization technique which is using the best of multiple estimates (BME)³⁰. It is an improvement of the simulated annealing technique in which the dependency on a cooling scheme is removed. The cost function is evaluated by using the combination of increasing the slope of the edges, reducing the intensity values in the dark area, and reducing the ringing artifacts caused by the deconvolution in Fourier space³¹ from the intensity profile as shown in Fig. 33. The calculation routine is as follows:

Figure 33 shows the results of determining the x-direction MTF of a detecting system using a scintillator made of GGG:Eu, 17.5 micron thickness. The calibration image (Fig. 32) was taken at 53 keV. Here we used n=1/2 and N=200. 200 iterations were performed. Note that the calculation routine is applied to the image without flat field correction. As can be seen in Fig. 33(c) the slope of the edges which defines the sharpness of the image is increased and the exponential decrease of the intensity in the dark area is corrected.

Figure 33.

Results of the determination of the MTF. (a) Convergence of the cost function. (b) 1D determined MTF. (c) Intensities profile before deconvolution (Blue) and after convolution (Red).

The MTF in y-direction is determined independently using the vertical intensity profile at the middle of the calibration image in Fig. 32(a). Then it is combined with the x-direction MTF to form the 2D-MTF (Fig. 34(a)). Deconvolution of the calibration image using this MTF and the difference is shown in Fig. 34(b) and (c).

Figure. 34.

Deconvolution results using the determined MTF. (a) 2D MTF combined from the x-direction MTF and the y-direction MTF. (b) Deconvolved image in Fig. 32(a). (c) Difference between (b) and Fig. 32(a).