# Mean distance between 2 points in adjacent voxels (cubes)

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?

Thanks.

• 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. – Wyllich May 24 '18 at 15:21
• 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}$? – Lee David Chung Lin May 24 '18 at 18:46
• Yes, the volume is equal to that fraction. I will correct that right away. – Wyllich May 24 '18 at 19:37
• 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. – Lee David Chung Lin May 24 '18 at 21:31
• 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. – Wyllich May 24 '18 at 21:41

## Notations

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$$. For a different configuration this approach needs modification since the space doesn't decompose nicely in general.

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 using $$E(a,b,c)$$, successively building up starting from $$I_f$$.

It will involve 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 (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 \left(\frac{1}{2} \left(\sqrt{5}+1\right)\right) \\ + 42 \log \left(\sqrt{5}+2\right) - 1344 \arcsin\bigl( \frac15 \bigr) - 168 \arcsin\bigl( \sqrt{\frac25 } \bigr) \Bigr] \end{multline*}$$ I coded $$E(a,b,c)$$ with Mathematica and would never consider doing this by hand. The $$E(1,\frac12,\frac12)$$ here as well as later $$E(a,b,c)$$ evaluations are done by the computer algebraic system.

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

Moving on, $$E(1,\frac12,\frac12) \approx 0.91455248459664313930734$$ 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 (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.1320874696118224667$$ 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 (spaces) 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 are the reminders that the exact solution is obtainable, 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 integral solely by using $$R_b$$. 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 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 from scratch for this particular $$\{ V_1, \, V_2 \}$$ pair. The transformation of random variables involved is pretty basic. It is just that due to the "hole" in $$V_2$$, inevitably one will be forced to decompose the space.

"Decomposing the space into disjoint spaces" is a standard technique in dealing with probability inquiries, and the decomposition is not necessarily into integer valued partitions. There are indeed limitations to this approach, but the concept is rather simple.

In terms of numerical calculation, once you have coded the $$E(a,b,c)$$ in whatever system you're using, this procedure is pretty straightforward.

I don't know your original motivation for doing this calculation, but I'd say this approach is practically useful across disciplines, from theoretical to applications.

• 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) – Wyllich May 28 '18 at 12:46
• ... was a good experience to have. Thank you very much again. (2) – Wyllich May 28 '18 at 12:47
• 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. – Lee David Chung Lin May 28 '18 at 13:12
• 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. – Wyllich May 28 '18 at 13:25
• Thank you for saying that. It's been a real pleasure. – Lee David Chung Lin May 28 '18 at 13:46