What is the 3d location of the tip of a triangular pyramid with faces identically angled from vertical but an arbitrarily tilted base face?

Question

How can I find the 3d location of the tip of a pyramid which has a 3 sided base arbitrarily angled relative to the xy plane and remaining three faces with identical anglation relative to the z-axis.

The known variables are:

The three xyz vertices for the base face

The angle relative to the z-axis for the other three faces

(this would be acceptible if fixed at pi/4, but much better if it could be determined on a case-by-case basis)

Why?

I am automatically generating supports for printing of a 3d modelled structure. In 3d printing, faces which "overhang" (ie are inverted and within approx 45 degrees of horizontal) often fail to print correctly without support.

The intention is to develop a more robust support structure for these faces which spreads force more evenly and thus requires less point adhesion to the printed object which therefore reduces surface damage when removed.

The option being currently tested is a series of inverted triangular pyramids. The pyramid faces will have an identical angle relative to horizontal (typically 45 deg, but it would be ideal if this angle could be arbitrarily set between 0 and PI/2). To minimise waste and optomise quality this angle should be identical for each of the 3 non-base faces irrespective of the orientation of the base.

Although many models to have faces with more than 3 sides, for practical purposes any faces with more than 3 sides will have previously been split into multiple 3-sided "base" faces which are each processed individually.

I have coded the following with outputs as per linked graphics which is probably "good enough". However although I don't believe it is robust (as evidenced when recursively looped), I haven't managed to get closer to a better solution.

triangular pyramidal supports underneath hand

triangular pyramidal supports underneath icosohedron

def support_interface(faceVertices, maxAngle):
    vertices = []
    for vert in faceVertices:
        vertices.append((vert[X], vert[Y], vert[Z]))
# find a point which is maxAngle (currently locked at 45 degrees) below and between each pair
#   vertices loops anticlockwise around the face
# then average this to get an approximate intersection point
#potentially repeat the process to successively approximate the solution

curV = 0                 # which vertex pair is currently being evaluated?
nextApproxVert = []      # temp list holding the most recent full iteration around all the vertex pairs
closeEnough = False      # not really needed - but set up for future successive approximation
loopCnt = 1              # the maximum number of loops to iterate over
#maxAngle = -3.14159/4   # locked at 45 degrees for testing
PI_2 = -3.14159/2        # 90 degrees

while not closeEnough:

    # load the next vertex pair in the loop around the edge of the face
    x0 = vertices[curV][0]
    y0 = vertices[curV][1]        
    z0 = vertices[curV][2]        

    if curV < vertexCnt -1:
        x1 = vertices[curV+1][0]        
        y1 = vertices[curV+1][1]        
        z1 = vertices[curV+1][2]        
    else:
        x1 = vertices[0][0]        
        y1 = vertices[0][1]        
        z1 = vertices[0][2]        

    # get the xyz axial distances of the triangle
    dx = x0 - x1
    dy = y0 - y1
    dz = z0 - z1

    #print("dxyz=",dx,dy,dz)


    # get the 2D distance on the z-plane between vertices
    xyLen = math.sqrt(dx*dx + dy*dy)
    #print("xyLen",xyLen)
    # get the 3D distance between vertices
    h0 = math.sqrt(dx*dx + dy*dy + dz*dz)
    #print("h0",h0)
    # get the angle of the hypotenuse relative to horizontal
    a0 = abs(math.asin(dz/h0))
    #print("a0",a0)
    #get the angle of the additional triangle needed to keep angle under max value
    a1 = PI_2 - a0 - maxAngle
    #print("a1",a1)
    #get the length of the opp and adj sides of the "safe" triangle
    opp = math.sin(a1) * h0   # the z drop from the lower vertex to the "safe" vertex
    adj = math.cos(a1) * h0
    #print("opAd",opp, adj)
    #now get the other side lengths for the triangle from the higher vertex to the safe vertex
    # at 90deg to horizontal
    zHigh = math.asin(PI_2 - maxAngle) * adj   # adj for the previous triangle is hyp of the new one
    zAdj = math.acos(PI_2 - maxAngle) * adj
    #print("zHigh/zAdj",zHigh,zAdj)

    xyProp = zAdj / xyLen
    #print("xyProp",xyProp)

    dx1 = dx * xyProp
    dy1 = dy * xyProp
    #print("dxy1",dx1,dy1)

    nextApproxVert.append((x0 - dx1, y0 - dy1, z0 + zHigh))

    curV +=1
    if curV >= vertexCnt:
        curV = 0
        loopCnt -= 1
        #print(nextApproxVert)
        vertices = nextApproxVert           
        nextApproxVert = []
        if loopCnt <= 0: 
            closeEnough = True

