Sharp edge near posterior of some brain masks

[psadil]

With the current version of MRIQC (22.0.1), some proportion of participants end up with brain masks that have a sharp boundary that excludes part of the occipital pole.

I haven't noticed a pattern in the scans that end up with masks that look like this. However, when I run synthstrip outside of MRIQC on the output of the pre_n4 node (which pulls docker image freesurfer/synthstrip:1.0)

synthstrip-docker -i clipped_corrected.nii.gz -o out.nii.gz -m mask.nii.gz -b 1

I get a mask that looks better. Here is a comparison, with the MRIQC mask in blue and the synthstrip-docker version in white.

This is how MRIQC is being called:

singularity exec -B $BIND:$BIND --cleanenv  \
  docker://nipreps/mriqc:22.0.1  mriqc  \
  $IN  $OUT  participant --participant-label 10357 \
  -w WORK --no-sub  --n_procs 10  --mem_gb 38 \
  --write-graph --modalities T1w --verbose-reports

Let me know if I should ask about this elsewhere, or if there are further suggestions for digging into this.

[oesteban]

SynthStrip inputs require some massaging, which are done with the freesurfer python bundle I didn't want to have in MRIQC.

I thought my implementation was similar to the original, but there must be some problem. It would be nice to look at the images that MRIQC and the SynthStrip container generate and pass along into the model.

Would you like to do that debugging?

cc/ @ahoopes

[ahoopes]

If it helps, I've updated the synthstrip code to no longer rely on the internal 'freesurfer' python package. Now it uses the 'surfa' library, which is available on pip and should be relatively lightweight. Might be worth updating the MRIQC implementation as well to keep things consistent.

[psadil]

If it helps, I've updated the synthstrip code to no longer rely on the internal 'freesurfer' python package. Now it uses the 'surfa' library, which is available on pip and should be relatively lightweight. Might be worth updating the MRIQC implementation as well to keep things consistent.

I'll let @oesteban , comment on whether depending on surfa is an option.

I have a bit of time now to try debugging. But when trying to run the model through either the MRIQC or current FS implementations, memory usage jumps way higher than what when using the SynthStrip container (i.e., to above 50GB). Is that expected? Are there settings in the SynthStrip container that keep memory usage lower?

environment.yml

