Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I'm trying to optimize the solution given by Michael E2 to my previous problem (Efficiently determining if 3D points are within a surface composed of polygons) which I figured would be relatively straight forward to do with Compile but I'm having trouble listing over two inputs to a compiled function.

Starting with Michael's code:

data = Flatten[
   Table[{x, y, z, x^2 + y^2 + z^2 + RandomReal[0.1]},
     {x, -2, 2, 0.2}, {y, -2, 2, 0.2}, {z, -2, 2, 0.2}], 2];

plot = ListContourPlot3D[data, Contours -> {1}, Mesh -> None];

polygonCoords = 
  Cases[Normal[plot[[1]]], Polygon[x_, ___] :> x, {0, Infinity}];

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

side[{P_, Q_, R_, ___}, X_] := Det@Differences[{X, P, Q, R}];
insideQ[polyhedron_, point_] := And @@ Positive[side[#, point] & /@ polyhedron];

points = RandomReal[{-1.5, 1.5}, {1000, 3}];
insideQ[polygonCoords, #] & /@ points; // AbsoluteTiming

Too slow: 0.5 sec for 1000 points and I have millions.

I inlined the formula for the det and replaced And + Positive with Min + UnitStep:

sideC = Compile[{{poly, _Real, 2}, {X, _Real, 1}}, 
   Block[{p1 = poly[[1, 1]], p2 = poly[[1, 2]], p3 = poly[[1, 3]], 
     q1 = poly[[2, 1]], q2 = poly[[2, 2]], q3 = poly[[2, 3]], 
     r1 = poly[[3, 1]], r2 = poly[[3, 2]], r3 = poly[[3, 3]], 
     x1 = X[[1]], x2 = X[[2]], x3 = X[[3]]}, 
    UnitStep[-p3 q2 r1 + p2 q3 r1 + p3 q1 r2 - p1 q3 r2 - p2 q1 r3 + 
      p1 q2 r3 + p3 q2 x1 - p2 q3 x1 - p3 r2 x1 + q3 r2 x1 + 
      p2 r3 x1 - q2 r3 x1 - p3 q1 x2 + p1 q3 x2 + p3 r1 x2 - 
      q3 r1 x2 - p1 r3 x2 + q1 r3 x2 + p2 q1 x3 - p1 q2 x3 - 
      p2 r1 x3 + q2 r1 x3 + p1 r2 x3 - q1 r2 x3]
    ], CompilationTarget -> "C", RuntimeAttributes -> {Listable}, 
   Parallelization -> True, RuntimeOptions -> "Speed"];

insideQC[polyhedron_, point_] := Min@sideC[polygonCoords[[All, ;; 3]], point]

This lists over all the polygons. Lets test it,

insideQC[polygonCoords, #] & /@ points; // AbsoluteTiming

Now this is over 10x faster, 0.04 seconds. Compiling this again, because of the CompiledFunctionCall provides only a tiny speed improvment.

I think it should be faster to list over both the list of polygon coordinates and the points. I suspect this question is related: Threading a compiled function over multiple arguments of different lengths, specifically ruebenko's last answer which works fine on MMA 9 on Win7 x64, It should be possible?

Why doesn't this work:

sideC[polygonCoords, points]

Shouldn't this list over both inputs? How can I make this work?

I tried to force inline-ing of the compiled function definition but it doesn't appear to be supported?

insideQC = 
  Compile[{{polyhedron, _Real, 3}, {points, _Real, 2}}, 
   Min@sideC[polyhedron, #] & /@ points, CompilationTarget -> "C", 
   RuntimeOptions -> "Speed", 
   CompilationOptions -> {"InlineCompiledFunctions" -> True, 
     "InlineExternalDefinitions" -> True}, 
   RuntimeAttributes -> {Listable}, Parallelization -> True];

Compile::cfinll: "The CompiledFunction could not be inlined because its use requires threading with the Listable runtime attribute."

share|improve this question
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.