Principal/mean curvature on surface mesh

OVITO
1

Hello everyone,

I would like to calculate the mean curvature along the surface of my surface mesh generated by ovito. I wrote this function:

def compute_mean_curvature(vertices, faces):
    """
    Discrete mean curvature
    """
    n_vertices = len(vertices)
    L = np.zeros((n_vertices, n_vertices))  # Laplace matrix
    A = np.zeros(n_vertices)  # Voronoi area per vertex

    for face in faces:
        i, j, k = face
        vi, vj, vk = vertices[i], vertices[j], vertices[k]

        # edge vectors
        e0 = vj - vk
        e1 = vk - vi
        e2 = vi - vj

        # compute angles
        l0 = np.linalg.norm(e0)
        l1 = np.linalg.norm(e1)
        l2 = np.linalg.norm(e2)

        angle_i = np.arccos(np.clip(np.dot(-e1, e2) / (l1 * l2), -1.0, 1.0))
        angle_j = np.arccos(np.clip(np.dot(-e2, e0) / (l2 * l0), -1.0, 1.0))
        angle_k = np.pi - angle_i - angle_j

        cot_i = 1 / np.tan(angle_i)
        cot_j = 1 / np.tan(angle_j)
        cot_k = 1 / np.tan(angle_k)

        # fill Laplace matrix (cotangent weights)
        for a, b, cot in [(j, k, cot_i), (k, i, cot_j), (i, j, cot_k)]:
            L[a, b] += cot
            L[b, a] += cot
            L[a, a] -= cot
            L[b, b] -= cot

        # Compute area for mean curvature normalization
        face_area = 0.5 * np.linalg.norm(np.cross(vj - vi, vk - vi))
        A[i] += face_area / 3
        A[j] += face_area / 3
        A[k] += face_area / 3

    # Laplace-Beltrami operator approximation
    mean_curv = np.zeros(n_vertices)
    for i in range(n_vertices):
        delta = np.sum([L[i, j] * (vertices[j] - vertices[i]) for j in range(n_vertices)], axis=0)
        mean_curv[i] = 0.25 * np.linalg.norm(delta) / A[i]     

    return mean_curv

and it seems to work on a mesh that I generated using python (I made a spherical mesh and then did the numerical integration of the mean curvature along the surface which gave me approx 1/R as it should for a sphere), but once I try the same thing with the surface mesh that I generated over my spherical particle consisting of many monomers it no longer gives me 1/R. So here is my question: am I treating the mesh data wrongly and how do I have to interpret ovito’s output:

pipeline = import_file(tmp_filename)

pipeline.modifiers.append(ConstructSurfaceModifier(
  method=ConstructSurfaceModifier.Method.AlphaShape,
  radius=Probe_radius,
  smoothing_level=Smoothing,
  identify_regions=True
))

data = pipeline.compute()
mesh = data.surfaces['surface']
mesh_volume = data.attributes['ConstructSurfaceMesh.filled_volume']
mesh_surface_area = data.attributes['ConstructSurfaceMesh.surface_area']
vertex_positions_mesh = mesh.vertices["Position"]
faces_mesh = mesh.get_face_vertices()
vertices = np.array(vertex_positions_mesh, dtype=float)
faces = np.array(faces_mesh, dtype=int)

I appreciate any input, thank you in advance!
Greetings
Leah

2

Hi Leah,

I have skimmed your code. It looks okay how you access the SurfaceMesh data structure. I couldn’t understand and check in detail if the calculation of the local curvature is really correct though.

I have converted your script into a Python modifier in OVITO Pro and tested it on a roughly spherical SurfaceMesh (see screenshot below). The mean of the local curvature values seems to fit to the expected value. For comparison, I estimated the spherical curvature 1/R from the volume and surface area calculated by OVITO, e.g. 1/R = A/(3V).

image
image1730×1170 564 KB

Was the deviation in your own test more significant?

3

Hi Alexander,

