ListPlot3D of data on "regular" lattice contained in an "irregular" region

Question

I have the following problem when trying to plot some 3D data in Mathematica. The data are initially computed on a regular lattice, that is I compute

data = Flatten[Table[{x,y,f[x,y]},{x,x0,x1,dx},{y,y0,y1,dy}],1]

The problem is that f takes real values only in a non-convex subset U on the plane. U is actually quite nasty: a "triangular" region where all corners are cuspy (imagine the region in between three identical circles with any two of them tangent to each other). When I try to plot data with ListPlot3D, the latter plots a surface over the convex hull of U.

I was wondering if there is a way to tell Mathematica to restrict itself only inside U. I was thinking that since my points already lie on a "regular" lattice this should not be too difficult but I have not found yet a solution.

I am aware of the RegionFunction option but it is very expensive to compute in my case so I am trying to figure out a way that uses only the already computed points in data.

Answer 1

This is probably not the optimal solution, so I will keep the question open in case someone has a better idea.

Here is how I solved the problem I described in my question. First, I replaced points {x,y,f[x,y]} in data for which f[x,y] was complex by {x,y,None}. Then the following function creates a 3D surface out of my data points. Note here that data is the result of

data = Table[{x,y,f[x,y]},{x,x0,x1,dx},{y,y0,y1,dy}]

that is, no flattening (which works better for the following function). The function is:

makeSurface[data_] := Module[{len1, len2, polys, a, b, c, d, p},
  len1 = Length[data];
  len2 = Length[data[[1]]];
  polys = Table[
    a = data[[i, j]];
    b = data[[i + 1, j]];
    c = data[[i + 1, j + 1]];
    d = data[[i, j + 1]];
    p = Select[{a, b, c, d}, #[[3]] =!= None &];
    If[Length[p] >= 2, Polygon[p], None],
    {i, 1, len1 - 1}, {j, 1, len2 - 1}];
  Graphics3D[Join[{EdgeForm[None],FaceForm[Directive[GrayLevel[0.5]]]},
    Select[Flatten[polys, 1], # =!= None &]]]]

The above code can be probably optimized but it worked well enough for me.

Answer 2

To generate a really nice picture, without using zillions of points, you might want to use a non-uniform distribution of points that includes the boundary of the region you want. Here's an example somewhat like what you describe. We start with three mutually tangent circles.

circPic = Graphics[{Circle[{0, Sqrt[3]}, 1],
  Circle[{-1, 0}, 1], Circle[{1, 0}, 1]}]

We write a Boolean function that determines whether a point in the rectangle {-1/2,1/2} by {0,Sqrt[3]/2} lies outside of all the circles and use this to generate some points in the region of interest.

inRegionQ[p:{x_, y_}] := Norm[p - {1, 0}] > 1 &&
  Norm[p + {1, 0}] > 1 && Norm[p - {0, Sqrt[3]}] > 1;
rectPoints = N[Flatten[Table[{x, y},
  {x, -1/2, 1/2, 0.02}, {y, 0.05, Sqrt[3]/2, 0.02}], 1]];
regionPoints = Select[rectPoints, inRegionQ];

Now we generate the boundary. The parameter n determines how many points we place on the boundary.

n = 120;
boundary = N[Join[
  Table[{1 - Cos[t], Sin[t]}, {t, Pi/n, Pi/3, Pi/n}],
  Table[{Cos[t], Sqrt[3] - Sin[t]}, {t, Pi/3 + Pi/n, 2 Pi/3, Pi/n}],
  Table[{Cos[t] - 1, Sin[t]}, {t, Pi/3 - Pi/n, 0, -Pi/n}]]];
points = Join[boundary, regionPoints];

Let's take a look.

Show[circPic, Graphics[Point[points]],
  PlotRange -> {{-3/4, 3/4}, {-0.3, 1.3}}]

Now, we define a function and use ListPlot3D to try to plot it.

f[x_, y_] := -(1 - Norm[{x - 1, y}]) (1 - Norm[{x + 1, y}])*
  (1 - Norm[{x, y - Sqrt[3]}]);
points3D = {#[[1]], #[[2]], f[#[[1]], #[[2]]]} & /@ points;
pic = ListPlot3D[points3D, Mesh -> All]

Somehow, we've got to delete that stuff that lies outside the region. In this particular example, we can use the fact that the function is zero on the boundary.

DeleteCases[Normal[pic], Polygon[{
  {x1_, y1_, z1_?(Abs[#] < 1/10.0^6 &)},
  {x2_, y2_, z2_?(Abs[#] < 1/10.0^6 &)},
  {x3_, y3_, z3_?(Abs[#] < 1/10.0^6 &)}}, ___],
  Infinity]

Pretty good, but there are a couple of problems near the cusps and it's definitely not very general since it used a specific property of the function. If you examine the structure of pic, you'll find that it contains a GraphicsComplex and the first n points in the first argument of that GraphicsComplex is exactly the boundary. Here's proof:

Most /@ pic[[1, 1, 1 ;; n]] == boundary

Now the boundary comes in three components and we want to delete any triangle that is formed by points chosen from just one of those components. The following code does this. Note that boundaryComponents contains the indices of those points that form the boundary and someSubsetQ[A,Bs] returns true if A is a subset of any one of the Bs. We want to delete those triangle indices in the multi-Polygon that are subsets of one of the boundaryComponents. That's accomplished in the following code by the DeleteCases command.

Oh, and let's add some decoration, too.

subsetQ[A_, B_] := Complement[A, B] == {};
someSubsetQ[A_, Bs_] := Or @@ Map[subsetQ[A, #] &, Bs];
boundaryComponents = Partition[Prepend[Range[n], n], 1 + n/3, n/3];
Show[pic /. Polygon[triangles_] :> {EdgeForm[Opacity[0.3]],
  Polygon[DeleteCases[triangles, _?(someSubsetQ[#, boundaryComponents] &)]]},
Graphics3D[{Thick, Line[Table[Append[pt, 0],
  {pt, Prepend[boundary, Last[boundary]]}]]}]]