Listing over multiple inputs (tensor rank > 1) in a compiled function

Question

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."