MorphologicalEulerNumber misbehaving in 3D?

Question

Context

I would like to compute the MorphologicalEulerNumber of 3D GaussianRandomField as a function of height above a given threshold.

Attempt

I proceed as follows:

  tt = Table[
  dat = GaussianRandomField[nn = 64, 3, 
     Function[k, Exp[-(k/10)^2] k^(-2)]] // Chop; 
  dat -= Mean[Flatten[dat]]; dat /= RootMeanSquare[dat // Flatten]; 
  ParallelTable[{x, 
    MorphologicalEulerNumber[
     Map[If[# < x, 1., 0.] &, dat, {3}] // Image3D]}, {x, -4, 4, 
    0.1}], {15}]; 

Then

pl1 = {Mean@tt, -Mean@tt} // ListLinePlot

which has clearly the wrong parity. It should behave like the cumulative alternate sum of critical points bellow the chosen threshold, so it should be even w.r.t. to height.

Indeed theory predicts that it should scale like

 euler[x_]=(-1+x^2)/(12*Sqrt[3]*E^(x^2/2)*Pi^2);
 Plot[euler[x] 10^5, {x, -4, 4}]

which looks quite different.

Question

It this a bug in MorphologicalEulerNumber in 3D? Or am I missing something obvious?

Comment

Note that the same test in 2D works as expected, up to edges effects which can be addressed by Henrik periodicMorphologicalEulerNumber function (which does not take Image but matrices).

tt = Table[
  dat = GaussianRandomField[nn = 1024, 2, 
     Function[k, Exp[-(k/10)^2] k^(-1)]] // Chop; 
  dat -= Mean[Flatten[dat]]; dat /= RootMeanSquare[dat // Flatten]; 
  ParallelTable[{x, 
    periodicMorphologicalEulerNumber[
     Map[If[# < x, 1., 0.] &, dat, {2}]]}, {x, -4, 4, 
    0.1}], {50}];

Then

  pl1 = {Mean@tt, -Mean@tt} // ListLinePlot

which matches qualitatively the theoretical expectation:

euler2D[x_]=-((E^(-(x^2/2)) x)/(4 Sqrt[2] \[Pi]^(3/2)))
Plot[euler[x] 2 10^3, {x, -4, 4}] 

I am using Mathematica 12.3 on MacOs X.15.

Answer 1

I thought a bit about this problem. My first idea was to create a hex-mesh for the periodic image and to compute MeshCellCounts from there. Unfortunately, MeshRegion does not allow me to do so. But even before I realized that, it came to my mind that building the whole mesh data structure is just unnecessarily heavy on memory when we only need the numbers of vertices, edges, quads, and hexes. So I came up with the following instead:

Each hex knows which vertices, edges, or quads belong to it. One can count the number of vertices as follows: Take the binary 3D image and add rotated copies to it; namely we use RotateLeft with all binary vectors from $\{0,1\}^3$. Each entry of the resulting image belongs to a vertex of the full grid. The vertices belonging to the white pixels (and only them) now have a value greater than 0. So we can simply apply Unitize and Total to count them.

For edges and quads, this is slightly more complicated, since there are three types of edges and three types of quads (according to each coordinate direction). But the counting works similarly: RotateLeft the image into certain directions, add them, and then count the number of nonzeroes. The code that realizes this for any dimension d and facets of dimension k is shown below:

ClearAll[rotations];

rotations[d_, k_] := rotations[d, k] = Module[{subsets, tuples, rot},
    subsets = Subsets[Range[d], {d - k}];
    tuples = Tuples[{0, 1}, d - k];
    rot = ConstantArray[0, {Length[subsets], Length[tuples], d}];
    Do[
     rot[[i, All, subsets[[i]]]] = tuples,
     {i, 1, Length[subsets]}];
    rot
    ];

cellcount[A_?TensorQ, k_Integer] := With[{d = TensorRank[A]},
   Sum[
    Total[Unitize[Sum[RotateLeft[A, tuple], {tuple, tuplelist}]], d],
    {tuplelist, rotations[d, k]}]
   ];

PeriodicArrayEulerCharacteristic[A_?TensorQ] := 
 Sum[(-1)^k cellcount[A, k], {k, 0, TensorRank[A]}]

rotations[d,k] is computed only once per session and memoized and cellcount uses vectorization. So repeated calls to cellcount should be really fast. In my first experiments, this turned out to be about 100 times faster than building a hex mesh first.

I have not checked the correctness very thoroughly though; please report if you experience any odd behavior.

Edit: Performance improvements

This should generate the random data a bit quicker. Notable changes are:

• Compiling cf only once (and via C into a library).

• Avoiding Outer in the computation of amplitude; it tends to be rather slow.

• Re instead of Chop to remove small imaginary numbers. This should leave packed arrays packed.

---

cf = With[{f = Function[k, Exp[-(k/10.)^2] k^(-2)]},
   Compile[{{vec, _Real, 1}},
    With[{nrm = Norm[vec]}, If[nrm <= 0., 0., Sqrt[f[nrm]]]],
    CompilationTarget -> "C",
    RuntimeAttributes -> Listable,
    Parallelization -> True
    ]
   ];

GaussianRandomField[size : (_Integer?Positive) : 256, 
   dim : (_Integer?Positive) : 2, cf_CompiledFunction] := 
  Module[{Pkn, fftIndgen, noise, amplitude, s2},
   s2 = Quotient[size, 2];
   fftIndgen = ArrayPad[Range[0, s2], {0, s2 - 1}, "ReflectedNegation"];
   noise = 
    Fourier[RandomVariate[NormalDistribution[], 
      ConstantArray[size, dim]]];
   amplitude = ArrayReshape[cf[Tuples[fftIndgen, dim]], ConstantArray[size, dim]];
   Re[InverseFourier[noise amplitude]]
   ];

Moreover, I suggest to use

ParallelTable[
 {x, PeriodicArrayEulerCharacteristic[UnitStep[Subtract[x,dat]]]},
 {x, -3.5, 3.5, 0.1}]

instead of

ParallelTable[
 {x, PeriodicArrayEulerCharacteristic[Map[If[# < x, 1., 0.] &, dat, {3}]]},
 {x, -3.5, 3.5, 0.1}]

Answer 2

I think there are three possibilities:

1. this is a bug

2. Mathematica is using a definition of "Euler number" that does not mean "Euler characteristic"

3. the image is not being interpreted in the expected manner

tl;dr I think this is a bug, because it calculates a quantity which is not the Euler number under any definition I've seen, and is not even an acceptable algebraic-topology invariant.

Here's my exploration of what could be going on, including some tests to support this conclusion.

---

The second one seems plausible at first glance. Note that Wolfram Mathworld defines the Euler number differently than the Euler characteristic: for a finite complex $K$, it defines it as

$$\chi'(K) = \sum_{n \in \mathbb{N}} (-1)^n\ \text{rank}(C_n(K)) $$

Contrast this with the definition of the Euler characteristic, which is

$$ \chi(K) = \sum_n (-1)^n\beta_n(K) = \sum_n (-1)^n\text{rank}(H_n(K)) $$

However, even though these are different definitions, Wikipedia implicitly says these agree by defining the Euler characteristic as the Euler number for finite CW-complexes. (Wikipedia uses the alternating sum of number of cells of a given dimension, which translates to the rank of the chain space. I'll come back with a better source (or just a proof!) if I have time, but hopefully the rest renders this immaterial anyway...)

Further, consider the following example.

MorphologicalEulerNumber[
 ReplacePart[ConstantArray[1, {3, 3, 1}], {2, 2, 1} -> 0] // Image3D //
   Echo]

This produces a 3x3x1 "solid torus" of cubes, as you can see from the Echo. But MorphologicalEulerNumber gives 1, despite the fact that the Euler characteristic is 0. And yet, even if you manually add up all the cells involved, with alternating signs,

-8 + (8 + 4 + 2*8 + 3*4) - (3*4 + 2*3*4 + 2*4*2 + 4*2 + 4) + 4*4*2

you get 0.

This makes me lean towards "bug", since both the Euler characteristic and number ought to be 0 for this shape.

However, I wondered if maybe Mathematica was interpreting voxels at the boundary of the image itself differently. However, padding with zeros,

MorphologicalEulerNumber[
 ArrayPad[ReplacePart[ConstantArray[1, {3, 3, 1}], {2, 2, 1} -> 0], 
    1] // Image3D // Echo]

still returns 1.

Interestingly, though, supplying a Padding -> 1 option to the original 3x3x1 torus changes things:

MorphologicalEulerNumber[
 ReplacePart[ConstantArray[1, {3, 3, 1}], {2, 2, 1} -> 0] // Image3D //
   Echo, Padding -> 1]

This now gives 0! But not for the manually padded version, unfortunately. The documentation on what exactly a "padding value" even is in this context is a bit scarce. If it means, as I'm guessing, "the value that should be padded around the boundary of the image", then that suggests that by padding with 1's all round, Mathematica is essentially creating a 1x1x1 3-dimensional hole (void) in the center, and is now subtracting that from the count of the single connected component.

Indeed, if we make one of those voids ourselves, we get 0:

MorphologicalEulerNumber[
 ReplacePart[ConstantArray[1, {3, 3, 3}], {2, 2, 2} -> 0] // Image3D //
   Echo]

---

My guess, then, is that Mathematica is completely ignoring the first Betti number, and is instead using the following formula in 3 dimensions:

$$ \text{MorphologicalEulerNumber}(K) = \beta_0(K) - \beta_2(K) $$

This explains the value of 1: the solid torus of cubes has a single connected component and no 3D void.

It also explains why the curve produced in your case is $\beta_0 - \beta_2$, and is entirely missing the $\beta_1$ component; see the figure on page 3 of this paper (though they might have the horizontal axis reversed compared to your convention).

This is also backed up by the description given, as long as we interpret "holes" as 3-dimensional holes (counted by $\beta_2$):

MorphologicalEulerNumber[image] by default gives the total number of connected white regions in image, minus the number of black holes that occur inside those regions.

In my opinion this is a serious issue, because this is not homotopy invariant. If you're doing any (computational) algebraic topology whatsoever, you'll generally want your invariants to not just be topological invariants, but homotopy ones as well. However, if you imagine retracting this solid-cube-torus down into its 2-dimensional version, the MorphologicalEulerNumber changes:

MorphologicalEulerNumber[
 ReplacePart[ConstantArray[1, {3, 3}], {2, 2} -> 0] // Image // Echo]

This gives 0. But since these two shapes are homotopy equivalent, their invariants should be the same.

My guess is that this might have been an ad-hoc generalization from 2D to 3D which was made without the knowledge that the Euler number is already well-defined in 3 (and higher) dimensions—or at least without the knowledge that its definition no longer corresponds to "number of connected components minus top-dimensional holes", since this is "just a coincidence" in 2D.

So, I would call this a bug, as this definition does not produce any sort of Euler number or characteristic, which in any definition I've seen is required to be homotopy invariant. (However, it is of course possible that some community uses this niche definition, though I would be surprised, and it would be out of step with a very widespread convention.) At the very least, since the quantity calculated is not an Euler number, the function should be called something different—and then, there should also be a function which calculates the Euler number! :)