# after it meets criteria for closeEnough then take the average of vertices           
(x,y,z)=(0,0,0)
for v in vertices:           
    x += v[0]
    y += v[1]
    z += v[2]
x=x/vertexCnt
y=y/vertexCnt
z=z/vertexCnt

return (x,y,z)

Answer 1

First, implement a function that computes a plane that passes through two given points and whose normal makes some given angle with the z-axis.

Here’s how to do that: suppose the two given points are $P$ and $Q$, and we want the normal of our plane to make an angle $A$ with the z-axis. Let the plane normal be given by the unit vector $N = (u,v,w)$. The equation of the plane can be written as $(X-P)\cdot N = 0$. All we need to do is figure out $N$. Since the vector $Q-P$ lies in the plane, we have $(Q-P)\cdot N = 0$. The angle requirement tells us that $N \cdot(0,0,1)= \cos A$, so $w= \cos A$. Also, since $N$ is a unit vector, $u^2 + v^2 +w^2 = 1$. These three equations give us enough information to calculate $N$.

Use this function to get the three planar sides of your “pyramid” (the correct term is “tetrahedron”, actually).

The point you want is the intersection of these three planes. So, you need to write a function that intersects three planes, if you don’t already have one.

Answer 2

The full code is here: https://github.com/dgm3333/3d-print-toolbox-modified/blob/master/supports.py

However the primary python routine which takes a list of face vertices and creates the support is as follows. Since it was a learning experience for me it's got lots of comments... Because bubba's advice got me there I've marked that answer as correct.

# This creates a triangular pyramid with one face in contact with the original face
#   and the other faces angled to ensure they do not exceed a 
#   predetermined maximum overhang angle
# It starts by calculating the length of a line perpendicular from an 
#   edge from the face to the centroid.
# It then calculates the height the line would intercept a vertical axis 
#   running through the centroid at maxAngle from vertical.
# The three vertices are then used to generate a plane from that edge
# This is repeated to generate three planes.
# The location at which all three intersect provides the 4th vertex to form the pyramid...
def support_interface(faceVerts, maxAngle):
if len(faceVerts) != 3:
    return

# will create with faces feeding under the centroid
# This isn't the ideal point as angle from some vertices will be much shallower than required
centroid =  Point((faceVerts[0].x + faceVerts[1].x + faceVerts[2].x)/3, 
             (faceVerts[0].y + faceVerts[1].y + faceVerts[2].y)/3, 
             (faceVerts[0].z + faceVerts[1].z + faceVerts[2].z)/3)

# two points forming a z-axis directly though the centroid
pZ1 = Point(centroid.x, centroid.y, centroid.z-1000)
pZ2 = Point(centroid.x, centroid.y, centroid.z+1000)


