6

RipleyRassonRegion and ConvexHull provide the smallest convex region enclosing all points in the set. For regions with real concavities, the volume of this region overestimates the actual volume. For example:

pts1 = RandomPointConfiguration[PoissonPointProcess[.6, 3], Ball[{0, 0, 0}, 5]];
pts2 = RandomPointConfiguration[PoissonPointProcess[.6, 3], Ball[{0, 10, 0}, 5]];
pts = Union[pts1["Points"], pts2["Points"]];
pltpts = Table[Point[pts[[i]]], {i, 1, Length[pts]}];
plt1 = Graphics3D[pltpts];
plt2 = RipleyRassonRegion[pts];
Show[plt1]
Show[Graphics3D[{Opacity[0.2], plt2}], plt1]
Volume[RipleyRassonRegion[pts]]
2 N[Volume[Ball[{0, 10, 0}, 5]]]

The regions I am analysing are brain nuclei and the points are neurons. I want the volume of the space filled by the points, shrunk down over the concave parts, like "boundary" function in MatLab: https://www.mathworks.com/help/matlab/ref/boundary.html

neurosci
  • 73
  • 4

2 Answers2

5

A crude method is to determine whether a given point is inside or outside based on how close it is to the nearest point in the given cloud. I'll illustrate with your example.

pts1 = RandomPointConfiguration[PoissonPointProcess[.6, 3], 
   Ball[{0, 0, 0}, 5]];
pts2 = RandomPointConfiguration[PoissonPointProcess[.6, 3], 
   Ball[{0, 10, 0}, 5]];
pts = Union[pts1["Points"], pts2["Points"]];
pltpts = Table[Point[pts[[i]]], {i, 1, Length[pts]}];
plt1 = Graphics3D[pltpts];

Define our "inside-outside" predicate based on a given epsilon.

nf = Nearest[pts -> "Distance"];
inside[pt : {_Real ..}, eps_] := nf[pt][[1]] <= eps

I'll choose epsilon as the mean of the all-closest-neighbors for the point cloud.

distances = nf[pts, 2][[All, 2]];
eps = Mean[distances]

(* Out[298]= 0.925333 *)

Now we can numerically integrate over the enclosing rectangle (if we did not know it, we could find it using MinMax on all point coordinates).

bnd = 5;
NIntegrate[
 Boole[inside[{x, y, z}, eps]], {x, -bnd, bnd}, {y, -bnd, 
  3*bnd}, {z, -bnd, bnd}]

(* During evaluation of In[302]:= NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.


During evaluation of In[302]:= NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 1007.5785551881089and 56.87577978459358 for the integral and error estimates. *)


(* Out[303]= 1007.58 *)

I will remark that his method is quite sensitive to the choice of the epsilon, and it will work poorly if the cloud is not in some sense dense. That it gave a viable result in this instance is mostly coincidence-- it really requires a thicker cloud.

Daniel Lichtblau
  • 58,970
  • 2
  • 101
  • 199
  • 1
    (+1) for a very nice solution and a thorough explanation. Regarding your last comment ...is not in some sense dense is there a way that you can determine -roughly- how dense it has to be to give reasonable results. I am asking just out of curiosity. – bmf Apr 14 '22 at 04:10
  • 1
    @DanielLichtblau I tried to get the points "inside" with your nice approach and get pts == Select[pts, (inside[#, 500]) &] (*True*) , no points outside? What's wrong here? Thanks! – Ulrich Neumann Apr 14 '22 at 06:25
  • 1
    @UlrichNeumann Recall I picked epsilon related to the average distance between closest points. If you use an epsilon larger than the max such distance then all of pts will be deemed on the inside. In this randomized example family the max seems to be typically in the range of 1 to 2. So an epsilon of 500 will indeed claim they are all inside points. – Daniel Lichtblau Apr 14 '22 at 16:00
  • 1
    @bmf I'm not certain but here is a vague idea. I played first with using spherical shells with points uniformly generated in the double-unit cube (bounds from -1 to 1 that is), then removing all that were not distance between 7/10 and 1 from the origin. When I started with 10^6 points this gave around 344000 points in the shell region (makes sense if you work out the volumes). At that density I could get consistently reasonable volume estimates. At a factor of 10 less they became notably less reliable though not horrible... – Daniel Lichtblau Apr 14 '22 at 16:05
  • @DanielLichtblau many thanks for your time and the explanation :) – bmf Apr 14 '22 at 16:08
  • 1
    ...The point is that the epsilon gives a trade between how much gets added as an outer layer, vs how much gets lost on the interior due to local paucity of points. So we want there to be sufficiently many points that nearest neighbor distances are small. That way we do gain nor too much if we make epsilon a bit on the large side, and at the same time we do not lose much volume on the actual interior. That at least is my view of what makes for a "good" situation for this method. – Daniel Lichtblau Apr 14 '22 at 16:09
  • @bmf You're quite welcome. It was a good exercise for me to try to sort out what I require for reasonable density, and why. – Daniel Lichtblau Apr 14 '22 at 16:11
  • @DanielLichtblau Thanks for your answer. I got your idea, but still pts == Select[pts, (inside[#, eps/100]) &] returns True! – Ulrich Neumann Apr 15 '22 at 08:26
  • @UlrichNeumann Yeah, that it will do though I'm not sure whether I should have coded it that way. Each member of pts is distance 0 from itself so the selection predicate is True. But I should maybe have "outliers" deemed outside, in which case the inside function should perhaps have compared against distance to the second nearest point rather than the first. – Daniel Lichtblau Apr 15 '22 at 15:18
5

Maybe another possible approach.

SeedRandom[2];
pts1 = RandomPointConfiguration[PoissonPointProcess[.6, 3], 
   Ball[{0, 0, 0}, 5]];
pts2 = RandomPointConfiguration[PoissonPointProcess[.6, 3], 
   Ball[{0, 10, 0}, 5]];
pts = Union[pts1["Points"], pts2["Points"]];

regs = GradientFittedMesh[pts, PerformanceGoal -> "Quality"] // 
   ConnectedMeshComponents;
reg = regs[[1]];
RepairMesh[reg] // Show // BoundaryDiscretizeGraphics // Volume
Graphics3D[{{Opacity[.2], EdgeForm[], reg}, Point[pts]}, 
 Boxed -> False, ViewPoint -> {3.09, -0.38, 1.31}]

652.469

enter image description here

cvgmt
  • 72,231
  • 4
  • 75
  • 133