``` name: mriqc channels: - conda-forge dependencies: - _libgcc_mutex=0.1=conda_forge - _openmp_mutex=4.5=2_kmp_llvm - alsa-lib=1.2.6.1=h7f98852_0 - atk-1.0=2.36.0=h3371d22_4 - attr=2.5.1=h166bdaf_0 - brotli=1.0.9=h166bdaf_7 - brotli-bin=1.0.9=h166bdaf_7 - brotlipy=0.7.0=py39hb9d737c_1004 - bzip2=1.0.8=h7f98852_4 - c-ares=1.18.1=h7f98852_0 - ca-certificates=2022.6.15=ha878542_0 - cached-property=1.5.2=hd8ed1ab_1 - cached_property=1.5.2=pyha770c72_1 - cairo=1.16.0=ha61ee94_1011 - certifi=2022.6.15=py39hf3d152e_0 - cffi=1.15.0=py39h4bc2ebd_0 - charset-normalizer=2.0.12=pyhd8ed1ab_0 - ci-info=0.2.0=pyh9f0ad1d_0 - click=8.1.3=py39hf3d152e_0 - cryptography=37.0.2=py39hd97740a_0 - cycler=0.11.0=pyhd8ed1ab_0 - dbus=1.13.6=h5008d03_3 - etelemetry=0.3.0=pyhd8ed1ab_0 - expat=2.4.8=h27087fc_0 - fftw=3.3.10=nompi_h77c792f_102 - filelock=3.7.1=pyhd8ed1ab_0 - font-ttf-dejavu-sans-mono=2.37=hab24e00_0 - font-ttf-inconsolata=3.000=h77eed37_0 - font-ttf-source-code-pro=2.038=h77eed37_0 - font-ttf-ubuntu=0.83=hab24e00_0 - fontconfig=2.14.0=h8e229c2_0 - fonts-conda-ecosystem=1=0 - fonts-conda-forge=1=0 - fonttools=4.33.3=py39hb9d737c_0 - freetype=2.10.4=h0708190_1 - fribidi=1.0.10=h36c2ea0_0 - future=0.18.2=py39hf3d152e_5 - gdk-pixbuf=2.42.8=hff1cb4f_0 - gettext=0.19.8.1=h73d1719_1008 - giflib=5.2.1=h36c2ea0_2 - glib=2.70.2=h780b84a_4 - glib-tools=2.70.2=h780b84a_4 - graphite2=1.3.13=h58526e2_1001 - graphviz=4.0.0=h5abf519_0 - gst-plugins-base=1.20.3=hf6a322e_0 - gstreamer=1.20.3=hd4edc92_0 - gtk2=2.24.33=h90689f9_2 - gts=0.7.6=h64030ff_2 - h5py=3.6.0=nompi_py39h7e08c79_100 - harfbuzz=4.3.0=hf9f4e7c_0 - hdf5=1.12.1=nompi_h2386368_104 - html5lib=1.1=pyh9f0ad1d_0 - icu=70.1=h27087fc_0 - idna=3.3=pyhd8ed1ab_0 - importlib-metadata=4.11.4=py39hf3d152e_0 - isodate=0.6.1=pyhd8ed1ab_0 - jack=1.9.18=h8c3723f_1002 - jinja2=3.0.3=pyhd8ed1ab_0 - joblib=1.1.0=pyhd8ed1ab_0 - jpeg=9e=h166bdaf_1 - keyutils=1.6.1=h166bdaf_0 - kiwisolver=1.4.3=py39hf939315_0 - krb5=1.19.3=h3790be6_0 - lcms2=2.12=hddcbb42_0 - ld_impl_linux-64=2.36.1=hea4e1c9_2 - lerc=3.0=h9c3ff4c_0 - libblas=3.9.0=15_linux64_mkl - libbrotlicommon=1.0.9=h166bdaf_7 - libbrotlidec=1.0.9=h166bdaf_7 - libbrotlienc=1.0.9=h166bdaf_7 - libcap=2.64=ha37c62d_0 - libcblas=3.9.0=15_linux64_mkl - libclang=14.0.5=default_h2e3cab8_0 - libclang13=14.0.5=default_h3a83d3e_0 - libcups=2.3.3=hf5a7f15_1 - libcurl=7.83.1=h7bff187_0 - libdb=6.2.32=h9c3ff4c_0 - libdeflate=1.12=h166bdaf_0 - libedit=3.1.20191231=he28a2e2_2 - libev=4.33=h516909a_1 - libevent=2.1.10=h9b69904_4 - libffi=3.4.2=h7f98852_5 - libflac=1.3.4=h27087fc_0 - libgcc-ng=12.1.0=h8d9b700_16 - libgd=2.3.3=h18fbbfe_3 - libgfortran-ng=12.1.0=h69a702a_16 - libgfortran5=12.1.0=hdcd56e2_16 - libglib=2.70.2=h174f98d_4 - libiconv=1.16=h516909a_0 - liblapack=3.9.0=15_linux64_mkl - libllvm14=14.0.6=he0ac6c6_0 - libnghttp2=1.47.0=h727a467_0 - libnsl=2.0.0=h7f98852_0 - libogg=1.3.4=h7f98852_1 - libopus=1.3.1=h7f98852_1 - libpng=1.6.37=h21135ba_2 - libpq=14.4=hd77ab85_0 - libprotobuf=3.19.4=h780b84a_0 - librsvg=2.54.3=h7abd40a_0 - libsndfile=1.0.31=h9c3ff4c_1 - libssh2=1.10.0=ha56f1ee_2 - libstdcxx-ng=12.1.0=ha89aaad_16 - libtiff=4.4.0=hc85c160_1 - libtool=2.4.6=h9c3ff4c_1008 - libudev1=249=h166bdaf_4 - libuuid=2.32.1=h7f98852_1000 - libvorbis=1.3.7=h9c3ff4c_0 - libwebp=1.2.2=h3452ae3_0 - libwebp-base=1.2.2=h7f98852_1 - libxcb=1.13=h7f98852_1004 - libxkbcommon=1.0.3=he3ba5ed_0 - libxml2=2.9.14=h22db469_0 - libxslt=1.1.35=h8affb1d_0 - libzlib=1.2.12=h166bdaf_1 - llvm-openmp=14.0.4=he0ac6c6_0 - looseversion=1.0.1=pyhd8ed1ab_0 - lxml=4.9.0=py39hb9d737c_0 - lz4-c=1.9.3=h9c3ff4c_1 - markupsafe=2.1.1=py39hb9d737c_1 - matplotlib=3.5.2=py39hf3d152e_0 - matplotlib-base=3.5.2=py39h700656a_0 - mkl=2022.1.0=h84fe81f_915 - munkres=1.1.4=pyh9f0ad1d_0 - mysql-common=8.0.29=haf5c9bc_1 - mysql-libs=8.0.29=h28c427c_1 - ncurses=6.3=h27087fc_1 - networkx=2.8.4=pyhd8ed1ab_0 - nibabel=4.0.1=pyhd8ed1ab_0 - nilearn=0.9.1=pyhd8ed1ab_0 - ninja=1.11.0=h924138e_0 - nipype=1.8.2=py39hf3d152e_0 - nspr=4.32=h9c3ff4c_1 - nss=3.78=h2350873_0 - numpy=1.23.0=py39hba7629e_0 - openjpeg=2.4.0=hb52868f_1 - openssl=1.1.1p=h166bdaf_0 - packaging=21.3=pyhd8ed1ab_0 - pandas=1.4.3=py39h1832856_0 - pango=1.50.7=hbd2fdc8_0 - pcre=8.45=h9c3ff4c_0 - pillow=9.1.1=py39hae2aec6_1 - pip=22.1.2=pyhd8ed1ab_0 - pixman=0.40.0=h36c2ea0_0 - portaudio=19.6.0=h57a0ea0_5 - prov=2.0.0=pyhd3deb0d_0 - psutil=5.9.1=py39hb9d737c_0 - pthread-stubs=0.4=h36c2ea0_1001 - pulseaudio=14.0=h7f54b18_8 - pycparser=2.21=pyhd8ed1ab_0 - pydicom=2.3.0=pyh6c4a22f_0 - pydot=1.4.2=py39hf3d152e_2 - pyopenssl=22.0.0=pyhd8ed1ab_0 - pyparsing=3.0.9=pyhd8ed1ab_0 - pyqt=5.15.4=py39h18e9c17_1 - pyqt5-sip=12.9.0=py39h5a03fae_1 - pysocks=1.7.1=py39hf3d152e_5 - python=3.9.13=h9a8a25e_0_cpython - python-dateutil=2.8.2=pyhd8ed1ab_0 - python_abi=3.9=2_cp39 - pytorch=1.10.2=cpu_py39h5e9ed0b_1 - pytz=2022.1=pyhd8ed1ab_0 - qt-main=5.15.4=ha5833f6_2 - rdflib=6.1.1=pyhd8ed1ab_0 - readline=8.1.2=h0f457ee_0 - requests=2.28.0=pyhd8ed1ab_1 - scikit-learn=1.1.1=py39h4037b75_0 - scipy=1.8.1=py39he49c0e8_0 - setuptools=59.5.0=py39hf3d152e_0 - simplejson=3.17.6=py39hb9d737c_1 - sip=6.5.1=py39he80948d_2 - six=1.16.0=pyh6c4a22f_0 - sleef=3.5.1=h9b69904_2 - sqlite=3.38.5=h4ff8645_0 - tbb=2021.5.0=h924138e_1 - threadpoolctl=3.1.0=pyh8a188c0_0 - tk=8.6.12=h27826a3_0 - toml=0.10.2=pyhd8ed1ab_0 - tornado=6.1=py39hb9d737c_3 - traits=6.3.2=py39hb9d737c_1 - typing_extensions=4.2.0=pyha770c72_1 - tzdata=2022a=h191b570_0 - unicodedata2=14.0.0=py39hb9d737c_1 - urllib3=1.26.9=pyhd8ed1ab_0 - webencodings=0.5.1=py_1 - wheel=0.37.1=pyhd8ed1ab_0 - xcb-util=0.4.0=h166bdaf_0 - xcb-util-image=0.4.0=h166bdaf_0 - xcb-util-keysyms=0.4.0=h166bdaf_0 - xcb-util-renderutil=0.3.9=h166bdaf_0 - xcb-util-wm=0.4.1=h166bdaf_0 - xorg-kbproto=1.0.7=h7f98852_1002 - xorg-libice=1.0.10=h7f98852_0 - xorg-libsm=1.2.3=hd9c2040_1000 - xorg-libx11=1.7.2=h7f98852_0 - xorg-libxau=1.0.9=h7f98852_0 - xorg-libxdmcp=1.1.3=h7f98852_0 - xorg-libxext=1.3.4=h7f98852_1 - xorg-libxrender=0.9.10=h7f98852_1003 - xorg-renderproto=0.11.1=h7f98852_1002 - xorg-xextproto=7.3.0=h7f98852_1002 - xorg-xproto=7.0.31=h7f98852_1007 - xvfbwrapper=0.2.9=pyhd8ed1ab_1005 - xz=5.2.5=h516909a_1 - zipp=3.8.0=pyhd8ed1ab_0 - zlib=1.2.12=h166bdaf_1 - zstd=1.5.2=h8a70e8d_1 - pip: - cython==0.29.30 - nitransforms==22.0.0 - surfa==0.0.11 ```

