Determine whether points lie within a cow

Question

I would like to determine whether randomly-generated points

lX = RandomVariate[UniformDistribution[{-0.4, 0.4}], {1000}];
lY = RandomVariate[UniformDistribution[{-0.2, 0.2}], {1000}];
lZ = RandomVariate[UniformDistribution[{-0.2, 0.2}], {1000}];
lPoints = Thread[{lX,lY,lZ}];

lie within a cow:

g = Graphics3D[{Opacity[.09], EdgeForm[Opacity[.1]], Polygon[#, 
        VertexColors -> Table[Hue[RandomReal[]], {Length[#]}]] & /@ 
        Cases[Normal[ExampleData[{"Geometry3D", "Cow"}]], 
        Polygon[x_, ___] :> x, {0, Infinity}]}, Lighting -> "Neutral", 
        ImageSize -> 400, Axes -> True];

Show[g, ListPointPlot3D[lPoints]]

which produces something like:

I have followed the advice given here. However, the two approaches that I tried didn't seem to work. The first does the following: it finds the coordinates of the three consecutive vertices, and calculates the determinant with the point. If all the signs are positive then the point lies within the polyhedron.

side[{P_, Q_, R_, ___}, X_] := Det@Differences[{X, P, Q, R}];
insideQ[polyhedron_, point_] := And @@ Positive[side[#, point] & /@ polyhedron];
pg = Graphics3D[Point[lPoints, 
         VertexColors -> (insideQ[polygonCoords, #] & /@ 
         lPoints /. {True -> Red, 
         False -> Directive[Opacity[0.5], Blue]})], Axes -> True];

However, the above method fails to detect whether points are in the cow. (I suspect that this is because the cow is non-convex?)

Another method on the page indicated that I could use the Mathematica function:

reg = BoundaryDiscretizeGraphics[g]

However, this failed because ``The boundary curves self-intersect or cross each other in BoundaryMeshRegion''. This made me think that any of the other methods on the page which (as far as I am aware) make this assumption, won't work.

Does anyone know how I can do this? For clarity I just want to determine whether the points lie within the boundaries of the cow, then graph them. Ideally I would be able to do this for a few thousand points.

Update: this seems to suggest that I should be able to use ``ray-casting'' to determine if a point lies within the cow. If I generate a large number of rays in random directions from each point then if all of them intersect the boundary an odd number of times. If instead a point is on or outside the boundary then it will intersect an even number of times. However am not sure whether this simplifies things, since I don't know how to tell if a line intersect's the cow's boundaries!

Best,

Ben

Answer 1

This is basically a rehash of code I posted in a prior thread on this topic. The underlying method is to shoot a ray from the point and see how many surface triangles it intersects.

elsie = ExampleData[{"Geometry3D", "Cow"}];

verts = 
 First[Cases[elsie, GraphicsComplex[a_, ___] :> a, Infinity]]; pgons =
  First[Cases[elsie, Polygon[x_, ___] :> x, Infinity]];

makeTriangles[tri : {aa_, bb_, cc_}] := {tri}
makeTriangles[{aa_, bb_, cc_, dd__}] := 
 Join[{{aa, bb, cc}}, makeTriangles[{aa, cc, dd}]]

trianglevnums = Flatten[Map[makeTriangles, pgons], 1];
triangles = trianglevnums /. j_Integer :> verts[[j]];

flats = Map[Most, triangles, {2}];
pts = verts;
xcoords = pts[[All, 1]];
ycoords = pts[[All, 2]];
zcoords = pts[[All, 3]];
xmin = Min[xcoords];
ymin = Min[ycoords];
xmax = Max[xcoords];
ymax = Max[ycoords];
zmin = Min[zcoords];
zmax = Max[zcoords];

n = 30;
xspan = xmax - xmin;
yspan = ymax - ymin;
dx = 1.05*xspan/n;
dy = 1.05*yspan/n;
midx = (xmax + xmin)/2;
midy = (ymax + ymin)/2;
xlo = midx - 1.05*xspan/2;
ylo = midy - 1.05*yspan/2;

edges[{a_, b_, c_}] := {{a, b}, {b, c}, {c, a}}

vertexBox[{x1_, y1_}, {xb_, yb_, dx_, dy_}] := {Ceiling[(x1 - xb)/dx],
   Ceiling[(y1 - yb)/dy]}

segmentBoxes[{{x1_, y1_}, {x2_, y2_}}, {xb_, yb_, dx_, dy_}] := 
 Module[{xmin, xmax, ymin, ymax, xlo, xhi, ylo, yhi, xtable, ytable, 
   xval, yval, index}, xmin = Min[x1, x2];
  xmax = Max[x1, x2];
  ymin = Min[y1, y2];
  ymax = Max[y1, y2];
  xlo = Ceiling[(xmin - xb)/dx];
  ylo = Ceiling[(ymin - yb)/dy];
  xhi = Ceiling[(xmax - xb)/dx];
  yhi = Ceiling[(ymax - yb)/dy];
  xtable = Flatten[Table[xval = xb + j*dx;
     yval = (((-x2)*y1 + xval*y1 + x1*y2 - xval*y2))/(x1 - x2);
     index = Ceiling[(yval - yb)/dy];
     {{j, index}, {j + 1, index}}, {j, xlo, xhi - 1}], 1];
  ytable = Flatten[Table[yval = yb + j*dy;
     xval = (((-y2)*x1 + yval*x1 + y1*x2 - yval*x2))/(y1 - y2);
     index = Ceiling[(xval - xb)/dx];
     {{index, j}, {index, j + 1}}, {j, ylo, yhi - 1}], 1];
  Union[Join[xtable, ytable]]]

pointInsideTriangle[
  p : {x_, y_}, {{x1_, y1_}, {x2_, y2_}, {x3_, y3_}}] := 
 With[{l1 = -((x1*y - x3*y - x*y1 + x3*y1 + x*y3 - x1*y3)/(x2*y1 - 
         x3*y1 - x1*y2 + x3*y2 + x1*y3 - x2*y3)), 
   l2 = -(((-x1)*y + x2*y + x*y1 - x2*y1 - x*y2 + x1*y2)/(x2*y1 - 
         x3*y1 - x1*y2 + x3*y2 + x1*y3 - x2*y3))}, 
  Min[x1, x2, x3] <= x <= Max[x1, x2, x3] && 
   Min[y1, y2, y3] <= y <= Max[y1, y2, y3] && 0 <= l1 <= 1 && 
   0 <= l2 <= 1 && l1 + l2 <= 1]

faceBoxes[
  t : {{x1_, y1_}, {x2_, y2_}, {x3_, y3_}}, {xb_, yb_, dx_, dy_}] := 
 Catch[Module[{xmin, xmax, ymin, ymax, xlo, xhi, ylo, yhi, xval, yval,
     res}, xmin = Min[x1, x2, x3];
   xmax = Max[x1, x2, x3];
   ymin = Min[y1, y2, y3];
   ymax = Max[y1, y2, y3];
   If[xmax - xmin < dx || ymax - ymin < dy, Throw[{}]];
   xlo = Ceiling[(xmin - xb)/dx];
   ylo = Ceiling[(ymin - yb)/dy];
   xhi = Ceiling[(xmax - xb)/dx];
   yhi = Ceiling[(ymax - yb)/dy];
   res = Table[xval = xb + j*dx;
     yval = yb + k*dy;
     If[pointInsideTriangle[{xval, yval}, 
       t], {{j, k}, {j + 1, k}, {j, k + 1}, {j + 1, k + 1}}, {}], {j, 
      xlo, xhi - 1}, {k, ylo, yhi - 1}];
   res = res /. {} :> Sequence[];
   Flatten[res, 2]]]

gridBoxes[pts : {a_, b_, c_}, {xb_, yb_, dx_, dy_}] := 
 Union[Join[Map[vertexBox[#, {xb, yb, dx, dy}] &, pts], 
   Flatten[Map[segmentBoxes[#, {xb, yb, dx, dy}] &, edges[pts]], 1], 
   faceBoxes[pts, {xb, yb, dx, dy}]]]

Here we do the preprocessing for our ray-shooting.

Timing[
 gb = DeleteCases[
   Map[gridBoxes[#, {xlo, ylo, dx, dy}] &, 
    flats], {a_, b_} /; (a > n || b > n), 2];
 grid = ConstantArray[{}, {n, n}];
 Do[Map[AppendTo[grid[[Sequence @@ #]], j] &, gb[[j]]], {j, 
   Length[gb]}];]

(* Out[163]= {3.875, Null} *)

planeTriangleParams[
  p : {x_, y_}, {p1 : {x1_, y1_}, p2 : {x2_, y2_}, p3 : {x3_, y3_}}] :=
  With[{den = 
    x2*y1 - x3*y1 - x1*y2 + x3*y2 + x1*y3 - 
     x2*y3}, {-((x1*y - x3*y - x*y1 + x3*y1 + x*y3 - x1*y3)/
      den), -(((-x1)*y + x2*y + x*y1 - x2*y1 - x*y2 + x1*y2)/den)}]

getTriangles[p : {x_, y_}] := 
 Module[{ix, iy, triangs, params, res}, {ix, iy} = 
   vertexBox[p, {xlo, ylo, dx, dy}];
  triangs = grid[[ix, iy]];
  params = Map[planeTriangleParams[p, flats[[#]]] &, triangs];
  res = Thread[{triangs, params}];
  Select[res, 
   0 <= #[[2, 1]] <= 1 && 
     0 <= #[[2, 2]] <= 1 && #[[2, 1]] + #[[2, 2]] <= 1.0000001 &]]

countAbove[p : {x_, y_, z_}] := 
 Module[{triangs = getTriangles[Most[p]], threeDtriangs, lambdas, 
   zcoords, zvals}, threeDtriangs = triangles[[triangs[[All, 1]]]];
  lambdas = triangs[[All, 2]];
  zcoords = threeDtriangs[[All, All, 3]];
  zvals = 
   Table[zcoords[[j, 1]] + 
     lambdas[[j, 1]]*(zcoords[[j, 2]] - zcoords[[j, 1]]) + 
     lambdas[[j, 2]]*(zcoords[[j, 3]] - zcoords[[j, 1]]), {j, 
     Length[zcoords]}];
  If[OddQ[Length[triangs]] && OddQ[Length[Select[zvals, z > # &]]], 
   Print[{p, triangs, Length[Select[zvals, z > # &]]}]];
  Length[Select[zvals, z > # &]]]

isInside[{x_, y_, 
    z_}] /; ! ((xmin <= x <= xmax) && (ymin <= y <= ymax) && (zmin <= 
       z <= zmax)) := False
isInside[p : {x_, y_, z_}] := OddQ[countAbove[p]]

We'll illustrate with 10K points, chosen so that they are all in roughly the bounding box of our cow. Could perhaps be made faster with Compile but I'll leave that order of magnitude for another day.

points = Map[{.5, .15, .25}*# &, RandomReal[{-1, 1}, {10000, 3}]];

Timing[
 pg = Graphics3D[
    Point[points, 
     VertexColors -> (isInside /@ points /. {True -> Red, 
         False -> Directive[Opacity[0.15], Blue]})], Axes -> True, 
    AxesLabel -> {"x", "y", "z"}];]

(* Out[172]= {9.39063, Null} *)

Here's how it comes out.

g = Graphics3D[{Opacity[.09], EdgeForm[Opacity[.1]], 
    Polygon[#, 
       VertexColors -> Table[Hue[RandomReal[]], {Length[#]}]] & /@ 
     Cases[Normal[elsie], Polygon[x_, ___] :> x, {0, Infinity}]}, 
   Lighting -> "Neutral", ImageSize -> 400, Axes -> True];

Show[g, pg, AxesLabel -> {"x", "y", "z"}]

Answer 2

One way is to compute the solid angle subtended by the cow viewed at the point by summing signed solid angles corresponding to the cow's polygonal faces. If the total is 4 pi, the point is inside the cow; if 0, outside.

Background

Quoting Wikipedia, "Solid angle is the two-dimensional angle in three-dimensional space that an object subtends at a point."

As a 2D analogy, think of a point and a line segment in a 2D plane. Center a unit-radius circle on the point. Project the line segment onto the circle by connecting the endpoints to the center point. The arc projected on the unit-radius circle has a certain length, the angle subtended by the line segment. It is measured in radians.

Now move up to 3D. Think of a point and a polygon in 3D. Center a unit-radius sphere on the point, and use a central projection to map the polygon onto the sphere. The area of the projected spherical polygon on the unit-radius sphere is the solid angle subtended by the polygon. The measure is steradians.

These angles are signed. In 2D, if the circular arc evolves counter-clockwise as the line is traversed from beginning to end, the angle is positive; if clockwise, negative. In 3D, if the polygon's normal vector is directed away from the point, the solid angle is positive; if toward, negative. (The normal vector is orthogonal to the polygon plane and by convention points away from the polyhedron/cow.)

Now think about composite objects. In 2D once again, suppose we have a point and a polygon whose boundary edges are a number of straight line segments. We sum the signed angles subtended at the point by all edges. If the point is within the polygon, the arcs entirely encircle the point and the angles sum to $2\pi$ radians. If the point is outside, the angles sum to 0 radians. This is because edges on the near side of the polygon have positive angles, while edges on the far side have negative angles. They entirely cancel. These results are true even if the polygon is non-convex -- appearing to double back on itself viewed from the point.

And in 3D, suppose we have a point and a polyhedron with polygonal faces. We sum the signed solid angles subtended at the point by all faces. If the point is within the polyhedron, the entire unit-sphere is covered and the solid angles sum to $4\pi$ (the area of a unit-radius sphere) steradians. If the point is outside, the sum is 0 steradians because solid angles of polygons facing away from the point are positive and those facing toward are negative. As in 2D, this is true even if the polyhedron is non-convex.

One reference for the 3D case is Spiegel, "Schaum's Outline Series: Theory and Problems of Vector Analysis", McGraw-Hill, Ch 6 Sec 27.

Code

So the problem boils down to summing the signed solid angle of the cow's polygonal faces subtended at the point. The total is $4\pi$ if the point is inside or 0 if outside.

Mathematica represents the cow using only triangular faces. Simple formulas for the signed solid angle subtended by a planar triangle appear in the following two papers:

Van Oosterom and Strackee, "The Solid Angle of a Plane Triangle", IEEE Transactions on Biomedical Engineering 30(2) 125-126.

Eriksson, "On the Measure of Solid Angles", Mathematics Magazine 63(3) 184-187.

This is adapted from the OP's question:

lX = RandomVariate[UniformDistribution[{-0.4, 0.4}], {1000}];
lY = RandomVariate[UniformDistribution[{-0.2, 0.2}], {1000}];
lZ = RandomVariate[UniformDistribution[{-0.2, 0.2}], {1000}];
lPoints = Thread[{lX, lY, lZ}];

cow = ExampleData[{"Geometry3D", "Cow"}];

polygonCoordList = Cases[Normal[cow], Polygon[x_, ___] :> x, {0, Infinity}];

I assume (but have not checked) that Mathematica winds faces consistently, meaning a normal vector computed from the three vertices points outward from the polyhedron. The Det[] function used below relies on this. If a face is wound incorrectly, its solid angle will have the wrong sign.

The Eriksson paper's formula for the solid angle of a single triangular face:

Clear[triangleSolidAngle];

triangleSolidAngle[pnt_, {vertex1_, vertex2_, vertex3_}] := Module[
{dir1, dir2, dir3},

(* unit-length direction vectors from pnt to the three vertices *)

dir1 = Normalize[vertex1 - pnt];
dir2 = Normalize[vertex2 - pnt];
dir3 = Normalize[vertex3 - pnt];

(* solid angle *)

2*ArcTan[1 + dir1 . dir2 + dir2 . dir3 + dir3 . dir1, Det[{dir1, dir2, dir3}]]
]

The solid angle of an entire triangular-faced polyhedron:

Clear[polyhedronSolidAngle];

polyhedronSolidAngle[pnt_, triangleList_] := Total[triangleSolidAngle[pnt, #] & /@ triangleList]

There can be numerical problems. We sum a great number (5804 faces for the cow) of small quantities. Results might not be exactly $4\pi$ or exactly 0.

Execution time might be another problem. On a 2008 MacBook Pro running Mathematica 8.0.4.0, evaluating one point takes about 0.35 seconds. So evaluating the OP's 1000 points would take about 6 minutes.

Trial run of only the first 100 points:

threshold = 0.05; (* empirical cutoff *)
pointInsidePolygonList = Select[Take[lPoints, 100], polyhedronSolidAngle[#, polygonCoordList] > threshold &];

Show[g, ListPointPlot3D[pointInsidePolygonList]]