Radial Basis Function (RBF) methods are important tools for scattered data interpolation and for the solution of PDEs in complexly shaped domains. The most straight forward approach that is used to evaluate the method while incorporating the “standard basis functions” involves solving a linear system which is typically poorly conditioned. Two variations of a method, dubbed the RBF-QR approach, use a different basis that spans the same space as the standard basis but which in some cases results in a better conditioned linear system. Software packages which implement the two RBF-QR approaches are freely available ( and ).
Extended precision floating point arithmetic can be used to accurately evaluate the ill-conditioned problem in the standard basis. This approach has been used and implemented in several different software environments that include: Mathematica , the Matlab Symbolic Toolbox , C++ , and Fortran . The extended precision approach is attractive because it retains one of the great strengths of the RBF method – simplicity. Whereas the RBF-QR approaches, as can be ascertained by browsing the software that implements the methods, is far more complex. A software package that implements the RBF method in extended precision is not currently available. A detailed comparison of the RBF-QR and extended precision with the standard basis approaches is made in .
In this work a Matlab  toolbox is described that features a regularization method for the ill-conditioned linear system, extended precision floating point arithmetic, and symmetry exploitation for the purpose of reducing flop counts. The toolbox is called the Matlab Radial Basis Function Toolbox (MRBFT). The toolbox uses an object oriented approach to organize its functionality via three main classes. Static methods in the class rbfx are used to implement functionality associated with RBF methods in general, while class methods are used to implement methods in subclasses of rbfx which apply to a particular RBF. The superclass rbfx has abstract methods which every subclass must implement. These include a definition of the RBF itself as well as a variety of derivative operators applied to the RBF. The complete list of abstract methods is as follows:
methods(Abstract = true) v = rbf(obj, r, s); % RBF definition d = D1(obj, r,s, x); % first derivative wrt x d = D2(obj, r, s, x); % second derivative wrt x d = D3(obj, r, s, x); % third derivative wrt x d = D4(obj, r, s, x); % fourth derivative wrt x d = G(obj, r, s, x, y); % Gradient d = L(obj, r, s); % Laplacian d = B(obj, r, s, x, y); % Biharmonic operator d = D12(obj, r, s, x, y); % mixed partial derivative d = D22(obj, r, s, x, y); % mixed partial derivative end
The rbfx class is subclassed by the gax and iqx classes which implement methods that are respectively particular to the Inverse Quadratic and Gaussian RBF. The class rbfCenters implements methods associated with scattered center locations and the class rbfCentro implements methods associated with algorithms that accurately and efficiently operate on structured matrices.
RBF interpolation uses a set of N distinct pointsin ℛd called centers. No restrictions are placed on the shape of problem domains or on the location of the centers. A RBF
is an infinitely differentiable (compactly supported and global RBFs without a shape parameter and with less smoothness exist but are not considered here) function of one variablethat is centered at and that contains a free parameter ɛ called the shape parameter. The RBF interpolant assumes the form
where a is a vector of expansion coefficients. The Gaussian (GA) RBF
and the inverse quadratic (IQ) RBF
are representative members of the class of global, infinitely differently RBFs containing a shape parameter that interpolate with exponential accuracy. The two RBFs and their various derivative are defined respectively in the classes gax and iqx. A particular RBF, for example the GA, is instantiated via
>> phi = iqx(); % phi is an instance of the inverse quadratic RBF class
Many other RBFs exist and may be added by the user to the toolbox by extending rbfx and using gax and iqx as examples.
The entire list of functions associated with an object is available as follows:
Methods for class iqx:
B D2 D4 iqx
D1 D22 G rbf
D12 D3 L
distanceMatrix1d distanceMatrix3d solve
distanceMatrix2d dm variableShape
First the methods that the class must define are listed and then the static methods inherited from the superclass are listed.
The RBF expansion coefficients are determined by enforcing the interpolation conditions
which result in a N × N linear system
The matrix B with entries
is called the system matrix. The solve method which is a static method of the rbfx class can be used to find the expansion coefficients. Static methods can be called in two ways. Either through the class
>> a = iqx.solve(B,f); % or a = rbfx.solve(B,f);
or from instances of the class
>> a = phi.solve(B,f);
Matlab is optimized for operations involving matrices and vectors. The process of revising loop-based, scalar-oriented code to use Matlab matrix and vector operations is called vectorization. At every opportunity, functions have been vectorized so that they execute as efficiently as possible.
The functions that define RBFs and their derivative matrices take distance matrices as their arguments. Taking for example 2d wherethe signed distance matrices
respectively contain the signed distance between the x and y coordinates of centers j and k. The distance matrix
stores the distance between centers j and k. With the x and y coordinates of the centers located in arrays xc and yc the distance matrices are formed via
>> [r, rx, ry] = rbfx.distanceMatrix2d(xc,yc);
and the RBF system matrix is constructed as
>> B = phi.rbf(r,s);
The evaluation of the interpolant (2) at M points xj can be accomplished by multiplying the expansion coefficients by the M × N evaluation matrix H that has entries entries
Continuing the 2d example, consider evaluation points x = (x1, x2). The distance matrices containing the distances between the centers and evaluation points
respectively contain the signed distances between the x and y coordinates of evaluation point j and center k. Then
With the x and y coordinates located in arrays x and y the distance matrices are formed via
>> [re, rxe, rye] = rbfx.distanceMatrix2d(xc,yc,x,y);
and the evaluation matrix is constructed via
>> H = phi.rbf(re,s);
The interpolant is then evaluated as
>> fa = H*a;
Derivatives are approximated by differentiating the RBF interpolant as
where 𝒟 is a linear differential operator. The operator 𝒟 may be a single differential operator or a linear differential operator such as the Laplacian. Evaluating (13) at the centers X can be accomplished by multiplying the expansion coefficients by the evaluation matrix H𝒟 with entries
That is, 𝒟f ≈ H𝒟a. For example, to approximate the first derivative with respect to x of a function of two variables using the MRBFT
>> Hd = phi.D1(r,s,rx); >> fa = Hd*a;
Alternatively, derivatives can be approximated by mutiplying the grid function valuesby the differentiation matrix D = H𝒟B–1 since
This is accomplished as
>> D = phi.dm(B,Hd); >> fa = D*f;
The shape parameter ɛk may take on different values at each center 8, 10, 22] that indicates that the use of a variable shape parameter may improve the conditioning of the system matrix as well as improve accuracy. A drawback of variable shape parameter strategies is that they cause the RBF system matrix to be non-symmetric. Reference  suggested the exponentially varying shape parameter(or equivalently in each column of the system or evaluation matrix). Such an approach is called a variable shape parameter strategy. Numerical evidence exists [
and well as the linearly varying shape parameter
Reference  gives examples of the benefits of using a random variable shape parameter
All three variable shape strategies are available in the MRBFT via
>> Bs, Hs = phi.variableShape(sMin,sMax,opt,N,M);
Both equations (6) for the expansion coefficients and (15) for the differentiation matrix assume that the system matrix is invertible. Both the GA and IQ system matrices are symmetric positive definite (SPD) if a constant shape parameter is used and therefore they are invertible. While it is invertible, the system matrix is typically very poorly conditioned and it may cease to be numerically SPD and the standard algorithm to factorize a SPD matrix may fail. The eigenvalues of B satisfy 0 < 𝜆min = 𝜆1 ≤ 𝜆2 ≤ … ≤ 𝜆N = 𝜆max and the matrix condition number in the 2-norm is κ(B) = 𝜆min/𝜆max. For a fixed set of centers, the shape parameter affects both the accuracy of the method and the conditioning of the system matrix. The RBF method is most accurate for smaller values of the shape parameter where the system matrix is ill-conditioned. Figure 1 illustrates a typical result. For fixed N, the error is reduced as the shape parameter is made smaller until the condition number of the system matrix reaches 𝒪(1018) in double precision and 𝒪(1036) in quadruple precision. A regularization technique to mitigate the effects of the ill-conditioning is discussed in the next section.
Reference  demonstrated that a simple regularization technique can mitigate the effects of the poor conditioning of RBF system matrices and in most cases ensure that the RBF system matrix remains numerically SPD so that Cholesky factorization can be used. Instead of solving the system
the regularized system
where C = B + µI is solved. The parameter µ is a small positive constant called the regularization parameter and I is the identity matrix. The technique is called the method of diagonal increments (MDI) and its first use dates back to the 1940’s . Matrix C is better conditioned than B as
For small µ, (B + µI)–1 is close to B–1 and MDI simply replaces B with (B + µI) in computing the solution of a system. Equation
quantifies how close that (B + µI)–1 and B–1 are . For very small µ the difference is negligible. A good choice of the parameter is µ = 5ɛM where ɛM is machine epsilon in the floating point number system being used.
All MRBFT functions that involve the factorization of a SPD matrix take two optional arguments: a regularization parameter µ and a logical variable safe. The default value of µ is 5e–15 which is the suggested value for double precision. The standard linear equation solver in Matlab is the mldivide function which may be evoked via the backslash operator. If a matrix is symmetric, Cholesky factorization is attempted. If Cholesky factorization fails, then the matrix is factorized with LU factorization. The default value of safe is true which causes all MRBFT routines to use the mldivide function. Setting safe to false forces the routines to use Cholesky factorization. The danger in directly calling Cholesky factorization is that if the regularization parameter is not large enough a matrix may fail to be numerically SPD and Choleksy factorization will fail due to taking the square root of a negative number if the matrix is severely ill-conditioned.
Figure 2 shows the results of a 2d interpolation problem with scattered centers with and without MDI regularization in double precision. Both the accuracy and condition number of the two approaches are virtually the same when the condition number of the system matrix can be accurately calculated in double precision which is for κ(B) ≈ 𝒪(1016) and smaller. When the condition number of the unregularized method reaches beyond that threshold the regularization keeps the condition number in the approximately 𝒪(1016) range and the regularized solution is approximately two decimal places more accurate in this example.
The MRBFT uses the Multiprecision Computing Toolbox for Matlab (MCT)  for its extended precision functionality. The MCT enables extended precision data types to be seamlessly used in place of the standard double type. As a result, existing Matlab programs can be converted to run with arbitrary precision with minimal changes. IEEE 754-2008 compliant quadruple precision is supported and the MCT is highly optimized for this case. Note that the MRBFT is in no way dependent on the MCT. The installation of the MCT is not necessary. However, without the MCT, the MRBFT is limited to double precision.
The following code that calculates the RBF expansion coefficients gives an example of how to change double precision computations to extended precision:
1 phi = iqx(); 2 mp.Digits(34); % digits of decimal precision 3 N = mp('30'); % replace with N = 30 to convert to double precision 4 xc = linspace(0,1,N); 5 r = phi.distanceMatrix1d(xc); 6 B = phi.rbf(r,2.5); 7 f = sin(xc); 8 a = phi.solve(B,f,0);
The only coding difference in the above code between double and extended precision is on lines 2 and 3. Once the number of digits is specified and N is changed to a mp object, all other operations involving the object are then done in extended precision. The script interpBenchExtended.m compares the execution speed of a 2d interpolation problem over a range of the shape parameter in both double and quadruple precision. In this particular example the quadruple precision calculation takes approximately 52 times longer to execute which is a typical result. In benchmarks, the MCT has been shown to be much more efficient in implementing extended precision than other available options. Additional comparisons of execution times can be found in reference .
The second major class of the MRBFT is rbfCenters. The methods of the class are the following:
Halton2d circleCenters squareCenters Hammersley2d circleUniformCenters
RBF methods place no restriction on the location of centers. However, randomly locating centers is unlikely to lead to good results. Theory dictates [2, 23] that centers should well-cover a domain in the sense that the centers are somewhat uniformly distributed with no large holes in their coverage and no centers clumped extremely close together. Centers located too close together hurt conditioning while large holes in the center coverage negatively affect convergence rates. Computational experience indicates that it is beneficial to locate centers more densely in boundary regions than in the interior of domains.
Quasirandom sequences  have become popular choices of centers for RBF methods. A quasirandom sequence is a n-tuples that fills n-space more uniformly than uncorrelated random points. Quasirandom sequences are also called low-discrepancy sequences. Halton points  are probably the most used quasirandom sequence in RBF methods due to the sequences being featured in the book . However, it is our experience that the Hammersley points  provide a superior coverage in many cases. The MRBFT implements functions that produce both Halton and Hammersley center distribution for circular and square domains. The functions
>> [x, y] = rbfCenters.squareCenters(N,a,b,clust er,ch,plt); % square [a,b] x [a,b] >> [x, y] = rbfCenters.circleCenters (N,cluster, ch,R,plt); % circle of radius R
provide the option to cluster the centers more densely around the boundaries of the domains. Figure 3 shows a set of Hammersley points and a set of clustered Hammersley points on the unit circle. An example of how the MRBFT can be used to distribute boundary clustered quasirandom centers in a more complexly shaped domain can be found in the script complexCentroCenters.m in the examples directory of the MRBFT.
The third major class of the MRBFT is rbfCentro which takes advantage of symmetry to reduce the computational expense and storage requirements of RBF methods. The methods of the class are:
centroCenters centroDecomposeMatrix hasSymmetry centroCircle centroEig isCentro centroCondition- centroMult isSkewCentro Number centroDM fullCentroMatrix solveCentro
A matrix B is centrosymmetric if B = JBJ and is skew-centrosymmetric if B = –JBJ where J, the contra-identity matrix, is a square matrix whose elements are all zero except those on its southwest-northeast diagonal which are all 1’s. RBF system and differentiation matrices (details in ) are (skew) centrosymmetric if the signed distance matrices (11) are (skew) centrosymmetric. Centrosymmetry does not depend on the location of centers, but rather on the distance between centers. Any center distribution in one dimension that is symmetrical about its mid-point causes the signed distance matrix to be skew-centrosymmetric. Center distributions in two dimensions that cause the signed distance matrix to be skew-centrosymmetric are easily generated in domains that are symmetric with respect to either the x-axis, y-axis, or the origin, or that can be made so by a linear transformation or rotation. Figure 4 shows two center distributions that result in centrosymmetric distance matrices. Centrosymmetry allows for significant flop count reductions in numerical linear algebra routines associated with RBF methods as well as reductions in computer memory requirements . The efficiency gains when working with centrosymmetric matrices come from the fact that a N × N (where N is even and M = N/2) centrosymmetric matrix has the block structure
where B11, B21, and J are M × M. The block structure allows many linear algebra operations to be performed while only forming and operating on half of the matrix. Additionally, the matrix (23) and
are orthogonally similar which allows for efficient eigenvalue calculation. Similar decompositions are possible for odd N as well.
The MRBFT functions for centrosymmetry include the following. Centrosymmetric center distributions in circular domains are produced by
>> [xc,yc] = rbfCentro.centroCircle (N,cluster,ch,R,plt);
and in domains with either x, y, or origin symmetry the function
>> [xc,yc] = rbfCentro. centroCenters(x,y,symType,plt);
uses half of the centers given as an argument to return a centrosymmetric center distribution.
Only half of the matrix needs to be formed and stored and then subsequently operated on in centrosymmetric linear algebra routines. The half-sized distance matrices are formed via calling the distanceMatrix functions with the arguments modified as follows:
>> r = rbfx.distanceMatrix1d(xc(1:N/2),xc); >> [r, rx, ry] = rbfx.distanceMatrix2d(xc(1:N/2) ,yc(1:N/2),xc,yc);
The half sized distance matrices can be used to form half-sized system and derivative evaluation matrices as
>> B = phi.rbf(r,s); >> H = phi.D2(r,s,rx);
Matrix condition numbers in the 2-norm are calculated via
>> [kappaB, kappaL, kappaM] = rbfCentro. centroConditionNumber(B,mu);
with a factor of four savings in the dominant term in the flop count compared to the standard algorithm. The RBF expansion coefficients for interpolation problems can be found from the half-sized system matrix as
>> a = rbfCentro.solveCentro(B,f,mu,safe);
Two quarter-sized matrices are factorized by a Cholesky factorization and the dominant term in the flop count is reduced by a factor of four fromto . In a similar manner, differentiation matrices are formed via
>> D = rbfCentro.centroDM(B,F,N,rho,mu,safe);
with a factor of four savings in computational effort. Centrosymmetric matrix-vector multiplication is asymptotically faster by a factor of two over the standard algorithm and can be accomplished via
>> [L,M] = rbfCentro.centroDecomposeMatrix(D,rho); % decompose into smaller matrices >> fp = rbfCentro.centroMult(f,L,M,rho); >> % multiply and reconstruct solution
Other miscellaneous functions concerning centrosymmetry are
>> c = rbfCentro.hasSymmetry(B); %tests a N x N matrix for (skew) centrosymmetry
which test matrices for symmetry and
>> D = rbfCentro.fullCentroMatrix(Dh,N,skew);
which constructs a full-size centrosymmetric matrix from a half-sized matrix.
The MRBFT distribution includes three directories that contains scripts that provide examples of using the functions in the toolbox. The examples directory contains scripts that use the functions for typical tasks associated with the RBF method. The tests directory contains scripts that verify that the various functions are working as claimed. The benchmarks directory contains scripts that measure the execution time of the centrosymmetric algorithms versus the standard algorithms and that measure the execution time of algorithms in double precision versus extended precision. Below a brief description of each script is given. More detailed information is contained in the comments of each script. Extensive comments are also contained in the source code the the classes rbfx, rbfCenters, and rbfCentro.
The following example scripts come with the toolbox:
The problem is solved two ways: 1) standard algorithms, 2) centrosymmetric algorithms. With N = 5000 the accuracy of the two approaches is the same but the centrosymmetric approach is approximately five times faster and requires only half the storage compared to the standard approach.
The following test scripts come with the toolbox:
The following benchmark scripts come with the toolbox:
Additional results demonstrating the benefits of using extending precision with RBF methods and the efficiency gains from exploiting centrosymmetry can be found in references  and  for which the MRBFT was used extensively to produce the examples and results.
The MRBFT is a freely available collection of Matlab files and scripts that implement RBF methods for scattered data interpolation and for the numerical solution of PDEs. The toolbox has been developed over a period of several years while it has been used in the author’s research. The class rbfx implements routines for common tasks associated with all RBFs and defines abstract methods for differential operators that must be implemented by all subclasses that define particular RBFs. The class rbfCenters implements methods for locating uniform and quasi-random centers in rectangular and circular domains. The class rbfCentro implements methods for locating centers in symmetric domains that cause RBF matrices to have a structure that allows for substantial saving in storage requirements and reductions in flop counts for associated linear algebra routines. All MRBFT functions can be implemented using extended precision floating point arithmetic provided that the MCT  is installed. Scripts that provide examples, tests, and benchmarks of the MRBFT are included with the distribution. Comments in the class files provide additional documentation.
The author uses the MRBFT in his own research and as a result the toolbox is constantly being improved and new features are being added. Future improvements to the MRBFT include: developing a class of methods for working with the local RBF method, implementing Mercer’s method for the GA RBF , developing a class that implements rational RBF methods, and adding methods to the rbfCenters class that distribute centers in more complexly shaped domains. Updates, bug fixes, and other improvements to the MRBFT are posted at  where the project is hosted.
A summary of the implementation of the MRBFT is as follows. The toolbox is implemented in Matlab which is widely used in Mathematics and other scientific disciplines. An object oriented approach is used to organize the functionality of the toolbox.
The MRBFT has been developed over a period of several years as it has been used in the author’s research. It was used extensively in the preparation of references  and  which are accompanied by scripts that use the MRBFT to produce the results in the manuscripts. Additionally, as discussed in the overview section, the toolbox comes with a collection of scripts that demonstrate its usage, benchmark its performance, and verify that its algorithms produce the correct results.
The MRBFT is programmed in Matlab which is available on Windows, OS X, and Linux.
The MRBFT was developed and tested using Matlab versions 2015b through 2016b . It is not clear as to which version of Matlab first featured object oriented programming (OOP). However, version R2008a made major changes in the way Matlab implements OOP. Thus, while not tested, it is possible that the MRBFT is compatible with Matlab versions as old as 2008a.
The MRBFT uses the Multiprecision Computing Toolbox for Matlab (MCT)  for its extended precision functionality. Note that the MRBFT is in no way dependent on the MCT. The installation of the MCT is not necessary. However, without the MCT, the MRBFT is limited to double precision floating point arithmetic.
The author is the only contributor to the software. However, the contributions of others are welcome.
Code repository GitHub
Licence: GNU GPL V3
Date published: 28/03/16
The potential for others to use the MRBFT is substantial. RBF methods have become very popular in applications. The MRBFT may potentially find use in industrial applications. Research in RBF methods is very active. The MRBF is a potentially useful tool for faculty research, graduate student research, and undergraduate research projects. The object oriented design of the toolbox makes it extremely easy for it to be extended by subclassing rbfx to define new RBFs. The MRBFT can be used to foster the idea of reproducible research in the area of RBF methods. The MRBFT provides core functionality that other research can be built on and for which the code is freely available.
The author has no competing interests to declare.
Advanpix (2016). Multiprecision computing toolbox for Matlab, version 3.9.9 for 64-bit Linux. http://www.advanpix.com/
Buhmann, M D (2003). Radial Basis Functions. Cambridge University Press, 0521633389DOI: https://doi.org/10.1017/CBO9780511543241
Fasshauer, G and McCourt, M (2012). Stable evaluation of Gaussian RBF interpolants. SIAM Journal on Scientific Computing 34: 737–762, DOI: https://doi.org/10.1137/110824784
Fasshauer, G E (2007). Meshfree Approximation Methods with Matlab. World Scientific, DOI: https://doi.org/10.1142/6437
Fornberg, B, Larsson, E and Flyer, N (2011). Stable computations with gaussian radial basis functions in 2-d. SIAM Journal on Scientific Computing 33: 869–892, DOI: https://doi.org/10.1137/09076756X
Henderson, H V and Searle, S R (1981). On deriving the inverse of a sum of matrices. SIAM Review 23(2): 53–60, DOI: https://doi.org/10.1137/1023004
Huang, C-S, Leebn, C-F and Cheng, A-D (2007). Error estimate, optimal shape factor, and high precision computation of multiquadric collocation method. Engineering Analysis with Boundary Elements 31: 614–623, DOI: https://doi.org/10.1016/j.enganabound.2006.11.011
Kansa, E J (1990). Multiquadrics – a scattered data approximation scheme with applications to computational fluid dynamics I: Surface approximations and partial derivative estimates. Computers and Mathematics with Applications 19(8/9): 127–145, DOI: https://doi.org/10.1016/0898-1221(90)90270-T
Kansa, E J (1990). Multiquadrics – a scattered data approximation scheme with applications to computational fluid dynamics II: Solutions to parabolic, hyperbolic, and elliptic partial differential equations. Computers and Mathematics with Applications 19(8/9): 147–161, DOI: https://doi.org/10.1016/0898-1221(90)90271-K
Kansa, E J and Carlson, R (1992). Improved accuracy of multiquadric interpolation using variable shape parameters. Computers and Mathematics with Applications 24: 99–120, DOI: https://doi.org/10.1016/0898-1221(92)90174-G
Larsson, E (2012). A MATLAB implementation of the RBF-QR method. http://www.it.uu.se/research/scicomp/software/rbf_qr
McCourt, M (2014). GaussQR: Stable Gaussian computation. http://math.iit.edu/~mccomic/gaussqr/
Niederreiter, H (1992). Random Number Generation and Quasi-Monte Carlo Methods, CBMS-NSF, SIAM. PhiladelphiaDOI: https://doi.org/10.1137/1.9781611970081
Piegorsch, W and Casella, G (1989). The early use of matrix diagonal increments in statistical problems. SIAM Review 31: 428–434, DOI: https://doi.org/10.1137/1031089
Sarra, S A (2011). Radial basis function approximation methods with extended precision floating point arithmetic. Engineering Analysis with Boundary Elements 35: 68–76, DOI: https://doi.org/10.1016/j.enganabound.2010.05.011
Sarra, S A (2014). Regularized symmetric positive definite matrix factorizations for linear systems arising from RBF interpolation and differentiation. Engineering Analysis with Boundary Elements 44: 76–86, DOI: https://doi.org/10.1016/j.enganabound.2014.04.019
Sarra, S A (2016). A Matlab radial basis function toolbox with symmetry exploitation, regularization, and extended precision. http://www.scottsarra.org/rbf/rbf.html
Sarra, S A and Cogar, S (2017). An examination of evaluation algorithms for the RBF method. Engineering Analysis with Boundary Elements 17: 36–45, DOI: https://doi.org/10.1016/j.enganabound.2016.11.006
Sarra, S A and Sturgill, D (2009). A random variable shape parameter strategy for radial basis function approximation methods. Engineering Analysis with Boundary Elements 33: 1239–1245, DOI: https://doi.org/10.1016/j.enganabound.2009.07.003
Wendland, H (2005). Scattered Data Approximation. Cambridge University Press, DOI: https://doi.org/10.1017/CBO9780511617539
Wong, T, Luk, W and Heng, P (1997). Sampling with Hammersley and Halton points. Journal of Graphics Tools 2: 9–24, DOI: https://doi.org/10.1080/10867651.1997.10487471