## How can I sample uniformly from a surface?

Given an equation of a parametric surface, is there a general way to sample of points uniformly distributed on that surface?

I'm interested in this problem for purposes of visualisation - rather than attempting to attempt to triangulate the surface and display with polygons, display a dense sample of points. This makes it easier to generalise to >3d.

Here's an example of a surface I'd like to display: the Klein bottle.

u = [-pi, pi]
v = [-pi, pi]

x1 = (r * cos(v) + a) * cos(u),
x2 = (r * cos(v) + a) * sin(u),
x3 = r * sin(v) * cos(u/2),
x4 = r * sin(v) * sin(u/2)


(where r and a are parameters that control the shape of the overall surface)

-

This paper may be of interest to you:

It directly answers your question, though you may need to do some of the computations numerically depending on how complicated your surfaces are. Moreover, stratifying your samples will help the sampling look uniform -- sampling uniformly over the entire domain tends to look blotchy.

-
 Perfect - exactly what I was looking for. – Hadley Dec 29 2009 at 8:31

just to develop what Fedja said above:

say you have a map $\phi(x,y)=(u(x,y),v(x,y),w(x,y))$ from $D = [a;b]^2$ to $\mathbb{R}^3$ that represents a surface (Klein bottle in your initial example). The infinitesimal surface $(x;x+dx) \times (y;y+dy)$ is mapped to a surface of area $|\partial_x \Phi \wedge \partial_y \Phi| dx dy = A(x,y) dx dy$. So it would suffice to sample points inside $D$ distributed according to a probability distribution $p(x,y) \propto A(x,y) dx dy$.

1: as Fedja mentioned it, you can use a rejection sampling approach - you just need to find a bound on $\|A\|_{\infty}$, which should be very easy since $D$ is simple.

2: you can use a MCMC approach, which is equally easily implemented (e.g. independence sampler) - might be a little bit faster.

-

The problem boils down to the simulation of the point in a 2D domain with the density proportional to some easily computable function $f(x)$ (the area distortion coefficient of the mapping from the parameter domain to the surface). The simplest way is to simulate a random point $x$ in some square containing the domain, look at whether the point is in the domain, and, if it is, keep it with the probability $f(x)/M$ where $M$ is some number greater than all values of $f$. If the parameterization is reasonable enough, this should work pretty well. What examples do you have in mind?

-
 For the fun, I just tried this approach with a spherical cap. There is no problem with the simulation speed or the distribution of points but the output looks quite terrible to the eye: some strange collection of dots that undergoes incomprehensible transformations when rotating the picture. Certainly not something you would want to "display"... – fedja Dec 29 2009 at 4:40 I'm not sure I understand your approach. What is the area distortion coefficient? – Hadley Dec 29 2009 at 8:15 If $\psi$ is the parameterization map, it is just the square root of the determinant of $(D\psi)^TD\psi$ where $D$ is the differential of $\psi$. I call it the area distortion coefficient because the image of a small neighborhood of a point $x$ in the parameterization domain has the area approximately equal to this coefficient times the area of the neighborhood. The advantage here is that you do not need to convert your parameterization into an area-preserving map (which asks for integration, etc), just to make sure that the area distortion coefficient doesn't vary too much. – fedja Dec 29 2009 at 14:05 Ah, I see. Thanks for the explanation. – Hadley Dec 29 2009 at 15:46

It's not quite the same as generating a single random point, but: if you generate (a Poisson sample of) uniformly random lines and then take all intersection points of the lines with the surface, you'll get (a Poisson sample of) uniformly random points on the surface.

-
But computing intersections of lines with a parametrized surface is "fun"... – Mariano Suárez-Alvarez Dec 29 2009 at 1:22

Triangulate it and sample each triangle with a density relative to the total area.

-
 Note that you want to triangulate the region of R^2 you are parameterizing by, but then map each triangle in by the parametric equations. So your surface will still be curved, and will not look like a bunch of triangles. – David Speyer Dec 29 2009 at 1:57 This is what I meant. Thanks for the details. – lhf Dec 29 2009 at 10:58