TractSeg icon indicating copy to clipboard operation
TractSeg copied to clipboard

IndexError: index 60 is out of bounds for axis 2 with size 60

Open mrsf23 opened this issue 3 years ago • 18 comments

Hello!

When we run tractometry, we get the following error for some, but not all of our subjects:

14%|███████████▊ | 7/50 [00:06<00:40, 1.07it/s] Traceback (most recent call last): File "/home/user/miniconda3/bin/Tractometry", line 138, in main() File "/home/user/miniconda3/bin/Tractometry", line 107, in main mean, std = tractometry.evaluate_along_streamlines(np.nan_to_num(scalar_image.get_fdata()), streamlines, File "/home/user/miniconda3/lib/python3.8/site-packages/tractseg/libs/tractometry.py", line 64, in evaluate_along_streamlines streamlines = _orient_to_same_start_region(streamlines, beginnings) File "/home/user/miniconda3/lib/python3.8/site-packages/tractseg/libs/tractometry.py", line 42, in _orient_to_same_start_region if beginnings[int(startpoint[0]), int(startpoint[1]), int(startpoint[2])] == 0: IndexError: index 60 is out of bounds for axis 2 with size 60

The data has isotropic spacing, is oriented correctly (we have also registered some of the data to MNI space, which did not help), and we have further tried using our own brain mask, which also did not help. The affine matrices in the headers appear to be the same for the different input images.

We have further noticed by looking at the data in MITK, that some of the bundles seem to be located outside the brain mask.

Do you have any idea what could cause the error? Thank you!

mrsf23 avatar Apr 16 '21 22:04 mrsf23

Hi, it might help to know a little more about how you are using TractSeg.

  1. For the subjects that go through tractometry successfully, do their outputs look correct?
  2. Did you try to rigidly align the data to MNI as outlined in the documentation (calculating FA and aligning to the FA template)?
  3. Are you using the latest version of TractSeg?
  4. How much work has been done on the input data before TractSeg? That is, is it raw DWI, preprocessed DWI, or pre-made peaks? If not raw, have you made sure preprocessing ran as intended?
  5. Have you looked at the peaks in MITK to make sure they look correct? I know for most of my data (Siemens scanner) I have to flip peaks across the x-axis.
  6. I am not sure what that "axis 2" is denoting (not a dev) here, but maybe it is the number of tracts created for that subject. Perhaps all 72 tracts were not created for those subjects.

Hope we can get to the bottom of this! Steven

smeisler avatar Apr 21 '21 03:04 smeisler

Hello Steven, thank you very much for your reply!

  1. Yes, the bundles look fine. However, they are a bit oblique in comparison to the brain mask and some streamlines pass the mask outlines. This is applies to all subjects, it seems worse for those who did not work (see pictures below).

  2. Yes, unfortunately it did not help. If we use the --preprocess command, the same error strangely already appears during the tracking command and the results look even worse (see third picture).

  3. Yes.

  4. We have tried pre-made peaks, as well as our preprocessed DWI data (processed in MRtrix with a standardised custom made pipeline) and both didn't work. As we've used the preprocessed DWI for other analyses, we are sure that there is no problem with the preprocessing.

  5. We looked at the peaks, they look correct in MITK after flipping them along the z-axis peaks (as decribed in the documentation). They look fine for both, the subjects for which the tractometry worked and also for those for which it didn't. We anyways tried to flip the peaks, what resulted in no bundles at all. So the peaks are most probably fine.

  6. it seems like in this issue somebody already had the same problem; see https://github.com/MIC-DKFZ/TractSeg/issues/123 Jakob Wasserthal replied that the dimensions of the dwi, peak, segmentation and FA image has to be the same. The "axis 2" is probably referring to some mismatching dimensions, we couldn't find any differences though.

In Conclusion, the only problem we could find what that the bundles are not located precisely within the brain mask. What we used in our bash script were the standard commands from the documentation.

Thank you for helping!

Subject which worked: subject_worked Subject which didn't work: subject_worked_not Results after --preprocess (aligning to MNI) MNI

mrsf23 avatar Apr 22 '21 18:04 mrsf23

Thanks for reporting all of that; I don't really know what may be the cause or what to suggest at this point, besides potentially sending the data to see if this error is replicable.

smeisler avatar Apr 22 '21 19:04 smeisler

File "/home/user/miniconda3/lib/python3.8/site-packages/tractseg/libs/tractometry.py", line 42, in _orient_to_same_start_region if beginnings[int(startpoint[0]), int(startpoint[1]), int(startpoint[2])] == 0: IndexError: index 60 is out of bounds for axis 2 with size 60

This error means that the endpoint of a streamline is not withing the start/endpoint mask generated by tractseg (endings_segmentation). Looking at the pictures it seems that things are not properly lined up. Therefore it makes sense that this error occurs.

However, it is strange that the tractogram ends up outside of the brainmask. If you look at the bundle masks generated by TractSeg (bundle_segmentations) are some of them also outside of the brainmask? Then I guess the affine matrix of your input image, of your brain mask and of the bundle_segmentations is not 100% identical.

If only the fibers (TOM_trackings) are outside you could check the --tracking_format argument. For propery visualization this should be trk_legacy and you have to use the same for both Tracking and Tractometry. But if you did not change the default this should be fine. But you could also try tck and look at the results with mrview just to make sure it is not an MITK issue (the visualization of trk and tck fibers in MITK is not so stable).

wasserth avatar Apr 23 '21 13:04 wasserth

The bundles (results from brundle_segmentations) are aligned with the brain mask, also in subjects which did not work: bundle_segmentation

We also checked the affine matrices, for which the only differences are in qform_code and sform_code (see pictures below). Could this be a problem?

We did not change the tracking_format default in our script. When changing to tck and opening the results in mrview, they look the same as in MITK: image

For the subjects which worked, the tracts are only slightly outside the mask: mrview_007

Is there something else we could try?

peaks.nii.gz image

Bundle_segmentations/AF_left.nii.gz image

Endings_segmentation/AF_left_e.nii.gz image

TOM/AF_left.nii.gz image

dwi/sub-004_dwi_preproc.nii.gz image

dwi/sub-004_dwi_b0mask.nii.gz image

mrsf23 avatar Apr 24 '21 00:04 mrsf23

As decribed previously, we had the issue that despite the bundles were nicely reconstructed, the bundles were misaligned to masks, peaks and FA maps. This caused the above mentioned error, despite peaks and masks were perfectly aligned and all files shared identical affine matrix.

I found out, that tracking on FODs using the options "--track_FODs SD_STREAM" or "--track_FODs iFOD2" worked perfectly. Further, I noticed that also tracking on TOMS worked perfectly when I changed the parameter "prob_tracking_software = 'mrtrix'" (instead of tractseg) in the code of bin/Tracking. I then found out, when I use default settings ("prob_tracking_software = 'tractseg'", tracking on TOMS) the misalignment is caused by the command to flip axes in the file tractseg_prob_tracking.py. When I commented out the code "streamlines = fiber_utils.invert_streamlines(streamlines, bundle_mask, affine, axis=axis)" in this file, the bundles were left-right flipped but aligned with mask and FA map. So the function "fiber_utils.invert_streamlines" flips this axis, what it actually does, but strangly it also applies a translation, which seems not to be necessary. When I comment out the line 417 in this function (file fiber_utils.py) " affine_invert[0, 3] = img_center_mm_space[1] * 2" everything works perfectly. So it seems to miscalculate a translation that is not necessary and that is why the bundles are all misaligned.

I have seen in previous GitHub issues that others had the same troubles, so I hope this helps for all struggling with this issue. I still would be very happy to hear any comment of @wasserth on that. Thank you!

mrsf23 avatar May 12 '21 14:05 mrsf23

Thanks a lot for this deep dive into the code. That really helps. I was not aware of this problem. Getting all the image and streamline transformations right is a big pain. Therefore, normally I rigidly align all images with MNI space before doing any analysis. In MNI space everything should work fine.

I remember there was a reason why I added this translation. I will try to figure that what is going on.

wasserth avatar May 12 '21 14:05 wasserth

Is it working when you add the argument --tracking_format tck to Tracking and Tractometry?

wasserth avatar Jun 08 '21 14:06 wasserth

Unfortunately no, at least if we do not comment out line 417.

mrsf23 avatar Jun 09 '21 15:06 mrsf23

Would it be possible to send one of the datasets where it is not working so I can try to reproduce the error? I would need peaks.nii.gz, dwi/sub-004_dwi_preproc.nii.gz, dwi/sub-004_dwi_preproc.bvals, dwi/sub-004_dwi_preproc.bvecs and dwi/sub-004_dwi_b0mask.nii.gz. You can find my email on the papers.

wasserth avatar Jun 10 '21 08:06 wasserth

I did some changes to the code. Could you try again with the newest master branch?

wasserth avatar Jun 21 '21 15:06 wasserth

One of our user ran into a similar issue with TractSeg 2.2.

+ Tractometry -i tractseg_output/TOM_trackings/ -o tractseg_output/Tractometry_peaks.csv -e tractseg_output/endings_segmentations/ -s tractseg_output/peaks.nii.gz --tracking_format tck --TOM tractseg_output/TOM --peak_length
Matplotlib created a temporary config/cache directory at /tmp/matplotlib-y_n6i7v7 because the default path (/.config/matplotlib) is not a writable directory; it is highly recommended to set the MPLCONFIGDIR environment variable to a writable directory, in particular to speed up the import of Matplotlib and to better support multiprocessing.

  0%|          | 0/50 [00:00<?, ?it/s]
  2%|▏         | 1/50 [00:00<00:36,  1.34it/s]
  2%|▏         | 1/50 [00:00<00:40,  1.21it/s]
Traceback (most recent call last):
  File "/usr/local/bin/Tractometry", line 138, in <module>
    main()
  File "/usr/local/bin/Tractometry", line 109, in main
    predicted_peaks=predicted_peaks, affine=scalar_image.affine)
  File "/usr/local/lib/python3.7/site-packages/tractseg/libs/tractometry.py", line 64, in evaluate_along_streamlines
    streamlines = _orient_to_same_start_region(streamlines, beginnings)
  File "/usr/local/lib/python3.7/site-packages/tractseg/libs/tractometry.py", line 42, in _orient_to_same_start_region
    if beginnings[int(startpoint[0]), int(startpoint[1]), int(startpoint[2])] == 0:
IndexError: index 66 is out of bounds for axis 0 with size 57

The image is very low res([ 4 57 73 83 125 1 1 ]) but affine is using RAS/neurological coordinate which is what MNI uses, right?

2021-09-09_11-26

I am trying to run it with the current master branch to see if it works out..

soichih avatar Sep 09 '21 15:09 soichih

It finished successfully with the current master branch!

image

Do you think I should try 2.3 branch also? Or is this fixed on latest master branch code?

soichih avatar Sep 09 '21 17:09 soichih

This issue is fixed only in the latest master branch. It is not yet fixed in 2.3. I am planning to do a few more changes and then I will tag a new version.

wasserth avatar Sep 10 '21 07:09 wasserth

Hello,

I am still seeing the same error message on the master branch ( IndexError: index 40 is out of bounds for axis 2 with size 40)

and I have no clue why it happing for some scanning but not others?

Any suggestions?

2021-11-15

Mano0olia avatar Nov 15 '21 15:11 Mano0olia

I'm running into the same "index X is out of bounds for axis 2 with size X". I managed to trace the error but maybe someone can elaborate on what's going on, or if I'm missing something.

  • Where can this error be traced back to (at least in my case)?
    • In the orient_to_same_start_region function under fiber_utils.py
  • Input variables:
    • streamlines, the array_sequence type that has an enumerable entry per streamline, each streamline with dimension (num_points_in_streamline, 3 coordinates). Loaded from the *.tck file for the tract under analysis.
    • beginnings, which is a 3D array of 0s/1s where 1 identifies voxels that belong to the 'beginning' segmentation of that tract denoted with [tract]_b.nii.gz.
  • What's happening? (as I understand it):
    • Line 468 of fiber_utils.py: for a streamline (sl), grab the first point coordinate. For me, this returns positional coordinates (e.g., "15, -0.4, 26") into the startpoint variable.
    • Subbing in the sample coordinates returned in startpoint, line 470 seems to evaluate: if beginnings[15, 0, 26] == 0. I think the idea here is to check whether the coordinates falls within the beginnings mask, but I'm not sure why it's using the startpoint coordinates as the indices of the 3D beginnings matrix - seems like there could be coordinates that are values > array index?
  • Additional notes:
    • I've overlaid the diffusion scan, tract (and ending) segmentations, and TOM tracking *.tck files in mrview and everything looks good and aligned.
    • Input data did have the same orientation as MNI space and is isotropic.
    • While manually debugging, I was fairly sure that the first coordinates of each streamline did indeed correspond to the *_b.nii.gz segmentation, so I tried commenting out the orient_to_same_start_region line. This fixed the 'index out of bounds' issue and the resulting figures look good.

