Uniform sampling (by volume) within a cone

Question

I'm looking for an algorithm that can generate points within a cone with a flat bottom (a disk).

I have the normalized axis along which the cone is being created (for our purposes let's just say it is the y-axis so (0, 1, 0) and the angle of the cone (let's say it is 45 degrees).

The only resources I could find online generate vectors within a cone, but they are based on sampling a sphere, so at the bottom you get a kind of "snow-cone" effect instead of a disk at the bottom.

That is done with the following pseudocode:

// Sample phi uniformly on [0, 2PI]
float phi = rand(0, 1) * 2 * PI

// Sample u uniformly from [cos(angle), 1]
float u = rand(0, 1) * (1 - cos(angle * PI/180)) + cos(angle * PI/180)

vec3 = vec3(sqrt(1 - u^2) * cos(phi), u, sqrt(1 - u^2) * sin(phi)))

The below picture is what I am going for. Having the ability to generate samples either on the surface or inside would be nice as well:

Answer 1

I could explain my solution in detail using integrals and probability distributions, but the lack of MathJax on this site makes that difficult. I'll keep my explanation at a simple level, but it should be clear. I'll also make the solution a little more general than you ask: we want a random point inside a right circular cone of height a and radius of base b, and we want the point with uniform sampling over the volume of that cone. This method directly chooses a random point in the cone without any rejection testing.

First let's consider the small cone of height h inside that larger cone, both cones with the same apex and parallel bases. The two cones are of course similar figures, and the square-cube law says that the volume of the smaller cone varies as the cube of its height. That height varies from 0 to a and we want its cube to be uniform over that range. Therefore we choose h to vary with the cube root of a uniform random variable, and we get (in Python 3 code),

h = a * (random()) ** (1/3)

We next consider the circular region that is the base of that smaller cone of height h. The radius of that base is (b / a) * h, by similar triangles. Now think of a smaller circular region of radius r inside that larger circular region, both circles in the same plane and with the same center. The area of the smaller circle varies with the square of its radius, so to get a uniform area over its range we take the square root of a uniform random variable. We get

r = (b / a) * h * sqrt(random())

We now want the angle t (for theta) of a point on the circumference of that smaller circle of radius r. The angle in radians obviously does not depend on the other factors, so we just use a uniform random variable to get

t = 2 * pi * random()

We now use those three random variables h, r, and t to choose our point inside the starting cone. If the apex of the cone is at the origin and the axis of the cone is along the positive y-axis, so that the center of the base is (0, a, 0) and a point on the circumference of the base is (b, a, 0), you can choose

x = r * cos(t)
y = h
z = r * sin(t)

---

When you asked about generating samples "on the surface" you did not clarify if you mean just the side (or is it "sides"?) of the cone, just the base, or the entire surface. Your second graphic appears to mean just the side, but I'll give code for all three.

The side only

Again we use a smaller cone of height h inside the larger cone. Its surface area varies as the square of its height, so we take the square root of a uniform random variable. The circle in its base is fixed, if our point is to be on the surface, and again the angle is just uniform. So we get

h = a * sqrt(random())
r = (b / a) * h
t = 2 * pi * random()

Use the same code for x, y, and z I used above for the interior of the cone to get the final random point on the side surface of the cone.

The base only

This is much like choosing a point in the interior, except the height is predetermined to equal the height of the entire cone. We get the following, somewhat simplified code:

h = a
r = b * sqrt(random())
t = 2 * pi * random()

Again, use the previous code for the final x, y, and z.

The entire surface

Here we can first decide, at random, whether to place our point on the base or on the surface, then place the point in one of the two ways above. The area of the base of a cone of height a and base radius b is pi * b * b while the surface area of the cone's side is pi * b * sqrt(a*a + b*b). We use the ratio of the base to the total of those areas to choose which subsurface to use for our point:

if random() < b / (b + sqrt(a*a + b*b)):
    return point_on_base(a, b)
else:
    return point_on_side(a, b)

Use my codes above for the side and base to complete that code.

---

Here are simple matplotlib 3D scatter plots of 10,000 random points, first inside the cone then on its side surface. Note that I made the apex angle 45°, as your text states but unlike your pictures. Viewing these from other angles seems to confirm that they are uniform in volume or area.