How can I generate a concave hull of 3D points?

I have 68 3d points (I'm guessing it's called a sparse point cloud). I want to connect them all together to create a mesh. I first tried this using Delaunay Triangulation. However, it didn't work well because Delaunay Triangulation only gives a convex hull, so it is giving a mesh that basically ignores the eyes which are concave for example.

Here is a picture illustrating what I mean:

So I tried using something else which is alphashape. I've been using this documentation: https://pypi.org/project/alphashape/ My problem is that it's simply just not working. Here are some pictures: The above picture shows the 3d points which I want to convert to a mesh Below pictures show the result of me using alpha shape!

I want to get a 3D convex hull.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import alphashape 

points_3d = np.array(sixtyEightLandmarks3D)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(points_3d[:, 0], points_3d[:, 1], points_3d[:, 2])
plt.show()


points_3d = [
    (0., 0., 0.), (0., 0., 1.), (0., 1., 0.),
    (1., 0., 0.), (1., 1., 0.), (1., 0., 1.),
    (0., 1., 1.), (1., 1., 1.), (.25, .5, .5),
    (.5, .25, .5), (.5, .5, .25), (.75, .5, .5),
    (.5, .75, .5), (.5, .5, .75)
]

points_3d = [
    (7, 191, 325.05537989702617), (6, 217, 330.15148038438355), (8, 244, 334.2528982671654), 
    (11, 270, 340.24843864447047), (19, 296, 349.17330940379736), (34, 320, 361.04985805287333), 
    (56, 340, 373.001340480165), (80, 356, 383.03263568526376), (110, 361, 387.06330231630074), 
    (140, 356, 383.08354180256816), (165, 341, 373.1621631409058), (187, 321, 359.4022815731698), 
    (205, 298, 344.64039229318433), (214, 272, 334.72376670920755), (216, 244, 328.54984401152893), 
    (218, 217, 324.34703636691364), (217, 190, 319.189598828032), (22, 166, 353.0056656769123), (33, 152, 359.0055709874152), 
    (52, 145, 364.0), (72, 147, 368.00135869314397), (91, 153, 372.0013440835933), (125, 153, 370.0013513488836), 
    (145, 146, 366.001366117669), (167, 144, 361.0), (186, 151, 358.00558654859003), (197, 166, 351.0056979594491), 
    (108, 179, 376.02127599379264), (108, 197, 381.02099679676445), (109, 214, 387.03229839381623), 
    (109, 233, 393.04579885809744), (87, 252, 383.03263568526376), (98, 255, 386.0323820614017), (109, 257, 387.03229839381623), 
    (120, 254, 385.0324661635691), (131, 251, 383.0469945058961), (44, 183, 360.01249978299364), (55, 176, 363.00550960006103), 
    (69, 176, 363.0123964825444), (81, 186, 364.0219773585106), (69, 188, 364.0219773585106), (54, 188, 364.0219773585106), 
    (136, 185, 361.01246515875323), (147, 175, 362.0013812128346), (162, 175, 361.00554012369395), 
    (174, 183, 357.0014005574768), (163, 188, 360.01249978299364), (149, 188, 362.01243072579706), 
    (73, 289, 384.04687213932624), (86, 282, 389.0462697417879), (100, 278, 391.0319680026174), 
    (109, 281, 391.0460330958492), (120, 277, 391.0319680026174), (134, 281, 387.03229839381623), 
    (147, 289, 380.0210520484359), (135, 299, 388.03221515745315), (121, 305, 392.04591567825315), 
    (110, 307, 392.04591567825315), (100, 306, 392.04591567825315), (86, 300, 391.0626548265636), 
    (78, 290, 386.0207248322297), (100, 290, 391.0319680026174), (109, 291, 391.0319680026174), 
    (120, 289, 391.0319680026174), (142, 289, 381.03280698648507), (120, 290, 391.0319680026174), 
    (109, 292, 392.03188645823184), (100, 291, 392.04591567825315), 
    (0., 1., 1.), (1., 1., 1.), (.25, .5, .5),
    (.5, .25, .5)
]



alpha_shape = alphashape.alphashape(points_3d, lambda ind, r: 0.3 + any(np.array(points_3d)[ind][:,0] == 0.0))
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.plot_trisurf(*zip(*alpha_shape.vertices), triangles=alpha_shape.faces)
plt.show()