Let me know if any of the above is confusing. Thanks for creating TractSeg - very neat tool!

j-tseng avatar Sep 29 '22 01:09 j-tseng

Thanks for this detailed analysis. This really helps. I pushed a change to the master branch which checks if the streamline startpoint is within the image dimensions of the beginnings mask. If it is not then the streamline is discarded. This should not happen very often. When it happens a warning is printed: "WARNING: steamline startpoint out of image. Discarding this streamline."

Could you check if this works for you and how often the warning is printed? It should only be printed a few times (meaning a few streamlines get discarded). If it is printed hundreds of times then something else is wrong.

wasserth avatar Sep 29 '22 09:09 wasserth

Hi @wasserth - glad the detailed analysis helps! I did a bit more digging and, unfortunately, I don't think your code addresses the root issue that I'm having:

  • The *.tck files generated by the Tracking call are storing the streamline point coordinates in a RAS+ space, not voxel space.
    • I ran the standard TractSeg/Tracking calls on a second participant to make sure this isn't a one-off.
    • Is this expected behaviour?
  • As a result, line 470 is using RAS+ coordinates (startpoint) to index the beginnings 3D binary mask image. This returns false when actually the coordinate does fall in the mask. Example from stepping through the left CST:
    • startpoint pulls the first point coordinate from a streamline as [-11.2, -4.2, 63.1].
    • Line 470 checks beginnings[-11, -4, 63], but beginnings is a 3D array with dim (x, y, z) = (122, 122, 70), so acceptable index values should be between [0, 121] or [0, 69].
    • Check returns false, but I open up the mask/images in fsleyes and can see that those coordinates do in fact fall in the mask.

I cobbled together a potential fix and opened a pull request - feel free to take a look and edit/adopt if you think it's okay! Gist of approach:

  1. plot_utils.py: feed inverted affine to orient_to_same_start_region as 3rd argument, since those matrices are not in existing input args.
  2. fiber_utils.py:
    1. Convert first and last point coordinates for each streamline to voxel space
    2. Check whether one of the ends is in the mask and flip/append accordingly
    3. Discard+warn if neither first nor last point coord is within the mask
    4. Verify that reoriented streamlines are all within the mask

Edit: Forgot to add a question: what's the +/- 0.5 shift for in the original code?

j-tseng avatar Sep 30 '22 16:09 j-tseng