Manipulating lists based on conditionals and 3D plotting the smooth surface that bounds a list of points

Question

First of all I want to create a list of evenly spaceed points in a cubic region in $\mathbb{R}^3$. Then I want to keep the points in this list which satisfy a certain condition, i.e. a function that takes these points in $\mathbb{R}^3$ and outputs true or false depending on the details of my problem. Then I need to plot these points using some sort of a interpolation function, which will be able to show the surface that bounds the volume that the points in my list are located in.

The selection criterion is the following:

There is a matrix $\mathcal{M}=\mathcal{M}(x,y,z)$ defined for every point in this list. Say for the point $P= (1,4,5)$, the matrix $\mathcal {M}(1,4,5)$ is Positive Semi-Definite, then I want to KEEP the point $(1,4,5)$. If at that point the matrix is nott Positive Semi-Definite, then I want to DISCARD that point.

My problem is

1. I could not create a list comprised of evenly spaced points in $\mathbb{R}^3$.
2. I could not discard or keep the elements of a given set based on the output that an element of the list gives when inputed into a function.
3. I could not create a smooth surface which bounds a list of points in $\mathbb{R}^3$.

Basically If at a point the matrix is positive semi-definite, then I want to keep that point. If at a point the matrix is not positive semi-definite, then I want to discard that point.

If needed I can also provide my function, which is the criteria for keeping or discarding certain points, but I am afraid it will only serve to clutter the discussion.

Answer 1

Here is a solution described in: Smooth convex hull of a large data set of 3D points.

mve[pts_?MatrixQ, tol_ : 1.*^-8] := 
 Module[{prec = Precision[pts], bet, bm, c, d, del, dp, h, j, kap, n, 
   qj, qm, sc, sig, u, zero}, 
  zero = SetPrecision[0, prec]; {n, d} = Dimensions[pts]; dp = d + 1;
  qm = PadRight[pts, {n, dp}, N[1, prec]];
  u = N[ConstantArray[1/n, n], prec];
  bm = Transpose[qm].DiagonalMatrix[u].qm;
  kap = Diagonal[qm.LinearSolve[bm, Transpose[qm]]];
  c = FixedPoint[(j = Ordering[kap, -1]; qj = Extract[qm, j];
       sc = 1/Extract[kap, j]; sig = (1 - sc dp)/(1 - sc);
       bet = sig/dp; h = 1 - bet; del = dp (1 - sc)/d;
       kap = del (kap - sc sig (qm.LinearSolve[bm, qj])^2);
       bm = h bm + bet Outer[Times, qj, qj];
       h # + SparseArray[{j -> bet}, n, zero]) &, u, 
     SameTest -> (SquaredEuclideanDistance[##] <= tol &)].pts;
  Ellipsoid[c, d (Drop[bm, -1, -1] - Outer[Times, c, c])]]

To test this, we first create some random points in the cube -1..1 and then plot the ellipsoide together with these points:

pts = RandomReal[{-1, 1}, {50, 3}];
Graphics3D[{{Opacity[0.2], mve[pts]}, Point[pts]}]