Let's imagine our central unit voxel $V_1=[0,1]^3$.

We allow the second one to be in $V_2=[-1,2]^3 \setminus [0,1]^3~$, the surrounding voxel (center removed).

My question:

What is the mean distance between a point in $V_1$ and a point in $V_2$?

Of course, the direct approach would be to write an integral:

$$ \frac{1}{1^3}\frac{1}{3^3-1^3} \times $$ $$\iiint_{(x_1,y_1,z_1)\in V_1} \iiint_{(x_2, y_2,z_2) \in V_2} \sqrt{(x_1-x_2)^2+(y_1-y_2)^2+(z_1-z_2)^2}\mathrm{d}x_{1,2} \mathrm{d}y_{1,2}\mathrm{d}z_{1,2} $$

but the one-variable function does not seem to admit a nice integral.

Has anyone gone through this calculation and would you be able to suggest any helpful manipulation before I dive into a numeric simulation?


  • $\begingroup$ Small update : succeeding a numeric simulation, I got $\approx 1.8288$. If this rings a bell to you or if it roughly seems odd, please tell me so that I can correct. I expected a smaller number but the distribution's boundaries I obtained seemed consistent with the problem. $\endgroup$
    – Wyllich
    May 24, 2018 at 15:21
  • $\begingroup$ An analytic solution is indeed available (using known results), but before I post an answer I need to make sure: if I understand your setup correctly, the volume of $V_2$ is $3^3 - 1$ so shouldn't the normalization constant be $\frac1{26}$? $\endgroup$ May 24, 2018 at 18:46
  • $\begingroup$ Yes, the volume is equal to that fraction. I will correct that right away. $\endgroup$
    – Wyllich
    May 24, 2018 at 19:37
  • 1
    $\begingroup$ Can you get access (by whatever means you feel appropriate) to the paper by Robbins and Bolis cited in the wiki page for Robbbins' constant? If not, in addition to me having to type up the equation (for box of general side lengths), you'll have to trust me that I quote their result correctly. $\endgroup$ May 24, 2018 at 21:31
  • $\begingroup$ Thank you very much for mentioning the name of Robbin's constant. I now possess a starting point from which I may derive the solution. Thank you again. I'll try to find the source material if possible, count on me. $\endgroup$
    – Wyllich
    May 24, 2018 at 21:41

1 Answer 1



Denote $E(a,b,c)$ as the expectation seen in the middle of p.278 (right before Eq.(2)) of the short paper (communication) by Robbins and Bolis. It is the average distance for a box of side lengths that are double of the inputs: $2a, 2b$, and $2c$.

Note that the inverse hyperbolic sine has a logarithmic form: $\sinh^{-1} x = \ln (x + \sqrt{x^2 + 1})$. That's how one ends up seeing those logs when $E(a,b,c)$ is evaluated numerically.

Denote $\mathcal{I}$ as the integral for your desired average but without the $1/26$ normalization.

Similarly, we will have the integrals $I_{f}$ for face-to-face, $I_{e}$ for edge-to-edge, and $I_{v}$ for vertex-to-vertex. All these are just integrals without normalization (thus not averages themselves). Details continue in the following sections.

Denote $R_b$ as the Robbins' constant, namely, $R_b = E(\frac12, \frac12, \frac12) \approx 0.6617$. Note that $R_b$ is the average, while it equals in value to its corresponding integral, which shall be denoted $I_0$. That is, numerically $R_b = I_0$ because the volume is $1$.

Strategy Outline

As you have noticed, there are obstacles in directly setting up the density. The symmetry and nice numbers (integers) strongly suggest that one should solve this via decomposing the space into blocks of unit cubes.

Any cube has $6$ faces, $12$ edges (sides), and $8$ vertices, therefore

$$\mathcal{I} = 6 I_f + 12 I_e + 8 I_v \tag{1} \label{Eq_26-decomp}$$

Your $V_1$ is the unit cube at the center, surrounded by $26$ unit cubes that are disjoint and which adds up exactly to your $V_2$.

The $26$ blocks consist of three groups: six of them are face-to-face with $V_1$, twelve of them are edge-touching-edge-only with $V_1$, and the remaining eight are vertex-touching-vertex-only with $V_1$.

The value of each $I_f$, $I_e$, and $I_v$ will be obtained successively using $E(a,b,c)$. Only after $I_f$ is derived can one calculate $I_e$. Then only after both $I_f$ and $I_e$ are known can one determine $I_v$.

