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)
 A: This paper may be of interest to you:
J. Arvo, Stratified Sampling of 2-Manifolds
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.
A: 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?
A: Triangulate it and sample each triangle with a density relative to the total area.
A: 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.
A: 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. 
