:bug: Surface initialization fails on non-isotropic voxels

[Hu-nie]

I am getting the following error while running the example, what is the reason?

[vanandrew]

I'll need more information to debug this. Can you please specify:

your OS
Python version
Version of brainextractor and dependencies (especially numba)

[Hu-nie]

@vanandrew Thanks for the reply, I'll send you an information to your question.

OS : mac os Python version : 3.7 Version of brainextractor and dependencies (especially numba) Package Version

---

aspose-cells 22.4.0 binvox 0.1.5 brainextractor 0.2.1 certifi 2021.10.8 charset-normalizer 2.0.12 ci-info 0.3.0 click 8.1.3 cycler 0.11.0 dicom2nifti 2.4.5 etelemetry 0.3.0 filelock 3.8.0 fonttools 4.30.0 freetype-py 2.3.0 geomdl 5.3.1 html-table-extractor 1.4.1 html-table-parser-python3 0.2.0 idna 3.3 imageio 2.22.1 importlib-metadata 5.0.0 isodate 0.6.1 JPype1 1.3.0 kiwisolver 1.3.2 llvmlite 0.39.1 looseversion 1.0.2 lxml 4.9.1 matplotlib 3.5.1 mpmath 1.2.1 networkx 2.6.3 nibabel 4.0.2 nipype 1.8.5 numba 0.56.3 numpy 1.21.5 opencv-python 4.5.5.64 packaging 21.3 pandas 1.1.5 Pillow 9.0.1 pip 22.2.2 prov 2.0.0 pydicom 2.2.2 pydot 1.4.2 pyglet 2.0.0 PyMCubes 0.1.2 PyOpenGL 3.1.0 pyparsing 3.0.7 PyQt5 5.15.6 PyQt5-Qt5 5.15.2 PyQt5-sip 12.9.1 pyrender 0.1.45 PySide6 6.3.0 PySide6-Addons 6.3.0 PySide6-Essentials 6.3.0 python-dateutil 2.8.2 python-gdcm 3.0.19 pytz 2021.3 PyWavefront 1.3.3 PyWavelets 1.3.0 rdflib 6.2.0 requests 2.27.1 scikit-image 0.19.3 scipy 1.7.3 seaborn 0.11.2 setuptools 58.0.4 Shapely 1.8.1.post1 shiboken6 6.3.0 SimpleITK 2.2.0 simplejson 3.17.6 six 1.16.0 soupsieve 2.3.2.post1 sympy 1.10.1 tifffile 2021.11.2 tqdm 4.63.0 traits 6.3.2 trimesh 3.15.8 typing_extensions 4.2.0 urllib3 1.26.9 voxlib 0.0.5 wheel 0.37.1 zipp 3.10.0

[vanandrew]

I think it might be easier to debug this without numba enabled. If you could follow these instructions, and then send me back what output you get:

1. Uninstall the current version of brainextractor you have installed (I see you are using anaconda, so it might be easier to just start a new fresh enviroment)
2. Clone this repo
3. Switch to the debug branch: git checkout debug
4. Install the debug version: pip install --user [path/to/the/git/repo]. When it is installed you should see version: 0.2.1d
5. Then re-run brainextractor. You should also see DEBUG Version at the top of the output to know you are using the right version.
6. If it errors, copy/paste back the error here.

[Hu-nie]

@vanandrew

When installing with the debug version, the following error occurs.

Also, the same error occurs when proceeding with the following command

pip install git+https://github.com/vanandrew/brainextractor@debug

[vanandrew]

That was an error on my part (forgot there is the PEP440 standard for versioning).

Just pushed a fix to the same branch, so try it again.

[Hu-nie]

While running, I get an error different from before.

I think it's my data problem

my data is 3d tof mra data.

I changed the original dicom format file to nifti format. The data format is (880,880,192).

[vanandrew]

Thanks! Seems like this is the original error that numba was "hiding".

What are the dimensions of your voxels? It looks like the brain extractor picked up that they are "non-isotropic", so they must be different in x/y vs. z (especially since you said the dimensions of your image is 880 x 880 x 192).

[vanandrew]

Ah. Ok I think I see the issue. The way the head radius estimation works, it assumes iso-tropic voxels. So my code for creating the initial sphere is off (and this issue is more severe when you have highly asymmetric dimensions in your data). Basically, this logic needs to be redone (but with some sort of ellipsoid, instead of a sphere) (L90-L101 in main.py):

    # compute 1/2 head radius with spherical formula
    self.r = 0.5 * np.cbrt(3 * np.sum(self.rdata > self.t) / (4 * np.pi))
    print("Head Radius: %f" % (2 * self.r))

    # get median value within estimated head sphere
    self.tm = np.median(self.data[sphere(self.shape, 2 * self.r, self.c)])
    print("Median within Head Radius: %f" % self.tm)

    # generate initial surface
    print("Initializing surface...")
    self.surface = trimesh.creation.icosphere(subdivisions=4, radius=self.r)
    self.surface = self.surface.apply_transform([[1, 0, 0, ci], [0, 1, 0, cj], [0, 0, 1, ck], [0, 0, 0, 1]])

I don't think I have the bandwidth at the moment to work on a fix for this edge case. But here's another alternative I can suggest: You can resample your data to have isotropic voxels (You'll lose resolution in x/y), but it'll allow you to get the brain mask with this code. Then resample that mask back to 880 x 880 x 192.

Or you can resample the z axis to be the same size as the x/y axes, but then the data would be quite large and I'm not sure you'll get much benefit.

[Hu-nie]

After resizing the voxel size to 1mm x 1mm x1mm, I think I need to find a way to put the data in the form of an array.

[vanandrew]

There are a bunch of ways to do resampling, but here's a solution in pure python (nibabel):

import nibabel as nib
from nibabel.processing import resample_to_output

img = nib.load("/path/to/nii")

img2 = resample_to_output(img, [1, 1, 1])  # 2nd argument is voxel dimensions, so this will resample to 1 mm x 1 mm x 1 mm

img2.to_filename("/path/to/new_file")

See resample_to_output for more details.

[vanandrew]

And to go back to the original space, you can use the resample_from_to function (e.g. so you can get the mask back in the original space):

# continuing from above...

from nibabel.processing import resample_from_to

mask = nib.load("/some/mask")
img3 = resample_from_to(mask, img)

# you may need to do a threshold here, resampling will create fractional values and the mask should be logical [0/1]
thresholded_data = img3.get_fdata() > 0.5

# save the resampled mask
nib.Nifti1Image(thresholded_data, img1.affine).to_filename("resampled_mask.nii.gz")