Each stage of the process will involve respectively a $2$-$1$-$1$ long box, a $2$-$2$-$1$ flat box, then finally the $2$-$2$-$2$ "double cube".

Calculation Details

Formally, the integrals can be defined as (with more shorthands in additional to yours) $$ D_{1,2} \equiv \sqrt{(x_1-x_2)^2+(y_1-y_2)^2+(z_1-z_2)^2} \\ \begin{align*} I_f = \iiint\limits_{ V_1} \iiint\limits_{(x_2, y_2,z_2) \in V_{\mathbf{F}}} D_{1,2} \,\mathrm{d}x_{1,2}\,\mathrm{d}y_{1,2}\, \mathrm{d}z_{1,2} \\ I_e = \iiint\limits_{ V_1} \iiint\limits_{(x_2, y_2,z_2) \in V_{\mathbf{E}}} D_{1,2} \,\mathrm{d}x_{1,2}\,\mathrm{d}y_{1,2}\, \mathrm{d}z_{1,2} \\ I_v = \iiint\limits_{ V_1} \iiint\limits_{(x_2, y_2,z_2) \in V_{\mathbf{V}}} D_{1,2} \,\mathrm{d}x_{1,2}\,\mathrm{d}y_{1,2}\, \mathrm{d}z_{1,2} \end{align*} $$ The volume $V_F$ can be any unit cube adjacent and face-to-face with $V_1$. For example, $V_F \equiv \bigl\{(x,y,z) \big|~~ x \in [1,2],\, y \in [0,1],\, z \in [0,1]\bigr\}$

Similarly, the edge-to-edge volume $V_E$ can be, for example, $\bigl\{x \in [1,2],\, y \in [0,1],\, z \in [1,2]\bigr\}$, and the vertex-to-vertex volume $V_V$ can be $\bigl\{x \in [1,2],\, y \in [1,2],\, z \in [1,2]\bigr\}$.

The formulation of $I_f$ (or $I_e$, or $I_v$) is not unique (up to $6$, $12$, $8$ multiplicity for any $V_1$, which itself can resides anywhere), but it doesn't matter. The value of each integral is unique and can be obtained using $E(a,b,c)$.

###Successive Build-up (1): $~$ starting from $I_f$

Consider a $2$-by-$1$-by-$1$ long box, which is equivalent to two $V_1$ joined face-to-face. Denote the box as $V_{_{2.1.1}}$ . We know that its volume is $|V_{_{2.1.1}}| = 2$ .

For the average distance, decompose the integral (decompose the space) into two within-cube and two across-cubes.

$$|V_{_{2.1.1}}|^2 \cdot E( 1,\tfrac12, \tfrac12) = 2 I_0 + 2 I_f \tag{2} \label{Eq_long-box_decomp} $$

There are two $I_0$ because each half of $V_{_{2.1.1}}$ are physically distinct so that ordering matters. For the same reason, there are two $I_f$.

Note that $I_0$ here emphasizes that it is the integral and not the average $R_b$, although numerically they are the same.

Now is the time we must make explicit use of $E(a,b,c)$ by Robbins and Bolis for the general box. \begin{multline*} E(1, \tfrac12, \tfrac12) = \frac1{1260} \Bigl[520 + 17\sqrt{2} - 20\sqrt{5} - 81 \sqrt{6} + 21 \log\left(\sqrt{2}+1\right) \\ + 168 \log \left(\sqrt{2}+\sqrt{3}\right) + 735 \log\Bigl( \frac{ \sqrt{6} + 1 }{ \sqrt{5} }\Bigr) + 1344 \log \Bigl(\frac{\sqrt5 +1}2 \Bigr) \\ + 42 \log \left(\sqrt{5}+2\right) - 1344 \arcsin( 1/5 ) - 168 \arcsin\bigl( \sqrt{2/5} \bigr) \Bigr] \end{multline*} In principle, the algebraic calculation of the $E(1,\frac12,\frac12)$ here as well as later $E(a,b,c)$ can be done by hand, but I opted to code $E(a,b,c)$ with Mathematica since numerical evaluations are also needed.

The above display is just to demonstrate what we are dealing with here. To save space, the exact algebraic expression for each use of $E(a,b,c)$ will be omitted.

