Shape recovery from sparse tomographic Xray data
Abstract
A twodimensional 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 computeraided design (CAD) software. However, the linear tomography task becomes a nonlinear inverse problem due to the NURBSbased 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 Xray 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.
Contents
 1 Introduction
 2 Materials and Methods
 3 Results
 4 Discussion
 5 Conclusions
1 Introduction
We propose a new reconstruction algorithm for twodimensional Xray tomography of objects with homogeneous attenuation coefficient. The basic idea is to represent the boundary of the object in terms of a nonuniform rational basis spline (NURBS) curve, thereby providing a direct connection to computeraided 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 backprojection, 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 illposed inverse problem. See Figure 1.
The proposed MCMCNURBS reconstruction method is first tested with two simulated phantoms, one of them convex and the other nonconvex. See [24, 25] for other examples. The Xray attenuation coefficient is constant inside the phantoms, and we simulate six projection images using fanbeam 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 nonnegative 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 Xray 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 nonconvex 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 Xray tomography that readily provide segmented images can be found in [19, 20, 15]. Nondestructive testing from very restricted data using Bayesian inversion in terms of a vectorgraphic 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 [8]. One of the applications of NURBSbased 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 [5].
The method proposed here is different from all the previous ones listed above.
2 Materials and Methods
2.1 Nurbs
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:
(1) 
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
(2) 
where by definition. A collection of breakpoints is then called knot vector:
(3) 
where .
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 rddegree 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
(4) 
and .
For computational reasons, we need to construct two discrete models: a pixelbased object model and a NURBSbased object model as in Figure 4. Pixelbased 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 pencilbeam model. NURBSbased 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 pixelbased model using an operator where NN is the resolution of the pixel image. Define
(5) 
Let consider as in Figure 5.
In the pixelbased 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 fanbeam projections as shown in Figure 3. Our measurement model can be written as follows.
(6) 
2.3 Tomographic Projection Data
The xray tomography data of the sugar phantom were acquired with the custombuilt 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 xray detector is a 12bit 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 xray 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 crosssection 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 illposed inverse problem. To compensate this illposedness, 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 [13].
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:
(7) 
where is a matrix model of the measurement and m is the xray measurement.
Prior distribution
In this part, we construct a prior density which will contribute to remove the illposedness. 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
(8) 
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:
(9) 
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 selfintersection can be avoided if holds.
Likelihood function
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.
Posterior distribution
The solution of the inverse problem is the posterior probability distribution which has the following form:
or
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:
(10) 
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:
(11) 
This CM estimate is the result of the recovered control points and attenuation coefficient.
A well known algorithm in MCMC is Metropolis Hastings [7, 9, 21, 34]. The Markov chains are generated as follow:

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 nsamples 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 edgepreserving reconstruction. Let us denote . The solution of TVregularized is defined by finding the vector that minimizes the penalty functional:
(12) 
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 [22] 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
(13) 
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
(14) 
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 nonconvex 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 plasticprinted using 3D printer machine. Both are then filled by white crystal sugar as shown in Figure 9.
3 Results
Throughout the chapter, we use the following measurement setup: sparseangle (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 nonadaptation 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 thresholdedTV reconstructions are presented. See Figures 10 and 11. The optimal TV parameter choice is calculated as discussed in Subsection 2.7. In the thresholdedTV 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, .
Optimal  Optimal  

Simulated  82.76  0.0249 
Simulated  82.76  0.0225 
Optimal  Optimal  

Physical  31.03  0.0254 
Physical  93.10  0.0215 
Thresholded Total Variation  MCMCNURBS  

Simulated  7.8%  0.37% 
Simulated  16.7%  0.74% 
Total Variation  MCMCNURBS  

Simulated  32.83%  7.48% 
Simulated  56.93%  3.9% 
Total Variation  MCMCNURBS  

Physical  50.95%  11.67% 
Physical  15.7%  5.95% 
The simulation time for each run is given in Table 6.
NURBSMCMC  10  12 
ThresholdedTV  0.67  0.67 
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 1d 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 burnin period. The TV regularizations for the 2D tomographic case using quadratic programming as well as the thresholdedTV reconstructions are presented in Figures 17 and 18.
4 Discussion
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 binaryimage 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 twodimensional chain plots and by Geweke’s convergence diagnostic [6].
Some of the histograms in Figures 12 and 15 are multimodal and clearly show the nonGaussian nature of the posterior distribution as a result of our nonlinear NURBSbased 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 vectorgraphic approach we use.
The reconstructions of the simulated and the physical phantoms using the MCMCNURBS 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 nonconvex example shapes we use.
The prior we ended up using is a tradeoff between several aspects: mathematical simplicity, not too many control points, and flexibility in representing nonconvex 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 upperleft 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 followup study would be investigating automatical choice of the number of control points. This could be based on a reversible jump MCMC strategy as in [10], but we do not discuss such possibilities further here.
In this sparsedata 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 overfitting 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, MCMCNURBS 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 [2] for such a study involving cracks in electrical conductivity.
5 Conclusions
Our NURBSbased 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 speedup 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 CADcompatible vectorgraphic format. Quantitative comparison to the baseline method, optimally thresholded TV regularization, is favorable to our method. In the case of recovering the attenuation value, NURBSMCMC delivers results with relative errors one order of magnitude smaller that the baseline method.
Acknowledgements
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.
References
 Alpers, Andreas, et al. "Geometric reconstruction methods for electron tomography." Ultramicroscopy 128 (2013): 4254.
 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): 343373.
 Anson, Oscar, Francisco J. Seron, and Diego Gutierrez. NURBSbased inverse reflector design. Proceedings of CEIG 2008 (2008): 6574.
 Bazilevs, Y., and T. J. R. Hughes. NURBSbased isogeometric analysis for the computation of flows about rotating components. Computational Mechanics 43.1 (2008): 143150.
 Brooks, Stephen P., and Gareth O. Roberts. Convergence assessment techniques for Markov chain Monte Carlo. Statistics and Computing 8.4 (1998): 319335.
 Chib, Siddhartha, and Edward Greenberg. Understanding the metropolishastings algorithm. The american statistician 49.4 (1995): 327335.
 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.599608 (1996): 42.
 Green, Peter J. "Reversible jump Markov chain Monte Carlo computation and Bayesian model determination." Biometrika 82.4 (1995): 711732.
 Haario, Heikki, Laine, M., Mira, A., & Saksman, E. DRAM: efficient adaptive MCMC. Statistics and Computing 16.4 (2006): 339354.
 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 magneticfield tomography (MFT) based on magnetoencephalography (MEG). Magnetics, IEEE Transactions on 45.3 (2009): 14161419.
 Klann, Esther, Ronny Ramlau, and Wolfgang Ring. A MumfordShah levelset approach for the inversion and segmentation of SPECT/CT data. Inverse Probl. Imaging 5.1 (2011): 137166.
 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 xray tomography with few radiographs: II. Application to dental radiology. Physics in Medicine and Biology 48.10 (2003): 1465.
 Ma, W., and JP. Kruth. NURBS curve and surface fitting for reverse engineering. The International Journal of Advanced Manufacturing Technology 14.12 (1998): 918927.
 MohammadDjafari, Ali, and Ken Sauer.: Shape reconstruction in Xray tomography from a small number of projections using deformable models. Maximum Entropy and Bayesian Methods. Springer Netherlands, 1998. 183198.
 MohammadDjafari, 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 MetropolisHastings algorithms with delayed rejection." Metron 59.34 (2001): 231241.
 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. 247256.
 Purisha, Zenith, and Samuli Siltanen. "Tomographic Inversion using NURBS and MCMC." Forging Connections between Computational Mathematics and Computational Geometry. Springer International Publishing, 2016. 153166.
 Renken, F. P., and G. Subbarayan.: NURBSbased solutions to inverse boundary problems in droplet shape prediction. Computer methods in applied mechanics and engineering 190.11 (2000): 13911406.
 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): 349367.
 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): 93112.
 Saini, Deepika, Sanjeev Kumar, and Tilak Raj Gulati. Reconstruction of freeform space curves using NURBSsnakes and a quadratic programming approach. Computer Aided Geometric Design (2015).
 Sarfraz, M. Computeraided reverse engineering using simulated evolution on NURBS. Virtual and Physical Prototyping 1.4 (2006): 243257.
 Siltanen, S., Kolehmainen, V., Järvenpää, S., Kaipio, J.P., Koistinen, P., Lassas, M., Pirttilä, J. and Somersalo, E. Statistical inversion for medical xray 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): 17011728.
 Tierney, Luke. "A note on MetropolisHastings kernels for general state spaces." Annals of Applied Probability (1998): 19.