Efficient barycentric point sampling on meshes
Jamie Portsmouth


The Study of the Bloch Transformed Fields Scattered by Locally Perturbed Periodic Surfaces
    Ruming ZhangCenter for Industrial Mathematics, University of Bremen
; 
    
==========================================================================================





We present an easy-to-implement and efficient analytical inversion algorithm for the unbiased random sampling of a set of points on a triangle mesh whose surface density is specified by barycentric interpolation of non-negative per-vertex weights. The correctness of the inversion algorithm is verified via statistical tests, and we show that it is faster on average than rejection sampling.  
 
Computer GraphicsG.3Probability and StatisticsProbabilistic algorithms








§ RELATED WORK



Point sampling on meshed surfaces is useful in a variety of computer graphics contexts <cit.>. 
Here we focus on the relatively simple problem of sampling a spatially inhomogeneous, unbiased random distribution of points on a mesh given a prescribed point density per unit area varying over the mesh, which has received surprisingly little attention. 

There has been much work on the generation of point distributions with desired local statistical characteristics for applications such as point stippling and blue noise sampling <cit.>, however these techniques are complex to implement, generally not very efficient, and deliberately introduce bias into the sampling.


Statistically unbiased inhomogeneous random sampling is clearly useful in a variety of contexts, for example in Monte Carlo sampling. For the homogeneous case, there is a standard inversion algorithm for sampling a random point on a triangle mesh with uniform probability <cit.>. 

For the inhomogeneous case,
rejection sampling <cit.> provides a general and relatively efficient solution, as was noted (for example) by Yan et. al <cit.>, however as a random algorithm this presents efficiency problems. Sik et. al <cit.> improved on rejection sampling by subdividing the mesh until each triangle can be considered to have uniform density without losing too much accuracy, then applying the homogeneous inversion algorithm. However this requires a relatively complex preprocessing stage.

Arvo et. al <cit.> extended the inversion algorithm to deal with a density which varies linearly within each face of a triangle mesh (corresponding to barycentric interpolation). However they expressed their solution in the form of a set of cubic equations, without simplifying further or developing a practical, tested algorithm. In this paper we revisit the approach of Arvo et al., and provide a more explicit algorithm than previously, with statistical validation and performance comparison to rejection sampling. 



§ METHOD


Consider the general problem of independently sampling random points on a three-dimensional triangle mesh, where the probability per unit area p_X(𝐱) of choosing a given point 𝐱 is proportional to a non-negative scalar weight on the mesh, ϕ(𝐱).
This weight can be interpreted as specifying the relative surface density (i.e. number per unit area) of sampled points. 



Normalization then implies:

    p_X(𝐱) = ϕ(𝐱)/∑_i ∫_T_iϕ(𝐲)  dA_i(𝐲)

where dA_i(𝐲) is the area element on the three-dimensional surface of triangle T_i.
This may be factored into the discrete probability of choosing a given triangle, multiplied by the conditional PDF (with area measure) of choosing a point within that triangle:

    p_X(𝐱)    =    p_T(T_i)   p_X|T(𝐱| T_i) 
       =   ∫_T_iϕ(𝐲)  dA_i(𝐲)/∑_k∫_T_kϕ(𝐲) dA_k(𝐲)·ϕ(𝐱)/∫_T_iϕ(𝐲)  dA_i(𝐲) .

Defining the area-weighted average of ϕ over triangle T_i (with area A_i) as ⟨ϕ⟩_i ≡∫_T_iϕ(𝐱)  dA_i(𝐱) / A_i, then the discrete probability of each triangle can be written as

    p_T(T_i)    =   A_i⟨ϕ⟩_i/∑_k A_k ⟨ϕ⟩_k  ,

and the conditional area-measure PDF of a point 𝐱 in a given triangle T_i is 

    p_X|T(𝐱| T_i) = ϕ(𝐱)/A_i ⟨ϕ⟩_i .

Sampling from the discrete probability distribution p_T(T_i) to choose a triangle is done via the CDF

    P_T(T_i) = ∑_j≤ i A_j⟨ϕ⟩_j/∑_k A_k ⟨ϕ⟩_k  ,

i.e. sample a uniform random deviate (a random variable drawn from the uniform distribution on the open interval [0,1)), then find index k such that P_T(T_k-1) ≤ξ < P_T(T_k) (with P_T(T_k)=0 by convention for k<k_min). This is usually done via bisection search, though note that more efficient algorithms exist such as that described by Sik et al. <cit.>.