Note: I am getting the points of the face and then putting it in the 3d_points variable manually. Also for some reason, I found that if I don't add: "(0., 1., 1.), (1., 1., 1.), (.25, .5, .5), (.5, .25, .5)" to the array, then alphashape won't work and will give me the following error: too many indices for array: array is 1-dimensional, but 2 were indexed

The reason for the error you got is that you used an alpha value (alpha=0.3 in your case) that was inappropriate for the input you had. This caused an empty 3D alpha shape with no triangles, and the alphashape package probably has a bug that doesn't handle this case gracefully - hence the cryptic error message. However, for the four points you added the alpha value of 0.3 was appropriate, so when you add them what you actually see is just the alpha-shape of these four points (notice the scale in your lower figure).

The alpha value in the package denotes 1 / radius of a circumscribing circle (for 2D) or sphere (for 3D) of a Delaunay simplex (triangle for 2D and tetrahedron for 3D). The simplices that remain in the alpha-shape are those with radius smaller than 1 / alpha (see here for example). In your case of alpha=0.3 the radius is approximately 3 and since your original points are much further apart all the Delaunay tetrahedra were purged away and the alpha-shape became empty as I explained above.

The following figures are the result of building the alpha-shape with your original points (without the additional four points) but with different alpha values.

Calling alpha_shape = alphashape.alphashape(points_3d, 0.02) with your drawing code gives the first figure, and calling with alpha=0.007 gives the second. Notice how a smaller alpha value (i.e., a larger radius) results in more triangles being part of the alpha-shape.

In some 3D Delaunay applications, people find it suitable to perform the Delaunay construction in 2D and elevate the 3D points back to 3D with the triangles from the 2D triangulations. For example, for modelling terrains these are usually the methods of choice. From your answer I understand that your application is similar.

Unfortunately I couldn't find an interface in the alphashape package that does this 2D filtering and returns the filtered triangles to be used on the elevated points. The following code I wrote builds a 2D alpha-shape for the 2D projections of your input. It follows more-or-less similar code I've written in past StackOverflow answers (e.g., here and here), but is modified and refactored for the needs of your specific problem.

The function alpha_shape_filter below basically builds the 2D Delaunay triangulation and filters out triangles with radius > alpha_radius.

from scipy.spatial import Delaunay

def circ_radius(p0,p1,p2):  
    """Computing the radius of the circumcircle of three points."""
    a = p1-p0
    b = p2-p0

    norm_a = np.linalg.norm(a, axis=1)
    norm_b = np.linalg.norm(b, axis=1)
    norm_a_b = np.linalg.norm(a-b, axis=1)
    cross_a_b = np.cross(a,b)
    return (norm_a*norm_b*norm_a_b) / np.abs(2.0*cross_a_b)


def alpha_shape_filter(points, alpha_radius):
    """
    Compute the triangles of that are part of the alpha shape (concave hull) of a set of points.
    :param points: np.array of shape (n,2) points.
    :param alpha_raius: alpha radius value.
    :return: scipy delaunay triangulation and boolean array corresponding to triangles in alpha shape.
    """
    assert points.shape[0] > 3, "Need at least four points"
    dt = Delaunay(points)

    p0 = points[dt.vertices[:,0],:]
    p1 = points[dt.vertices[:,1],:]
    p2 = points[dt.vertices[:,2],:]

    rad1 = (circ_radius(p0, p1, p2))

    is_in_shape_filter = (rad1 < alpha_radius)
    return dt, is_in_shape_filter

Calling the function (on your projected input points) with the following code results in the figure below.

dt, is_in_shape = alpha_shape_filter(points_3d[:, [0,1]], 1 / 0.02)
simplices = dt.simplices[is_in_shape]

plt.figure()
triang = mtri.Triangulation(dt.points[:, 0], dt.points[:, 1], simplices)
plt.triplot(triang, 'b-')
plt.plot(points_3d[:, 0], points_3d[:, 1], "ok")

When elevating the points back and plotting the results with the following code we get the following 3D figure below.

fig = plt.figure()
ax = plt.axes(projection='3d')
ax.plot_trisurf(points_3d[:, 0], points_3d[:, 1], triangles=simplices, Z=points_3d[:, 2])
ax.scatter(points_3d[:, 0], points_3d[:, 1], points_3d[:, 2], c="k")