Adding genus g surfaces

[Sai271828]

Currently, there are only spherical and torus point cloud generators. It would be great to extend this functionality to include genus-g surfaces.
I’m not aware of a straightforward way to sample genus-g surfaces directly, but one possible approach could be to construct them by gluing together multiple tori. While this may not represent a smooth manifold in the strictest sense, I believe the persistent homology of such constructions would closely resemble that of a true genus-g surface.
If there are better alternatives for constructing or sampling genus-g surfaces, I’d be very interested to learn more!

[catanzaromj]

This seems like a great idea. Other approaches that might work:

• sampling from explicit parameterizations for low genus
• for genus >1, could we work in the universal cover/hyperbolic plane, sample there, and then project down via a covering map?
• attempt to generate a mesh of the surface and then sample from it?

It isn't clear to me yet how much of a problem it would be that your sampling method doesn't produce a smooth manifold in a strict sense. It might not be a problem at all.

[ctralie]

Here's another idea: use implicit functions. One approach could be to make the union of signed distance functions of individual tori, then use "walk on spheres" to do sampling. Have a look at Figure 2 in this paper, for example:
https://www.cs.cmu.edu/~kmcrane/Projects/MonteCarloGeometryProcessing/paper.pdf

[Sai271828]

Thank you for the suggestions!

This seems like a great idea. Other approaches that might work:

• sampling from explicit parameterizations for low genus
  As for explicit parametrizations, I am not aware of any that are both good and easy to compute, even for genus 2. If you have any references on this topic, I would be very interested.
• for genus >1, could we work in the universal cover/hyperbolic plane, sample there, and then project down via a covering map?
  This is quite appealing to me as well. I have primarily used the Poincaré disk in the past. However, I am curious if you have any suggestions for visualizing this quotient. For example, I would like to use persistent homology to confirm that it is indeed a genus g surface, and ideally, I would also like to plot the point cloud.
• attempt to generate a mesh of the surface and then sample from it?
  I don't quite follow this suggestion. Do you have any references for this?

It isn't clear to me yet how much of a problem it would be that your sampling method doesn't produce a smooth manifold in a strict sense. It might not be a problem at all.

Here is the current implementation I have, it is not perfect but it does produce decent (at least visually speaking) genus g point clouds:

def genusg_surface(
n: int = 100,
genus: int = 2,
c: float = 2.0,
a: float = 1.0,
noise: Optional[float] = None,
ambient: Optional[int] = None,
seed: Optional[int] = None,
uniform: bool = False,
) -> np.ndarray:
"""
Sample n data points from a genus-g surface. The genus is assumed to be at least 2.

Parameters
----------
n : int, default=100
    Number of data points in shape.
genus : int, default=2
    Genus of the surface.
c : float, default=2.0
    Distance from center to center of tube.
a : float, default=1.0
    Radius of tube.
noise: float, optional
    Standard deviation of normally distributed noise added to data.
ambient : int, optional
    Embed the surface into a space with ambient dimension equal to `ambient`. The surface is randomly rotated in this high dimensional space.
seed : int, optional
    Seed for random state.
uniform : bool, default=False
    If True, sample points uniformly on the surface.

Returns
-------
data : np.ndarray
    An ``(n,3)`` np.ndarray.
"""
rng = np.random.default_rng(seed)
# Sample points on the surface

assert genus >= 2, "Genus should be at least 2"
assert a <= c, "The major radius should be larger than the minor radius"
# Generate the left most donut
# For uniform sampling, we can use rejection sampling as done for genus 1

n_samples_per_genus = n // genus

left_donut = []
while len(left_donut) < n_samples_per_genus + n - n_samples_per_genus * genus:
    u = np.random.uniform(0, 2 * np.pi)
    v = np.random.uniform(0, 2 * np.pi)
    x = (c + a * np.cos(v)) * np.cos(u)
    y = (c + a * np.cos(v)) * np.sin(u)
    z = a * np.sin(v)
    if x < c + 3 * a / 4:
        if uniform:
            # REJECTION SAMPLING TO ENSURE UNIFORMITY
            # The map from u,v to x,y,z is not area-preserving, so we use rejection sampling to ensure uniformity
            jacobian = a * (c + a * np.cos(v))
            acceptance_prob = jacobian / (a * (c + a))
            if np.random.uniform(0, 1) < acceptance_prob:
                left_donut.append((x, y, z))

        else:
            left_donut.append((x, y, z))