# loop through each edge 
#   identify the z intercept (using the centroid as the z-axis)
#   generate the 3 points to define the plane
#   add the plane generated to the list of planes
planes = []
for v in range(3):
    p1 = faceVerts[v]
    if v == 2:
        p2 = faceVerts[0]   # it's a loop so vertex[0] is the next in line from vertex[2]
    else:
        p2 = faceVerts[v+1]


    # identify the intersection of the lines
    # This is calculated as the dot product of the line defined by the two points with the centroid z-axis
    # Blender could also do this step:  mathutils.geometry.intersect_line_line(v1, v2, v3, v4)
    lli = LineLineIntersect3D(p1,p2,pZ1,pZ2)


    if (lli.uv.x == 0.0 and lli.uv.x == 0.0 and lli.uv.x == 0.0):
        #print("ERROR: unit vector has length zero. This means both lines are in contact.")
        #print("  This may be a zero-area face (but then it shouldn't be overhanging so seems odd), as it should never happen using the centroid z-axis (or the face wouldn't be a triangle)")
        #print("  aborting calculation on this face")
        return
        #print(" NB it would be a risk if using the true z-axis")
        #print("Thus the z-axis intercept cannot be used to provide the 3rd point of the plane")
        #print("And in this situation could try using (1,1,0)->(1,1,1) line instead)")
        #lli = LineLineIntersect3D(p1,p2,Point(1,1,0),Point(1,1,1))

    #print("Calculating plane from p1p2 edge with coordinates of:")
    #print(p1)
    #print(p2)
    #print("UV (unit vector):", lli.uv)

    # the point on p1p2 will be exactly maxAngle from the point the plane intersects the centroid z-axis
    #print("Nearest point on p1p2 to z-axis:",lli.Pmem1)
    # the hypotenuse of a horizontal xy triangle becomes the adjacent side for the triangle down to the intercept
    # thus the opposite side (the difference in z-axis height from the current point is Opp = Adj*tan(maxAngle)
    #formula if using the z-axis at x=0,y=0
    #xyHyp = sqrt(lli.Pmem1.x*lli.Pmem1.x + lli.Pmem1.y*lli.Pmem1.y) 
    dx = lli.Pmem1.x - centroid.x
    dy = lli.Pmem1.y - centroid.y
    #dz = lli.Pmem1.z - centroid.z
    print("dx=",dx,"dy=",dy)
    xyHyp = sqrt(dx*dx + dy*dy)

    print("xy hypotenuse = xyz adjacent side:",xyHyp)
    opp = abs(xyHyp * math.tan(maxAngle))
    if opp==0.0:
        #print("ERROR: this function will not work for horizontal planes")
        #print("  This should never happen if using centroid and non-zero face size")
        #print("Function will ignore this face until script is updated to handle this")
        return            

    # Because we are using the centroid for our z-axis, we know the correct point must be below the height calculated for the intersection
    # and the normal N.z needs to point up to be towards the centroid
    p3 = Point(centroid.x, centroid.y, lli.Pmem1.z-opp)
    planeC = Plane3D(p1, p2, p3)
    if planeC.N.z < 0:
        planeC.N = -planeC.N

    planes.append(planeC)
    print("Verts of plane", p1, p2, p3)

    '''
    # If using the true z-axis the check is a bit more fiddly we need to test which of the two points produces a plane in the correct orientation
    #print("z intercept should be +/-",opp)
    # now we have two potential intercepts (one above and one below the level of the intersection)
    p3a = Point(0,0,lli.Pmem1.z+opp)
    p3b = Point(0,0,lli.Pmem1.z-opp)

    #print("zIntercept candidates are:-")     
    print(p3a)
    print(p3b)

    # Now there are four potential normals for the plane.

    # Selecting the normal which is on the same side as one of the other face points and which has uv.z >0 (or it's negative) 
    #   will identify the correct one. 
    # Use the face centron to identify whether the normal is on the same or opposide side of the plane as the face
    #   and whether the normal points up or down (it should point up towards the face or down away)
    #print ("Testing plane candidate generated using point p3a")
    planeC = Plane3D(p1, p2, p3a)
    correctPlane = False
    lieCheckA = planeC.lie_check(centroid)

    if lieCheckA == 0.0:
        #print("Warning: The face centroid lies on current plane.")
        #print("  The face  may be at exactly maxAngle (but if this is the case it shouldn't have been selected for support)")
    elif lieCheckA < 0.0:
        #print ("The centroid lies on OPPOSITE side of current plane as normal N.")
        if planeC.N.z < 0:
            #print ("The normal N.z points DOWN so this is the correct plane")
            correctPlane = True
            planeC.N = -planeC.N
    else:
        #print ("The centroid lies on SAME side of current plane as normal N.")
        if planeC.N.z > 0:
            #print ("The normal N.z points UP so this is the correct plane")
            correctPlane = True
            planeC.N = -planeC.N

    if not correctPlane:        
        #print ("p3a plane didn't meet selection criteria")
        #print ("Testing plane candidate generated using  point p3b")
        planeC = Plane3D(p1, p2, p3b)
        correctPlane = False
        lieCheckA = planeC.lie_check(centroid)
        if lieCheckA == 0.0:
            #print ("ERROR: The face centroid lies on current plane - THIS SHOULDN'T HAPPEN IF THE PLANE IS CORRECT")
        elif lieCheckA < 0.0:
            #print ("The centroid lies on OPPOSITE side of current plane as normal N.")
            if planeC.N.z < 0:
                #print ("The normal N.z points DOWN so this is the correct plane")
                correctPlane = True
                planeC.N = -planeC.N
        else:
            #print ("The centroid lies on SAME side of current plane as normal N.")
            if planeC.N.z > 0:
                #print ("The normal N.z points UP so this is the correct plane")
                correctPlane = True
                planeC.N = -planeC.N

    if not correctPlane:
        #print ("ERROR: Neither plane matches - THIS SHOULDN'T HAPPEN IF THE PLANE IS CORRECT")
        return
    else: 
        planes.append(planeC)
    '''