Thank you so much for your fast and detailed reply!! I printed out np.mean(mean_curv) and I get: Average mean curvature H: 0.087145 (expected: 0.050000, since my radius of my spherical test particle is 20). Further I want to calculate:

    for face in faces:
        if len(face) != 3:
            print("Non-triangular face detected:", face)
        # Get vertex indices for the triangle
        v0, v1, v2 = vertices[face[0]], vertices[face[1]], vertices[face[2]]

        # Compute vectors
        a = v1 - v0
        b = v2 - v0

        # Cross product and area
        cross_product = np.cross(a, b)
        area = 0.5 * np.linalg.norm(cross_product)

        triangle_areas.append(area)

    print("calculated surface area:")
    print(np.sum(triangle_areas))
    Delta_A = np.sum((2*H_faces) * triangle_areas)
    R_A = np.sqrt( mesh_surface_area/(4*np.pi) )
    Delta_a = Delta_A/(8*np.pi*R_A)

and Delta_a should be almost 1 but I get 7.4932. This I think might have two issues: The deviating mean curvature and also how I calculate the triangle areas maybe, because I summed them up and expected them to be the same value as the surface mesh surface area but there is a big difference: surface mesh surface area:
4910.381599263607 and calculated surface area:
13769.115021524864 maybe I am using some areas not only once?

Thank you again for your help!
Best,
Leah

4

Oh, I think I have an idea of what is causing this problem: Once you activate the option ConstructSurfaceModifier.identify_regions, the modifier outputs a two-sided mesh. That means every triangular face has an opposite face with reserve winding order, which is connected to the same three vertices. This likely breaks your surface curvature calculation method.

A solution is to exclude all faces of the mesh from the calculation that belong to one side of the surface. I’ll get back to you here with more instructions soon.

5

Wow, would have never thought of these details. Thank you so much!! <3 Looking forward to hearing from you soon!

6

Another critical question: Is your model using periodic boundary conditions? If it is, then calculations like the following will likely be wrong if the three vertices of the face are on different sides of the periodic simulation cell:

# Compute vectors
a = v1 - v0
b = v2 - v0

You would have to modify the calculation to correctly take into account PBCs and wrap around the computed vectors. See SimulationCell.delta_vector() method for more information.

7

It is but I have a lot of objects through which I loop which is why I generate a new temporary file/data set for each of their coordinates and put it in the middle of the box to avoid these things.

8

I see, then PBCs may not be a concern here.

To get rid of one side of the mesh, you can add the following modifiers to your data pipeline:

pipeline.modifiers.append(ExpressionSelectionModifier(expression = 'Filled == 1', operate_on = 'surface_regions:surface/regions'))
pipeline.modifiers.append(DeleteSelectedModifier(operate_on = {'surface_regions'}))

They select and delete all mesh regions that fulfill the criterion Filled==1. With these regions also all faces that belong to these regions get deleted. What remains are only faces that are facing toward the unfilled spatial regions (Filled==0).

Can you try this, please? In my own test, if I don’t do the face deletion, the total area computed from the sum of the faces (your Python code) is exactly twice as large as the area reported by data.attributes['ConstructSurfaceMesh.surface_area']. When I include the deletion step, both area values match exactly.

9

So I did:

    pipeline = import_file(tmp_filename)

    pipeline.modifiers.append(ConstructSurfaceModifier(
        method=ConstructSurfaceModifier.Method.AlphaShape,
        radius=Probe_radius,
        smoothing_level=Smoothing,
        identify_regions=True
    ))
    pipeline.modifiers.append(ExpressionSelectionModifier(expression = 'Filled == 1', operate_on = 'surface_regions:surface/regions'))
    pipeline.modifiers.append(DeleteSelectedModifier(operate_on = {'surface_regions'}))

    data = pipeline.compute()
    mesh = data.surfaces['surface']
    mesh_volume = data.attributes['ConstructSurfaceMesh.filled_volume']
    mesh_surface_area = data.attributes['ConstructSurfaceMesh.surface_area']

