2
$\begingroup$

I get the random discrete data as follows

data={{65, 5, 0.1}, {70, 5, 0.1}, {75, 5, 0.1}, {80, 5, 0.1}, {70, 10, 
  0.1}, {75, 10, 0.1}, {80, 10, 0.1}, {70, 15, 0.1}, {75, 15, 
  0.1}, {80, 15, 0.1}, {75, 20, 0.1}, {80, 20, 0.1}, {75, 25, 
  0.1}, {80, 25, 0.1}, {80, 30, 0.1}, {80, 35, 0.1}, {55, 5, 
  0.2}, {60, 5, 0.2}, {65, 5, 0.2}, {70, 5, 0.2}, {75, 5, 0.2}, {80, 
  5, 0.2}, {60, 10, 0.2}, {80, 10, 0.2}, {65, 15, 0.2}, {80, 15, 
  0.2}, {65, 20, 0.2}, {80, 20, 0.2}, {70, 25, 0.2}, {80, 25, 
  0.2}, {75, 30, 0.2}, {80, 30, 0.2}, {75, 35, 0.2}, {80, 35, 
  0.2}, {80, 40, 0.2}, {45, 5, 0.3}, {50, 5, 0.3}, {55, 5, 0.3}, {60, 
  5, 0.3}, {65, 5, 0.3}, {70, 5, 0.3}, {75, 5, 0.3}, {80, 5, 
  0.3}, {50, 10, 0.3}, {80, 10, 0.3}, {55, 15, 0.3}, {80, 15, 
  0.3}, {60, 20, 0.3}, {80, 20, 0.3}, {65, 25, 0.3}, {80, 25, 
  0.3}, {70, 30, 0.3}, {80, 30, 0.3}, {75, 35, 0.3}, {80, 35, 
  0.3}, {75, 40, 0.3}, {80, 40, 0.3}, {40, 5, 0.4}, {45, 5, 0.4}, {50,
   5, 0.4}, {55, 5, 0.4}, {60, 5, 0.4}, {65, 5, 0.4}, {70, 5, 
  0.4}, {75, 5, 0.4}, {80, 5, 0.4}, {45, 10, 0.4}, {80, 10, 0.4}, {50,
   15, 0.4}, {80, 15, 0.4}, {55, 20, 0.4}, {80, 20, 0.4}, {60, 25, 
  0.4}, {80, 25, 0.4}, {65, 30, 0.4}, {80, 30, 0.4}, {70, 35, 
  0.4}, {80, 35, 0.4}, {75, 40, 0.4}, {80, 40, 0.4}, {35, 5, 
  0.5}, {40, 5, 0.5}, {45, 5, 0.5}, {50, 5, 0.5}, {55, 5, 0.5}, {60, 
  5, 0.5}, {65, 5, 0.5}, {70, 5, 0.5}, {75, 5, 0.5}, {80, 5, 
  0.5}, {40, 10, 0.5}, {80, 10, 0.5}, {45, 15, 0.5}, {80, 15, 
  0.5}, {50, 20, 0.5}, {80, 20, 0.5}, {55, 25, 0.5}, {80, 25, 
  0.5}, {60, 30, 0.5}, {80, 30, 0.5}, {65, 35, 0.5}, {80, 35, 
  0.5}, {70, 40, 0.5}, {75, 40, 0.5}, {80, 40, 0.5}, {30, 5, 
  0.6}, {35, 5, 0.6}, {40, 5, 0.6}, {45, 5, 0.6}, {50, 5, 0.6}, {55, 
  5, 0.6}, {60, 5, 0.6}, {65, 5, 0.6}, {70, 5, 0.6}, {75, 5, 
  0.6}, {80, 5, 0.6}, {30, 10, 0.6}, {80, 10, 0.6}, {40, 15, 
  0.6}, {80, 15, 0.6}, {45, 20, 0.6}, {80, 20, 0.6}, {50, 25, 
  0.6}, {80, 25, 0.6}, {55, 30, 0.6}, {80, 30, 0.6}, {60, 35, 
  0.6}, {80, 35, 0.6}, {65, 40, 0.6}, {70, 40, 0.6}, {75, 40, 
  0.6}, {80, 40, 0.6}, {20, 5, 0.7}, {25, 5, 0.7}, {30, 5, 0.7}, {35, 
  5, 0.7}, {40, 5, 0.7}, {45, 5, 0.7}, {50, 5, 0.7}, {55, 5, 
  0.7}, {60, 5, 0.7}, {65, 5, 0.7}, {70, 5, 0.7}, {75, 5, 0.7}, {80, 
  5, 0.7}, {25, 10, 0.7}, {80, 10, 0.7}, {30, 15, 0.7}, {80, 15, 
  0.7}, {40, 20, 0.7}, {80, 20, 0.7}, {45, 25, 0.7}, {80, 25, 
  0.7}, {50, 30, 0.7}, {80, 30, 0.7}, {55, 35, 0.7}, {80, 35, 
  0.7}, {65, 40, 0.7}, {70, 40, 0.7}, {75, 40, 0.7}, {80, 40, 
  0.7}, {15, 5, 0.8}, {20, 5, 0.8}, {25, 5, 0.8}, {30, 5, 0.8}, {35, 
  5, 0.8}, {40, 5, 0.8}, {45, 5, 0.8}, {50, 5, 0.8}, {55, 5, 
  0.8}, {60, 5, 0.8}, {65, 5, 0.8}, {70, 5, 0.8}, {75, 5, 0.8}, {80, 
  5, 0.8}, {20, 10, 0.8}, {80, 10, 0.8}, {25, 15, 0.8}, {80, 15, 
  0.8}, {30, 20, 0.8}, {80, 20, 0.8}, {35, 25, 0.8}, {80, 25, 
  0.8}, {45, 30, 0.8}, {80, 30, 0.8}, {50, 35, 0.8}, {80, 35, 
  0.8}, {55, 40, 0.8}, {60, 40, 0.8}, {65, 40, 0.8}, {70, 40, 
  0.8}, {75, 40, 0.8}, {80, 40, 0.8}, {10, 5, 0.9}, {15, 5, 0.9}, {20,
   5, 0.9}, {25, 5, 0.9}, {30, 5, 0.9}, {35, 5, 0.9}, {40, 5, 
  0.9}, {45, 5, 0.9}, {50, 5, 0.9}, {55, 5, 0.9}, {60, 5, 0.9}, {65, 
  5, 0.9}, {70, 5, 0.9}, {75, 5, 0.9}, {80, 5, 0.9}, {10, 10, 
  0.9}, {15, 10, 0.9}, {20, 10, 0.9}, {25, 10, 0.9}, {30, 10, 
  0.9}, {35, 10, 0.9}, {40, 10, 0.9}, {45, 10, 0.9}, {50, 10, 
  0.9}, {55, 10, 0.9}, {60, 10, 0.9}, {65, 10, 0.9}, {70, 10, 
  0.9}, {75, 10, 0.9}, {80, 10, 0.9}, {15, 15, 0.9}, {20, 15, 
  0.9}, {25, 15, 0.9}, {30, 15, 0.9}, {35, 15, 0.9}, {40, 15, 
  0.9}, {45, 15, 0.9}, {50, 15, 0.9}, {55, 15, 0.9}, {60, 15, 
  0.9}, {65, 15, 0.9}, {70, 15, 0.9}, {75, 15, 0.9}, {80, 15, 
  0.9}, {20, 20, 0.9}, {25, 20, 0.9}, {30, 20, 0.9}, {35, 20, 
  0.9}, {40, 20, 0.9}, {45, 20, 0.9}, {50, 20, 0.9}, {55, 20, 
  0.9}, {60, 20, 0.9}, {65, 20, 0.9}, {70, 20, 0.9}, {75, 20, 
  0.9}, {80, 20, 0.9}, {25, 25, 0.9}, {30, 25, 0.9}, {35, 25, 
  0.9}, {40, 25, 0.9}, {45, 25, 0.9}, {50, 25, 0.9}, {55, 25, 
  0.9}, {60, 25, 0.9}, {65, 25, 0.9}, {70, 25, 0.9}, {75, 25, 
  0.9}, {80, 25, 0.9}, {35, 30, 0.9}, {40, 30, 0.9}, {45, 30, 
  0.9}, {50, 30, 0.9}, {55, 30, 0.9}, {60, 30, 0.9}, {65, 30, 
  0.9}, {70, 30, 0.9}, {75, 30, 0.9}, {80, 30, 0.9}, {40, 35, 
  0.9}, {45, 35, 0.9}, {50, 35, 0.9}, {55, 35, 0.9}, {60, 35, 
  0.9}, {65, 35, 0.9}, {70, 35, 0.9}, {75, 35, 0.9}, {80, 35, 
  0.9}, {45, 40, 0.9}, {50, 40, 0.9}, {55, 40, 0.9}, {60, 40, 
  0.9}, {65, 40, 0.9}, {70, 40, 0.9}, {75, 40, 0.9}, {80, 40, 0.9}}
  ListPointPlot[data]

