Facet Subdivision

[timtitan]

Hello, I am looking at numpy-stl to use in to import stl files for a python based geometrical optics model and I need a way to set a maximum facet area. I had assumed a median subdivision algorithm was the best way to go but do you have any advice on the most efficient way to implement it? I'm still working out the structure for adding and removing facets.

[wolph]

If I understand correctly, you are trying to calculate the convex hull, right?
The qhull library is pretty good at that, the scipy module has an interface to calculate it: http://scipy.github.io/devdocs/generated/scipy.spatial.ConvexHull.html

Just note that it is a very heavy operation and could take a long time.

Alternatively, you can get a pretty nice approximation using an octree or kd-tree

[timtitan]

Not the Convex Hull no, I will be passing the facet vertices to a ray
casting algorithm to work out which facets are 'visible' from a set angle.
Due to demands of accuracy I need to ensure that no facet has an area
larger than a tenth of a wavelength or similar. Happily if I calculate the
midpoint between each vertices then those give the central point and 6
equal area smaller triangles. I'm still digging into your excellent code
and I thought getting your advice on how its meant to be manipulated would
be most helpful.
On 7 Jan 2016 6:13 p.m., "Rick van Hattem" notifications@github.com wrote:

If I understand correctly, you are trying to calculate the convex hull,
right?
The qhull library is pretty good at that, the scipy module has an
interface to calculate it:
http://scipy.github.io/devdocs/generated/scipy.spatial.ConvexHull.html

Just note that it is a very heavy operation and could take a long time.

Alternatively, you can get a pretty nice approximation using an octree or
kd-tree

—
Reply to this email directly or view it on GitHub
#20 (comment).

[wolph]

Ah, yes... I see your point. It's more like a ray tracer in that case... I would recommend the octree/kd-tree approach in that case, I've written something like that before in C++ https://github.com/WoLpH/computer_graphics

It's not too difficult to write but it takes quite a bit of time.

(oops, wrong button there)

[timtitan]

I think i'm going about this the wrong way. I can access the facet points through either mesh.mesh.vectors or mesh.mesh.points, but when I have calculated the six new triangles to replace the single facet which is to large, how do I remove the existing facet, then add new ones on? I assumed that pop() would work but it isn't an included method. Is the most robust way to alter the facets of an object to take the mesh.Mesh.vectors, and alter that array and add the new vertices to that, then pass it back in as data['vectors'] as shown in your modifying mesh objects example?

[wolph]

The basic data store within the mesh objects is the Mesh.data Numpy array. That limits modification somewhat however as numpy arrays have a fixed size, the only way to add/remove rows is by creating a new array. So if you wish to remove rows I would recommend creating a new array with the correct amount of rows and write the data to that. After that you can create a new Mesh object.

[timtitan]

Ok thank you, I had come to that to conclusion, will it calculate the normals and areas for the new mesh object automatically?

[wolph]

When you create a mesh object it automatically takes care of the normals: https://github.com/WoLpH/numpy-stl/blob/master/stl/base.py#L154-L155

The areas isn't updated automatically but can be updated through the update_areas function: https://github.com/WoLpH/numpy-stl/blob/master/stl/base.py#L204

[timtitan]

I've written a simple function to split any stl facets larger than a specified area into six smaller triangles with identical areas, but it seems to generating some odd results, i'm not sure if the normals are updating correctly either. I've included a rendering of the original test object and the new one with smaller facets. Both rectangles should be identical in size, just defined by different sized facets.

This is the rough code

#Facet Update Module
from stl import mesh
import numpy as np

def perp(a):
    b = np.empty_like(a)
    b[0] = -a[1]
    b[1] = a[0]
    return b

def line_intersect(a1,a2, b1,b2) :
    da = a2-a1
    db = b2-b1
    dp = a1-b1
    dap = perp(da)
    denom = np.dot( dap, db)
    num = np.dot( dap, dp )
    return (num / denom.astype(float))*db + b1

def mid(a,b): 
    midpoint=np.zeros([2])
    midpoint=[(a[0]+b[0])/2.,(a[1]+b[1])/2.,(a[2]+b[2])/2.] 
    return midpoint

def triangle_split(vertices):
    six_triangle_vertices=np.zeros([6,3,3])
    centroid=line_intersect(vertices[0],mid(vertices[1],vertices[2]),vertices[1],mid(vertices[0],vertices[2]))
    six_triangle_vertices[0]=[vertices[0],mid(vertices[0],vertices[1]),centroid]
    six_triangle_vertices[1]=[mid(vertices[0],vertices[1]),vertices[1],centroid]
    six_triangle_vertices[2]=[vertices[1],mid(vertices[1],vertices[2]),centroid]
    six_triangle_vertices[3]=[mid(vertices[1],vertices[2]),vertices[2],centroid]
    six_triangle_vertices[4]=[vertices[2],mid(vertices[0],vertices[2]),centroid]
    six_triangle_vertices[5]=[mid(vertices[0],vertices[2]),vertices[0],centroid]
    return six_triangle_vertices

def area_fix(object_list,area_limit):
    areas=[]
    large_facets=[]
    fixed_list=[None]*len(object_list)
    #large_facets becomes an index of the index of facets in each object which are larger than requested
    for count in range(len(object_list)):
        areas.append(object_list[count].areas)
        large_facets.append([v for v,x in enumerate(areas[count]) if x >= area_limit])
        object_vectors=np.zeros([len(object_list[count].data['vectors']),3,3])
        object_vectors=object_list[count].data['vectors']
        trimmed_object_vectors=np.delete(object_vectors,[large_facets[count]],0)
        for count2 in range(len(large_facets[count])):
            facet_vectors=object_list[count].data['vectors'][large_facets[count][count2]]
            new_object_vectors=np.zeros([6,3,3])
            new_object_vectors=triangle_split(facet_vectors)
            trimmed_object_vectors=np.append(trimmed_object_vectors,new_object_vectors,axis=0)
        new_object=np.zeros(len(trimmed_object_vectors), dtype=mesh.Mesh.dtype)
        new_object['vectors']=trimmed_object_vectors         
        fixed_list[count]=mesh.Mesh(new_object.copy())
        fixed_list[count].update_areas

    return fixed_list```

[wolph]

Well, that indeed looks a tad broken, I can't really tell you what the problem is straight up but I see at least 1 small error in the code. You're not calling update_areas, it's a function that needs to be called for it to work.

[timtitan]

I have found the problem with the code, the algorithm I was using for generating the central point was returning NaN instead of the coordinates in many circumstances, and this has been fixed now.