Mixed one-bit compressive sensing with applications to overexposure correction for CT reconstruction

01/03/2017 ∙ by Xiaolin Huang, et al. ∙ Shanghai Jiao Tong University FAU Michigan State University Stanford University FUDAN University 0

When a measurement falls outside the quantization or measurable range, it becomes saturated and cannot be used in classical reconstruction methods. For example, in C-arm angiography systems, which provide projection radiography, fluoroscopy, digital subtraction angiography, and are widely used for medical diagnoses and interventions, the limited dynamic range of C-arm flat detectors leads to overexposure in some projections during an acquisition, such as imaging relatively thin body parts (e.g., the knee). Aiming at overexposure correction for computed tomography (CT) reconstruction, we in this paper propose a mixed one-bit compressive sensing (M1bit-CS) to acquire information from both regular and saturated measurements. This method is inspired by the recent progress on one-bit compressive sensing, which deals with only sign observations. Its successful applications imply that information carried by saturated measurements is useful to improve recovery quality. For the proposed M1bit-CS model, alternating direction methods of multipliers is developed and an iterative saturation detection scheme is established. Then we evaluate M1bit-CS on one-dimensional signal recovery tasks. In some experiments, the performance of the proposed algorithms on mixed measurements is almost the same as recovery on unsaturated ones with the same amount of measurements. Finally, we apply the proposed method to overexposure correction for CT reconstruction on a phantom and a simulated clinical image. The results are promising, as the typical streaking artifacts and capping artifacts introduced by saturated projection data are effectively reduced, yielding significant error reduction compared with existing algorithms based on extrapolation.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 2

page 7

page 8

page 9

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

I Introduction

In a real-world measuring system, the measurable range could be limited and the final output for a linear measurement can be described as

(1)

where

is the signal to be estimated,

is a sensing vector,

stands for the noise, and and are the lower and upper saturation thresholds, respectively. When an observation equals or , we say that this measurement is saturated. The saturation phenomenon may happen when the detector has a narrow range [1, 2, 3]. In this paper, we discuss how to recover from a sensing system with saturation and then apply the proposed method to correct overexposure for C-arm computed tomography (C-arm CT).

C-arm CT is characterized with high spatial resolution, large field of view, and 3D imaging capability. It has become a valuable tool for medical diagnoses, therapy guidance, and interventions. The key technology that enables these advances is the flat-panel detector, which provides good soft-tissue imaging performance, distortion-free images as well as high detective quantum efficiency [4]. However, although typical C-arm flat-panel detectors have a dynamic range of 14-bit that is sufficient for conventional fluoroscopy, it may not be high enough to eliminate overexposure in all projections acquired during a 3D acquisition, such as imaging relatively thin body parts (e.g., the knee). If no counter-measures are carried out, strong direct radiation may hit certain detector regions in some views (the intensity range of the traveled X-rays are greater than the detector’s inherent detectable range) and thus overexpose them. Consequently, reconstruction from overexposed projections typically leads to incorrect Hounsfield unit (HU) values, streaking artifacts, and loss of clear outer boundaries. Moreover, HU values of a homogeneous object tend to become smaller when they are far away from the object center, resulting in so-called capping artifacts. Since all these overexposure artifacts impair the final visualization of low contrast objects, it is of practical significance to develop an appropriate overexposure correction scheme on the saturated images to compensate for these artifacts.

To give an intuitive impression of saturation phenomenon happens in CT scan, a knee phantom is designed in Fig.1. Its projection data without and with saturation are shown in Fig.1 and 1, respectively. If there is no saturation, then the object can be reconstructed well by classical analytical algorithms, such as the standard filtered back projection (FBP, [5]) or iterative reconstruction methods, such as the simultaneous algebraic reconstruction technique (SART, [6]). When there is saturation in the sinogram data, e.g. Fig.1(d), CT reconstruction from the saturated data is quite hard. For example, directly applying FBP on the saturated data in Fig.1 outputs Fig.1, which is far from satisfactory.

Fig. 1: (a) A knee phantom; (b) reconstructed by FBP from saturated projections shown by (d); (c) full projection data; (d) saturated projection data (the measurable range ; see Section V-A for details).

Due to the loss of information, signal recovery from saturated measurements are not straightforward. In CT reconstruction, primarily targeted to correct the overexposure artifacts for knee imaging, an optimization-based multiple shape fitting approach was suggested to extrapolate the missing data using a pair of cylinder shapes that are fitted in the sinogram domain [7]. The parameters describing the cylinders are estimated from the overexposed data by minimizing the least square error. However, the method would face challenges when it comes to complicated shapes that cannot be modeled with a composition of cylindrical or elliptical shapes. Furthermore, additional hardware, such as a shaped pre-filteration (bowtie) can be deployed to correct the saturation during the data acquisition [8]. However, that method requires the imaged object to be exactly centered in the rotation center and the use of a bowtie filter also appears not flexible enough for the applications of C-arm CT systems. Another way to prevent saturation at the object is to apply modeling clay wrapped around the imaged object as an additional absorber during the data acquisition [9, 10]. However, the method is not ideal since the modeling clay may cause great discomfort for the patient. Alternatively, a Kinect camera was used to obtain the object depth information, which can be combined into a C-arm system for correcting overexposured projections [11]. More specifically, the Kinect depth data is used to find the points of intersection between the X-ray beam path and object surface. Overexposed pixels on the projection images are corrected by extrapolating the absorption along the corresponding line integrals.

