mincGetWorldVoxel is really mincGetWorldVoxelNeighbour in disguise, but only sometimes

[yohanyee]

Issue

Basically, mincGetWorldVoxel() sometimes returns the value of a neighbouring voxel, as compared to what Display renders. Given world coordinates that lie between image sampling points, the output of mincGetWorldVoxel() is also sometimes inconsistent with the output of mincConvertWorldToVoxel() %>% mincGetVoxel(...).

It appears that mincGetWorldVoxel() follows the convention where voxels fall between sampling points. For example, consider an image that contains three voxels, starting at the origin and with 1 mm separations, with values: 21, 10, and 28. mincGetWorldVoxel() sees this as:

# World coordinates : ... 0.8---0.9---1.0---1.1---1.2---1.3---1.4---1.5---1.6---1.7---1.8---1.9---2.0---2.1 ...
# Sampling          :                  |                                                           |
# Voxel             : ...  0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2  ...
# Value             : ...  21 21 21 21 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 28 28 28 ...

On the other hand, mincConvertWorldToVoxel() and Display follow the convention where voxels are centered at sampling points:

# World coordinates : ... 0.8---0.9---1.0---1.1---1.2---1.3---1.4---1.5---1.6---1.7---1.8---1.9---2.0---2.1 ...
# Sampling          :                  |                                                           |
# Voxel             : ...  1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2  ...
# Value             : ...  10 10 10 10 10 10 10 10 10 10 10 10 10 10 28 28 28 28 28 28 28 28 28 28 28 28 28 ...

---

Test data

Below is an R script to reproduce the issue, and a Python script to generate a test minc file for this:

Generate test minc file (python)

#!/usr/bin/env python3

from pyminc.volumes.factory import volumeFromDescription
import numpy as np

img = np.array([[[21, 10, 28]]], dtype=np.int16)

vol = volumeFromDescription("test_file.mnc", dimnames=["zspace", "yspace", "xspace"], sizes=[1,1,3], starts=[0,0,0], steps=[1,1,1])
vol.data = img
vol.writeFile()
vol.closeVolume()

R code that demonstrates the issue

# R version: 4.1.3
# RMINC version: 1.5.3.0

library(RMINC)

# Test file
# This image contains three voxels along the xspace dimension, starting at the origin and with 1 mm separations
# Values are: 21, 10, 28
test_file <- "test_file.mnc"

# Confirm shape
# z, y, x
minc.dimensions.sizes(test_file)

# Confirm intensities given voxel coordinates
# z, y, x; starts at 0
as.integer(round(mincGetVoxel(test_file, 0, 0, 0)))==21
as.integer(round(mincGetVoxel(test_file, 0, 0, 1)))==10
as.integer(round(mincGetVoxel(test_file, 0, 0, 2)))==28

# Confirm intensities given world coordinates
# x, y, z; starts at 0
as.integer(round(mincGetWorldVoxel(test_file, 0, 0, 0)))==21
as.integer(round(mincGetWorldVoxel(test_file, 1, 0, 0)))==10
as.integer(round(mincGetWorldVoxel(test_file, 2, 0, 0)))==28

# mincGetWorldVoxel follows the convention where voxels fall between sampling points
# World coordinates : ... 0.8---0.9---1.0---1.1---1.2---1.3---1.4---1.5---1.6---1.7---1.8---1.9---2.0---2.1 ...
# Sampling          :                  |                                                           |
# Voxel             : ...  0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2  ...
# Value             : ...  21 21 21 21 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 28 28 28 ...
as.integer(round(mincGetWorldVoxel(test_file, 1, 0, 0)))   # Returns: 10
as.integer(round(mincGetWorldVoxel(test_file, 0.9, 0, 0))) # Returns: 21
as.integer(round(mincGetWorldVoxel(test_file, 1.9, 0, 0))) # Returns: 10
as.integer(round(mincGetWorldVoxel(test_file, 1.6, 0, 0))) # Returns: 10
as.integer(round(mincGetWorldVoxel(test_file, 1.5, 0, 0))) # Returns: 10

# But this is inconsistent with the conversion to voxel coordinates, which centers voxels at sampling points
# World coordinates : ... 0.8---0.9---1.0---1.1---1.2---1.3---1.4---1.5---1.6---1.7---1.8---1.9---2.0---2.1 ...
# Sampling          :                  |                                                           |
# Voxel             : ...  1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2  ...
# Value             : ...  10 10 10 10 10 10 10 10 10 10 10 10 10 10 28 28 28 28 28 28 28 28 28 28 28 28 28 ...
attr(mincGetWorldVoxel(test_file, 1, 0, 0), "voxelCoord")
attr(mincGetWorldVoxel(test_file, 0.9, 0, 0), "voxelCoord")
attr(mincGetWorldVoxel(test_file, 1.9, 0, 0), "voxelCoord")
attr(mincGetWorldVoxel(test_file, 1.6, 0, 0), "voxelCoord")
attr(mincGetWorldVoxel(test_file, 1.5, 0, 0), "voxelCoord")

# Using mincGetVoxel on voxel coordinates taken from the mincGetWorldVoxel attributes (which in turn relies on
# mincConvertWorldToVoxel) returns different values as compared to a direct mincGetWorldVoxel call, due to this
# mismatch in conventions
mincConvertWorldToVoxel(test_file, 1, 0, 0)   # Voxel coord: 0 0 1
mincConvertWorldToVoxel(test_file, 0.9, 0, 0) # Voxel coord: 0 0 1
mincConvertWorldToVoxel(test_file, 1.9, 0, 0) # Voxel coord: 0 0 2
mincConvertWorldToVoxel(test_file, 1.6, 0, 0) # Voxel coord: 0 0 2
mincConvertWorldToVoxel(test_file, 1.5, 0, 0) # Voxel coord: 0 0 2
mincGetVoxel(filenames = test_file, 0, 0, 1) # Returns: 10
mincGetVoxel(filenames = test_file, 0, 0, 1) # Returns: 10
mincGetVoxel(filenames = test_file, 0, 0, 2) # Returns: 28
mincGetVoxel(filenames = test_file, 0, 0, 2) # Returns: 28
mincGetVoxel(filenames = test_file, 0, 0, 2) # Returns: 28

# Note: Display renders the voxels based on this second convention that centers voxels at sampling points

[gdevenyi]

The official MINC standard is coordinates of a voxel are the center, see slide 14
https://www.bic.mni.mcgill.ca/software/workshop2003/02-2003MINC-meeting_js.pdf