# Now we should have three planes we just need to establish the intersect point
if len(planes) == 3:
    #print("3 PLANES IDENTIFIED")
    pass
else:
    #print("ERROR: Number of planes != 3")
    #print("  Not sure why as zero area faces should already have been filtered, but function will ignore this face")
    return

# The intersection point of the three planes "M" is given by:
# M = (d1*(N1 cross N2) + d2*(N2 cross N0) + d3*(N0 cross N1)) / (N0 dot (N1 cross N2))
# 'cross' indicates the cross product and 'dot' indicates the dot product
# blender could do this step with:
# blender could do this with
# mathutils.geometry.intersect_plane_plane(plane_a_co, plane_a_no, plane_b_co, plane_b_no)
#print("Normals and D coefficients for the three planes are:-")   
N0 = planes[0].N
N1 = planes[1].N
N2 = planes[2].N
D0 = planes[0].D
D1 = planes[1].D
D2 = planes[2].D
#k0 = round(N0.dot(planes[0].p1), 6)
#k1 = round(N1.dot(planes[1].p1), 6)
#k2 = round(N2.dot(planes[2].p1), 6)
#print(N0.x, N0.y, N0.z, "   ", D0)#, "   ", k0)
#print(N1.x, N1.y, N1.z, "   ", D1)#, "   ", k1)
#print(N2.x, N2.y, N2.z, "   ", D2)#, "   ", k2)

N12 = N1.cross(N2)
N20 = N2.cross(N0)
N01 = N0.cross(N1)
N0dN12 = round(N0.dot(N12), 6)


#print("Original Face Vertices and Centroid are:-")
#for vert in faceVerts:
    #print(vert.x, vert.y, vert.z)
#print(centroid.x, centroid.y, centroid.z, "(=centroid)")


#print("Intersection Point Is:-")    
# TODO: figure out why Point overloads aren't working
# M = (D0*N12 + D1*N20 + D2*N01) / N0dN12
# point overloads not working, so calculate individually first
D0N12 = Point(N12.x*D0, N12.y*D0, N12.z*D0)    
D1N20 = Point(N20.x*D1, N20.y*D1, N20.z*D1)    
D2N01 = Point(N01.x*D2, N01.y*D2, N01.z*D2)    

CPs = D0N12 + D1N20 + D2N01


pyramidTip = Point(CPs.x / N0dN12, CPs.y / N0dN12, CPs.z / N0dN12)   # = M
#print(pyramidTip.x, pyramidTip.y, pyramidTip.z)

# now create a list with the vertices from the original face + the new vertex
verts = []
for vert in faceVerts:
    verts.append(vert)
faceVertCnt = len(verts)
verts.append(pyramidTip)
for v in verts:
    print(v)
##print("faceVertCnt=",faceVertCnt)

faces = []
faces.append((range(faceVertCnt)))
faces.append((0,faceVertCnt-1,faceVertCnt))
for i in range(0,faceVertCnt-1):
    faces.append((i+1,i,faceVertCnt))

##print("Verts=",verts1)
##print("Faces=",faces1)
obj = createMesh('interface', (0,0,0), verts, [], faces)

bpy.ops.collection.objects_remove_all()
bpy.data.collections['AutoSupports'].objects.link(obj)
print()
return pyramidTip

```