In contrast to the existing saturation correction methods that require prior-knowledge, in this paper we propose a data-driven method to find directly from measurements observed by (1). The key is to acquire information from the saturated measurements, which do not provide analog information but do contain one-bit information: when a measurement is upper-saturated (lower-saturated), i.e., , we know that is greater (smaller) than or equal to ().

Data mining from inequalities is theoretically possible when is sparse, which has been verified by the recent progress on one-bit compressive sensing (1bit-CS). The concept of 1bit-CS was firstly introduced by [12] to recover sparse signals from measurements containing only one-bit information, i.e., sign measurements. It is rooted in compressive sensing (CS, [13, 14, 15]), and there are some interesting discussions on both theory and algorithms; see, e.g., [16] [17] [18] [19] [20] [21]. Unlike the pure 1bit-CS, (1) has analog measurements as well. Therefore, we refer to the method dealing with (1) as mixed one-bit compressive sensing (M1bit-CS). There is few study on this topic. To the best of our knowledge, it has been considered only by [22, 23]. But they put major emphasis on quantization error and assumed that there was no noise on the saturation measurements.

In this paper, we are going to recover sparse signals from noise-corrupted measurements with saturation. The major contributions include:

  • propose M1bit-CS models;

  • establish the corresponding alternating direction method of multipliers (ADMM);

  • design an iterative saturation detection scheme;

  • evaluate M1bit-CS on one-dimensional signal recovery problems;

  • apply M1bit-CS on overexpose correction for C-arm CT and achieve very promising performance.

Before discussing M1bit-CS, we introduce some notations to describe (1) that is a linear sensing system with saturation. The dimension of the signal is denoted by , the number of total measurements by , and the number of saturated measurements by . A boolean vector is used to indicate a measurement is saturated or not. Mathematically,

Furthermore, for saturated measurements, i.e., , we use to distinguish the upper-saturation () and lower-saturation ():

(2)

where and . Then our task of signal recovery from (1) is to find a sparse signal such that is close to when and coincides with the saturation inequality (2) when .

The rest of this paper is organized as follows. In Section II, M1bit-CS models are established. Fast algorithms are designed in Section III. Section IV evaluates the proposed algorithms in numerical experiments. Section V applies the proposed method to overexposure correction for C-arm CT reconstruction. In Section VI, a conclusion is given to end this paper.

Ii Mixed One-bit Compressive Sensing

Ii-a (One-bit) Compressive Sensing

A classical compressive sensing (CS) model for unsaturated data, i.e., , can be formulated as the following well-known model:

(3)

The loss function in the above formulation could be the least squares loss. The related theory and algorithms have been insightfully investigated in the last decade; see

[24] and references therein.

If we only have sign measurements, which could be regarded as an extreme quantization process, or equivalently and being very close to zero in (1)111We assume that the negative measurements are upper bounded by ., we obtain the so-called one-bit compressive sensing (1bit-CS, [12] [16]). Many 1bit-CS models can be expressed as follows:

(4)
s.t.

where is a given constant. There are two differences between (3) and (4). First, a different loss function is used for the inconsistency of the sign. Second, because multiplying a signal by a positive constant does not change the signs of linear measurements, the norm constraint ( is always assumed) is needed.

In recent works on 1bit-CS ([19] [20] [21]), the nonconvex constraint is relaxed to for computational efficiency. But this relaxation makes it inconvenient to use the hinge loss, which is a natural choice for sign inconsistency: minimizing in a norm-ball leads to a trivial and unreasonable solution . Therefore, in convex 1bit-CS models, usually takes the linear loss , which could be explained as to pursue high correlation between and .

Ii-B Mixed One-bit Compressive Sensing

To consider both analog measurements and one-bit information, we propose the following mixed one-bit compressive sensing (M1bit-CS):

(5)
s.t.

Instead of putting the norm in the constraint, we can also place it in the objective function and come to the following problem,

(6)

where . This modification is denoted as M1bit-CS with -norm regularization (M1bit-CSR) and the original one (5) is M1bit-CS with -norm constraint (M1bit-CSC). When is known or can be estimated accurately, one can choose M1bit-CSC, otherwise, M1bit-CSR is applicable.

For the analog part, we use the least squares loss for in this paper. Its robust modifications, such as Huber loss, are also applicable. For the saturation part, directly inheriting the idea of 1bit-CS leads to the linear loss. But here we have also analog measurements that can provide magnitude information, and it makes using the natural hinge loss for sign inconsistency possible. The linear and the hinge loss are covered by the pinball loss defined below:

(7)

where describes the slope in the negative part: when , it is the hinge loss and when , it reduces to the linear loss. The pinball loss has been applied for both regression ([25] [26]) and classification ([27] [28]). For M1bit-CS, it is possible to take any value from to for . The best one is problem-dependent and we will investigate its choice in numerical experiments.

