# Why won't DiscretizeRegion work with a custom distance function?

I am trying (in 10.0) to create a region from a list of data points. I have created myDistance to use as the boolean test for whether or not a point should be included in the region

myDistance[x_?NumericQ, y_?NumericQ, z_?NumericQ] :=
Norm[{x, y, z} - First[Nearest[pts, {x, y, z}]]]


where pts is a list of (lattice) points that I want to help define the region. I can define a region based on this function:

myRegion = ImplicitRegion[myDistance[x, y, z] <= 1.01 Sqrt, {x, y, z}]


but if I try to discretize this region I get the vague error: DiscretizeRegion::drf: DiscretizeRegion was unable to discretize the region ImplicitRegion[<<2>>]. >>

Any idea what is going on?

On the other hand, RegionPlot3D[myDistance[x,y,z]<= 1.01 Sqrt,{x,-10,10},{y,-10,10},{z,-10,10}] produces a perfectly fine representation of the region.

UPDATE: The problem only seems to occur for large point sets (>20 on my machine). Smaller point sets work fine in both cases.

• Try increasing PlotPoints. This works for me on several point sets: DiscretizeGraphics@Normal@RegionPlot[myDistance[x, y, z] <= 1.01 Sqrt, , {x, -10, 10}, {y, -10, 10}, {z, -10, 10}, PlotPoints -> 30]. – Michael E2 Nov 30 '15 at 4:36
• This is working so far for me for quite large point sets (assuming you mean RegionPlot3D...). Why won't DiscretizeRegion do the same thing? – djphd Nov 30 '15 at 5:24
• This works, too, on your specific example: myRegion2 = RegionUnion @@ (Sphere[#, 1.01 Sqrt] & /@ pts); DiscretizeRegion[myRegion2], but it is done without a custom distance function. – Michael E2 Dec 1 '15 at 2:44
• Nice! I like that much better, and it's probably faster. We'll see how it goes. – djphd Dec 1 '15 at 13:48

Another method from another of my comments, which may work on earlier V10 versions:

Clear[myDistance2];
pts = RandomReal[{-8, 8}, {111, 3}];
With[{d = Nearest[pts]},
myDistance2[dist_][x_?NumericQ, y_?NumericQ, z_?NumericQ] :=
Length@d[{x, y, z}, {1, dist}]];

DiscretizeGraphics@
Normal@RegionPlot3D[
myDistance2[1.01 Sqrt][x, y, z] >= 1, {x, -10, 10}, {y, -10,
10}, {z, -10, 10}, PlotPoints -> 30] // AbsoluteTiming (I changed the distance function slightly so that it would be several times faster than the OP's. Note also that RegionPlot3D does not always produce graphics with closed boundaries, even if the surface is mathematically closed. Sometimes increasing PlotPoints or MaxRecursion can help.)

Response to comment: @djphd Nearest computes a data structure that makes searching for the nearest point very efficient. The data structure depends only on the points pts. So you can save time over many queries by constructing the data structure only once. In the first argument to With, the data structure is returned in the form of a NearestFunction returned by Nearest; With injects this function d into the definition of myDistance2. The form d[point, {n, radius}] finds the nearest n points within the radius. It will return an empty list of length zero if there are none; otherwise, the length will be at least 1 and at most n. In the definition of myDistance2, n is 1, which gives the fastest performance. In fact, using Nearest[pts -> Automatic] is a bit faster, but it does not seem to matter much in this case.

• Works for 10.0 and 10.1! Thank you! Can you discuss a little about why using With and Length make this so much faster? Essentially the algorithm is the same as my original effort, but this is significantly faster than I would have imagined. Tested so far on my system with thousands of points and no problems whatsoever (especially when reducing PlotPoints a little because the original pts are so close). – djphd Dec 2 '15 at 0:57
• @djphd To answer your question was too long for a comment, so I included it in my answer. – Michael E2 Dec 2 '15 at 2:13

Your code produces a perfectly good 3D region. To visualize it you just have run DiscretizeRegion on it.

SeedRandom;
With[{pts = RandomReal[1., {5, 3}]},
inRegion[x_?NumericQ, y_?NumericQ, z_?NumericQ] :=
Norm[{x, y, z} - First[Nearest[pts, {x, y, z}]]]]
r = ImplicitRegion[inRegion[x, y, z] <= 1.01 Sqrt, {x, y, z}];
dr = DiscretizeRegion[r] Now further computations can be done, such as

Volume[dr]


7.8467

I don't know what you were expecting to see as the result of ImplicitRegion. The documentation makes it clear that it is a wrapper function that encapsulates region information for various region savvy functions to use.

• I agree that this works, and my question has been reworded. Your code and mine both break for more than ~20 points in pts on my machine. Could it be a memory problem? – djphd Nov 30 '15 at 4:30

Using regions directly is faster and works:

pts=RandomReal[{-8,8},{111,3}];

myRegion2=RegionUnion@@(Sphere[#,1.01 Sqrt]&/@pts);
DiscretizeRegion[myRegion2]//AbsoluteTiming • Nice code, but it does not work on my system (10.0 or 10.1 on Windows 7). Instead I get {463.711, EmptyRegion} with the errors DiscretizeRegion::drcd: Available methods not able to resolve all components of dimension 2 and DiscretizeRegion::dre: DiscretizeRegion did not find any sample points in the region RegionUnion[<<111>>] with bounds {{-10.8139,10.8925},{-10.9602,10.9442},{-10.9636,10.9659}}. If this is not correct, different bounds may help. >> I have tried different bounds, different accuracy goals and different precision goals to no avail. I don't think it's a memory issue. – djphd Dec 1 '15 at 21:59
• @djphd The region functionality in 10.0 was not completely developed. It's been filled in more and more with each release, 10.1, 10.2, and (I think) 10.3. I have 10.2 currently. – Michael E2 Dec 1 '15 at 22:49