[ahoopes]

PyTorch tends to use a lot more memory than it necessarily needs to (if space is available). Docker has a default setting that might limit container memory on your machine - maybe it's related to that? I have not set anything specific in the SynthStrip configuration.

[psadil]

PyTorch tends to use a lot more memory than it necessarily needs to (if space is available). Docker has a default setting that might limit container memory on your machine - maybe it's related to that? I have not set anything specific in the SynthStrip configuration.

Something to do with Docker seems plausible, though I'm not sure. I suppose that's not critical. FWIW, bringing in the cpuonly package results in memory usage comparable to when the container is run.

@oesteban, I think there might be an issue in how mriqc prepares the input for SynthStrip. A clue was that, when testing on an image with shape (256, 256, 176), the variable input_tensor as given to the model ends up with shape (1,1,256,256,192). At the analogous step in mri_synthstrip, the shape is (1,1,256,192,256). I think this difference comes from how the target shape is calculated. For example, mriqc uses the coordinates of the corner voxels e.g.,, but I think that result is in the input space when it should be in the output--accounting for rotations to get there.

The image that I'm working with has affine

nb.load('test.nii.gz').affine
array([[   1. ,    0. ,    0. , -122. ],
       [   0. ,    1. ,    0. , -122. ],
       [   0. ,    0. ,    1. , -102.5],
       [   0. ,    0. ,    0. ,    1. ]])