In M1bit-CS model, the parameter controls the signal sparsity, and it exists also in regular CS models. This parameter can be selected according to prior-knowledge or by cross-validation.

is a trade-off between the analog and the saturated measurements. Heuristically, we set it as

such that is smaller when there are more saturated measurements. For M1bit-CSR (6), there is an additional parameter which is used to adjust the norm of the reconstructed signal. In 1bit-CS, the requirement on the signal norm is crucial. For M1bit-CS, it becomes less important since there are also analog measurements but still it helps the signal recovery.

Ii-C Iterative Saturation Detection

In the sensing system (1), if the observation is or , we know that it is saturated and use it in M1bit-CS as a one-bit measurement. But in some systems, itself has a lower/upper bound and then it is interesting and useful to distinguish whether a measurement reaches the bounds. In this paper, we focus on the lower-saturation and assume that the lower bound is zero. Discussions on upper-saturation and other bounding values are similar.

With prior-knowledge on non-negativeness of the measurements, we have

For example, in a CT scan system, is related to the integral of the tissue density along a X-ray. It is certainly non-negative and means that the corresponding X-ray hits the detector directly without crossing any tissue. If those observations with are found, we can set and use them as analog measurements in M1bit-CS. They will be helpful, especially for reconstructing the boundary.

The detection criterion is:

Accordingly, we can use the reconstructed result to update and then to solve M1bit-CS with the new guess of , iteratively. Note that incorrectly detecting a saturated measurement as a zero one is more harmful than regarding a zero measurement as a saturated one. Therefore, the initial guess is for all . For a lower-saturated and non-negative sensing system, the iterative saturation detection (ISD) scheme is described below:

  1. set saturation indicator ;

  2. signal recovery with , denote the result as ;

  3. calculate the measurements of , i.e., ;

  4. update

  5. set for ;

  6. go to 2) until there is no change on .

Iii Alternating Direction Method of Multipliers

In this paper, we use the least squares loss and the pinball loss for M1bit-CS (5) and (6). Then they are convex problems and many existing methods can be applied. We tried several methods including primal-dual algorithms [29] and different Alternating Direction Methods of Multipliers (ADMM, [30]) implementations [31]. In this paper, we introduce an implementation of ADMM, which is efficient in our experiments, but of course is not necessarily optimal for other tasks.

In order to apply ADMM, we reformulate (5) by introducing two auxiliary variables and , and the equivalent problem is

(8)
s.t.

where is an indicator function that returns 0 if and otherwise. Note here is -dimensional for consistency but only the components with have effect.

The corresponding augmented Lagrangian with the dual variables (the Lagrange multipliers: , ) is

where and . In every iteration, the ADMM updates and in a Gauss-Seidel style with fixed and updates the dual variable . Updating can be decoupled into updating and independently. We choose the updating order in our ADMM implementation.

  1. -subproblem: we have to solve

    of which the optimal solution is analytically given by

    Here is the shrinkage operator for the pinball loss:

  2. -subproblem: we have to solve

    of which the analytical solution is

    Here is the projection of a vector onto the ball :

  3. -subproblem: we have to solve

    It is essentially a quadratic programming plus an item:

    (9)

    It does not have analytical solutions. We approximately solve this problem using FISTA [32].

The ADMM for M1bit-CSC (5) is summarized in Algorithm 1. When the -subproblem is exactly solved, the convergence of Algorithm 1 follows the classical analysis of ADMM; see, e.g., [30]. Otherwise, the convergence results for inexact ADMM in [33] and [34] can be applied.

Input ; give parameters
Initialize   , repeat
       update by solving (3)
until stopping criteria is satisfied;
Algorithm 1 ADMM for M1bit-CSC (5)

We can establish ADMM in the same way for M1bit-CSR (6), and the algorithm is similar to Algorithm 1 except that the update of becomes solving the following problem:

of which the optimal solution is analytically given by

(10)

Iv Experimental Validation

In this section, we evaluate the proposed algorithms on one-dimensional sparse signal recovery problems. The parameters will be discussed numerically, and then the algorithms will be compared with two existing methods: i) simply dropping those saturated measurements and using the remaining measurements to recover the sparse signal by lasso; ii) Robust Dequantized Compressive Sensing (RDCS, [23]), which is designed for both saturation and quantization error. When there is no quantization error, RDCS takes the following formulation:

s.t. (11)

It assumes that there is no noise in the saturated measurements. Thus model (IV) is not suitable when there is noise in measurements before quantization. In fact, model (IV) is equivalent to model (5) with , , and . The ADMM algorithms for lasso and RDCS can be found in [30] and [23], respectively. All the experiments are done with Matlab 2014b on Windows 7 with Core i5-3.10 GHz and 8.0 GB RAM.

Let us consider an experiment with . Denote the true signal as

, of which the non-zero components are generated firstly following the standard Gaussian distribution and then normalized such that

. The sensing vectors are drawn from the standard Gaussian distribution. There is Gaussian noise in the measurements. The noise level is measured by the ratio of the summed squared magnitude of noise-free measurements to that of noise and denoted by . For the saturation part, we choose the saturation levels such that there are upper saturated measurements and lower saturated ones. Suppose that the reconstructed signal is , its quality is measured by the signal-to-noise ratio (SNR) in dB:

