Quantifying a manifold folding unto itself I have a dataset of ~7k scattered points in 3D which represents a manifold that may or may not "fold unto itself". Here's an example where this does happen (look at the top-right yellow triangles):

My ultimate goal is to have a metric to compare many such datasets, in order to optimize the hyperparameters leading to the formation of this hypersurface/manifold/volume, while minimizing the unwanted folding. 
One way to detect whether such folding occurs is by examining the errors when performing leave-one-out cross-validation. Specifically, we can quantify the folding, by a) interpolating the value (color) of each point based on k nearest neighbors; b) comparing it the interpolated and correct value; and c) aggregating the errors.
When I tried implementing this in MATLAB, using a combination of fitcknn, crossval and kfoldLoss, processing took on the order of minutes - which is unacceptably slow - as I must be able to perform this error computation on the order of a second or less. I know about CUDA-based algorithms for accelerated knn search1. However, as I don't come from this field I can't tell if I can leverage a fast knn search algorithm to speed up the leave-one-out error computation. (In fact, I don't even know if this is the most suitable approach to the problem.)
Is there something "costly" I can compute once for each dataset (kdtree? pairwise distances?) and use that for fast computation of the aggregate error? How else can I speed this up?
I would also be happy to hear about any other algorithm that could quantify the folding effect.
 A: I've ended up approximating cross-validation using inverse-distance-weighted (IDW) interpolation using the nearest neighbors. This produced useful results and was also fast to compute. Below is the outline of the implementation:
Part 1: Find K nearest neighbors of each point


*

*Decide on the amount of nearest neighbors, K, to use.

*Compute squared Euclidean distances between all points. (Note: this is slightly faster than the actual distances, since we avoid taking the root.)

*Establish the indices of the K nearest neighbors, for each point.


Part 2: Find internal and boundary points


*Find a non-convex hull of the set of points. Points on the boundary will be known as "boundary points" and rest will be "internal points". (Note: I know this is easier said than done. I used MATLAB's boundary with a shrink factor of 1.)


Part 3: Perform interpolation and compute error


*Compute the value at each internal point as a weighted sum of the values of its K neighbors, where the weights are the inverse distances in some power P. Note: using only the internal points is crucial for correctness, since IDW, which doesn't use/store direction information, cannot extrapolate correctly (as opposed to trilinear interpolation, for example).

*Subtract the interpolated value from the known value, at each point, to obtain the error.
For visualization purposes, a histogram appeared useful. For computational purposes, a mean absolute error appeared useful.