and it gave me the same result. Do I have to write it differently into one pipeline modifier?

10

The code looks correct. When you do a print(mesh.faces.count), you should see the effect of the face deletion. The face count should drop to 1/2 of the original value.

I would suggest you send me the entire package (Python script and the test data file with the spherical particle). Then I can investigate in more detail what’s happening. In particular, I can perform the analysis in OVITO Pro, which provides more options to visualize the script output and check intermediate calculation results.

I don’t know if you already have the permission to attach files to MatSci forum posts. If not, please send them via email to [email protected].

11

Thank you! I sent you an email (regards: surface mesh - mean curvature)
See you there
Leah

12

The problem lies in your file temporary_debug.lammps: When I look at that file in OVITO, I see that some of the points are located outside the primary simulation box (see screenshot). They don’t form a contiguous object.

This does not represent a problem for the ConstructSurfaceModifier as it still builds the correct spherical surface mesh. However, the mesh vertices inherit their positions from the original input points, i.e., some of the vertices will have coordinates outside the simulation box. You just don’t see this immediately in the visualization, because OVITO wraps the mesh around when rendering it in the viewports and PBCs are turned on.

You can, however, make the problem visible by color coding the surface based on the z-coordinate values of the mesh vertices. That’s what I have done in the following screenshot, where the discontinuity in the mesh’s vertex coordinates can be seen.

image
image1836×1766 654 KB

So this is again where the PBC problem comes in, which I mentioned earlier. Your face area calculation routine is not prepared for this situation and yields incorrect results.

An easy way to fix this is by wrapping all particle coordinates back into the primary simulation box prior to constructing the surface mesh:

pipeline.modifiers.append(WrapPeriodicImagesModifier())
pipeline.modifiers.append(ConstructSurfaceModifier(...))

Then your script outputs the correct values.

By the way, it is not necessary to write the input points to a temporary file and then load that file again. You can directly hand over the point coordinates and build an ad-hoc pipeline in Python:

    # Transfer data to OVITO pipeline
    input_data = DataCollection()
    input_data.create_cell([
        [box_max-box_min, 0, 0, box_min],
        [0, box_max-box_min, 0, box_min],
        [0, 0, box_max-box_min, box_min]])
    particles = input_data.create_particles(count=len(positions_array))
    particles.create_property('Position', data=positions_array)
    pipeline = Pipeline(source=StaticSource(data=input_data))
13

ok I changed my function that calculates my observables to this: def calculate_surface_mesh_volume(positions, box_bounds):
“”"
Computes surface mesh volume and ΔA using curvature,
returns volume, surface area, ΔA.
“”"

positions_array = np.array(positions)
num_atoms = len(positions_array)

box_min = box_bounds[0][0]
box_max = box_bounds[0][1]
box_size = box_max - box_min

# Build OVITO data directly in memory
input_data = DataCollection()
input_data.create_cell([
    [box_size, 0, 0, box_min],
    [0, box_size, 0, box_min],
    [0, 0, box_size, box_min]
])
particles = input_data.create_particles(count=num_atoms)
particles.create_property('Position', data=positions_array)

pipeline = Pipeline(source=StaticSource(data=input_data))

# Wrap all points back into the box
pipeline.modifiers.append(WrapPeriodicImagesModifier())

pipeline.modifiers.append(ConstructSurfaceModifier(
    method=ConstructSurfaceModifier.Method.AlphaShape,
    radius=Probe_radius,
    smoothing_level=Smoothing,
    identify_regions=True
))

pipeline.modifiers.append(ExpressionSelectionModifier(
    expression='Filled == 1',
    operate_on='surface_regions:surface/regions'
))
pipeline.modifiers.append(DeleteSelectedModifier(operate_on={'surface_regions'}))

# Run pipeline
data = pipeline.compute()
mesh = data.surfaces['surface']

mesh_volume = data.attributes['ConstructSurfaceMesh.filled_volume']
mesh_surface_area = data.attributes['ConstructSurfaceMesh.surface_area']