Iv-a Parameters Selection

The parameter in the pinball loss is important to characterize the M1bit-CS model. As discussed previously, is for 1bit-CS, is for classification problems. To numerically investigate the performance of different values, the following situations are considered:

  1. signal sparsity , number of measurements , number of saturated observations , noise level ;

  2. are as the same as 1), but ;

  3. are as the same as 1), but ;

  4. are as the same as 1), but .

We then apply Algorithm 1 to recover the signal and report the average SNR over 100 trials in Fig.2. From the result, one can observe that a small and negative performs well. The trend is that when the saturation ratio is large, the absolute value of the best is also large, which coincides with the experiments reported in 1bit-CS literatures: is better than in 1bit-CS where . In the rest of this paper, we always use . Additional tuning can improve the performance but requires more computational burden.

snr[c]SNR tau[c]

Fig. 2: Reconstructed SNRs for different values. CASE 1): blue solid curve; CASE 2): red dashed curve; CASE 3): black dotted curve; CASE 4): green dot-dashed curve.

For M1bit-CSR (6), there is an additional parameter which is used to adjust the norm of the reconstructed signal. In 1bit-CS, the requirement on the signal norm is crucial. For M1bit-CS, it becomes less important since there are also analog measurements but still it helps the signal recovery. We in the following consider different values for , and . The average SNR over 100 trials for different values are displayed by the blue solid curve in Fig.3, which also shows the norm of the reconstructed signals by the red dashed curve. As expected, the highest SNR is obtained when is near one. We in this paper use . Note that when there are fewer saturation data, the influence of the norm constraint is weaken. In that case, could be simply regarded as a regularization term.

SNR[c]SNR gamma[c] norm[c]

Fig. 3: Reconstruction performance for different values. SNR: blue solid curve; The norm of the reconstructed signal: red dashed curve.

Iv-B Synthetic Experiments

In this subsection, we compare M1bit-CS with lasso and RDCS on synthetic data. The sparsity parameter is tuned based on lasso, and the same value is then applied to other methods. We first set , , , and consider the recovery performance with different saturation ratios varying from to . The average results over trials are displayed in Fig.4. When the saturation ratio is low, there are plenty of analog measurements, then lasso can reconstruct the signal well. When increases, lasso’s recovery quality dramatically decreases. Saturated measurements are helpful to improve the recovery performance, which is the reason why RDCS, M1bit-CSC, and M1bit-CSR have significantly better results than lasso. Among the three algorithms, M1bit-CSC performs the best. Its advantage over RDCS is due to the facts that noise on saturated measurements are suppressed and the signal’s norm information is useful for recovery. If the norm of the true signal is unknown or cannot be well estimated, then we use M1bit-CSR, whose performance is slightly worse than M1bit-CSC. Similar performance could be observed in Fig.4, which is the average recovery results over 100 trails for different sparse levels with a fixed saturation ratio .

snr[c]SNR (dB) ratio[c]saturation ratio data1[l]lasso data2[l]RDCS data3[l]M1bit-CSC data4[l]M1bit-CSR

snr[c]SNR (dB) K[c]number of non-zero components data1[l]lasso data2[l]RDCS data3[l]M1bit-CSC data4[l]M1bit-CSR

Fig. 4: Recovery performance of lasso (black dotted), RDCS (blue dashed), M1bit-CSC (red solid line with triangle), and M1bit-CSR (red solid) for (a) different saturation ratios and (b) different numbers of non-zero components.

We then keep the signal sparsity as , the saturation ratio as , and vary the number of measurements from to . The average SNRs over 100 trials are reported in Fig.5. We can observe the effect of the proposed algorithms in mining knowledge from saturated measurements. For example, M1bit-CSC with measurements (640 analog measurements and 160 saturated ones) has a similar SNR as lasso with 1000 measurements (800 analog measurements are used and the left 200 ones are dropped). It means that saturated measurements have almost the same power as analog ones in recovering the sparse signal in that case and our algorithms are able to explore all the power with a longer computational time. Generally, the computational time of M1bit-CS is about 10 times of that of lasso, as reported in Fig.5. We think that when there are saturated measurements and if the computational time requirement is not critical, it is worthy to use M1bit-CS to improve the recovery accuracy.

snr[c]SNR (dB) m[c]number of measurements data1[l]lasso data2[l]RDCS data3[l]M1bit-CSC data4[l]M1bit-CSR

time[c]time (s) m[c]number of measurements data1[l]lasso data2[l]RDCS data3[l]M1bit-CSC data4[l]M1bit-CSR

Fig. 5: Recovery performance of lasso (black dotted), RDCS (blue dashed), M1bit-CSC (red solid line with triangle), and M1bit-CSR (red solid) for different numbers of measurements. (a) SNR; (b) computational time.

V Overexposure correction for CT Reconstruction

V-a Problem Formulation

Now we are at the stage to apply the proposed M1bit-CS to CT reconstruction to correct overexposure. In a CT system, the object density image can be represented in its discrete form as a vector , and the acquired projection data as . When there is no saturation, the X-ray transform of the object can be written as