Conformation swaps the second and third axes. The original z-axis was rather narrow, and so the brain in the resulting image gets cropped. I guess that's the source of the sharp edge.

Perhaps it would work to use the target_affine when calculating the corners in mm?

    # ...calculate target_affine above here...
    # Get corner voxel centers in mm
    corners_xyz = np.abs(
        affine
        @ target_affine
        @ np.hstack((corner_centers_ijk, np.ones((len(corner_centers_ijk), 1)))).T
    )

Using that gives a better mask, see below. As above, the mask from synthstrip-docker is in white and the one from the (updated) script is in blue.

[psadil]

Hi Maintainers,

Following up on this. Did that patch sound feasible? Or, how about adding a dependence on surfa, following the official SynthStrip implementation? If either of those sound okay, I'd be happy to try submitting a pull request.

[oesteban]

Hi @psadil - sorry you caught me on vacation. Yes, your patch sounds reasonable. We'd love to see your PR adding it. Let's fix this bug first and then consider the dependency on surfa, if we believe that will be the best option.

[oesteban]

Actually, I'm not completely sure of your fix, the target_affine has an offset which happens to work with you image but may not work anywhere else. Could you share this particular image with me?

[psadil]

I can't share that image, but I may be able to share another that shows the same issue. If this is the problem, then it might also be that the sharp edge can be elicited by manually clipping the z-axis. Let me look in to this, and I'll get back to you in the next few days.

[oesteban]

Perhaps it would be better if you could just report the original shape and affine of that failing image. I would work from there.

[psadil]

Sure, they are

dcmmeta_shape:  [256, 256, 176]
dcmmeta_affine: [
  [1.0, 0.0, 0.0, -122.0], 
  [0.0, 1.0, 0.0, -117.0], 
  [0.0, 0.0, 1.0, -102.5], 
  [0.0, 0.0, 0.0, 1.0]
]

[oesteban]

Yes, it seems our reorientation is not producing the right ordering at the output.

[psadil]

Just to confirm: was it still going to be helpful to see another example? It's my impression that the reorientation you mention is the cause of the sharp masks (when combined with a relatively small axis)

[oesteban]

Thanks for your patience, really hard for me to focus on MRIQC these days.

Yes, I'm pretty sure this is a result of the wrong reorientation happening. Let me first try to figure it out without more data. Thanks!

[psadil]

Thanks for your patience, really hard for me to focus on MRIQC these days.

Thanks for the fix, and in general for maintaining this package!