Moving on, $E(1,\frac12,\frac12) \approx 0.9145524845966431393$ and we have \begin{align*} I_f = \frac{ |V_{_{2.1.1}}|^2 }2 E( 1,\tfrac12, \tfrac12) - I_0 &= 2 E( 1,\tfrac12, \tfrac12) - Rb \tag{3.a} \label{Eq_Int-face_in_Rb} \\ &\approx 1.16739778692611 \tag{3.b} \label{Eq_Int-face_value} \end{align*}

###Successive Build-up (2): $~$ from $I_f$ to $I_e$

Consider a $2$-by-$\mathbf{2}$-by-$1$ flat box, which is equivalent to four $V_1$ joined together with pairwise face-to-face. Denote this box as $V_{_{2.2.1}}$. We know that its volume is $|V_{_{2.2.1}}| = 4$.

The decomposition of the integral (the decomposition of space) goes like this: for each of the four unit cubes, there is one within-cube integral, two face-to-face integral, and one edge-to-edge integral. $$|V_{_{2.2.1}}|^2 \cdot E( 1, 1, \tfrac12) = 4 \left( I_0 + 2 I_f + I_e \right) \tag{4} \label{Eq_flat-box_decomp} $$ Again, the cubes are physically distinct and order matters. Note that the decomposition counts the spaces as $4^2 = 4\cdot (1 + 2 + 1)$.

Plugin the $E(a,b,c)$ formula to get $E(1, 1,\frac12) \approx 1.13208746961$ and then replacing $I_f$ with Eq\eqref{Eq_Int-face_in_Rb}, we have: \begin{align*} I_e = \frac{ |V_{_{2.2.1}}|^2 }4 E( 1, 1, \tfrac12) - I_0 - 2I_f &= 4E( 1,1, \tfrac12) - R_b - 2 \bigl( 2 E( 1,\tfrac12, \tfrac12) - R_b \bigr) \\ &= 4 \bigl( E( 1, 1, \tfrac12) - E( 1,\tfrac12, \tfrac12) \bigr) + R_b \tag{5.a} \label{Eq_Int-edge_in_Rb} \\ &\approx 1.5318471223279 \tag{5.b} \label{Eq_Int-edge_value} \end{align*}

###Successive Build-up (3): $~$ get $I_v$ from $I_f$ and $I_e$

Consider a $2$-by-$2$-by-$2$ "double cube", which is equivalent to eight $V_1$ joined together with pairwise face-to-face. It is also just $V_1$ scaled up by a factor of two. Denote this box as $V_{_{2.2.2}}$.

We know that its volume is $|V_{_{2.2.2}}| = 2^3 = 8$, and this time we also know that $E(1,1,1) = 2 E( \frac12, \frac12, \frac12) = 2R_b$ simply due to the linear scalability of the geometry.

The integrals (the space) decompose similarly. For each of the eight unit cubes, there is one within-cube integral, three face-to-face integral, three edge-to-edge integral, and one vertex-to-vertex integral. $$|V_{_{2.2.2}}|^2 \cdot E( 1, 1, 1) = 8 \left( I_0 + 3 I_f + 3I_e + I_v \right) \tag{6} \label{Eq_double-cube_decomp} $$ As before, the order matters and the decomposition says $8^2 = 8\cdot (1 + 3 + 3 + 1)$.

Take $E( 1, 1, 1) = 2R_b\,,\,$ take $I_f$ as Eq\eqref{Eq_Int-face_in_Rb}, and take $I_e$ as Eq\eqref{Eq_Int-edge_in_Rb}, we end up with

\begin{align*} I_v &= \frac{ |V_{_{2.2.2}}|^2 }8 E( 1, 1, 1) - I_0 - 3I_f - 3 I_e \\ &= 8 \cdot 2R_b - R_b - 3 \Bigl[ 2 E( 1,\tfrac12, \tfrac12) - R_b + 4 \bigl( E( 1, 1, \tfrac12) - E( 1,\tfrac12, \tfrac12) \bigr) + R_b \Bigr] \\ &= 15 R_b - 6\bigl[ 2E( 1, 1, \tfrac12) - E( 1,\tfrac12, \tfrac12) \bigr] \tag{6.a} \label{Eq_Int-vertex_in_Rb} \\ &\approx 1.82787300624563276 \tag{6.b} \label{Eq_Int-vertex_value} \end{align*}

Finally, putting things together into Eq\eqref{Eq_26-decomp} we get:

