Solution with matplotlib:

#!/usr/bin/python3

import sys

import matplotlib
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

import numpy
from numpy.random import randn
from scipy import array, newaxis


# ======
## data:

DATA = array([
    [-0.807237702464, 0.904373229492, 111.428744443],
    [-0.802470821517, 0.832159465335, 98.572957317],
    [-0.801052795982, 0.744231916692, 86.485869328],
    [-0.802505546206, 0.642324228721, 75.279804677],
    [-0.804158144115, 0.52882485495, 65.112895758],
    [-0.806418040943, 0.405733109371, 56.1627277595],
    [-0.808515314192, 0.275100227689, 48.508994388],
    [-0.809879521648, 0.139140394575, 42.1027499025],
    [-0.810645106092, -7.48279012695e-06, 36.8668106345],
    [-0.810676720161, -0.139773175337, 32.714580273],
    [-0.811308686707, -0.277276065449, 29.5977405865],
    [-0.812331692291, -0.40975978382, 27.6210856615],
    [-0.816075037319, -0.535615685086, 27.2420699235],
    [-0.823691366944, -0.654350489595, 29.1823292975],
    [-0.836688691603, -0.765630198427, 34.2275056775],
    [-0.854984518665, -0.86845932028, 43.029581434],
    [-0.879261949054, -0.961799684483, 55.9594146815],
    [-0.740499820944, 0.901631050387, 97.0261463995],
    [-0.735011699497, 0.82881933383, 84.971061395],
    [-0.733021568161, 0.740454485354, 73.733621269],
    [-0.732821755233, 0.638770044767, 63.3815970475],
    [-0.733876941678, 0.525818698874, 54.0655910105],
    [-0.735055978521, 0.403303715698, 45.90859502],
    [-0.736448900325, 0.273425879041, 38.935709456],
    [-0.737556181137, 0.13826504904, 33.096106049],
    [-0.738278724065, -9.73058423274e-06, 28.359664343],
    [-0.738507612286, -0.138781586244, 24.627237837],
    [-0.738539663773, -0.275090412979, 21.857410904],
    [-0.739099040189, -0.406068448513, 20.1110519655],
    [-0.741152200369, -0.529726022182, 19.7019157715],
])

Xs = DATA[:,0]
Ys = DATA[:,1]
Zs = DATA[:,2]


# ======
## plot:

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

surf = ax.plot_trisurf(Xs, Ys, Zs, cmap=cm.jet, linewidth=0)
fig.colorbar(surf)

ax.xaxis.set_major_locator(MaxNLocator(5))
ax.yaxis.set_major_locator(MaxNLocator(6))
ax.zaxis.set_major_locator(MaxNLocator(5))

fig.tight_layout()

plt.show() # or:
# fig.savefig('3D.png')

Result:

Probably not very beautiful. But it will be, if You provide more points.

For surfaces it's a bit different than a list of 3-tuples, you should pass in a grid for the domain in 2d arrays.

If all you have is a list of 3d points, rather than some function f(x, y) -> z, then you will have a problem because there are multiple ways to triangulate that 3d point cloud into a surface.

Here's a smooth surface example:

import numpy as np
from mpl_toolkits.mplot3d import Axes3D  
# Axes3D import has side effects, it enables using projection='3d' in add_subplot
import matplotlib.pyplot as plt
import random

def fun(x, y):
    return x**2 + y

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
x = y = np.arange(-3.0, 3.0, 0.05)
X, Y = np.meshgrid(x, y)
zs = np.array(fun(np.ravel(X), np.ravel(Y)))
Z = zs.reshape(X.shape)

ax.plot_surface(X, Y, Z)

ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')

plt.show()

It is possible to plot the 3D surface over your scatter plot using the plt.plot_trisurf(...) function as long as you find the right ordering of vertices for the triangles. There is a function from SciPy called ConvexHull which finds the simplices of the points on the outside of the data set. This is very handy, but does not immediately work on this example because your data set is not convex!

The solution is to make the data convex by expanding the points away from the center until they form a sphere. See below for a visualization of this.

After turning the head into a sphere it is now possible to call ConvexHull(...) to get the desired triangulation. This triangulation can be applied to the spherical head first (see below). Then, the head can be shrunk back into its original form, and the triangulation's vertices are still valid!

This is the final product!

Code

import numpy as np
import matplotlib.pyplot as plt
import csv
from scipy.spatial import KDTree
from scipy.spatial import ConvexHull
from matplotlib import cm
from matplotlib import animation

plt.style.use('dark_background')

# Data reader from a .csv file
def getData(file):
    lstX = []
    lstY = []
    lstZ = []
    with open(file, newline='\n') as f:
        reader = csv.reader(f, quoting=csv.QUOTE_NONNUMERIC)
        for row in reader:
            lstX.append(row[0])
            lstY.append(row[1])
            lstZ.append(row[2])
    return lstX, lstY, lstZ

# This function gets rid of the triangles at the base of the neck
# It just filters out any triangles which have at least one side longer than toler
def removeBigTriangs(points, inds, toler=35):
    newInds = []
    for ind in inds:
        if ((np.sqrt(np.sum((points[ind[0]]-points[ind[1]])**2, axis=0))<toler) and
            (np.sqrt(np.sum((points[ind[0]]-points[ind[2]])**2, axis=0))<toler) and
            (np.sqrt(np.sum((points[ind[1]]-points[ind[2]])**2, axis=0))<toler)):
            newInds.append(ind)
    return np.array(newInds)

# this calculates the location of each point when it is expanded out to the sphere
def calcSpherePts(points, center):
    kdtree = KDTree(points) # tree of nearest points
    # d is an array of distances, i is array of indices
    d, i = kdtree.query(center, points.shape[0])
    spherePts = np.zeros(points.shape, dtype=float)
    
    radius = np.amax(d)
    for p in range(points.shape[0]):
        spherePts[p] = points[i[p]] *radius /d[p]
    return spherePts, i # points and the indices for where they were in the original lists
    

x,y,z = getData(".\coords3Ddetailed.csv")

pts = np.stack((x,y,z), axis=1)

# generating data
spherePts, sphereInd = calcSpherePts(pts, [0,0,0])
hull = ConvexHull(spherePts)
triangInds = hull.simplices # returns the list of indices for each triangle
triangInds = removeBigTriangs(pts[sphereInd], triangInds)

# plotting!
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter3D(pts[:,0], pts[:,1], pts[:,2], s=2, c='r', alpha=1.0)
ax.plot_trisurf(pts[sphereInd,0], pts[sphereInd,1], pts[sphereInd,2], triangles=triangInds, cmap=cm.Blues, alpha=1.0)
plt.show()

You can do this using ax.plot_trisurf to construct the surface from triangles:

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

X = np.array([-4.3, -3, 1.2, 1.1])
Y = np.array([-2.3, 2, -4.5, 0.8])
Z = np.array([1., 1., 1., 1.])

fig = plt.figure(figsize=(8,8))
ax = plt.axes(projection='3d')
surf1 = ax.plot_trisurf(X, Y, Z, antialiased=True)

plt.show()

The output:

If you want to also remove the line dividing the triangles, pass the same color and edgecolor to plot_trisurf.