It remains to sample from the conditional PDF p_X|T(𝐱| T_i).

[caption=Triangle selection code, label=lst:trianglesample, float]
void <@TRIANGLE_CDF@>( vector<double>    tri_cdf )

	double cdf = 0.0;
	for (size_t ti=0; ti<<@NUM_TRIS@>(); i++) 
	
		double phi_avg = <@WEIGHT_AVG@>(ti);
		cdf += phi_avg * <@AREA@>(ti);
		tri_cdf.push_back(cdf);
	
	for (size_t ti=0; ti<<@NUM_TRIS@>(); i++) tri_cdf[i] /= cdf;


static inline size_t <@CHOOSE_TRI@>(const vector<double>    tri_cdf)

	double xi = <@RAND@>();
	size_t jl = 0; size_t ju = tri_cdf.size()+1;
	while (ju-jl > 1)
	
		size_t jm = (ju+jl) >> 1;
		if (xi >= tri_cdf[jm-1]) jl=jm;
		else                     ju=jm;
	
	return jl;




A completely general method for sampling from p_X|T(𝐱| T_i)  given any function ϕ(𝐱) is provided by rejection sampling. 
We first find ϕ_max,i = max(ϕ(𝐱): 𝐱∈ T_i) (this can be precomputed). We then choose a random point in the triangle drawn from a uniform distribution. This is most easily done by parameterizing points 𝐱∈ T_i as 𝐱 = u 𝐯^i_u + v 𝐯^i_v + w 𝐯^i_w, where 𝐯^i_u, 𝐯^i_v, 𝐯^i_w are the triangle vertices and the per-triangle barycentric coordinates (u, v, w) are each in the range [0,1] with u+v+w=1. Then uniform random sampling on the triangle is done via the formulas <cit.>

    u = 1 - √(ξ_1) ,    v = (1 - u)  ξ_2

where (ξ_1, ξ_2) are uniform random deviates.
We then decide to accept this trial sample by drawing another uniform random deviate ξ_3 and testing whether ξ_3  ϕ_max,i < ϕ(𝐱). If this is true, we accept 𝐱 as the sample, otherwise we draw another trial sample of 𝐱 and continue, until acceptance.


This procedure is outlined in Algorithm <ref>.
*EndIf







While the rejection sampling method is very general, as a random algorithm it does not provide any strict guarantees about the number of random samples which will be taken. 



However, in the common case of a weight defined by per-vertex barycentric weighting, the inversion method provides an analytical formula for the sampled points<cit.>. Per-vertex weighting means we associate with each of the three vertices of triangle T_i a non-negative real weight. Let us denote the weights of the three vertices 𝐯^i_u, 𝐯^i_v, 𝐯^i_w as ϕ_u, ϕ_v, ϕ_w respectively. Assuming the triangle is not degenerate, we may express the barycentric coordinates in terms of position: u(𝐱), v(𝐱), w(𝐱). 
Barycentric interpolation then defines the weight at all points 𝐱∈ T_i via (equivalent to linear interpolation within the triangle):

    ϕ(𝐱) = u(𝐱) (ϕ_u-ϕ_w) + v(𝐱)  (ϕ_v-ϕ_w) + ϕ_w .

This is a common scheme for interpolating per-vertex weights to produce a C^0 continuous function on the mesh. 
Integrals over the triangle area elements dA_i may be completed by a change of variables to barycoordinates, via the standard identity dA_i = 2A_i  du  dv.  





It follows that the area-averaged barycentric weight is equal to the mean vertex weight:

    ⟨ϕ⟩_i = ϕ_u + ϕ_v + ϕ_w/3 .

Using Eqn. (<ref>), the normalized conditional PDF in barycentric coordinate measure is then given by

    p_U,V(u, v) = p_X|T(𝐱| T_i) dA_i/du  dv = 2 ϕ(u, v)/⟨ϕ⟩_i .


We now introduce the normalized relative weights

    Φ_u    ≡   ϕ_u-ϕ_w/⟨ϕ⟩_i ,   Φ_v ≡ϕ_v-ϕ_w/⟨ϕ⟩_i .