(12)

where is the system matrix. Suppose that the intensity of one X-ray before entering the object is and the intensity departing the object is , then Lambert-Beer’s law tells us:

Due to the physical limits of the detector or the analog-to-digital converter, there is a maximal intensity and any intensity exceeding this threshold cannot be distinguished. In other words, we cannot observe the analog value when

In practice, the maximal and minimal detectable value of the detector can be adjusted at each projection angle for exposure control but the measurable dynamic range of the detector is fixed and denoted by , where and are the maximal and minimal intensities that are measured by the detector at projection angle . At each projection, is determined by the object, and is adjusted correspondingly afterwards by the fixed as . Then we can choose a dynamic threshold , where is the maximum projection value after the minus log transform at projection . Thus, the observation is a truncation of , i.e., , where takes the value of with being the angle containing the -th projection. Since the observation stands for the integral of the tissue intensity, it should be non-negative. In fact, there are amount of X-rays that hit the detector directly, i.e., the corresponding measurements is zero. Thus, when a measurement get saturated, people usually set its value as zero. Therefore, the observations in CT reconstruction are

When overexposure happens, traditional CT reconstruction methods are not suitable, as shown by Fig.1. To correct overexposure, the researchers typically considered extrapolation to compensate the missing information; see e.g., [35] [36] [37]. Among them, the method given by [38] that assumes the missing data as line integrals of a partial water cylinder is efficient. It is based on FBP algorithm and in the following, we refer this method as FBP-WCE (FBP with water cylinder extrapolation). Its reconstruction result is displayed in Fig.6. One can see significant improvement from the standard FBP method, of which the result is given by Fig.1.

V-B Knee Phantom

The overexposure phenomenon is a kind of one-sided saturation. If is known, we can directly apply M1bit-CSR (6) to do overexposure correction (we do not know the norm of the image and M1bit-CSC is not suitable). Note that the sparsity for CT reconstruction is on the gradient domain, i.e., total variation (TV) is used as the regularization term; see [39] [40] [41] for TV regularization CT reconstruction. When applying M1bit-CS on CT reconstruction, we do not use FISTA but a TV minimization algorithm to solve the -subproblem (3). The reconstruction image is shown in Fig.6, which preliminarily shows that acquiring information from saturated projections largely helps CT reconstruction.

Fig. 6: Reconstructed results from saturated projection Fig.1 (a) FBP-WCE; (b) M1bit-CSR.

In practice, when , we do not know whether this projection ray is saturated or does not hit the object. Thus, the ISD scheme proposed in Section II-C should be used together with M1bit-CSR (abbreviated as M1bit-CSR-ISD). Again consider the projections shown in Fig.1. Fig.7 displays the detected . The difference between the ideal and the detected is shown in Fig.7, where the false detections (zero measurements that are detected as saturated ones) are marked in blue and missing detection (saturated measurements that are regarded as zero ones) are in green.

Fig. 7: Saturation indicator for saturated data shown in Fig.1: (a) the true ; (b) indicator detected by M1bit-CSR-ISD; (c) difference between the true and the detected saturation indicator (blue: false detection; green: missing detection).

The final reconstruction image is given by Fig.8, which shows that a small number of incorrect detections can be tolerated by M1bit-CSR (6) and the reconstruction result is still satisfactory. As plotted in Fig.8, root of mean square error (RMSE, in HU) is decreasing with the increasing of ISD iterations. Since ISD requires to solving M1bit-CSR multiple times, the detection process is time-consuming. Utilizing prior-knowledge may reduce the detection time, which could be an interesting study in the future.

RMSE (HU)[c]RMSE (HU) ISD iteration[c]ISD iteration

Fig. 8: (a) Reconstructed result of M1bit-CS-ISD from saturated projections shown in Fig.1; (b) RMSE (HU) in different ISD iterations.

V-C Simulated Clinical Data

Finally, we apply the proposed method on a clinical head dataset. The data are acquired with a Siemens Artis zee angiographic C-arm system (Siemens Healthcare GmbH, Forchheim, Germany). In this experiment, we choose one slice of a 3D clinical head dataset as the ground truth image (Fig.9) and reproject it to simulate the acquired sinogram data in the fan-beam system with the following trajectory parameters: the source-to-isocenter distance is 750 mm and isocenter-to-detector distance is 450 mm. The angular step is 1 degree and the total scan range is 360 degrees. The equal-spaced detector length is 620 mm with the pixel length mm.

The full projections are shown in Fig.10. When the detection range is limited, there could be saturation for the projection. In Fig.10 and Fig.10, observations for and are displayed, respectively. Our task is to recover the image from the saturated projections via M1bit-CS-ISD. The results are compared with FBP and SART, two standard CT reconduction frameworks. For FBP, we apply the modification given by [38] that utilizes water cylinder extrapolation to remedy missing projections caused by truncation or overexposure. For SART, one can remove those saturated projections when they are found, for which the ISD can be used as well. We denote this method as SART-ISD, of which the detection scheme is as the same as M1bit-CSR-ISD but SART is used as the reconstruction method.

