Create surface or mesh from the point cloud of the surface of a color volume in CIELAB

I would like to be able to calculate the coverage of color spaces by other color spaces or display gamuts in CIELAB, so with color volumes, rather than the more common and far easier color triangles in the chromaticity diagrams. And I'd also like to make nice 3D plots of those color spaces in CIELAB.

Using the math on brucelindbloom.com I've managed to get a point cloud of the surface of a color volume (code below) in the CIELAB space, but that's where I get stuck.

tic
format long
N = 7;        %color depth per channel in bits
Gamma = 2.2;
Rx = 0.64;        %Rx-By are the primaries of Adobe RGB color space in CIE 1931 xy chromaticity diagram
Ry = 0.33;
Gx = 0.21;
Gy = 0.71;
Bx = 0.15;
By = 0.06;
Xw = 0.95047;
Yw = 1;
Zw = 1.08883;
max_lvl = 2^N-1;
R = 0:1:max_lvl;
R = transpose(R);
G = R;
B = R;
RGB = allcomb(R,G,B);
RGB_norm = RGB/max_lvl;
surface_mask = (RGB_norm(:,1)>0 & RGB_norm(:,1)<1 & RGB_norm(:,2)>0 & RGB_norm(:,2)<1 & RGB_norm(:,3)>0 & RGB_norm(:,3)<1);
RGB_gamma = RGB_norm.^Gamma;
clear RGB RGB_norm
Xr = Rx/Ry;
Yr = 1;
Zr = (1-Rx-Ry)/Ry;
Xg = Gx/Gy;
Yg = 1;
Zg = (1-Gx-Gy)/Gy;
Xb = Bx/By;
Yb = 1;
Zb = (1-Bx-By)/By;
XYZ_rgb = [Xr Xg Xb; Yr Yg Yb; Zr Zg Zb];
XYZ_rgb_inv = inv(XYZ_rgb);
XYZ_w = [Xw; Yw; Zw];
S_rgb = XYZ_rgb_inv * XYZ_w;
M= zeros(size(XYZ_rgb));
M(:,1) = XYZ_rgb(:,1) * S_rgb(1,1);
M(:,2) = XYZ_rgb(:,2) * S_rgb(2,1);
M(:,3) = XYZ_rgb(:,3) * S_rgb(3,1);
XYZ = transpose(M * transpose(RGB_gamma));
epsilon = 216/24389 * ones(size(XYZ,1),1);
kappa = 24389/27;
x_r = XYZ(:,1)/XYZ_w(1,1);
y_r = XYZ(:,2)/XYZ_w(2,1);
z_r = XYZ(:,3)/XYZ_w(3,1);
clear XYZ

mask = x_r > epsilon;
x_r( mask ) = x_r( mask ) .^ (1/3);
x_r( ~mask ) = ( kappa * x_r( ~mask ) + 16 ) / 116;

mask = y_r > epsilon;
y_r( mask ) = y_r( mask ) .^ (1/3);
y_r( ~mask ) = ( kappa * y_r( ~mask ) + 16 ) / 116;

mask = z_r > epsilon;
z_r( mask ) = z_r( mask ) .^ (1/3);
z_r( ~mask ) = ( kappa * z_r( ~mask ) + 16 ) / 116;

L = 116*y_r-16;
a = 500*(x_r-y_r);
b = 200*(y_r-z_r);
clear epsilon x_r y_r z_r mask
toc
tic
figure(1);
abL = horzcat(a,b,L);
abL(surface_mask,:) = [];
RGB_gamma(surface_mask,:) = [];
patchinfo.Vertices = abL;
patchinfo.Faces = [1:size(abL,1)]';
p = patch(patchinfo,'Marker','.','MarkerSize',1,'EdgeColor','interp','FaceVertexCData',RGB_gamma);
axis equal;
axis([min(a)-10 max(a)+10 min(b)-10 max(b)+10 min(L)-10 max(L)+10]);
grid on;
whitebg([0.5,0.5,0.5]);
view(3);
toc

I've asked Bruce Lindbloom what method he uses for calculating the values in the first table on the RGB Working Space Information page of his site and he said that he tessellates the gamut surfaces into tiny triangles to form a closed polyhedron (and from there he calculates the inscribed volume, etc.).

I've tried doing that with Delaunay triangulation, but the result is not what I want, because it doesn't follow the surface (or rather the point cloud of it) properly. Concave areas turn out straight, convhull has the same problem.

tri = delaunay(a,b,L);
figure(2);
trisurf(tri,a,b,L,'EdgeColor','none','LineStyle','none','FaceLighting','gouraud');
axis equal;
axis([min(a)-10 max(a)+10 min(b)-10 max(b)+10 min(L)-10 max(L)+10]);
grid on;
whitebg([0.5,0.5,0.5]);
view(3);

And how can it be that I end up with just over 6 times as many rows in tri (that's the number of faces / triangles, right?) as the number of points (vertices), showing only minor variations with a different color bit depth N (tested from 3 to 8)? That doesn't sound right...

So how do I make a proper (triangulated) surface from the point cloud of the surface? If possible also in the right colors. RGB_gamma describes the color of all points / vertices, these should be easy to interpolate across the surface.

Calculating the inscribed volume and subsequent steps are of later concern. I might even use a completely different method for that, like the Monte Carlo method used by Mark Meyer for a similar purpose here.