0

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.

ghost
  • 381
  • 2
  • 10
  • 3
    "I could not create a list comprised of evenly spaced points in $\mathbb R^3$" - look up CoordinateBoundsArray[]. "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." - look up Select[] and Cases[]. – J. M.'s missing motivation Feb 27 '21 at 21:50
  • @J.M. What about the 3rd problem. Do you also have a suggestion regarding that? – ghost Feb 27 '21 at 21:51
  • 3
    Edit your question to include your selection criterion, and other people might be able to give concrete suggestions. – J. M.'s missing motivation Feb 27 '21 at 21:54
  • 2
    Considering question 3. A convex hull is a surface bounding your points, but it is not smooth. You would have to define more specifically what the conditions for smooth are. If an ellipsoid will do, then look e.g. at: https://mathematica.stackexchange.com/questions/57838/smooth-convex-hull-of-a-large-data-set-of-3d-points – Daniel Huber Feb 28 '21 at 10:17
  • @DanielHuber Yes an ellipsoid will do. Let me check that. – ghost Feb 28 '21 at 10:28
  • 1
    I did not commit to answering number 3 because you did not really indicate a concrete selection criterion, and things will be different depending on whether the region implied by your selected points is convex or not. If you are sure of the convexity, you can use what Daniel linked to. – J. M.'s missing motivation Feb 28 '21 at 10:31
  • @J.M. Yes, it is convex. – ghost Feb 28 '21 at 10:40

1 Answers1

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]}]

enter image description here

Daniel Huber
  • 51,463
  • 1
  • 23
  • 57