\begin{align*} \mathcal{I} &= 6I_f + 12I_e + 8 I_v \\ &= 2 \Bigg\{ 3\bigl( 2 E(1,\tfrac12,\tfrac12 ) - R_b\bigr) + 6 \Bigl[ R_b + 4 \bigl( E(1, 1,\tfrac12 ) - E(1,\tfrac12,\tfrac12 )\bigr) \Bigr] + 4 \Bigl[ 15R_b - 6 \bigl( 2E(1, 1,\tfrac12 ) - E(1,\tfrac12,\tfrac12 )\bigr) \Bigr]\Bigg\} \\ &= 2 \bigl( 63R_b + 6 E(1,\tfrac12,\tfrac12 ) - 24 E(1,1,\tfrac12 ) \bigr) \tag{7.a} \label{Eq_Int-desired_in_Rb} \\ &\approx 40.0095362394564449 \tag{7.b} \label{Eq_Int-desired_value} \\ \end{align*}

The desired average is $\frac1{26} \mathcal{I}\approx 1.53882831690217095767637627$

The respective sub-equations Eq\eqref{Eq_Int-face_in_Rb}, Eq\eqref{Eq_Int-edge_in_Rb}, Eq\eqref{Eq_Int-vertex_in_Rb}, and Eq\eqref{Eq_Int-desired_in_Rb} showcase the expression in minimal terms of $E(a,b,c)$ where terms are converted to $R_b$ whenever possible.

These equations are reminders that the exact solution exists and practically available, that the functional form still exhibits some recognizable geometry, and that one can calculate any intermediate or final results to as many digits as one like.

Concluding Remarks

It should be clear that one cannot arrive at the desired integrals by using solely $R_b$ and nothing else. After carving out the contribution of $V_1$ at the center from a triple-cube, one still has to deal with the $V_2$-to-$V_2$ cross terms. In the end, one will have to do something similar to the derivation above.

When the space is in a nice configuration (smooth boundary with no zig-zags, no holes within ,etc), the exact integral in general is manageable, yet already requiring quite some attention to details. It is very telling that for the derivation of $E(a,b,c)$, Robbins and Bolis wrote that "By tedious yet routine successive integration one finds that ..."

For the sake of completeness, for related in-site posts one can start with rectangle in $2$-dim and more generally for hypercubes.

It is totally feasible to derive the density (of the distance between two points) from scratch for this particular $\{ V_1, \, V_2 \}$ pair. The transformations of random variables involved are all standard (squaring, addition, square-root). It is just that due to the "hole" in $V_2$, inevitably one will be forced to decompose the space.

Examining the whole process of "decomposing the space into disjoint spaces", one can see that this simple concept can be modified and applied to more complicated configurations, while having fundamental limitations.

In terms of numerical evaluation, once you have coded the $E(a,b,c)$ in the computer algebraic system of your choosing, this procedure is pretty straightforward. Regardless of your original motivation for doing this calculation, but I'd say this approach is practically useful across disciplines, from theoretical to applications.

  • $\begingroup$ What a ride! Thank you for posting such a detailed and comprehensive answer. I will definitely remember the idea of dividing the integral space as you masterfully did to obtain intermediate integrals through multiple iterations. I guess I owe you an explanation about the whys. In LiDAR image processing, I asked myself the possibility of predicting the resolution (here understand mean distance between a point and its closest neighbour) of a downsampled cloud using a voxel grid. In spite of my hypothesis based on the distance you calculated being disproved by experience, I still believe it (1) $\endgroup$
    – Wyllich
    May 28, 2018 at 12:46
  • $\begingroup$ ... was a good experience to have. Thank you very much again. (2) $\endgroup$
    – Wyllich
    May 28, 2018 at 12:47
  • $\begingroup$ Very glad I could help. I know only a little bit about image process or data compression in general, but what you described sounds no lesser than some innovative ideas that eventually became successful. $\endgroup$ May 28, 2018 at 13:12
  • $\begingroup$ Now that I think about it, your method gave me new ideas to precise my hypothesis which is based on isotropy (hence V2 containing V1). However, I may add some coefficients in front of each term of equations (2), (4), (6) to take into account anisotropy and capture an idea about the cloud's global structure. I must repeat myself : your dedication gave birth to many new ideas I may now explore. $\endgroup$
    – Wyllich
    May 28, 2018 at 13:25
  • $\begingroup$ Thank you for saying that. It's been a real pleasure. $\endgroup$ May 28, 2018 at 13:46

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .