Shape recovery from sparse tomographic X-ray data
A two-dimensional tomographic problem is studied. The target is assumed to be a homogeneous object bounded by a smooth curve. A Non Uniform Rational Basis Splines (NURBS) curve is used as computational representation of the boundary. This approach conveniently provides the result in a format readily compatible with computer-aided design (CAD) software. However, the linear tomography task becomes a nonlinear inverse problem due to the NURBS-based parameterization. Therefore, Bayesian inversion with Markov chain Monte Carlo (MCMC) sampling is used for calculating an estimate of the NURBS control points. The reconstruction method is tested with both simulated data and measured X-ray projection data. The proposed method recovers the shape and the attenuation coefficient significantly better than the baseline algorithm (optimally thresholded total variation regularization), but at the cost of heavier computation.
- 1 Introduction
2 Materials and Methods
- 2.1 NURBS
- 2.2 Tomographic Measurement Model
- 2.3 Tomographic Projection Data
- 2.4 Bayesian Inversion
- 2.5 Markov chain Monte Carlo
- 2.6 Total Variation regularization
- 2.7 Optimal TV parameter choice
- 2.8 Simulated phantoms
- 2.9 Physical phantoms
- 3 Results
- 4 Discussion
- 5 Conclusions
We propose a new reconstruction algorithm for two-dimensional X-ray tomography of objects with homogeneous attenuation coefficient. The basic idea is to represent the boundary of the object in terms of a non-uniform rational basis spline (NURBS) curve, thereby providing a direct connection to computer-aided design (CAD) software [5, 8, 23, 29]. The control points of the parametric curve are the degrees of freedom in the inverse problem, which then becomes nonlinear. Therefore, we resort to Bayesian inversion and Markov chain Monte Carlo (MCMC) sampling [13, 16, 27, 33] for estimating the unknown curve from projection data.
Why not reconstruct the object with a traditional reconstruction method, such as filtered back-projection, and fit a NURBS curve to the boundary of the reconstructed domain? This indeed is quicker and more reliable in the case of a comprehensive tomographic dataset with dense angular sampling. However, often there are time constraints, geometric obstructions, or radiation dose issues preventing the collection of a detailed dataset. With such cases in mind, we test our algorithm with projection data collected from only six directions around the object. It turns out that the combination of NURBS and MCMC outperforms many traditional methods in this extremely ill-posed inverse problem. See Figure 1.
The proposed MCMC-NURBS reconstruction method is first tested with two simulated phantoms, one of them convex and the other non-convex. See [24, 25] for other examples. The X-ray attenuation coefficient is constant inside the phantoms, and we simulate six projection images using fan-beam geometry. The projection directions are uniformly distributed over the full circle. We use Delayed Rejection Adaptive MCMC (DRAM) approach for effective Monte Carlo computation [3, 11, 27]. The degrees of freedom in the inverse problem are the coordinates of the control points of the NURBS curve and the unknown attenuation coefficient.
We use Total Variation (TV) regularization as a baseline method for comparison and quantitative assessment of tomographic reconstruction quality. The degrees of freedom in the inverse problem are the non-negative pixel values. To be as fair as possible to the baseline method, we choose the regularization parameter optimally by comparing the reconstruction to the true simulated phantom. We force the reconstruction to be piecewise constant (with only one possible positive attenuation value) by optimally thresholding the TV regularized solution.
Furthermore, we test both methods using X-ray projection data measured from physical objects, similarly shaped than the simulated phantoms.
The proposed method is found to outperform the optimized TV regularization algorithm in all the cases studied. See Figure 1 for the example concerning the physical non-convex phantom.
NURBS curves have been used as building blocks in CAD software for decades. Some works in optimizing the NURBS representation in certain inverse problems, including reverse engineering, fitting strategies, and recovering shapes from photographs, have been done in [1, 4, 12, 14, 18, 26, 32, 31]. Other works in X-ray tomography that readily provide segmented images can be found in [19, 20, 15]. Non-destructive testing from very restricted data using Bayesian inversion in terms of a vector-graphic format are discussed in [19, 20, 16].
Also, there are various branches of computational science using NURBS curves representation. For example, isogeometric analysis provides a computational approach for integrating finite element analysis (FEA) and CAD . One of the applications of NURBS-based isogeometric analysis is computing flows from spinning propellers in opposite direction. In this case, a discretization using NURBS gives more accurate results in computing the flows compared to standard finite elements .
The method proposed here is different from all the previous ones listed above.
2 Materials and Methods
We consider an unknown physical domain and model the object boundary defined by a parametrized curve . In our computational problem, a NURBS curve as a piecewise rational function defined on , is introduced. We divide the interval [0,1] into pieces and called breakpoints which map into the endpoints of the polynomial segment.
We introduce the curve representation as follows:
where the points configure the curve shape. The points are called control points and is the following rational function of degree p:
where are nonnegative weights for all and are basis functions that describe how strongly the control points attract the NURBS curve. They are defined recursively as
and for as
where by definition. A collection of breakpoints is then called knot vector:
If are evenly spaced, then is called a periodic uniform knot vector. In our discussion, we implement a closed NURBS curve using periodic uniform knot vector to recover the boundary of the object. By repeating the first control points after the last point, an unclamped closed NURBS curve is obtained.
Two examples of closed NURBS curves of rd-degree basis functions with seven control points are given. Figure 2 illustrates the closed curve with the basis functions defined on the same periodic uniform knot vector:
In those examples and in our computations, we assume that the weights corresponding to all of the control points are the same.
2.2 Tomographic Measurement Model
Consider a physical domain and a continuous tomography model where as
For computational reasons, we need to construct two discrete models: a pixel-based object model and a NURBS-based object model as in Figure 4. Pixel-based object model is a discretization of the continuous model into pixels in two dimensional image, so then it discretizes a line integral to a standard pencil-beam model. NURBS-based object model is a model where the boundary of the object is represented in terms of a NURBS curve. We represent the control points in polar coordinates where is the number of control points:
Consider a vector and construct , where is an attenuation parameter.
The line integral is discretized by switching to pixel-based model using an operator where NN is the resolution of the pixel image. Define
Let consider as in Figure 5.
In the pixel-based model, the measurement is defined as follows:
where is a distance that the line travels in the jth pixel. We use the collection of lines arising from several fan-beam projections as shown in Figure 3. Our measurement model can be written as follows.
2.3 Tomographic Projection Data
The x-ray tomography data of the sugar phantom were acquired with the custom-built CT device nanotom 180 supplied by Phoenix| Xray Systems + Services GmbH (Wunstorf, Germany). The chosen geometry resulted in a magnification with resolution of 62.2 m/pixel for physical phantom and 63.8 m/pixel for physical phantom . The phantoms can be seen in Figure 9. The x-ray detector is a 12-bit CMOS flat panel detector with pixels of 100 m size (Hamamatsu Photonics, Japan). A set of 120 projection images were acquired over a full 360 degree rotation with uniform angular step of 3 degrees between projections. Each projection image was composed of an average of six 250 ms exposures. The x-ray tube acceleration voltage was 80 kV and tube current 300 A, and the full polychromatic beam was used for image acquisition. For this work we choose the projections corresponding to the middle cross-section of a sugar phantom as Figure 9, thus the task is only 2D problem. We picked 6 projections from the measured data with uniform angular sampling from a total opening angle of 360 degrees.
2.4 Bayesian Inversion
In our problem, we have a sparse tomographic model which leads to ill-posed inverse problem. To compensate this ill-posedness, a priori knowledge of the target should be explored well. For computational reasons, the information needs to be transformed to quantitative form. Bayesian inversion provides a flexible way to synthesize this extra information of the target .
The main idea of this approach is to recast the inverse problem as a problem of Bayesian inference. We use probability theory to model our lack of information in the inverse problem. All the variables in the model are considered as random variables.
There are three main steps in constructing the problem to Bayesian framework. Firstly, gathering information prior to the measurement and construct it as a prior density. Secondly, constructing a likelihood function that expresses how likely the observation outcomes with unknown parameter given. Thirdly, constructing and exploring the posterior probability density as a solution of the inverse problem.
In our case, the tomographic model leads to the nonlinear inverse problem with additive Gaussian errors as follows:
where is a matrix model of the measurement and m is the x-ray measurement.
In this part, we construct a prior density which will contribute to remove the ill-posedness. In this case, we have a homogeneous object with diameter not more than , where is determined as a radius from the origin. This condition means that the control points satisfy . They are shown as the dots in the Figure 7. The dashed circle is the NURBS curve produced by the control points with a fixed uniform periodic knot vector.
Assume that the prior information is Gaussian distributed with variance . Then a priori information has the following quantitative form
where , .
In addition, to avoid uncontrollable movement of control points, where they can rotate only to one side and get trapped, we condition the range of the control points movement by setting as an angle of each control point:
where is a lower bound for and is an upper bound for . The condition of the radius length is provided as well as: where is a constant. To minimize oscillations in the curve, the following condition should be satisfied
where are the control points as in Figure 6 and is a constant.
To avoid a control polygon intersection, another hard prior setting is added. Assume that is another polygon in the Figure 6, then the self-intersection can be avoided if holds.
We have , so then
From Equation 7 we then model the measurement process as:
and it is called likelihood function, where is a normalization constant.
The solution of the inverse problem is the posterior probability distribution which has the following form:
where denotes equality up to a normalization constant.
To get an estimation of , we use the conditional mean method that is averaging the values of the generated samples:
Since we face integration problem, Markov chain Monte Carlo (MCMC) technique is proposed to solve this problem. In the following section, details of MCMC method are discussed.
2.5 Markov chain Monte Carlo
Let us consider the Equation 10. To handle the integration problem, MCMC approach is applied by generating samples from the posterior distribution. In probability theory, the law of large numbers ensures that:
Hence, the conditional mean (CM) estimate for can be written as follows:
This CM estimate is the result of the recovered control points and attenuation coefficient.
For , give an initial state .
Set the proposal state , where
If then set .
Draw a random number from uniform distribution on . If then set , else set .
If then stop; else set and go to step.
To improve the efficiency of the Metropolis Hastings algorithm: Delayed Rejection Adaptive Metropolis (DRAM) is proposed [3, 11, 30, 35]. DRAM is a method that combines two powerful ideas: Delayed rejection (DR) and Adaptive Metropolis (AM). In DR, if the proposal state is rejected, instead of remaining in the same position a second stage is proposed. The strategy of AM combining to DR method allows the covariance matrix calibrated using n-samples path of the MCMC chain regardless at which stage these samples are accepted in DR algorithm.
2.6 Total Variation regularization
In this tomographic model, we have only a homogeneous object and a sharp boundary that divides the background and the domain of the object. One of the well known methods to solve this problem is total variation (TV) regularization which produces edge-preserving reconstruction. Let us denote . The solution of TV-regularized is defined by finding the vector that minimizes the penalty functional:
where is a regularization parameter and TV is a total variation defined as follows.
To apply that in the computation, 2D adaptation of the quadratic programming approach is implemented. See  Section 6.2.
2.7 Optimal TV parameter choice
In this subsection, we discuss about how to choose an optimal regularization parameter for TV in terms of a piecewise constant image. Because TV is the baseline method in this study, we want it to perform optimally. Therefore, we choose the regularization parameter in an unrealistically effective way, using the knowledge of the true (simulated) attenuation coefficient.
In the case of simulated data, we implement the thresholding method to yield the result as a binary image by choosing the optimal regularization parameter and the attenuation threshold , where
Let us write , where is the optimal regularization parameter. Define as thresholded TV reconstruction, where
The error of the reconstructions is computed as follows. Denote as the image of the target 2D object and as the image of the reconstruction. Set for points that belong to the original object but not to the reconstruction and for points that belong to the reconstruction but not to the original object. The relative error in the reconstruction is written as
For every , the relative error is measured for all . Eventually, we pick the regularization parameter that contains the minimum error from each thresholding value .
2.8 Simulated phantoms
In this part, two simulated phantoms are presented. We build images with one homogeneous convex shape and one homogeneous non-convex shape as in Figure 8.
Each phantom has the image resolution and they are built without using NURBS to avoid an inverse crime. The images have only two values: 0 for the background and 0.027 for the inner shape.
2.9 Physical phantoms
In this paper, we also present the reconstruction of real objects with the same shape as the simulated phantoms. One of the physical phantom is simply made by cardboard while another phantom is plastic-printed using 3D printer machine. Both are then filled by white crystal sugar as shown in Figure 9.
Throughout the chapter, we use the following measurement setup: sparse-angle (6 angles) from full angles of fan beam geometry are presented.
3.1 Simulated phantom reconstruction
The projection data of both simulated phantoms, and , is corrupted by Gaussian noise each. In the real projection data, we normalize it as . Based on the a priori information that the target is solid homogeneous object, we set the small values (less than ) in the to be zero. This intends to be an air background.
Prior to implementing the MCMC algorithm, some parameters need to be fixed. Based on the experiments, to create the perfectly target shapes of and in terms of NURBS curve, the minimum number of control points that are needed are 6 and 12, respectively. Those numbers of control points are then proceeded as our parameters of interest for both simulated and real problems. The periodic uniform knot vector is used and the weights are set equally the same for all control points. As discussed in Subsection 2.2, the coordinates of the control points are defined in the polar coordinates. Since the parameter of the attenuation value needs to be recovered as well, therefore the total number of parameters of interests for and are 13 and 25, respectively. In DRAM scheme, as a rule of thumb, the length of non-adaptation period in low dimensional problems, , is fixed to be .
The simulations are performed with a MATLAB code on the machine equipped with 2.9 GHz Intel Core i5 CPU and 8 GB memory. The image resolution in simulation run is set to have a size of pixels for the sake of computing speed. A unit disc with the radius over the resolution width and the attenuation value are fixed as starting points. A Gaussian distribution with centre at is set for the prior where n is the number of control points. For both targets and , for and . The angles are set to be for phantom and for phantom .
The total number of evaluations is 6 000 000. The histograms and the chains for simulated data reconstructions show relatively the same behaviour as those of the physical data. Therefore the figures are omitted and the histograms and the MCMC chains of physical data reconstruction are given in the Subsection 3.2.
For comparison, the TV regularizations for the 2D tomographic case using quadratic programming as well as the thresholded-TV reconstructions are presented. See Figures 10 and 11. The optimal TV parameter choice is calculated as discussed in Subsection 2.7. In the thresholded-TV reconstructions, we choose the range of regularization parameter range of and the attenuation threshold value of . The optimal parameters choices for each data are given in Table 1 (Table 2 is the optimal parameters choices for physical problem). The absolute relative errors are presented as well as a comparison to CM estimate of the attenuation value, .
|Thresholded Total Variation||MCMC-NURBS|
The simulation time for each run is given in Table 6.
3.2 Physical phantom reconstruction
The geometry of tomography projection of the real data is given in Subsection 2.3. The construction of the prior knowledge is similar with simulated case. The total number of evaluations is 4 500 000. Figure 12 and Figure 15 present the histograms of the 1-d marginal posterior distribution of radii from physical data phantom and angles from physical data , respectively. The histograms show the distribution of the values of the samples in the MCMC chain. The chains of the angles from physical data phantom and the chains of the radii from physical data phantom are presented in Figure 13 and Figure 14. All the figures are presented after omitting the burn-in period. The TV regularizations for the 2D tomographic case using quadratic programming as well as the thresholded-TV reconstructions are presented in Figures 17 and 18.
We reconstruct both simulated shapes and outlines of real objects from extremely sparsely collected tomographic data: we only use 6 projection directions (sparsely) spanning the full 360. In our new method we represent the boundaries of homogeneous attenuating objects as a NURBS curve, whose control points are degrees of freedom in the inverse problem. The Bayesian inversion is done by sampling the posterior distribution with an MCMC method. By referring to Section 2.2, another degree of freedom, the attenuation coefficient is set as well as our parameter of interest. The final result, in the form of the NURBS curve, is obtained by computing the conditional mean (CM estimate) of the chains of control points and the attenuation coefficient.
Our baseline method for comparison is Total Variation regularization followed by a thresholding step to yield a homogeneous object. To be as fair as possible to the comparison method, we chose both the regularization parameter and the binary-image threshold optimally.
Let us discuss the features of our Bayesian inversion computation. The MCMC chains of all radii and angles seem to be mixing very well. The convergence is determined by visually looking at one- and two-dimensional chain plots and by Geweke’s convergence diagnostic .
Some of the histograms in Figures 12 and 15 are multimodal and clearly show the non-Gaussian nature of the posterior distribution as a result of our nonlinear NURBS-based parametrization of the unknown. Estimating the full posterior distribution by MCMC allows us to observe nonlinearities such as the multimodality of marginal distributions as shown in Figure 12, 15 and 16. We can quantify the uncertainty in the results as well. All of this comes in addition to getting useful reconstructions.
After computing the NURBS curve from the CM estimate of the control points, the result is given as a binary image as an implementation of our tomographic measurement model. Although in the reconstruction process we used a grid, the final image can be conveniently presented at higher resolution, in this case as grid. This is an advantage of the vector-graphic approach we use.
The reconstructions of the simulated and the physical phantoms using the MCMC-NURBS approach produce comparably similar shapes and the attenuation value as the targets. The retrieval of the attenuation parameters of the proposed method for reconstructing the simulated and physical phantoms yields the relative errors of and , respectively. See Table 3. Curiously, the histograms of attenuation value chains, shown in Figure 16, are close to Gaussian distributions, in constrast to those of the control points.
As always in Bayesian inversion, the design of the prior distribution is crucial. In this research we spent considerable time for specifying a prior that
is simple to write down as a mathematical formula, and
enables good enough recovery of cavities the non-convex example shapes we use.
The prior we ended up using is a trade-off between several aspects: mathematical simplicity, not too many control points, and flexibility in representing non-convex shapes. For example, in Figure 19, the reconstructions from several highest posterior are presented and depicted the change of MCMC chain. In particular, if we look closer to the upper-left of the reconstruction, some of the chain form a spiky shape. This behaviour is allowed by the prior construction which gives the flexibility of the control points to form the cavities. The spiky appearance, of course, can be removed by applying a stricter prior (reduce the possibility of the oscillation more), but consequently, the control points will not form the cavities well. In the case of , this is how the cavities are represented: the chains of , , and , combined with , , and , form the bottom cavity. The same holds for , and , and for , and , in the case of the upper cavity. The prior settings allow the chain to create such cavities.
A natural follow-up study would be investigating automatical choice of the number of control points. This could be based on a reversible jump MCMC strategy as in , but we do not discuss such possibilities further here.
In this sparse-data study involving 6 projections only, our proposed method is significantly better than the baseline TV method, as shown in Table 4 and Table 5. Clearly, even our new method still has a room for improvement as the relative errors in shape reconstruction are in the interval 3%-12%. However, such improvement must be based on using more a priori information since the measurement data is very sparse. In this study we do not want to impose more or stricter prior information in the fear of over-fitting our model. In a practical application one might have better and more accurate a priori information about the target, which could then be used for lowering the reconstruction error further.
There are several interesting avenues for further investigation. Namely, MCMC-NURBS tomography can be used for imaging corrosion inside pipelines; in that case the inner boundary of the pipeline is represented as a NURBS curve. Also, one can extend the new method to targets comprising several components. In that case the MCMC algorithm can have random “birth” and “death” events varying the number of components. See  for such a study involving cracks in electrical conductivity.
Our NURBS-based nonlinear Bayesian inversion performs very well in the tomographic task of recovering homogeneous 2D objects from extremely sparse projection data. By using advanced MCMC techniques, we have shown that the MCMC approach is computationally feasible. To tackle the heavy computation, further speed-up strategies are available, such as parallellization, improving the choice of initial value and optimizing the covariances of the sampling strategy.
The results are conveniently in CAD-compatible vector-graphic format. Quantitative comparison to the baseline method, optimally thresholded TV regularization, is favorable to our method. In the case of recovering the attenuation value, NURBS-MCMC delivers results with relative errors one order of magnitude smaller that the baseline method.
This work was supported by the Academy of Finland through the Finnish Centre of Excellence in Inverse Problems Research 2012–2017, decision number 250215, and Väisälä Foundation. We warmly thank Stephan Anzengruber for contributing the integrated code for NURBS.
- Alpers, Andreas, et al. "Geometric reconstruction methods for electron tomography." Ultramicroscopy 128 (2013): 42-54.
- Andersen, Kim E., Stephen P. Brooks, and Martin B. Hansen. "A Bayesian approach to crack detection in electrically conducting media." Inverse Problems 17.1 (2001): 121.
- Andrieu, Christophe, and Johannes Thoms. A tutorial on adaptive MCMC. Statistics and Computing 18.4 (2008): 343-373.
- Anson, Oscar, Francisco J. Seron, and Diego Gutierrez. NURBS-based inverse reflector design. Proceedings of CEIG 2008 (2008): 65-74.
- Bazilevs, Y., and T. J. R. Hughes. NURBS-based isogeometric analysis for the computation of flows about rotating components. Computational Mechanics 43.1 (2008): 143-150.
- Brooks, Stephen P., and Gareth O. Roberts. Convergence assessment techniques for Markov chain Monte Carlo. Statistics and Computing 8.4 (1998): 319-335.
- Chib, Siddhartha, and Edward Greenberg. Understanding the metropolis-hastings algorithm. The american statistician 49.4 (1995): 327-335.
- Cottrell, J. A., Hughes, T.J.R. and Bazilevs, Y.: Isogeometric analysis : toward integration of CAD and FEA. John Wiley & Sons, Ltd 2009.
- Gelman, Andrew, G. Roberts, and W. Gilks. "Efficient metropolis jumping rules." Bayesian statistics 5.599-608 (1996): 42.
- Green, Peter J. "Reversible jump Markov chain Monte Carlo computation and Bayesian model determination." Biometrika 82.4 (1995): 711-732.
- Haario, Heikki, Laine, M., Mira, A., & Saksman, E. DRAM: efficient adaptive MCMC. Statistics and Computing 16.4 (2006): 339-354.
- Jing, Zhang, Feng Shaowei, and Cui Hanguo. Optimized NURBS curve and surface fitting using simulated annealing. Computational Intelligence and Design, 2009. ISCID’09. Second International Symposium on. Vol. 2. IEEE, 2009.
- Kaipio, J., Somersalo, E.: Statistical and Computational Inverse Problems. Springer (2004)
- Khan, Sanowar H., Kirill Y. Aristovich, and Alexey I. Borovkov.: Solution of the forward problem in magnetic-field tomography (MFT) based on magnetoencephalography (MEG). Magnetics, IEEE Transactions on 45.3 (2009): 1416-1419.
- Klann, Esther, Ronny Ramlau, and Wolfgang Ring. A Mumford-Shah level-set approach for the inversion and segmentation of SPECT/CT data. Inverse Probl. Imaging 5.1 (2011): 137-166.
- Kolehmainen, V : Novel approaches to image reconstruction in diffusion tomography. Ph.D thesis.
- Kolehmainen, V., Siltanen, S., Järvenpää, S., Kaipio, J.P., Koistinen, P., Lassas, M., Pirttilä, J. and Somersalo, E.. Statistical inversion for medical x-ray tomography with few radiographs: II. Application to dental radiology. Physics in Medicine and Biology 48.10 (2003): 1465.
- Ma, W., and J-P. Kruth. NURBS curve and surface fitting for reverse engineering. The International Journal of Advanced Manufacturing Technology 14.12 (1998): 918-927.
- Mohammad-Djafari, Ali, and Ken Sauer.: Shape reconstruction in X-ray tomography from a small number of projections using deformable models. Maximum Entropy and Bayesian Methods. Springer Netherlands, 1998. 183-198.
- Mohammad-Djafari, Ali.: A Bayesian approach to shape reconstruction of a compact object from a few number of projections. arXiv preprint physics/0111120 (2001).
- Mira, Antonietta. "On Metropolis-Hastings algorithms with delayed rejection." Metron 59.3-4 (2001): 231-241.
- Mueller J. and Siltanen S.: Linear and Nonlinear Inverse Problems with Practical Applications. SIAM, Computational Science and Engineering (2012).
- Piegl, Les, and Wayne Tiller. The NURBS book. 1997. Monographs in Visual Communication (1997).
- Purisha, Zenith, and Siltanen, Samuli . "Tomographic Reconstruction of Homogeneous 2D Geometric Models with Unknown Attenuation." System Modeling and Optimization. Springer Berlin Heidelberg, 2014. 247-256.
- Purisha, Zenith, and Samuli Siltanen. "Tomographic Inversion using NURBS and MCMC." Forging Connections between Computational Mathematics and Computational Geometry. Springer International Publishing, 2016. 153-166.
- Renken, F. P., and G. Subbarayan.: NURBS-based solutions to inverse boundary problems in droplet shape prediction. Computer methods in applied mechanics and engineering 190.11 (2000): 1391-1406.
- Robert, Christian, and George Casella. Monte Carlo statistical methods. Springer Science & Business Media, 2013.
- Roberts, Gareth O., and Jeffrey S. Rosenthal. Examples of adaptive MCMC. Journal of Computational and Graphical Statistics 18.2 (2009): 349-367.
- Rogers, D. F.: An Introduction to NURBS : with historical perspective, vol. 1 of Academic Press, Morgan Kaufmann (2001)
- Rosenthal, Jeffrey S. Optimal proposal distributions and adaptive MCMC. Handbook of Markov Chain Monte Carlo (2011): 93-112.
- Saini, Deepika, Sanjeev Kumar, and Tilak Raj Gulati. Reconstruction of free-form space curves using NURBS-snakes and a quadratic programming approach. Computer Aided Geometric Design (2015).
- Sarfraz, M. Computer-aided reverse engineering using simulated evolution on NURBS. Virtual and Physical Prototyping 1.4 (2006): 243-257.
- Siltanen, S., Kolehmainen, V., Järvenpää, S., Kaipio, J.P., Koistinen, P., Lassas, M., Pirttilä, J. and Somersalo, E. Statistical inversion for medical x-ray tomography with few radiographs: I. General theory. Physics in medicine and biology 48.10 (2003): 1437.
- Tierney, Luke. Markov chains for exploring posterior distributions. the Annals of Statistics (1994): 1701-1728.
- Tierney, Luke. "A note on Metropolis-Hastings kernels for general state spaces." Annals of Applied Probability (1998): 1-9.