Each possible set of relative weights maps into a point in the (Φ_u, Φ_v) plane (i.e. each such point represents all triangles which have the same relative weights). Expressing the problem in terms of this two-dimensional vector of weights leads to some simplification in the analytical expressions compared to Arvo's approach <cit.>. Note that the case where ⟨ϕ⟩_i=0, i.e. all zero weights, can be ignored as the probability of sampling a point in such a zero-weighted triangle is zero. It follows from the non-negativity of the weights that (Φ_u, Φ_v) satisfy various inequalities





and thus lie within the triangular region in the (Φ_u, Φ_v) plane depicted in Figure <ref>. The uniform weighting corresponds to the origin (Φ_u, Φ_v)= (0,0).  

[style=bottom, caption=Routine to sample the u barycoordinate., label=lst:U_NEWTON, float]
double <@U@>( double Phi_u, double Phi_v, const double tol = 5.0e-3 )

    double r = <@RAND@>(); double l = (2.0*Phi_u - Phi_v)/3.0;
    const int maxIter = 20; double u = 0.5; int n=0;
    while (n++ < maxIter) 
        double u1 = 1.0-u;
        double P  = u*(2.0-u) - l*u*u1*u1 - r;
        double Pd = max(u1*(2.0 + l*(3.0*u-1.0)), DBL_EPSILON);
        double du = max(min(P/Pd, 0.25), -0.25); u -= du;
        u = max(min(u, 1.0-DBL_EPSILON), DBL_EPSILON);
        if (fabs(du) < tol) break;
    
    return u;


In these coordinates, the barycentric-measure PDF Eqn. (<ref>) reduces to

    p_U,V(u, v) = 2 ( u Φ_u + v Φ_v + 1 - Φ_u + Φ_v/3)  .


We now describe how to sample points from this PDF.
Integration gives the marginal PDF of u:

    p_U(u)    =   ∫_0^1-u p_U,V(u, v)  dv  .

