Hello, I'm using marching square function on a planar slice of sdf values, but am getting weird zigzag lines in the output. Here is the link to the sdf value and a minimal script that reproduces the image below.

Thank you in advance for looking into it!

Our doublearea_instrinsic currently contains

dblA = 0.5 * np.sqrt((a+(b+c)) * (c-(a-b)) * (c+(a-b)) * (a+(b-c)))

If the areas are very small, this can NaN. Ideally, we would do

arg = (a+(b+c)) * (c-(a-b)) * (c+(a-b)) * (a+(b-c))
dblA = 0.5 * np.sqrt(np.maximum(arg, 0.))

However, if I make this change, the unit test for reach_for_the_spheres fails.
I speculate that it is because the optimization in reach_for_the_spheres does not terminate as early anymore, so it reaches a later stage where it can crash.

What should be done about it?

I ran the function
poisson_surface_reconstruction(nerf_pts, nerf_normals, gs=np.array([50,50,50])
and it results in the error
File "<__array_function__ internals>", line 200, in ravel_multi_index ValueError: invalid entry in coordinates array

A possibly fix for it is using "mode: clip" in the np.ravel_multi_index function, or exposing that option to the user.

This crashes, and it shouldn't. My current thinking is that the low resolution means there are labels that cause empty regions (e.g., labels with no associated tets). For example, I up-ed the resolution to 10, and that seems to run without issue.

import gpytoolbox, time
import numpy as np
 
labels = 4
resolution = 5
v, t = gpytoolbox.regular_cube_mesh(resolution)
w = np.zeros((v.shape[0],labels))
for i in range(labels):
    p = np.random.rand(3)
    print( "Center: ", p )
    for j in range(v.shape[0]):
        w[j,i] = -np.linalg.norm( v[j] - p )
 
print("Number of nodes: ", v.shape[0])
print("Number of labels: ", w.shape[1])
time_0 = time.time()
u, g, l = gpytoolbox.upper_envelope(v,t,w)
time_1 = time.time()
print("UE time: ", time_1-time_0)
print("Output tets: " , g.shape[0] )

The error occurs at the line

     46 try:
---> 47     from gpytoolbox_bindings_copyleft import _mesh_union_cpp_impl
     48     from gpytoolbox_bindings_copyleft import _mesh_intersection_cpp_impl

ImportError: DLL load failed while importing gpytoolbox_bindings_copyleft: The specified module could not be found.

During handling of the above exception, another exception occurred:

Any idea how to fix this?

Should catch this and return an empty line!

I tried to import some models from the Engineering Shape Benchmark (ESB)

For example backdoor.obj and backdoor1.obj from the Flat-Thin Wallcomponents.tar.gz
But the imported shape using the .obj type is (0, 0).

Importing the STL file works fine.

When the grinding wheel swept along the helix, I set the align_rotations_with_velocity to True and False respectively, and the result of True is a bit strange. Following is my code:

import polyscope

from test.context import gpytoolbox
from test.context import numpy as np
from test.context import unittest
from gpytoolbox import read_mesh
from gpytoolbox.copyleft import swept_volume
from gpytoolbox.copyleft import mesh_boolean
from gpytoolbox.copyleft import do_meshes_intersect

import polyscope as ps
ps.init()

# Read sample mesh
v, f = gpytoolbox.read_mesh("D:/pycharmpro/gpytoolbox/test/unit_tests_data/wheel.obj")
ps.register_surface_mesh("Wheel",v,f)

# Translation vectors to make Catmull-Rom spline
def generate_spiral_translations(radius, start_y, end_y, num_steps, pitch):
    """
    Generate translation vectors along a helical path around the y-axis,
    with specified start and end y-coordinates.

    Parameters:
    radius (float): Radius of the helix in the xz-plane.
    start_y (float): Starting y-coordinate of the helix.
    end_y (float): Ending y-coordinate of the helix.
    num_steps (int): Number of steps along the helix.
    pitch (float): Vertical distance (height) for one full rotation around the y-axis.

    Returns:
    translations (np.ndarray): Array of translation vectors (x, y, z).
    """
    # Calculate the height per step
    total_height = end_y - start_y
    h_step = total_height / num_steps

    # Initialize an empty array to store the translations
    translations = np.zeros((num_steps, 3))

    # Generate the translations
    for i in range(num_steps):
        # Calculate the angle for the current step
        # Use pitch to determine the angle per step, considering the full height of the helix
        theta = 2 * np.pi * i / (num_steps * pitch)

        # Calculate the x and z coordinates based on the angle and radius
        x = radius * np.cos(theta)
        z = radius * np.sin(theta)
        y = start_y + i * h_step  # Adjust y based on start_y and step height

        # Store the translation vector
        translations[i] = [x, y, z]

    return translations

# Example usage
radius = 0.025
start_y = -0.07
end_y = 0.07
num_steps = 50
pitch = 0.2  # Controls the number of rotations over the height

translations = generate_spiral_translations(radius, start_y, end_y, num_steps, pitch)

# Call swept volume function
u1,g1 = swept_volume(v,f,rotations=None,translations=translations,eps=0.001,verbose=True,align_rotations_with_velocity=True)
ps.register_surface_mesh("Swept_Volume1",u1,g1)

u2,g2 = swept_volume(v,f,rotations=None,translations=translations,eps=0.001,verbose=True,align_rotations_with_velocity=False)
ps.register_surface_mesh("Swept_Volume2",u2,g2)
ps.show()

The result:

align_rotations_with_velocity=True

align_rotations_with_velocity=False

It looks like one side is normal, and the other side is wrong. I don't know if it's due to my helix generation function. You can try to set a sweep with a path as a spiral and see if a similar problem occurs?

The write_mesh function can optionally take in texture coordinates UV and normals N, which require corresponding index matrices Ft and Fn, respectively. However, the function doesn't seem to check for cases where, for example, UV is provided but Ft is not provided. This results in situations where a written .obj file can contain vt rows but the f rows only contain vertex indices and no texture indices, so the file does not actually have a valid UV map.

I think a reasonable default in these cases would be to fall back on F, if there are as many texture coordinates/normals as vertices and Ft/Fn are not provided. If this isn't satisfied, I think the function should return an error. If it's preferable to avoid defaults and make the API explicit, I still think the function should return an error when UV/N is provided but Ft/Fn is not, so the semantics are clear.

I tried to remesh some models from the Engineering Shape Benchmark (ESB)

For example backdoor.stl from the Flat-Thin Wallcomponents.tar.gz

The result using remesh_botsch is exactly the same as the input.

import gpytoolbox as gpy
v, f = gpy.read_mesh(r"C:\Users\Michael\Desktop\Flat-Thin-Wallcomponents\Flat-Thin Wallcomponents\Back Doors\backdoor.stl")
u, g = gpy.remesh_botsch(v, f, h=None)

I discovered all vertices of the model are detected as boundary vertices by boundary_vertices(f)

But the model has clearly some "inner" vertices (at least on the rounded part in the middle.

Are you interested in me implementing the FEM for 3D Tets?

I already have a working version for every element type using gmsh.
You probably want to stay with numpy / scipy only, correct?
If you don't need it next week, I can adapt it (for linear elastic tet elements) to only use numpy/scipy 😀.

• Can you give me some hints where to find existing code/functionality for 3D Tets in gpytoolbox?

• Do you use voigt notation for stresses/straits?
  

• Or diagonal notation for stresses/straits?
  
  (only relevant for anisotropic materials)

Here's a png:

Running

from gpytoolbox import png2poly
poly = png2poly("octopus.png")[0]
import matplotlib.pyplot as plt
plt.plot(poly[:,0],poly[:,1])
plt.show()

for some reason produces a rotated polyline, as if the x and y coordinates were flipped:

Sometimes, you'll see that a polyline is given with the last point being equal to the first to denote that it's closed. Currently, signed_distance_polygon does not support this and instead returns nan since it finds an edge with zero length. It'd be nice if the function just ignored zero-length edges instead of giving nans.

This issue was discovered by Jack Zhang.

Our PLY reader fails to read in the face list for a triangle mesh. libigl reads the ply just fine. If I open the triangle mesh in blender and reexport as PLY, it works though.
Simple example (the mesh files are in meshes.zip):

import gpytoolbox as gpy
import igl

V_ply,F_ply = gpy.read_mesh("mesh.ply")
print(f"V_ply: {V_ply.shape}, F_ply: {F_ply.shape}")

V_ply_igl,F_ply_igl = igl.read_triangle_mesh("mesh.ply")
print(f"V_ply_igl: {V_ply_igl.shape}, F_ply_igl: {F_ply_igl.shape}")

V_ply_blender,F_ply_blender = gpy.read_mesh("mesh_exported_by_blender.ply")
print(f"V_ply_blender: {V_ply_blender.shape}, F_ply_blender: {F_ply_blender.shape}")

V_obj,F_obj = gpy.read_mesh("mesh.obj")
print(f"V_obj: {V_obj.shape}, F_obj: {F_obj.shape}")

This produces:

V_ply: (2802, 3), F_ply: (0, 0)
V_ply_igl: (2802, 3), F_ply_igl: (934, 3)
V_ply_blender: (2802, 3), F_ply_blender: (934, 3)
V_obj: (2802, 3), F_obj: (934, 3)

I.e., gpytoolbox returns an empty triangle list on the provided mesh.ply, but not on the reexported mesh_exported_by_blender.ply.

What could cause this?
Looking at the two PLY binary files in detail, one difference that I can see is that for mesh.ply, the vertex indices are given as int, and for mesh_exported_by_blender.py the vertex indices are given as uint:

header of mesh.ply:
element face 934.property list uchar int vertex_index

header of mesh_exported_by_blender.py:
element face 934.property list uchar uint vertex_indices

Indeed, our PLY binding just seems to assume vertex indices in the face list will always be uint. Libigl's does not, they have a switch statement supporting various types.

I tried to remesh a non-closed/ watertight mesh with the remesh_botsch algorithm.
Compared to a closed mesh (with almost the same size) it takes a lot longer (>100 times) and returns an error by exiting my python console.

• Is your implementation not useable for non-closed meshes? The original paper has options for the boundary.

The SSL certificate in gmplib.org, from which libigl (and therefore us) downloads gmp, has expired, which causes curl (called to download gmp during our build) to fail. I merged #41 , a hotfix that changes our version of libigl to my PR branch sgsellan/libigl, but I am opening this issue so we have a reminder to revert this in the future when the libigl folks have converged on a solution.

My understanding is that the swept volume computation does not support computing the swept area of a 2D mesh.
I would greatly appreciate if it if your implementation also supported 2D meshes. Thanks!

hello @sgsellan
I execute pip install gpytoolbox==0.1.0
I checked the previous version and the latest version seems to be missing something

The following problems occurred：
from gpytoolbox_bindings_copyleft import _swept_volume_impl
ImportError: DLL load failed while importing gpytoolbox_bindings_copyleft: 找不到指定的模块。

thank you!

Hey!

I'm trying to run reach for the sphere, but the code crashes on some models:

V,F = gpy.read_mesh("53159.obj")
j = 32
sdf = lambda x: gpy.signed_distance(x, V, F)[0]
gx, gy, gz = np.meshgrid(np.linspace(-1.0, 1.0, j+1), np.linspace(-1.0, 1.0, j+1), np.linspace(-1.0, 1.0, j+1))
U = np.vstack((gx.flatten(), gy.flatten(), gz.flatten())).T
V0, F0 = gpy.icosphere(2)
Vr,Fr = gpy.reach_for_the_spheres(U, sdf, V0, F0, min_h=.01)

Output:

 Segmentation fault (core dumped)

Is there some requirements on the input? The one I used (53159 of the Thingi10k dataset) is watertight and normalized to -0.8/0.8, but its genus is not zero.

libigl has removed mat_max in recent versions. We use it in src/cpp/upper_envelope.cpp.

We should update our code so it will work with newer versions of libigl when we eventually upgrade.

This dumb assert assumes 3D, we should remove it

converged = gpy.reach_for_the_spheres_iteration(state)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/silviasellan/miniconda3/lib/python3.11/site-packages/gpytoolbox/reach_for_the_spheres.py", line 667, in reach_for_the_spheres_iteration
if np.any((np.isnan(state.V))) or np.any(doublearea(state.V, state.F)==0) or len(non_manifold_edges(state.F))>0 :

EEEEEEEE

ERROR: test_bunny (main.TestRemeshBotsch)

Traceback (most recent call last):
File "D:\gpytoolbox\test\test_remesh_botsch.py", line 17, in test_bunny
v,f = gpytoolbox.read_mesh("test/unit_tests_data/bunny_oded.obj")
File "D:\gpytoolbox\src\gpytoolbox\read_mesh.py", line 89, in read_mesh
V,F,UV,Ft,N,Fn = _read_obj(file,return_UV,return_N,reader)
File "D:\gpytoolbox\src\gpytoolbox\read_mesh.py", line 135, in _read_obj
raise Exception(f"The file {file} could not be opened.")
Exception: The file test/unit_tests_data/bunny_oded.obj could not be opened.

Excuse me，Does Gpytoolbox have a visual function, for example, I want to display the swept body calculated by the swept_volume function, how can I do it?

Thank you for recommending Polyscope, it's easy to use. #128
But I want to demonstrate the process of spawning a swept body, as demonstrated in your swept_volumes' talk.
Probably because the swept_volume function in Gpytoolbox gets the result after the computation is completed, and Polyscope can't show the process of sweeping.
Do you have any suggestions?

The documentation for the torus function states that the first parameter (nR) is the number of vertices along the large perimeter and the second parameter (nr) is the number of vertices along the small perimeter, but in the code the meaning of these parameters are switched.

V,F = gpy.torus(10, 3, R=1., r=0.1) produces:

But according to the documentation should produce: