21
$\begingroup$

I've fallen down a rabbit hole recently, and need a reality check.

Suppose you have an arbitrary unit vector, $u \in \mathbf{R}^3$, and you cut up the unit cube into slices orthogonal to $u$ with width (measured along $u$) $\delta$. What are the volumes of the slices?

Here's what I mean:

enter image description here

(Each image shows every other slice. Code for creating the images is at the end of the question).

Seems like a simple enough question -- canonical, almost -- and here's a simple enough solution using ImplicitRegion:

cubevertices = Tuples[{0., 1.}, 3];
slicevolumeIR[v_, width_] := Module[{slicerange, slicepoints, slices},

  (* Limits of the projection of the cube onto u *)
  slicerange = MinMax[cubevertices.v]; 

  (* Where the slices occur along u *)
  slicepoints = Union[Range[Sequence @@ slicerange, width], {slicerange[[2]]}];

  (* Find the slices as ImplicitRegions *)
  slices = 
   ImplicitRegion[#1 <= {x, y, z}.v <= #2, {{x, 0, 1}, {y, 0, 1}, {z, 0, 1}}]
    & @@@ Partition[slicepoints, 2, 1];

  (* Find the volumes of the slices *)
  RegionMeasure /@ slices
 ]

Then, for the example that generated the above figure,

u = Normalize@RandomReal[{-1, 1}, 3]
δ = 0.1;
AbsoluteTiming[
 volumesIR = slicevolumeIR[u, δ];
 ]

(* {1.37398, {0.000903296, 0.00632307, 0.0171626, 0.0334219, 0.0550925, 
    0.0806039, 0.105649, 0.122993, 0.129516, 0.125199, 0.110043, 0.085887, 
    0.0601134, 0.0373954, 0.0200329, 0.00809016, 0.00156718, 7.61824*10^-6}} *)

Which is all well and good, but that took over a second and there were only 18 slices (as shown in the image... there's a tiny slice of leftovers at the top right of the plot on the right). I would like to be able to use this as part of a function to be fed into NMinimize as $u$ runs across a range, with on the order of several hundred slices. Clearly, RegionMeasure[ImplicitRegion[...]] isn't up to the task.

So the question is...

What can I do to make it fast enough?

And, while we're on the subject,

How fast can it be done?

It feels like this should be straightforward -- it's just unit cubes and unit vectors, after all. How hard could it be? On the other hand, it's not like simple things can't be complicated. Which is why I need a reality check.

I've tried a number of improvements to the slicevolumesIR approach, the most successful of which I'll describe below. But it seems rather convoulted and verbose to me. I can't help feeling that there's still a lot more that could be done to simplify and speed up the process.

  1. Are there ways to restructure the code to speed things up? (For example, I couldn't find a clean and effective way to use the Listable attribute, but that doesn't mean there isn't one).
  2. Or ways to simplify the algebra or exploit mathematical structures that I've missed?
  3. Or are there completely different, subtler ways of approaching the question? Can it be translated to some other mathematical framework that makes a solution more forthcoming?
  4. Or is it just a more computationally complex problem than it first appears, and I'm waay overthinking it?

(Compiled and/or Parallelized code is fine. I haven't done either so far because I didn't want to use Compile to cover up inefficient code and I might want to save parallelization for the outer code. But if you can make use of it here, then it's fair game.)

Cross section area approximation

The main observation here is that the volume of each slice of the cube is approximately proportional to the area of the cross section (orthogonal to $u$) down the middle of the slice (in a Riemann sum sort of way). Since the volume has to sum to 1, we can normalize these areas to get pretty good estimates of the slice volumes.

A couple of points about how I've gone about it:

  1. I've hard-coded as much of the algebra as I can (cubeedges parameterizes the edges of the cube, and tvalues determines the values of those parameters for a particular cross section). So I can read off the points where a cross section cuts the edges of the cube and don't need to Solve anything.

  2. After finding the vertices where the cross section intersects the cube edges, I project it onto a 2D subspace and find the area by direct calculation (in polyarea) rather than using RegionMeasure in 3D (which is slow).

  3. I should also note that the function should to be able to cope with degeneracies in $u$, such as the possibility of it lying along an axis, and therefore not intersecting with some of the edges. (Hence the safedivide function.)

Set up some global variables for describing the cube:

cubevertices = Tuples[{0., 1.}, 3]; 
cubeedges = Function[t, #1 + t (#2 - #1)] & @@@ 
    Select[Subsets[cubevertices, {2}], HammingDistance @@ # == 1 &];

Now a few helper functions:

safedivide[num_, denom_] := If[denom == 0, ∞, num/denom]

(* Calculate the values of t in `cubeedgeeqns` for which each edge 
  intersects the plane passing through projection points *)
tvalues[v_, projectionpoints_] :=
 safedivide @@@ {{#, v[[3]]}, {#, v[[2]]}, {#, v[[1]]}, {# - v[[3]], v[[2]]}, 
      {# - v[[3]], v[[1]]}, {# - v[[2]], v[[3]]}, {# - v[[2]], v[[1]]}, 
      {# - (v[[2]] + v[[3]]), v[[1]]}, {# - v[[1]], v[[3]]}, 
      {# - v[[1]], v[[2]]}, {# - (v[[1]] + v[[3]]), v[[2]]}, 
      {# - (v[[1]] + v[[2]]), v[[3]]}} & /@ projectionpoints

(* Calculate the area of a 2D polygon embedded in R^3 by projecting
  onto a 2D subspace orthogonal to u *)
polyarea[vertexlist_, orthobasis_] := Block[{
   len = Length@vertexlist, 
   list2D = Transpose[vertexlist.# & /@ orthobasis], 
   slist2D},

  (* Sort the vertices *)
  slist2D = SortBy[# - Mean[list2D] & /@ list2D, ArcTan @@ # &];

  (* Calculate the area of the polygon. Using Sum with Mod is faster than the 
    Det formula, and much faster than Area or RegionMeasure on Polygon*)
  0.5 Sum[With[{j = Mod[i, len] + 1}, 
    slist2D[[i, 1]] slist2D[[j, 2]] - slist2D[[i, 2]] slist2D[[j, 1]]],
   {i, len}]
 ]

Putting it together:

slicevolume2D[v_, width_] := Block[{
   vrange = MinMax[cubevertices.v], 
   orthobasis = Rest@Orthogonalize[Join[{v}, RandomReal[1, {2, 3}]]], 
   vpoints, tvals, crosssecvertices}, 

  (* Find where the planes occur along v -- at the slice midpoints *)
  vpoints = Mean /@ Partition[
     Union[Range[Sequence @@ vrange, width], {vrange[[2]]}], 2, 1];

  (* Find where the planes cut the cube edges *)
  tvals = tvalues[v, vpoints];

  (* Find the vertices of the cross sections *)
  crosssecvertices = Union[#1[#2] & @@@ 
         Pick[Transpose[{cubeedges, #}], 0 <= # <= 1 & /@ #] &@#, 
      SameTest -> (Chop[#1 - #2] == {0, 0, 0} &)] & /@ tvals;

  (* Find the areas and Normalize *)
  Normalize[polyarea[#, orthobasis] & /@ crosssecvertices, Total]
 ]

Try it out on the example from above:

AbsoluteTiming[
 volumes2D = slicevolume2D[u, δ];
 ]

(* {0.00351377, Null} *)

Which is a pretty good improvement over the ImplicitRegion version, and fairly accurate 

Max[Abs[volumes2D - volumesIR]]

(* 0.000443282 *)

(and accuracy will increase with the number of slices).

Plot the planes (code in Appendix):

enter image description here

Benchmarks

Here's a comparison with a larger number of slices (288, as it turns out).

SeedRandom[12]
u = Normalize@RandomReal[{-1, 1}, 3];
δ = 0.005;
AbsoluteTiming[
 volumesIR = slicevolumeIR[u, δ];
 ]

AbsoluteTiming[
 volumes2D = slicevolume2D[u, δ];
 ]

(* {31.9965, Null}
   {0.0572544, Null} *)

Compare the results:

Max[Abs[volumes2D - volumesIR]]
ListLinePlot[{volumes2D, volumesIR}]

(* 2.95272*10^-7 *)

enter image description here

So there is very good agreement, and slicevolume2D is over 500 times faster than slicevolumeIR. Still, I'm left wondering if it's overly complicated and there isn't a simpler/faster way to do this.

Appendix 1: Code for the figures

Plotting the slices (for slicevolumesIR) or polygons (for slicecolumes2D) is simply a case of returning the regions instead of their computed measures. In the first case:

slicesIR[v_, width_] := Module[{slicerange, slicepoints},
  slicerange = N@MinMax[cubevertices.v]; 
  slicepoints = 
   Union[Range[Sequence @@ slicerange, width], {slicerange[[2]]}];
  ImplicitRegion[
   #1 <= {x, y, z}.v <= #2, {{x, 0, 1}, {y, 0, 1}, {z, 0, 1}}
   ] & @@@ Partition[slicepoints, 2, 1]
  ]

Then plot with 

sliceregionsIR = slicesIR[u, δ];

GraphicsRow[{Show[
   RegionPlot3D[#, Mesh -> None] & /@ sliceregionsIR[[1 ;; ;; 2]], 
   Graphics3D[{Thick, Arrow[{{0, 0, 0}, u}]}]], 
  Show[RegionPlot3D[#, Mesh -> None] & /@ sliceregionsIR[[2 ;; ;; 2]], 
   Graphics3D[{Thick, Arrow[{{0, 0, 0}, u}]}]]}]

Getting the slices from slicevolume2D is a little more complicated because we have to sort the vertices in 3D -- in slicevolume2D the vertices are sorted after being projected onto a 2D subspace.

slices2D[v_, width_] := 
 Block[{vrange = MinMax[cubevertices.v], vpoints, tvals, crosssecvertices}, 
  vpoints = 
   Mean /@ Partition[
     Union[Range[Sequence @@ vrange, width], {vrange[[2]]}], 2, 1];
  tvals = tvalues[v, vpoints];
  crosssecvertices = Union[#1[#2] & @@@ 
         Pick[Transpose[{cubeedges, #}], 0 <= # <= 1 & /@ #] &@#, 
      SameTest -> (Chop[#1 - #2] == {0, 0, 0} &)] & /@ tvals;
  Polygon /@ (With[{centre = Mean[#], mp = Mean[#[[{1, 2}]]]}, 
       SortBy[#, 
        ArcTan[Cross[# - centre, mp - centre].v, 
         (mp - centre).(# - centre)] &]
       ] & /@ crosssecvertices)
  ]

Then plot with 

sliceregions2D = slices2D[u, δ];

Graphics3D[{sliceregions2D[u, δ], Thick, Arrow[{{0, 0, 0}, u}]}]

Appendix 2: Update for @ybeltukov's answer

Here's the benchmark from earlier for @ybeltukov's sliceMeasure function. I am restricting sliceMeasure to a 3-dimensional unit cube fixed with origin {0., 0., 0.} and $u$ being a unit vector. Note that it works just as well on arbitrary $n$-dimensional parallelepipeds floating around in $\mathbf{R}^n$ with $u$ being whatever you like.

SeedRandom[12]
u = Normalize@RandomReal[{-1, 1}, 3];
δ = 0.005;

origin = {0., 0., 0.};
base = {{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}};
urange = MinMax[cubevertices.u];
dist = Union[Range[Sequence @@ urange, δ], {urange[[2]]}];

AbsoluteTiming[
 volumes = sliceMeasure[origin, base, u, dist];
 ]

(* {0.000222371, Null} *)

It seems absurd to compare accuracy, since sliceMeasure doesn't make approximations in the same way that slicevolume2D does, and I have no reason to suppose that slicevolumeIR is more accurate than sliceMeasure (but Max[Abs[Differences[volumes] - volumesIR]] == 7.62888*10^-15, so they're very close).

To put those times into their proper context; sliceMeasure is about 143,888 times faster than volumesIR, whereas New Horizons (58,536 km/h) was merely 58,536 times faster than a sloth (1 km/h).

But that's not the half of it, since sliceMeasure really comes into it's own on longer lists:

timings = Table[
  u = Normalize@RandomReal[{-1, 1}, 3];
  urange = MinMax[cubevertices.u];
  dist = Union[Range[Sequence @@ urange, δ], {urange[[2]]}];
  {Length@dist,
   First@AbsoluteTiming[slicevolume2D[u, δ]],
   First@AbsoluteTiming[sliceMeasure[origin, base, u, dist]]},
  {δ, 10^-Range[1, 6, 0.5]}
  ];

ListLogLogPlot[{timings[[;; , {1, 2}]], timings[[;; , {1, 3}]]}, 
 Joined -> True, AxesLabel -> {"Number of slices", "Time"}, 
 PlotLegends -> {"slicevolume2D", "sliceMeasure"}]

enter image description here

In summary: sliceMeasure is clearly "fast enough", and I'd be very surprised if any kind of significant speed increases were possible. But, more importantly, it is also an elegant mathematical formulation of the problem, and very concisely coded.

Improve this question
$\endgroup$
2
  • $\begingroup$ Great work! I think it is already self-contained problem. However, I think little more info on the context of this question would be helpful. $\endgroup$
    – yarchik
    Nov 24, 2017 at 11:07
  • $\begingroup$ @yarchik The original context was that I had this question rattling round in the back of my head, and eventually the question of slice volumes cropped up and took on a life of its own. So, although I might, in principle, use this code somewhere else, and that's my nominal reason for wanting the speed increase, my real motivation is that I feel like cutting up a cube should be a simpler problem than I've made it. Hope that clarifies things. $\endgroup$
    – aardvark2012
    Nov 24, 2017 at 21:18

1 Answer 1

Reset to default
12
$\begingroup$

There is a general problem of slicing $n$-dimensional parallelepiped by a hyperplane. Fortunately, it can be solved using the following piece of code:

sliceMeasure[origin_, base_, vec_, dist_] := 
  With[{d = dist - vec.origin, c = base.vec, n = Length@base}, 
   Total[Array[Ramp[d - {##}.c]^n (-1)^+## &, Table[2, {n}], 0], 
     n] Abs@Det[base/c]/n!];

The main idea is to sum up volumes of $2^n$ simplexes, which origins are located at corners of the parallelepiped with a proper choice of signs. In the above code, origin is the origin of the parallelepiped, base is the list of basis vectors of the parallelepiped, vec and dist specify the halfspace x.vec <= dist. The function is listable over dist, which greatly improves the performance.

In three dimensions, we can denote base as $({\bf b}_1,{\bf b}_2,{\bf b}_3)$, vec as $\bf v$, and dist - vec.origin as $d$. In this case, the volume of the intersection of parallelepiped and halfspace has a form $$ V= V_0\sum_{i,j,k=0}^1(-1)^{i+j+k}\max(d-i {\bf b}_1\cdot{\bf v}-j{\bf b}_2\cdot{\bf v}-k{\bf b}_3\cdot{\bf v},0)^3,\\ V_0=\frac{1}{3!}\left|\det\left(\frac{{\bf b}_1}{{\bf b}_1\cdot{\bf v}},\frac{{\bf b}_2}{{\bf b}_2\cdot{\bf v}},\frac{{\bf b}_3}{{\bf b}_3\cdot{\bf v}}\right)\right|. $$

Let me demonstrate some examples:

SeedRandom[5];
dims = 2;
origin = RandomReal[1, {dims}];
base = RandomReal[1, {dims, dims}];
vec = RandomReal[1, {dims}];
dist = RandomReal[3];

RegionMeasure@RegionIntersection[Parallelepiped[origin, base], 
   HalfSpace[vec, dist]];
(* 0.144488 *)

sliceMeasure[origin, base, vec, dist]
(* 0.144488 *)

Graphics[{Red, Parallelepiped[origin, base], Blue, 
  Opacity[0.5], HalfSpace[vec, dist]}, Frame -> True]

enter image description here

SeedRandom[1];
dims = 3;
origin = RandomReal[1, {dims}];
base = RandomReal[1, {dims, dims}];
vec = RandomReal[1, {dims}];
dist = RandomReal[3];

RegionMeasure@RegionIntersection[Parallelepiped[origin, base], 
   HalfSpace[vec, dist]];
(* 0.0140496 *)

sliceMeasure[origin, base, vec, dist]
(* 0.0140496 *)

Graphics3D[{Parallelepiped[origin, base], Opacity[0.5], 
  HalfSpace[vec, dist]}]

enter image description here

Let us calculate volumes for $10^7$ different slices:

dist = 5 Rescale@N@Range[10^7];

(volumes = sliceMeasure[origin, base, vec, dist];) // AbsoluteTiming
(* {0.743143, Null} *)

Less than a second!

Here is the dependence of the volume as a function of the slice distance (we have too many points so we plot each 100-th point)

ListLinePlot[Transpose[{dist, volumes}][[;; ;; 100]], 
 AxesLabel -> {"distance", "volume"}]

enter image description here

P.S. the volume between two planes can be easily calculated as a difference for two different slices.

Improve this answer
$\endgroup$
1
  • $\begingroup$ I've added an update to the question with some comparisons for sliceMeasure. Very nicely done. $\endgroup$
    – aardvark2012
    Nov 27, 2017 at 12:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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