synthstrip fails with large borders

[neurolabusc]

Recent commits to this repository attempt to fix an issue with the current stable release (7.4.1) when mri_synthstrip is run with a large border (-b). However the latest solution still fails with many images even though the results look fine for a lower b value. It is easy to replicate with the images from the CQ500 dataset (e.g. subject 289). To elicit this bug you would run:

mri_synthstrip -i sub-289_CT.nii.gz -o 6.nii.gz -b 6

and the masking looks wrong:

The issue can be fixed by fixing the edt code from:

    if args.border >= max_dist:
        print(f'specified border {args.border} greater than max dtrans {max_dist} - computing sdt')
        dif = args.border - (max_dist-1)
        mask1 = (sdt >= (max_dist-1))  # region that original sdt has real distances
        dtrans = scipy.ndimage.morphology.distance_transform_edt(mask1) + (max_dist-2)
        sdt = dtrans  # negative interior distances are lost, but doesn't matter

to read:

    if args.border >= max_dist:
        print(f'specified border {args.border} greater than max dtrans {max_dist} - computing sdt')
        mask1 = np.where(sdt > 0.0, 1, 0)
        sdt = scipy.ndimage.distance_transform_edt(mask1)

with this change, the images look fine:

as an aside, one could simply always use the Python edt method instead of using the tensor compute for small distances. This would simplify the code, and also allow non-integer values for the border. This would generalize to sub-millimeter conformed images. While one assumes that EDT should be computationally expensive, the Python implementation is extremely fast, leveraging the unexpected and elegant fact that the EDT is separable.

I do think supporting large borders is useful. The robustness of synthstrip allows it to be used across modalities to crop images - both to remove facial features and speed us subsequent processing. A larger border aids tools that compute intensity homogeneity correction and prefer neighborhood information. Here is the proposed code with a 18mm border:

[mu40]

Thanks for spotting this issue and the proposed solution, @neurolabusc! We'll look into it in a couple of days.

[mu40]

I'm afraid I couldn't reproduce the issue for any of the subject-289 scans from CQ500 with the latest SynthStrip script in this repository. Is it possible that you obtained the above results with version 1.5 from the Docker Hub, and could you share the image file you used?

The web page on the EDT provided an interesting read. While I do like the idea of always extending the narrow-band SDT estimated by the model, the SciPy implementation is a little slow for 256-by-256-by-256 images, for example, and it would need to run twice in the general case. Until we get around to perhaps adding the marching parabolas algorithm to surfa, recomputing the EDT only when needed appears more efficient.

[buildqa]

Please re-post this to the freesurfer help list if you still have an issue, but closing the issue here.