The cumulative density function (CDF) for the marginal PDF of u, ∫_0^u p_U(u')  du' is given by (and similarly for the CDF of v, P_V(v))

    P_U(u)    =     u(2-u) - (2Φ_u - Φ_v)/3 u(u-1)^2	 .






Inversion of P_U(u)=ξ_u where ξ_u is a uniform random deviate yields the sampled value of u. It is efficient to solve this by Newton's method via the update rule: 

    u_n+1 = u_n - P_U(u_n)-ξ_u/P'_U(u_n) .

An example implementation of this sampling routine in C is provided in Listing <ref> (where  generates a uniform random deviate). This includes tolerances to keep the solution within bounds, and limits the step size to 1/2 to aid convergence (as statistically verified in Section <ref>).









In order to sample v, we require the CDF for v conditional on the sampled value of u, given by 

    P_V(v | u) = ∫_0^v p_V|U(v' | u)  dv' = 1/p_U(u)∫_0^v p_U,V(u, v')  dv'  .

Evaluating this gives

    P_V(v | u) = 2v/p_U(u)[ 1 + (u-1/3) Φ_u + 
      (v/2-1/3) Φ_v ]  .

Inversion of P_V(v | u)=ξ_v (where ξ_v is again a uniform random deviate) gives the sampled value of v. As P_V(v | u) is quadratic in v, this inversion has two solutions

    v_± = τ±√( τ^2 (1-ξ_v) + (τ + u - 1)^2ξ_v )
               
where the square root here denotes the principal square root, and 

    τ(u, Φ_u, Φ_v)    ≡   1/3 - 1+(u-1/3) Φ_u/Φ_v .
 
Here τ ranges over the entire real line, i.e. τ∈ [-∞, ∞].
As Φ_v approaches zero (uniform weighting), τ diverges; however, in this limit the inversion can be simplified to:

    v → (1-u)  ξ_v  + O(|τ|^-1)   as  |Φ_v| → 0, |τ| →∞ .
 

If τ≫ 0, then clearly the v_- branch must be chosen (in order that v ∈ [0,1]). Similarly, if τ≪ 0, then the v_+ branch must be chosen. Figure <ref> shows how the solution for v varies as a function of the random deviate ξ_v and τ. 
From this figure it is intuitively clear that for general τ, the correct choice of branch is given by

    v = {[              v_+  τ≤ (1-u)/2    ,;              v_- τ > (1-u)/2    . ].


[style=bottom, caption=Routine to sample the v barycoordinate., label=lst:V, float]
double <@V@>( double u, double Phi_u, double Phi_v )

	double r = <@RAND@>();
	const double epsilon = 1.0e-6 
	if (fabs(Phi_v) < epsilon) return (1.0 - u)*r;
	double tau = 1.0/3.0 - (1.0 + (u-1.0/3.0)*Phi_u)/Phi_v;
	double tmp = tau + u - 1.0;
	double q = sqrt(tau*tau*(1.0-r) + tmp*tmp*r);
	return tau <= 0.5*(1.0-u) ? tau + q : tau - q;



Let τ = γ(1 - u), which implies:

    v_±/(1-u) = γ±√(γ^2 + ξ_v(1-2γ)) .

First consider the case γ>1/2. Then the term inside the square root satisfies the following inequality

    γ^2 + ξ_v(1-2γ) > |γ-1|^2

since (1-2γ)<0 and ξ_v<1. Thus the following inequalities are satisfied:

    v_+/(1-u) > γ + |γ-1| > 1  ,   v_-/(1-u) < γ - |γ-1| < 1  .

Therefore since v≤(1-u), if γ>1/2 the v_- branch must be chosen. The analogous argument for γ<1/2 shows that the v_+ branch must be chosen in that case.



Algorithm <ref> summarises the resulting method for inversion sampling (u, v).
A C implementation for sampling v via this method is provided in routine V (Listing <ref>).

In Figure <ref> we show a simple example of a point distribution sampled via this inversion method, in which the per-vertex weight is taken to be the magnitude of the local discrete vertex curvature. In Figure <ref> we show another example in which the per-vertex weight function has a periodic 3d variation of the form |cos(x/L) cos(y/L) cos(z/L)| with some length-scale L. 











§ VALIDATION AND PERFORMANCE
 


In Figure <ref>, we plot the calculated CDF P_V(v) and empirically measured CDF obtained via the inversion sampling of Algorithm <ref>, for a representative set of triangle weightings (here those which fit in a 16x16 grid covering the valid region in the (Φ_u, Φ_v) plane in Figure <ref>). Only 0.1% of the CDF points are shown for clarity. 

The Kolmogorov-Smirnov statistic for each empirical CDF is at most D=0.005, which is less than the critical statistic at 99% confidence level (D_crit=1.63/√(N)=0.007),  giving confidence that the sampling algorithm is correct over the whole (Φ_u, Φ_v) region. 

We performed an alternative verification with the chi-squared test. 

We covered the valid region in the (Φ_u, Φ_v) plane with a 16x16 grid, and for each of these triangle weightings, drew samples and binned them in a grid in the uv plane (where at each bin resolution enough samples were drawn to ensure approximately 100 per uv bin). In Figure <ref>, we plot the logarithm of the chi-squared p-values for all the triangle weightings as a function of the sample binning resolution N in the uv plane (again using Newton's method).


We repeated this for increasingly high bin resolution. The exponential drop in p-value logarithm to <-50 for all the samples as the uv bin resolution is increased to 64x64 gives high confidence that the sampling algorithm is correct.



We focus here only on performance of the sampling within a given triangle, as we do not deal with optimization of the triangle selection itself. In Figure <ref> we show the relative performance of the rejection and inversion methods applied to independent samples from a single triangle (with Φ_u=-3, Φ_v=-3, i.e. the maximal per-triangle “inhomogeneity”) where each data point indicates the average sample time in nanoseconds averaged over 10^7 samples. These timings were taken running single-threaded on a 2.3 Intel Core i7 processor. The inversion method runs approximately 20-80% faster than rejection, where the efficiency depends strongly on the required tolerance for the sampled u barycoordinate, so some trade-off between accuracy and speed is involved. Of course, both inversion and rejection algorithms can also be easily multi-threaded. 



§ CONCLUSION


We derived a simple and efficient inversion method for sampling points on a triangle mesh with density defined by barycentrically interpolated per-vertex weights, and verified that it produces statistically correct point samples. We showed that weighted point sampling on a triangle mesh via the inversion method is faster on average than rejection sampling.

We note that the method presented here can overall be regarded as complementary to that of Sik et al. <cit.>. 
Their method involves usage of a more sophisticated method for triangle sampling, and is more general as it can deal with arbitrarily varying weight functions (as can rejection sampling), however it is also considerably more complex to implement. While an algorithm such as their fast triangle selection method is required for optimal performance, using our analytical inversion sampling of linearly varying weights should allow for further improvement as it would allow less subdivision to achieve the same sampling accuracy. We suggest that in future work it would be interesting to explore the combination of these methods.





eg-alpha-doi