left_donut = np.array(left_donut)

# Generate the right most donut
right_donut = []
while len(right_donut) < n_samples_per_genus:
    u = np.random.uniform(0, 2 * np.pi)
    v = np.random.uniform(0, 2 * np.pi)
    x = (c + a * np.cos(v)) * np.cos(u)
    y = (c + a * np.cos(v)) * np.sin(u)
    z = a * np.sin(v)
    if x > a / 4 - c - a:
        if uniform:
            # REJECTION SAMPLING TO ENSURE UNIFORMITY
            # The map from u,v to x,y,z is not area-preserving, so we use rejection sampling to ensure uniformity
            jacobian = a * (c + a * np.cos(v))
            acceptance_prob = jacobian / (a * (c + a))
            if np.random.uniform(0, 1) < acceptance_prob:
                right_donut.append((x, y, z))

        else:
            right_donut.append((x, y, z))
right_donut = np.array(right_donut)

# Generate middle donuts
middle_donuts = []
while len(middle_donuts) < n_samples_per_genus:
    u = np.random.uniform(0, 2 * np.pi)
    v = np.random.uniform(0, 2 * np.pi)
    x = (c + a * np.cos(v)) * np.cos(u)
    y = (c + a * np.cos(v)) * np.sin(u)
    z = a * np.sin(v)
    if (x > a / 4 - c - a) & (x < c + 3 * a / 4):
        if uniform:
            # REJECTION SAMPLING TO ENSURE UNIFORMITY
            # The map from u,v to x,y,z is not area-preserving, so we use rejection sampling to ensure uniformity
            jacobian = a * (c + a * np.cos(v))
            acceptance_prob = jacobian / (a * (c + a))
            if np.random.uniform(0, 1) < acceptance_prob:
                middle_donuts.append((x, y, z))

        else:
            middle_donuts.append((x, y, z))
middle_donuts = np.array(middle_donuts)

# Combine all points
# If genus is 2, we have only left and right donuts
if genus == 2:
    # MOVE RIGHT DONUT TO THE RIGHT
    right_donut[:, 0] += 2 * c + 3 * a / 2
    data = np.vstack((left_donut, right_donut))
else:
    # Move middle donuts to the right positions
    move_distances = [(2 * c + 3 * a / 2) * (i + 1) for i in range(genus - 1)]
    for i in range(0, genus - 2):
        temp_donut = middle_donuts.copy()
        temp_donut[:, 0] += move_distances[i]
        left_donut = np.vstack((left_donut, temp_donut))
    # Move right donut to the rightmost position
    right_donut[:, 0] += move_distances[-1]
    data = np.vstack((left_donut, right_donut))

if noise:
    data += noise * rng.standard_normal(data.shape)

if ambient:
    data = embed(data, ambient)

return data

[Sai271828]

Here's another idea: use implicit functions. One approach could be to make the union of signed distance functions of individual tori, then use "walk on spheres" to do sampling. Have a look at Figure 2 in this paper, for example: https://www.cs.cmu.edu/~kmcrane/Projects/MonteCarloGeometryProcessing/paper.pdf

Very interesting! However, I am not sure how the paper is related. Could you please expand a bit?

[ctralie]

@Sai271828 Sure, so I was coming at this by thinking about what it would take to glue torii together and sample from them. It's a bit messy to do that directly from mesh representations of the torii, so an easier way could be to switch to an implicit representation, like a signed distance function, where the union of two surfaces is much easier to perform.

At that point, we just need to sample from the surface. This walk on spheres gives a good method to sample points a levelset of an implicit function. We just have to uniformly sample initial points inside of the volume, which can be done with rejection sampling (by, e.g. drawing a point in the bounding box of the shape and only keeping it if it's inside). Then, we follow walk on spheres starting at that point (which is a bit like ray marching)

All that said, I haven't looked carefully at your technique yet, and it could be totally fine! I just thought this was an interesting problem.