# Geometry for further analysis
vertex_positions_mesh = mesh.vertices["Position"]
faces_mesh = mesh.get_face_vertices()
vertices = np.array(vertex_positions_mesh, dtype=float)
faces = np.array(faces_mesh, dtype=int)

# Compute curvature
H = compute_mean_curvature(vertices, faces)
H_faces = np.mean(H[faces], axis=1)

# Surface area via triangle areas
triangle_areas = []
for face in faces:
    v0, v1, v2 = vertices[face[0]], vertices[face[1]], vertices[face[2]]
    a = v1 - v0
    b = v2 - v0
    cross_product = np.cross(a, b)
    area = 0.5 * np.linalg.norm(cross_product)
    triangle_areas.append(area)

triangle_areas = np.array(triangle_areas)
Delta_A = np.sum((2 * H_faces) * triangle_areas)

# Diagnostics
print("surface mesh volume:", mesh_volume)
print("surface mesh surface area:", mesh_surface_area)
print("1/R from OVITO (S/(3V)):", mesh_surface_area / (3 * mesh_volume))
print("calculated triangle surface area:", np.sum(triangle_areas))

return mesh_volume, mesh_surface_area, Delta_A

and then I get
Average mean curvature H: 0.072944 (expected: 0.050000)
Average mean curvature H_faces: 0.058708 (expected: 0.050000)
surface mesh volume: 32080.61708786852
surface mesh surface area: 4910.381599263605
1/R from OVITO (S/(3V)): 0.05102127125354347
calculated triangle surface area: 4910.381599263608

and for the H_faces it seems to work, also for my real data files: when I output the surface area and the volume for each molecule and check it with the visualization on my laptop I get the same values.
I think it works now :smiley:
Thank you!
Leah

14

Excellent. Looks like we were able to get this problem solved. Let me know if you have any further questions.

The OVITO team would be delighted if you (or your supervisor/principal investigator) would consider licensing OVITO Pro to support our work on the OVITO Python package. But this is optional of course.

Finally, if you want to post Python code in this forum with correct formatting, you can find instructions here.

15

Thank you so much!! I will propose this to my supervisor. :slight_smile:

1 Like
16

Dear Alexander,

I am sorry to bother you again, but I would like to add a feature to my code where I only consider vertices and faces from the outside region of the surface mesh. Meaning my objects sometimes have holes in them which I would like to ignore. Is that possible with the filter you added to my code?

Best,
Leah

17

Please take a look at the documentation of the Construct Surface modifier, where the meanings of the properties “Filled” and “Exterior” are specified. What you are interested in is the part of the surface that bounds the exterior region. Thus, you should select and delete all mesh parts that belong to regions that are not exterior regions. This can be done by changing the selection criterion from Filled == 1 to Exterior == 0.

But keep this note from the documentation in mind:

Note that, for 3d periodic systems which do not have any open boundaries, the distinction between internal and external empty regions cannot be made. Empty regions will never be marked as exterior in 3d periodic models.

That means the (infinite) exterior region will only be identified if you turn off PBCs:

input_data.create_cell([
    [box_size, 0, 0, box_min],
    [0, box_size, 0, box_min],
    [0, 0, box_size, box_min]
], 
pbc=(False, False, False))
18

Ok! I get this error though: RuntimeError: Data pipeline failure: Modifier ‘Expression selection’ reported: Unexpected token “Exterior” found at position 0. I probably will not be able to solve this without upgrading ovito.

Thank you!

19

That’s odd. You only changed the expression string from Filled==1 to Exterior==0. Like this, right?

pipeline.modifiers.append(ConstructSurfaceModifier(
    method=ConstructSurfaceModifier.Method.AlphaShape,
    radius=Probe_radius,
    smoothing_level=Smoothing,
    identify_regions=True
))

pipeline.modifiers.append(ExpressionSelectionModifier(
    expression='Exterior==0',
    operate_on='surface_regions:surface/regions'
))
20

exactly, any imports I might be missing?