Integer overflow in `get_transformation_matrix()` due to dtype-dependent behavior of `np.indices` (NumPy 1.x vs 2.x)

[Abd-zak]

Title:

Description

The function SpaceConverter.get_transformation_matrix() can produce incorrect transformation matrices and voxel sizes depending on the NumPy version and build, due to implicit use of integer dtypes.

The issue arises because np.indices(self._shape) returns an array with a default integer dtype that is platform- and version-dependent. In some environments (e.g., NumPy 1.26.x on Windows/linux), this dtype is int32, which leads to integer overflow during subsequent matrix computations.

---

Root cause

The problematic code is:

ortho_grid = []
for i in np.indices(self._shape):
    ortho_grid.append(i - i[shift_voxel])
ortho_grid = np.array(ortho_grid).reshape(3, np.prod(self._shape))

G = np.dot(ortho_grid, ortho_grid.T)

If ortho_grid is int32, then:

• G = ortho_grid @ ortho_grid.T is computed in int32
• Values exceed the 32-bit integer range (~2.1e9)
• This causes silent overflow
• The resulting G is incorrect (e.g., negative diagonal values)
• This propagates through inversion and produces incorrect transformation matrices and voxel sizes

---

Observed behavior

Good environment

• NumPy: 2.4.2
• Python: 3.12
• np.indices(...).dtype: int64

Results:

• Correct Gram matrix (G)
• Correct transformation matrix
• Correct voxel size (~8–10 nm)

Bad environment

• NumPy: 1.26.3
• Python: 3.10 / 3.12
• np.indices(...).dtype: int32

Results:

• Overflowed Gram matrix (negative / incorrect values)
• Incorrect transformation matrix
• Voxel size underestimated by ~10–80×

---

Important note

This issue is:

• ❌ Not OS-specific
• ❌ Not strictly a NumPy version bug
• ✅ A dtype safety issue in cdiutils

NumPy 2.0 changed the default integer type on 64-bit Windows from int32 to int64, which explains why the issue disappears in newer environments. However, relying on this implicit behavior is unsafe.

---

Suggested fix

Force a safe dtype explicitly when creating ortho_grid.

Recommended fix

ortho_grid = np.asarray(
    [i - i[shift_voxel] for i in np.indices(self._shape, dtype=np.int64)],
    dtype=np.float64,
).reshape(3, np.prod(self._shape))

Minimal fix

ortho_grid = np.array(ortho_grid, dtype=np.float64).reshape(3, np.prod(self._shape))

This ensures:

• No integer overflow
• Consistent behavior across environments
• Numerical stability

---

Expected behavior

get_transformation_matrix() should produce identical results for identical input data, regardless of NumPy version or platform.

---

Actual behavior

Results depend on the dtype returned by np.indices(), leading to silent numerical errors in some environments.

---

Conclusion

The function should not rely on NumPy’s default integer dtype. Explicitly casting to int64 or float64 is necessary to ensure correctness and reproducibility.

---