Fig. 9: Reconstruction results for the clinical data (): (a) ground truth; (b) FBP-WCE; (c) SART-ISD; (d) M1bit-CSR-ISD.
Fig. 10: Projections of Fig.9: (a) full projection data; (b) saturated projection data with ; (c) saturated projection data with .

For , the reconstruction results of FBP-WCE and SART-ISD are given in Fig.9 and 9, respectively. As shown before, the traditional FBP method cannot handle the saturated data. With water cylinder extrapolation, the reconstruction quality has been improved but loss of clear patient boundaries still happens. The overall performance of SART-ISD is slightly better than FBP-WCE but capping artifact can be identified at the object border. Further improvement is obtained using the proposed M1Bit-CSR-ISD to acquire information from the saturated data. As shown in Fig.9, most of outer boundaries are nicely restored and streaking artifacts are effectively eliminated.

Next, we artificially add Gaussian noise on the projections and the standard deviation of the noise

. The reconstructed images of FBP-WCE (with a smooth filter), SART-ISD, and M1bit-CSR-ISD are displayed in Fig.11. Comparing the images in Fig.9 and 11, we can find that the proposed method is not sensitive to additive noise, which is rooted in the loss functions of (5) and (6

). If the projections contain outliers, then some robust losses can be applied.

Fig. 11: Reconstruction results from the noise-corrupted and saturated projections ( and ): (a) ground truth; (b) FBP-WCE; (c) SART-ISD; (d) M1bit-CSR-ISD.

When the detectable range is smaller such that , saturation becomes worse. Then the performance of FBP and SART, even with water cylinder extrapolation and the ISD scheme, dramatically drops. In practice, this heavy saturation rarely happens and is out of the scope of FBP-WCE. But even with this heavy saturation, the proposed M1bit-CSR-ISD can output a good result. The reconstructed images and enlargements of a region are illustrated in Fig.12 and Fig.13, respectively. In this case, both FBP-WCE and SART-ISD fail to restore the clear outer boundaries of the patient, while M1bit-CSR-ISD is still able to achieve this in a more accurate manner. We further report the saturation detection result in Fig.14, from which one can observe that most of the saturations have been properly detected.

Fig. 12: Reconstruction results for the clinical data (): (a) ground truth (yellow rectangle is enlarged in Fig.13); (b) FBP-WCE; (c) SART-ISD; (d) M1bit-CSR-ISD.
Fig. 13: An enlargement of a region of the clinical data, which is marked by yellow rectangle in Fig.12: (a) ground truth; (b) FBP-WCE; (c) SART-ISD; (d) M1bit-CSR-ISD.
Fig. 14: Saturation indicator for the clinical phantom with . (a) the true saturation indicator; (b) detected by M1bit-CSR-ISD; (c) difference between the true and the detected saturation indicator (blue: false detection; green: missing detection).

In Table I, we list the RMSE (in HU) for FBP, FBP-WCE, SART-ISD, M1bit-CSR with ISD and with an ideal saturation matrix. Compared with classical FBP, FBP-WCE and SART-ISD have significant improvements. But the schemes are basically different: FBP-WCE applies water cylinder extrapolation as the complements of the missing projections, and SART-ISD is to iteratively detect and drop incorrect projections. The proposed M1bit-CSR can adequately use the information carried by overexposed measurements and then outputs good reconstruction results. The gap between M1bit-CSR-ISD and M1bit-CSR with the ideal is not very big, which on the one hand shows the accuracy of the proposed ISD method, and on the other hand confirms the robustness of M1bit-CSR against a part of incorrect detection.

FBP FBP SART M1bit-CSR M1bit-CSR
data WCE ISD ISD (ideal )
knee phantom
() 182.4 97.93 54.91 8.547 5.461
head
() 43.22 35.34 31.97 11.62 6.930
head with noise
() 43.84 36.71 34.22 11.72 7.001
head
() 132.2 87.93 70.63 14.81 12.39
TABLE I: Root of Mean Square Error (HU) of Reconstruction Results

In M1bit-CS models, the saturation indicator is a key input. In this paper, we designed a simple but efficient iterative method to detect . Further improvement may come from the use of prior-knowledge and reasonable assumptions on tissue structures. At least, those prior-knowledge provides a good initial guess. An accurate initial guess can reduce both the reconstruction error and the detection time.

Vi Conclusion

Aiming at signal reconstruction when saturation happens, we established mixed one-bit compressive sensing that could be regarded as a bridge between the regular and one-bit CS. We developed the corresponding ADMM and an iterative saturation detection method. They are then successfully applied to overexposure correction for C-arm CT reconstruction. The improvement from the existing methods is quite significant. Generally, the reconstruction performance of the proposed method is satisfactory even in the presence of severe detector saturation, showing the benefit of our method on considering saturated projections.

CT scanning that has overexposure is a typical example of sensing systems with saturation. The promising performance of M1bit-CS implies its potential use on other tasks with measurements beyond the reliable sensing region. We can even take advantage of saturation phenomena for specific purposes. For example, one can sample signs for a part of the measurements to reduce the expense of the analog-to-digital conversion without significantly sacrificing the recovery accuracy. In low-dose computed tomography, when we largely reduce the radiation dose to the level such that some projections are below the threshold, then those projection values become unreliable but they do provide one-bit information. Then M1bit-CS can also be applied to improve the performance.