I get this

Image Missing

Then,I use ListSurfacePlot[data], but the result is wrong. I want to get a 3D plot like RegionPlot3D

Edit graphic representation:

 Graphics3D[{PointSize[.01], Hue[#[[3]]], Point[# {1, 1, 10} ]} & /@ data]

enter image description here

$\endgroup$
3
  • $\begingroup$ Try ListPlot3D. $\endgroup$ Commented Apr 8, 2015 at 16:49
  • 3
    $\begingroup$ I pasted in a Graphic3D figure just so folks can see what the data looks like... $\endgroup$
    – george2079
    Commented Apr 8, 2015 at 18:30
  • $\begingroup$ Do you get what you need with ListPointPlot3D[data] or ListPlot3D[data] or ListSurfacePlot3D[data, BoxRatios -> {1, 1, 1/3}]? $\endgroup$
    – kglr
    Commented Apr 8, 2015 at 18:46

2 Answers 2

3
$\begingroup$

Motivation

The following gives some idea of what the data looks like:

Graphics3D[
  Polygon /@ (Module[
    {center, pts},
    pts = #;
    center = Mean[pts];
    SortBy[pts, (N[ArcTan @@ Most[# - center]] &)]
  ] &) /@ GatherBy[data, #[[3]] &],
  Axes -> True,
  Boxed -> True,
  BoxRatios -> {1, 1, 1},
  ViewPoint -> 1000 {1, 1, 1}
]

zlayers

For each $z$-layer, I have connected all points to created a non-self-intersecting polygon (with codes adapted from #48091).

The result makes clear that the top and bottom $z$-layers have points on both the boundary and inside. The inner $z$-layers have points on the boundary only.

Likewise for the $x$-layers:

Graphics3D[
  Polygon /@ (Module[
    {center, pts},
    pts = #;
    center = Mean[pts];
    SortBy[pts, (N[ArcTan @@ Rest[# - center]] &)]
  ] &) /@ GatherBy[data, #[[1]] &],
  Axes -> True,
  Boxed -> True,
  BoxRatios -> {1, 1, 1},
  ViewPoint -> 1000 {1, 1, 1}
]

xlayers

and $y$-layers:

Graphics3D[
  Polygon /@ (Module[
    {center, pts},
    pts = #;
    center = Mean[pts];
    SortBy[pts, (N[ArcTan @@ Drop[# - center, {2}]] &)]
  ] &) /@ GatherBy[data, #[[2]] &],
  Axes -> True,
  Boxed -> True,
  BoxRatios -> {1, 1, 1},
  ViewPoint -> 1000 {1, 1, 1}
]

ylayers

All three layers together:

beehive

What you want is a concave hull. Let me present two approaches.

Approach 1

Because, as shown by the layers above, your points come in discrete cubic cells, each cell has a convex hull constructed from the data points it contains. The concave hull can be constructed my merging these cellular convex hulls.

My procedure is as follows.

  1. Your data is separated into layers in each direction.
  2. Each consecutive layers are formed into a 3D region (DelaunayMesh) bounded by a convex hull (xstrips and so on).
  3. Cubic regions (cuboids) are generated by partitioning points (pts) in the bounding box of your data.
  4. The points in the bounding box are chosen for the construction of the convex hulls only if they are in a region in which xstrips, ystrips, and zstrips intersect (soliddata).
  5. The chosen points (soliddata) within each of the cuboids are made into convex hulls.
  6. Polygon faces from all the hulls are collected together into the final result.

Below is the code:

data = Rationalize /@ data;
xlayers = (Module[
  {center, pts},
  pts = #;
  center = Mean[pts];
  SortBy[pts, (N[ArcTan @@ Rest[# - center]] &)]
] &) /@ Select[GatherBy[data, #[[1]] &], Length@# > 2 &];
ylayers = (Module[
  {center, pts},
  pts = #;
  center = Mean[pts];
  SortBy[pts, (N[ArcTan @@ Drop[# - center, {2}]] &)]
] &) /@ Select[GatherBy[data, #[[2]] &], Length@# > 2 &];
zlayers = (Module[
  {center, pts},
  pts = #;
  center = Mean[pts];
  SortBy[pts, (N[ArcTan @@ Most[# - center]] &)]
] &) /@ Select[GatherBy[data, #[[3]] &], Length@# > 2 &];
xstrips = DelaunayMesh[Level[#, {-2}]] & /@ Partition[xlayers, 2, 1];
ystrips = DelaunayMesh[Level[#, {-2}]] & /@ Partition[ylayers, 2, 1];
zstrips = DelaunayMesh[Level[#, {-2}]] & /@ Partition[zlayers, 2, 1];
xs = Range[Min[#[[1]] & /@ data], Max[#[[1]] & /@ data], 5];
ys = Range[Min[#[[2]] & /@ data], Max[#[[2]] & /@ data], 5];
zs = Range[Min[#[[3]] & /@ data], Max[#[[3]] & /@ data], 1/10];
pts = Table[
  {x, y, z},
  {x, xs}, {y, ys}, {z, zs}
];
partitions = Partition[pts, {2, 2, 2}, 1];
cellpoints = Flatten[#, 2] & /@ Level[partitions, {-5}];
cuboids = Cuboid[First@Sort[#], Last@Sort[#]] & /@ cellpoints;
ptsflat = Level[pts, {-2}];
soliddata = Select[
  ptsflat,
  Function[
    p,
    AllTrue[
      {xstrips, ystrips, zstrips}, 
      Function[strip, AnyTrue[strip, RegionMember[#, p] &]]
    ]
  ]
];
hullpoints = Function[p, Select[soliddata, RegionMember[p, #] &]] /@ cuboids;
realhullpoints = Select[hullpoints, Length@# > 3 &];
hulls = DeleteCases[ConvexHullMesh /@ realhullpoints, _EmptyRegion];
polygons = MeshPrimitives[#, 2] & /@ hulls;
Graphics3D[
  {
    EdgeForm[{Thin, Opacity[1/30, Black]}],
    polygons
  },
  AxesLabel -> {"x", "y", "z"},
  Axes -> True,
  Boxed -> False,
  BoxRatios -> {1, 1, 1},
  ViewPoint -> 1000 {-1, 1, -1}
]

reconstruction

Approach 2

You can also try what's described here: #8753. We don't actually need to explicitly load TetGenLink in Mathematica 10 as many of the needed functions have been incorporated into the Wolfram Language.

concave hull

tetrahedra = Level[MeshPrimitives[DelaunayMesh[data], 3], {-3}];
radius[p_] := Circumsphere[p][[2]];
radii = radius /@ tetrahedra;
alphashape[rmax_] := Pick[tetrahedra, radii, r_ /; r < rmax]
faces[tetras_] := Flatten[
  tetras /. {a_, b_, c_, d_} :> {{a, b, c}, {a, b, d}, {a, c, d}, {b, c, d}},
  1
];
externalfaces[faces_] := Cases[Tally[Sort /@ faces], {face_, 1} :> face];
radiusmax = 150;
polys = externalfaces@faces@alphashape[radiusmax];
Graphics3D[
  {EdgeForm[{Thin, Opacity[1/30, Black]}], Polygon@polys},
  AxesLabel -> {"x", "y", "z"},
  Axes -> True,
  Boxed -> False,
  BoxRatios -> {1, 1, 1},
  ViewPoint -> 1000 {-1, 1, -1}
]
$\endgroup$
2
  • 1
    $\begingroup$ Why re-compute the radius of the circum-sphere when you can get it directly from Circumsphere as follows: Circumsphere[p][[2]] $\endgroup$
    – RunnyKine
    Commented Jul 9, 2015 at 5:45
  • $\begingroup$ @RunnyKine Thanks! Didn't notice that before. $\endgroup$
    – Taiki
    Commented Jul 9, 2015 at 10:25
1
$\begingroup$
Show[ConvexHullMesh[data],
 BoxRatios -> {1, 1, 1},
 Axes -> True,
 ImageSize -> 500]

enter image description here

$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.