scikit-tda
Adding genus g surfaces
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!
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.
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
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
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?
@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.