Acknowledgment

The authors would like to thank Prof. Ji Liu from the University of Rochester for sharing the code of RDCS.

References

  • [1] J. S. Feng, “Sensor saturation effects on NLTS measurements [of magnetoresistive sensors],” IEEE Transactions on Magnetics, vol. 34, no. 4, pp. 1952–1954, 1998.
  • [2] Y. Fang, Y. Zhang, N. Qi, and X. Dong, “AM-AFM system analysis and output feedback control design with sensor saturation,” IEEE Transactions on Nanotechnology, vol. 12, no. 2, pp. 190–202, 2013.
  • [3] E. dos Santos, G. Cardoso, P. Farias, and A. de Morais, “CT saturation detection based on the distance between consecutive points in the plans formed by the secondary current samples and their difference-functions,” IEEE Transactions on Power Delivery, vol. 28, no. 1, pp. 29–37, 2013.
  • [4] N. Strobel, O. Meissner, J. Boese, T. Brunner, B. HeiGl, M. Hoheisel, G. Lauritsch, M. NaGel, M. Pfister, E.-P. Rührnschopf, B. Scholz, B. Schreiber, M. Spahn, M. Zellerhoff, and K. Klingenbeck-Regn, “3D imaging with flat-detector C-arm systems,” in Multislice CT.   Springer, 2009, pp. 33–51.
  • [5] L. Feldkamp, L. Davis, and J. Kress, “Practical cone-beam algorithm,” Journal of the Optical Society of America A, vol. 1, no. 6, pp. 612–619, 1984.
  • [6] A. H. Andersen and A. C. Kak, “Simultaneous algebraic reconstruction technique (SART): A superior implementation of the ART algorithm,” Ultrasonic Imaging, vol. 6, no. 1, pp. 81–94, 1984.
  • [7] A. Preuhs, M. Berger, Y. Xia, A. Maier, J. Hornegger, and R. Fahrig, “Over-exposure correction in CT using optimization-based multiple cylinder fitting,” in Bildverarbeitung für die Medizin 2015, H. Handels, Ed., 2015, pp. 35–40. [Online]. Available: https://www5.informatik.uni-erlangen.de/Forschung/Publikationen/2015/Preuhs15-OCI.pdf
  • [8] N. Mail, D. Moseley, J. Siewerdsen, and D. Jaffray, “The influence of bowtie filtration on cone-beam CT image quality,” Medical Physics, vol. 36, no. 1, pp. 22–32, 2009.
  • [9] J.-H. Choi, R. Fahrig, A. Keil, T. F. Besier, S. Pal, E. J. McWalter, G. S. Beaupré, and A. Maier, “Fiducial marker-based correction for involuntary motion in weight-bearing C-arm CT scanning of knees (Part I): Numerical model-based optimization,” Medical Physics, vol. 40, no. 9, p. 091905, 2013.
  • [10] J.-H. Choi, A. Maier, A. Keil, S. Pal, E. J. McWalter, G. S. Beaupré, G. E. Gold, and R. Fahrig, “Fiducial marker-based correction for involuntary motion in weight-bearing C-arm CT scanning of knees (Part II): Experiment,” Medical physics, vol. 41, no. 6, p. 061902, 2014.
  • [11] J. Rausch, A. Maier, R. Fahrig, J.-H. Choi, W. Hinshaw, F. Schebesch, S. Haase, J. Wasza, J. Hornegger, and C. Riess, “Kinect-based correction of overexposure artifacts in knee imaging with C-Arm CT systems,” International Journal of Biomedical Imaging, vol. 2016, 2016.
  • [12] P. T. Boufounos and R. G. Baraniuk, “1-bit compressive sensing,” in the 42nd Annual Conference on Information Sciences and Systems.   IEEE, 2008, pp. 16–21.
  • [13]

    E. J. Candès and T. Tao, “Decoding by linear programming,”

    IEEE Transactions on Information Theory, vol. 51, no. 12, pp. 4203–4215, 2005.
  • [14] D. L. Donoho, “Compressed sensing,” IEEE Transactions on Information Theory, vol. 52, no. 4, pp. 1289–1306, 2006.
  • [15] E. J. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information,” IEEE Transactions on Information Theory, vol. 52, no. 2, pp. 489–509, 2006.
  • [16] J. N. Laska, Z. Wen, W. Yin, and R. G. Baraniuk, “Trust, but verify: Fast and accurate signal recovery from 1-bit compressive measurements,” IEEE Transactions on Signal Processing, vol. 59, no. 11, pp. 5289–5301, 2011.
  • [17] M. Yan, Y. Yang, and S. Osher, “Robust 1-bit compressive sensing using adaptive outlier pursuit,” IEEE Transactions on Signal Processing, vol. 60, no. 7, pp. 3868–3875, 2012.
  • [18] L. Jacques, J. N. Laska, P. T. Boufounos, and R. G. Baraniuk, “Robust 1-bit compressive sensing via binary stable embeddings of sparse vectors,” IEEE Transactions on Information Theory, vol. 59, no. 4, pp. 2082–2102, 2013.
  • [19]

    Y. Plan and R. Vershynin, “Robust 1-bit compressed sensing and sparse logistic regression: A convex programming approach,”

    IEEE Transactions on Information Theory, vol. 59, no. 1, pp. 482–494, 2013.
  • [20] L. Zhang, J. Yi, and R. Jin, “Efficient algorithms for robust one-bit compressive sensing,” in

    Proceedings of the 31st International Conference on Machine Learning

    , 2014, pp. 820–828.
  • [21] R. Zhu and Q. Gu, “Towards a lower sample complexity for robust one-bit compressed sensing,” in Proceedings of the 32nd International Conference on Machine Learning, 2015, pp. 739–747.
  • [22] J. N. Laska, P. T. Boufounos, M. A. Davenport, and R. G. Baraniuk, “Democracy in action: Quantization, saturation, and compressive sensing,” Applied and Computational Harmonic Analysis, vol. 31, no. 3, pp. 429–443, 2011.
  • [23] J. Liu and S. J. Wright, “Robust dequantized compressive sensing,” Applied and Computational Harmonic Analysis, vol. 37, no. 2, pp. 325–346, 2014.
  • [24] Y. C. Eldar and G. Kutyniok, Compressed Sensing: Theory and Applications.   Cambridge University Press, 2012.
  • [25]

    A. Christmann and I. Steinwart, “How SVMs can estimate quantiles and the median,” in

    Advances in Neural Information Processing Systems, 2007, pp. 305–312.
  • [26] I. Steinwart and A. Christmann, “Estimating conditional quantiles with the help of the pinball loss,” Bernoulli, vol. 17, no. 1, pp. 211–225, 2011.
  • [27]

    X. Huang, L. Shi, and J. A.K. Suykens, “Support vector machine classifier with pinball loss,”

    IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 36, no. 5, pp. 984–997, 2014.
  • [28] Y. Xu, Z. Yang, and X. Pan, “A novel twin support-vector machine with pinball loss,”

    IEEE Transactions on Neural Networks and Learning Systems

    , 2016.
  • [29] A. Chambolle and T. Pock, “A first-order primal-dual algorithm for convex problems with applications to imaging,” Journal of Mathematical Imaging and Vision, vol. 40, no. 1, pp. 120–145, 2011.
  • [30] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Foundations and Trends in Machine Learning, vol. 3, no. 1, pp. 1–122, 2011.
  • [31] M. Yan and W. Yin, “Self equivalence of the alternating direction method of multipliers,” in Spllitiing Methods in Communication and Imaging, Science and Engineering, R. Glowinski, S. Osher, and W. Yin, Eds.   Springer-Verlag, 2016.
  • [32] A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM Journal on Imaging Sciences, vol. 2, no. 1, pp. 183–202, 2009.
  • [33] M. K. Ng, F. Wang, and X. Yuan, “Inexact alternating direction methods for image recovery,” SIAM Journal on Scientific Computing, vol. 33, no. 4, pp. 1643–1668, 2011.
  • [34] L. Shen and S. Pan, “Inexact indefinite proximal ADMMs for 2-block separable convex programs and applications to 4-block DNNSDPs,” arXiv preprint arXiv:1505.04519, 2015.
  • [35] B. Ohnesorge, T. Flohr, K. Schwarz, J. P. Heiken, and J. P. Bae, “Efficient correction for CT image artifacts caused by objects extending outside the scan field of view,” Medical Physics, vol. 27, no. 1, pp. 39–46, 2000.
  • [36] A. Maier, B. Scholz, and F. Dennerlein, “Optimization-based extrapolation for truncation correction,” in 2nd CT Meeting, 2012, pp. 390–394.
  • [37] Y. Xia, H. Hofmann, F. Dennerlein, C. Schwemmer, S. Bauer, G. Chintalapani, P. Chinnadurai, J. Hornegger, and A. Maier, “Towards clinical application of a Laplace operator-based region of interest reconstruction algorithm in C-arm CT,” IEEE Transactions on Medical Imaging, vol. 33, pp. 593–606, 2014.
  • [38] J. Hsieh, E. Chao, J. Thibault, B. Grekowicz, A. Horst, S. McOlash, and T. J. Myers, “A novel reconstruction algorithm to extend the CT scan field-of-view,” Medical Physics, vol. 31, no. 9, pp. 2385–2391, 2004.
  • [39] H. Yu and G. Wang, “A soft-threshold filtering approach for reconstruction from a limited number of projections,” Physics in Medicine and Biology, vol. 55, no. 13, p. 3905, 2010.
  • [40] L. Ritschl, F. Bergner, C. Fleischmann, and M. Kachelrieß, “Improved total variation-based CT image reconstruction applied to clinical data,” Physics in Medicine and Biology, vol. 56, no. 6, p. 1545, 2011.
  • [41] Z. Chen, X. Jin, L. Li, and G. Wang, “A limited-angle CT reconstruction method based on anisotropic TV minimization,” Physics in Medicine and Biology, vol. 58, no. 7, p. 2119, 2013.