poldracklab / pydeface

defacing utility for MRI images
MIT License
109 stars 42 forks source link

defacing issue crop image #27

Closed alexfoias closed 4 years ago

alexfoias commented 5 years ago

I tried to deface some T1w images (brain and spine up to C7 level). The problem is that the images are cropped to C2 level in the spine.

Is there a way to enlarge the FOV of the pydeface ?

muschellij2 commented 4 years ago

Can you provide an MWE https://stackoverflow.com/help/minimal-reproducible-example ? Also, you can either find the mask by looking at voxels that are different from the non-masked image or you can dig into the temporary files (currently) that are generated https://github.com/poldracklab/pydeface/blob/master/pydeface/utils.py#L116. You can then cut the mask to where it needs to stop (aka the neck).

muschellij2 commented 4 years ago

Looking a bit more I'd recommend setting --nocleanup so that the mask is output from the operation so you can inspect the mask and you may be able to manipulate these or figure out some rules from these.

alexfoias commented 4 years ago

@muschellij2 I followed your advice and I saved the mask. The mask doesn't include the complete VOI. Do I have to manually modify the mask to include the complete neck region or is there a parameter that I could change to expand the coverage?

muschellij2 commented 4 years ago

Again, almost impossible to diagnose without an MWE https://stackoverflow.com/help/minimal-reproducible-example ?

alexfoias commented 4 years ago

Here is an example image : https://openneuro.org/datasets/ds002393/versions/1.0.2/file-display/sub-chiba750:anat:sub-chiba750_T1w.nii.gz

muschellij2 commented 4 years ago

Is this what you're looking for? This relies on my implementation of pydeface in the R package fslr, which is almost identical. The extension to the rest of the image

library(RNifti)
library(fslr)
#> Loading required package: oro.nifti
#> oro.nifti 0.10.1
#> 
#> Attaching package: 'oro.nifti'
#> The following object is masked from 'package:RNifti':
#> 
#>     origin
#> Loading required package: neurobase
fname = "sub-chiba750_anat_sub-chiba750_T1w.nii.gz"
download.file(url = paste0("https://openneuro.org/crn/datasets/ds002393/", 
                           "snapshots/1.0.2/files/", 
                           "sub-chiba750:anat:sub-chiba750_T1w.nii.gz"),
              destfile = fname)
img = RNifti::readNifti(fname)
out = face_removal_mask(fname)
#> Warning in get.fsl(): Setting fsl.path to /usr/local/fsl
#> Warning in get.fsloutput(): Can't find FSLOUTPUTTYPE, setting to NIFTI_GZ
#> FSLDIR='/usr/local/fsl'; PATH=${FSLDIR}/bin:${PATH};export PATH FSLDIR; sh "${FSLDIR}/etc/fslconf/fsl.sh"; FSLOUTPUTTYPE=NIFTI_GZ; export FSLOUTPUTTYPE; ${FSLDIR}/bin/flirt -in "/usr/local/fsl/data/standard/MNI152_T1_1mm.nii.gz" -ref "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpK4U4ug/reprex17ffa74bf4484/sub-chiba750_anat_sub-chiba750_T1w.nii.gz" -out "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpKjcGt4/file58d7925ddcf" -dof 12 -omat "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpKjcGt4/file58d19f2c8d.mat" -cost mutualinfo
#> FSLDIR='/usr/local/fsl'; PATH=${FSLDIR}/bin:${PATH};export PATH FSLDIR; sh "${FSLDIR}/etc/fslconf/fsl.sh"; FSLOUTPUTTYPE=NIFTI_GZ; export FSLOUTPUTTYPE; ${FSLDIR}/bin/flirt -in "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpKjcGt4/file58dbe034b1.nii.gz" -ref "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpK4U4ug/reprex17ffa74bf4484/sub-chiba750_anat_sub-chiba750_T1w.nii.gz" -applyxfm -init "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpKjcGt4/file58d19f2c8d.mat" -interp nearestneighbour -out "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpKjcGt4/file58d25057f2c"
out_img = readnii(out)
ortho2(img)

Here is my extension. The idea is to find all the indices in x, y, and z-planes that are in the mask. We then extend down in the z-plane (dim3) and "forward" in the y-plane (dim2) and show the masked output:


xx = which(out_img == 0, arr.ind = TRUE)
dim3 = 1:(min(xx[,3]) + floor(  (max(xx[,3]) - min(xx[,3]))/2))
dim2 = min(xx[,2]):dim(out_img)[2]
# eg = expand.grid(dim1 = unique(xx[,1]), 
                 # dim2 = unique(xx[,2]), dim3 = dim3)
eg = expand.grid(dim1 = unique(xx[,1]),
                 dim2 = dim2, dim3 = dim3)
eg = as.matrix(eg)
new_mask = out_img
new_mask[eg] = 0
final_img = img * new_mask
ortho2(final_img)

Created on 2020-01-23 by the reprex package (v0.3.0.9001)

Session info ``` r sessioninfo::session_info() #> ─ Session info ─────────────────────────────────────────────────────────────── #> setting value #> version R version 3.6.0 (2019-04-26) #> os macOS Mojave 10.14.6 #> system x86_64, darwin15.6.0 #> ui X11 #> language (EN) #> collate en_US.UTF-8 #> ctype en_US.UTF-8 #> tz America/New_York #> date 2020-01-23 #> #> ─ Packages ─────────────────────────────────────────────────────────────────── #> package * version date lib source #> abind 1.4-5 2016-07-21 [1] CRAN (R 3.6.0) #> assertthat 0.2.1 2019-03-21 [1] CRAN (R 3.6.0) #> backports 1.1.5 2019-10-02 [1] CRAN (R 3.6.0) #> bitops 1.0-6 2013-08-17 [1] CRAN (R 3.6.0) #> cli 2.0.1 2020-01-08 [1] CRAN (R 3.6.0) #> crayon 1.3.4 2017-09-16 [1] CRAN (R 3.6.0) #> curl 4.3 2019-12-02 [1] CRAN (R 3.6.0) #> digest 0.6.23 2019-11-23 [1] CRAN (R 3.6.0) #> evaluate 0.14 2019-05-28 [1] CRAN (R 3.6.0) #> fansi 0.4.1 2020-01-08 [1] CRAN (R 3.6.0) #> fs 1.3.1 2019-05-06 [1] CRAN (R 3.6.0) #> fslr * 2.24.4 2019-12-03 [1] local #> glue 1.3.1 2019-03-12 [1] CRAN (R 3.6.0) #> highr 0.8 2019-03-20 [1] CRAN (R 3.6.0) #> htmltools 0.4.0 2019-10-04 [1] CRAN (R 3.6.0) #> httr 1.4.1 2019-08-05 [1] CRAN (R 3.6.0) #> knitr 1.26.1 2020-01-05 [1] Github (muschellij2/knitr@f5af631) #> magrittr 1.5 2014-11-22 [1] CRAN (R 3.6.0) #> matrixStats 0.55.0 2019-09-07 [1] CRAN (R 3.6.0) #> mime 0.8 2019-12-19 [1] CRAN (R 3.6.0) #> neurobase * 1.29.1 2020-01-14 [1] local #> oro.nifti * 0.10.1 2019-12-14 [1] local #> pillar 1.4.3 2019-12-20 [1] CRAN (R 3.6.0) #> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 3.6.0) #> purrr 0.3.3 2019-10-18 [1] CRAN (R 3.6.0) #> R.methodsS3 1.7.1 2016-02-16 [1] CRAN (R 3.6.0) #> R.oo 1.23.0 2019-11-03 [1] CRAN (R 3.6.0) #> R.utils 2.9.0 2019-06-13 [1] CRAN (R 3.6.0) #> R6 2.4.1 2019-11-12 [1] CRAN (R 3.6.0) #> Rcpp 1.0.3 2019-11-08 [1] CRAN (R 3.6.0) #> reprex 0.3.0.9001 2020-01-05 [1] Github (tidyverse/reprex@5ae0b29) #> rlang 0.4.2 2019-11-23 [1] CRAN (R 3.6.0) #> rmarkdown 2.0.7 2020-01-17 [1] Github (rstudio/rmarkdown@2faf16a) #> RNifti * 1.0.1 2019-11-27 [1] CRAN (R 3.6.0) #> rstudioapi 0.10.0-9003 2020-01-05 [1] Github (rstudio/rstudioapi@abe596d) #> sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.6.0) #> stringi 1.4.5 2020-01-11 [1] CRAN (R 3.6.0) #> stringr 1.4.0 2019-02-10 [1] CRAN (R 3.6.0) #> styler 1.1.1 2019-05-06 [1] CRAN (R 3.6.0) #> tibble 2.1.3 2019-06-06 [1] CRAN (R 3.6.0) #> withr 2.1.2 2018-03-15 [1] CRAN (R 3.6.0) #> xfun 0.11 2019-11-12 [1] CRAN (R 3.6.0) #> xml2 1.2.2 2019-08-09 [1] CRAN (R 3.6.0) #> yaml 2.2.0 2018-07-25 [1] CRAN (R 3.6.0) #> #> [1] /Library/Frameworks/R.framework/Versions/3.6/Resources/library ```
muschellij2 commented 4 years ago

This may be more clear as to what is getting added.

library(RNifti)
library(fslr)
#> Loading required package: oro.nifti
#> oro.nifti 0.10.1
#> 
#> Attaching package: 'oro.nifti'
#> The following object is masked from 'package:RNifti':
#> 
#>     origin
#> Loading required package: neurobase
fname = "sub-chiba750_anat_sub-chiba750_T1w.nii.gz"
download.file(url = paste0("https://openneuro.org/crn/datasets/ds002393/", 
                           "snapshots/1.0.2/files/", 
                           "sub-chiba750:anat:sub-chiba750_T1w.nii.gz"),
              destfile = fname)
img = RNifti::readNifti(fname)
out = face_removal_mask(fname)
#> Warning in get.fsl(): Setting fsl.path to /usr/local/fsl
#> Warning in get.fsloutput(): Can't find FSLOUTPUTTYPE, setting to NIFTI_GZ
#> FSLDIR='/usr/local/fsl'; PATH=${FSLDIR}/bin:${PATH};export PATH FSLDIR; sh "${FSLDIR}/etc/fslconf/fsl.sh"; FSLOUTPUTTYPE=NIFTI_GZ; export FSLOUTPUTTYPE; ${FSLDIR}/bin/flirt -in "/usr/local/fsl/data/standard/MNI152_T1_1mm.nii.gz" -ref "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpK4U4ug/reprex17ffa365ce2a8/sub-chiba750_anat_sub-chiba750_T1w.nii.gz" -out "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpwhfPxR/file5f0382ef3d6" -dof 12 -omat "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpwhfPxR/file5f042b7e932.mat" -cost mutualinfo
#> FSLDIR='/usr/local/fsl'; PATH=${FSLDIR}/bin:${PATH};export PATH FSLDIR; sh "${FSLDIR}/etc/fslconf/fsl.sh"; FSLOUTPUTTYPE=NIFTI_GZ; export FSLOUTPUTTYPE; ${FSLDIR}/bin/flirt -in "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpwhfPxR/file5f043b9ef4f.nii.gz" -ref "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpK4U4ug/reprex17ffa365ce2a8/sub-chiba750_anat_sub-chiba750_T1w.nii.gz" -applyxfm -init "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpwhfPxR/file5f042b7e932.mat" -interp nearestneighbour -out "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpwhfPxR/file5f0128a836b"
out_img = readnii(out)
ortho2(img)

ortho2(img, out_img)


xx = which(out_img == 0, arr.ind = TRUE)
dim3 = 1:(min(xx[,3]) + floor(  (max(xx[,3]) - min(xx[,3]))/2))
dim2 = min(xx[,2]):dim(out_img)[2]
# eg = expand.grid(dim1 = unique(xx[,1]), 
                 # dim2 = unique(xx[,2]), dim3 = dim3)
eg = expand.grid(dim1 = unique(xx[,1]),
                 dim2 = dim2, dim3 = dim3)
eg = as.matrix(eg)
new_mask = out_img
new_mask[eg] = 0
final_img = img * new_mask
ortho2(final_img)

ortho2(img, new_mask == 0)

Created on 2020-01-23 by the reprex package (v0.3.0.9001)

Session info ``` r sessioninfo::session_info() #> ─ Session info ─────────────────────────────────────────────────────────────── #> setting value #> version R version 3.6.0 (2019-04-26) #> os macOS Mojave 10.14.6 #> system x86_64, darwin15.6.0 #> ui X11 #> language (EN) #> collate en_US.UTF-8 #> ctype en_US.UTF-8 #> tz America/New_York #> date 2020-01-23 #> #> ─ Packages ─────────────────────────────────────────────────────────────────── #> package * version date lib source #> abind 1.4-5 2016-07-21 [1] CRAN (R 3.6.0) #> assertthat 0.2.1 2019-03-21 [1] CRAN (R 3.6.0) #> backports 1.1.5 2019-10-02 [1] CRAN (R 3.6.0) #> bitops 1.0-6 2013-08-17 [1] CRAN (R 3.6.0) #> cli 2.0.1 2020-01-08 [1] CRAN (R 3.6.0) #> crayon 1.3.4 2017-09-16 [1] CRAN (R 3.6.0) #> curl 4.3 2019-12-02 [1] CRAN (R 3.6.0) #> digest 0.6.23 2019-11-23 [1] CRAN (R 3.6.0) #> evaluate 0.14 2019-05-28 [1] CRAN (R 3.6.0) #> fansi 0.4.1 2020-01-08 [1] CRAN (R 3.6.0) #> fs 1.3.1 2019-05-06 [1] CRAN (R 3.6.0) #> fslr * 2.24.4 2019-12-03 [1] local #> glue 1.3.1 2019-03-12 [1] CRAN (R 3.6.0) #> highr 0.8 2019-03-20 [1] CRAN (R 3.6.0) #> htmltools 0.4.0 2019-10-04 [1] CRAN (R 3.6.0) #> httr 1.4.1 2019-08-05 [1] CRAN (R 3.6.0) #> knitr 1.26.1 2020-01-05 [1] Github (muschellij2/knitr@f5af631) #> magrittr 1.5 2014-11-22 [1] CRAN (R 3.6.0) #> matrixStats 0.55.0 2019-09-07 [1] CRAN (R 3.6.0) #> mime 0.8 2019-12-19 [1] CRAN (R 3.6.0) #> neurobase * 1.29.1 2020-01-14 [1] local #> oro.nifti * 0.10.1 2019-12-14 [1] local #> pillar 1.4.3 2019-12-20 [1] CRAN (R 3.6.0) #> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 3.6.0) #> purrr 0.3.3 2019-10-18 [1] CRAN (R 3.6.0) #> R.methodsS3 1.7.1 2016-02-16 [1] CRAN (R 3.6.0) #> R.oo 1.23.0 2019-11-03 [1] CRAN (R 3.6.0) #> R.utils 2.9.0 2019-06-13 [1] CRAN (R 3.6.0) #> R6 2.4.1 2019-11-12 [1] CRAN (R 3.6.0) #> Rcpp 1.0.3 2019-11-08 [1] CRAN (R 3.6.0) #> reprex 0.3.0.9001 2020-01-05 [1] Github (tidyverse/reprex@5ae0b29) #> rlang 0.4.2 2019-11-23 [1] CRAN (R 3.6.0) #> rmarkdown 2.0.7 2020-01-17 [1] Github (rstudio/rmarkdown@2faf16a) #> RNifti * 1.0.1 2019-11-27 [1] CRAN (R 3.6.0) #> rstudioapi 0.10.0-9003 2020-01-05 [1] Github (rstudio/rstudioapi@abe596d) #> sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.6.0) #> stringi 1.4.5 2020-01-11 [1] CRAN (R 3.6.0) #> stringr 1.4.0 2019-02-10 [1] CRAN (R 3.6.0) #> styler 1.1.1 2019-05-06 [1] CRAN (R 3.6.0) #> tibble 2.1.3 2019-06-06 [1] CRAN (R 3.6.0) #> withr 2.1.2 2018-03-15 [1] CRAN (R 3.6.0) #> xfun 0.11 2019-11-12 [1] CRAN (R 3.6.0) #> xml2 1.2.2 2019-08-09 [1] CRAN (R 3.6.0) #> yaml 2.2.0 2018-07-25 [1] CRAN (R 3.6.0) #> #> [1] /Library/Frameworks/R.framework/Versions/3.6/Resources/library ```
alexfoias commented 4 years ago

@muschellij2 Thanks for your help. I tried to follow the instructions above, but I ran into trouble .

``` > library(RNifti) > library(fslr) Loading required package: oro.nifti oro.nifti 0.9.1 Attaching package: ‘oro.nifti’ The following objects are masked from ‘package:RNifti’: origin, pixdim, pixdim<- Loading required package: neurobase Please cite the fslr package using: Jenkinson M, Beckmann CF, Behrens TE, Woolrich MW, Smith SM (2012). “FSL.” _Neuroimage_, *62*(2), 782-790. Muschelli J, Sweeney E, Lindquist M, Crainiceanu C (2015). “fslr: Connecting the FSL Software with R.” _The R Journal_, *7*(1), 163-175. > fname = "sub-chiba750_anat_sub-chiba750_T1w.nii.gz" > download.file(url = paste0("https://openneuro.org/crn/datasets/ds002393/", + "snapshots/1.0.2/files/", + "sub-chiba750:anat:sub-chiba750_T1w.nii.gz"), + destfile = fname) trying URL 'https://openneuro.org/crn/datasets/ds002393/snapshots/1.0.2/files/sub-chiba750:anat:sub-chiba750_T1w.nii.gz' downloaded 51.8 MB > img = RNifti::readNifti(fname) > ortho2(img) Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘pixdim’ for signature ‘"niftiImage"’ > sessioninfo::session_info() ─ Session info ───────────────────────────────────────────────────────────────────────────────────── setting value version R version 3.6.2 (2019-12-12) os macOS Mojave 10.14.6 system x86_64, darwin15.6.0 ui RStudio language (EN) collate en_CA.UTF-8 ctype en_CA.UTF-8 tz America/Montreal date 2020-01-24 ─ Packages ───────────────────────────────────────────────────────────────────────────────────────── package * version date lib source abind 1.4-5 2016-07-21 [1] CRAN (R 3.6.0) assertthat 0.2.1 2019-03-21 [1] CRAN (R 3.6.0) bitops 1.0-6 2013-08-17 [1] CRAN (R 3.6.0) cli 2.0.1 2020-01-08 [1] CRAN (R 3.6.0) crayon 1.3.4 2017-09-16 [1] CRAN (R 3.6.0) fansi 0.4.1 2020-01-08 [1] CRAN (R 3.6.0) fslr * 2.24.1 2019-08-05 [1] CRAN (R 3.6.0) glue 1.3.1 2019-03-12 [1] CRAN (R 3.6.0) matrixStats 0.55.0 2019-09-07 [1] CRAN (R 3.6.0) neurobase * 1.29.0 2020-01-14 [1] CRAN (R 3.6.0) oro.nifti * 0.9.1 2017-10-26 [1] CRAN (R 3.6.0) R.methodsS3 1.7.1 2016-02-16 [1] CRAN (R 3.6.0) R.oo 1.23.0 2019-11-03 [1] CRAN (R 3.6.0) R.utils 2.9.2 2019-12-08 [1] CRAN (R 3.6.0) Rcpp 1.0.3 2019-11-08 [1] CRAN (R 3.6.0) RNifti * 1.0.1 2019-11-27 [1] CRAN (R 3.6.0) sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.6.0) withr 2.1.2 2018-03-15 [1] CRAN (R 3.6.0) [1] /Library/Frameworks/R.framework/Versions/3.6/Resources/library ```

I noticed that I have a different version of fslr & oro.nifti than you.

muschellij2 commented 4 years ago

Yeah, I'm trying to update the packages on my end, I thikn oyu need to install:

remotes::install_github("muschellij2/oro.nifti")
alexfoias commented 4 years ago

Thanks. I managed to reproduce the results for the give case. I'll try to create a script to do a batch processing for the complete dataset.

alexfoias commented 4 years ago

@muschellij2 I tried to run the script for T2w, but I ran into some issues. Here is the output: https://gist.github.com/alexfoias/53b5a3b80502737edabb68e8e33d44b9

muschellij2 commented 4 years ago

So the overall process works on the fact you can register a brain to a brain, or really a head to a head. That image you are downloading really isn't a brain/head, but more like half a brain a lot of the spinal cord. Not only that, the image is very narrow and not a full head/next. The affine transformation won't handle this well and you need to use a different approach. Best, John

On Wed, Feb 12, 2020 at 11:56 AM alexfoias notifications@github.com wrote:

@muschellij2 https://github.com/muschellij2 I tried to run the script for T2w, but I ran into some issues. Here is the output: https://gist.github.com/alexfoias/53b5a3b80502737edabb68e8e33d44b9

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/poldracklab/pydeface/issues/27?email_source=notifications&email_token=AAIGPLVSASGG2JTJCIAS5ATRCQS5JA5CNFSM4HKWF56KYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOELRQ4XY#issuecomment-585305695, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAIGPLSFON7OE2S4PQRKAI3RCQS5JANCNFSM4HKWF56A .

muschellij2 commented 4 years ago

Also, someone from @poldrack @poldracklab probably should weigh in.

ofgulban commented 4 years ago

Agree with @muschellij2 's previous comment. Defacing does not work by default because your coverage is not a conventional brain image. I would recommend using --template and --facemask flags in addition to setup a custom spine-imaging-focused face mask. For instance, if you have 10 subjects, take one of the subjects, draw manually a face mask and use the anatomical image of this subject with the --template flag, and the face mask with the --facemask flag when defacing other subjects. Your pydeface command should look sth like this at the end:

pydeface /path/to/sub-02.nii --template /path/to/sub-01.nii --facemask /path/to/sub-01_facemask.nii
pydeface /path/to/sub-03.nii --template /path/to/sub-01.nii --facemask /path/to/sub-01_facemask.nii

Here --facemask /path/to/sub-01_facemask.nii should be created by you manually using sth like ITK-SNAP or fsleyes/fslview.

(Note: I am not from @poldrack 's lab, but I put these extra flags to the command line interface to cover for such use cases.)

ofgulban commented 4 years ago

@alexfoias with regards to T2w image defacing, see the --applyto flag sections of this pull request: #15
If you have higher signal-to-noise ratio T1w images together with T2w images, you might want to consider using this --applyto flag.

muschellij2 commented 4 years ago

You can make some adaptations, but you really need to see the above comments, but here's something that works for this one file (need to extend mask):

library(RNifti)
library(fslr)
#> Loading required package: oro.nifti
#> oro.nifti 0.10.2
#> 
#> Attaching package: 'oro.nifti'
#> The following object is masked from 'package:RNifti':
#> 
#>     origin
#> Loading required package: neurobase
library(extrantsr)
#> 
#> Attaching package: 'extrantsr'
#> The following object is masked from 'package:neurobase':
#> 
#>     zero_pad
#> The following objects are masked from 'package:oro.nifti':
#> 
#>     origin, origin<-
#> The following object is masked from 'package:RNifti':
#> 
#>     origin
fname = file.path(tempdir(), 
                  "sub-chibaIngeniasub-chibaIngenia_T2w.nii.gz")
if (!file.exists(fname)) {
  download.file(url = paste0(
    "https://openneuro.org/crn/datasets/ds002393/", 
    "snapshots/1.0.2/files/", 
    "sub-chibaIngenia:anat:sub-chibaIngenia_T2w.nii.gz"),
    destfile = fname)
}
img = RNifti::readNifti(fname)
tfname = mni_fname(mm = "1")
#> Warning in get.fsl(): Setting fsl.path to /usr/local/fsl
#> Warning in get.fsloutput(): Can't find FSLOUTPUTTYPE, setting to NIFTI_GZ
face_fname = mni_face_fname(mm = "1")
face = readnii(face_fname)
timg = readnii(tfname)
timg = copyNIfTIHeader(timg, timg[(91-30):(91+30), , 1:70])
face = copyNIfTIHeader(face, face[(91-30):(91+30), , 1:70])

noneck <- double_remove_neck(
  fname, 
  template.file = file.path(fsldir(), 
                            "data/standard", "MNI152_T1_1mm_brain.nii.gz"), 
  template.mask = file.path(fsldir(),
                            "data/standard", "MNI152_T1_1mm_brain_mask.nii.gz"))
#> # Registration to template
#> # Swapping Dimensions
#> FSLDIR='/usr/local/fsl'; PATH=${FSLDIR}/bin:${PATH};export PATH FSLDIR; sh "${FSLDIR}/etc/fslconf/fsl.sh"; FSLOUTPUTTYPE=NIFTI_GZ; export FSLOUTPUTTYPE; ${FSLDIR}/bin/fslhd "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmppwGYy3/sub-chibaIngeniasub-chibaIngenia_T2w.nii.gz"
#> FSLDIR='/usr/local/fsl'; PATH=${FSLDIR}/bin:${PATH};export PATH FSLDIR; sh "${FSLDIR}/etc/fslconf/fsl.sh"; FSLOUTPUTTYPE=NIFTI_GZ; export FSLOUTPUTTYPE; ${FSLDIR}/bin/fslorient -getorient "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmppwGYy3/sub-chibaIngeniasub-chibaIngenia_T2w.nii.gz"
#> FSLDIR='/usr/local/fsl'; PATH=${FSLDIR}/bin:${PATH};export PATH FSLDIR; sh "${FSLDIR}/etc/fslconf/fsl.sh"; FSLOUTPUTTYPE=NIFTI_GZ; export FSLOUTPUTTYPE; ${FSLDIR}/bin/fslorient -swaporient "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmppwGYy3/file6aaa35f6f8ac/sub-chibaIngeniasub-chibaIngenia_T2w.nii.gz"
#> FSLDIR='/usr/local/fsl'; PATH=${FSLDIR}/bin:${PATH};export PATH FSLDIR; sh "${FSLDIR}/etc/fslconf/fsl.sh"; FSLOUTPUTTYPE=NIFTI_GZ; export FSLOUTPUTTYPE; ${FSLDIR}/bin/fslswapdim "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmppwGYy3/file6aaa35f6f8ac/sub-chibaIngeniasub-chibaIngenia_T2w.nii.gz"  RL PA IS "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmppwGYy3/file6aaa693fffe1";
#> # Running Registration of file to template
#> # Applying Registration output is
#> $warpedmovout
#> antsImage
#>   Pixel Type          : float 
#>   Components Per Pixel: 1 
#>   Dimensions          : 63x320x320 
#>   Voxel Spacing       : 0.799999237060547x0.800000011920929x0.800000011920929 
#>   Origin              : -13.47086 120.9473 -144.1494 
#>   Direction           : 0.9991751 -0.02116131 0.03466073 -0.02100112 -0.9997671 -0.004979392 -0.03475802 -0.00424737 0.9993867 
#> 
#> 
#> $warpedfixout
#> antsImage
#>   Pixel Type          : float 
#>   Components Per Pixel: 1 
#>   Dimensions          : 182x218x182 
#>   Voxel Spacing       : 1x1x1 
#>   Origin              : -90 126 -72 
#>   Direction           : 1 0 0 0 -1 0 0 0 1 
#> 
#> 
#> $fwdtransforms
#> [1] "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmppwGYy3/file6aaa192704440GenericAffine.mat"
#> 
#> $invtransforms
#> [1] "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmppwGYy3/file6aaa192704440GenericAffine.mat"
#> 
#> $prev_transforms
#> character(0)
#> # Applying Transformations to file
#> # Applying Transforms to other.files
#> # Writing out file
#> # Writing out other.files
#> # Removing Warping images
#> # Reading data back into R
#> # Reading in Transformed data
#> # Dropping slices not in mask
#> # Swapping Dimensions Back
#> FSLDIR='/usr/local/fsl'; PATH=${FSLDIR}/bin:${PATH};export PATH FSLDIR; sh "${FSLDIR}/etc/fslconf/fsl.sh"; FSLOUTPUTTYPE=NIFTI_GZ; export FSLOUTPUTTYPE; ${FSLDIR}/bin/fslorient -swaporient "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmppwGYy3/file6aaa6c3724a6.nii.gz"
#> FSLDIR='/usr/local/fsl'; PATH=${FSLDIR}/bin:${PATH};export PATH FSLDIR; sh "${FSLDIR}/etc/fslconf/fsl.sh"; FSLOUTPUTTYPE=NIFTI_GZ; export FSLOUTPUTTYPE; ${FSLDIR}/bin/fslswapdim "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmppwGYy3/file6aaa6c3724a6.nii.gz"  LR PA IS "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmppwGYy3/file6aaa184745cb";
#> # Registration to template
#> # Swapping Dimensions
#> FSLDIR='/usr/local/fsl'; PATH=${FSLDIR}/bin:${PATH};export PATH FSLDIR; sh "${FSLDIR}/etc/fslconf/fsl.sh"; FSLOUTPUTTYPE=NIFTI_GZ; export FSLOUTPUTTYPE; ${FSLDIR}/bin/fslhd "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmppwGYy3/file6aaa6f371ee0.nii.gz"
#> FSLDIR='/usr/local/fsl'; PATH=${FSLDIR}/bin:${PATH};export PATH FSLDIR; sh "${FSLDIR}/etc/fslconf/fsl.sh"; FSLOUTPUTTYPE=NIFTI_GZ; export FSLOUTPUTTYPE; ${FSLDIR}/bin/fslorient -getorient "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmppwGYy3/file6aaa6f371ee0.nii.gz"
#> FSLDIR='/usr/local/fsl'; PATH=${FSLDIR}/bin:${PATH};export PATH FSLDIR; sh "${FSLDIR}/etc/fslconf/fsl.sh"; FSLOUTPUTTYPE=NIFTI_GZ; export FSLOUTPUTTYPE; ${FSLDIR}/bin/fslorient -swaporient "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmppwGYy3/file6aaadc3d81a/file6aaa6f371ee0.nii.gz"
#> FSLDIR='/usr/local/fsl'; PATH=${FSLDIR}/bin:${PATH};export PATH FSLDIR; sh "${FSLDIR}/etc/fslconf/fsl.sh"; FSLOUTPUTTYPE=NIFTI_GZ; export FSLOUTPUTTYPE; ${FSLDIR}/bin/fslswapdim "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmppwGYy3/file6aaadc3d81a/file6aaa6f371ee0.nii.gz"  RL PA IS "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmppwGYy3/file6aaa34a09a05";
#> # Running Registration of file to template
#> # Applying Registration output is
#> $warpedmovout
#> antsImage
#>   Pixel Type          : float 
#>   Components Per Pixel: 1 
#>   Dimensions          : 63x320x320 
#>   Voxel Spacing       : 0.799999237060547x0.800000011920929x0.800000011920929 
#>   Origin              : -13.47086 120.9473 -144.1494 
#>   Direction           : 0.9991751 -0.02116131 0.03466073 -0.02100112 -0.9997671 -0.004979392 -0.03475802 -0.00424737 0.9993867 
#> 
#> 
#> $warpedfixout
#> antsImage
#>   Pixel Type          : float 
#>   Components Per Pixel: 1 
#>   Dimensions          : 182x218x182 
#>   Voxel Spacing       : 1x1x1 
#>   Origin              : -90 126 -72 
#>   Direction           : 1 0 0 0 -1 0 0 0 1 
#> 
#> 
#> $fwdtransforms
#> [1] "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmppwGYy3/file6aaa53048b30GenericAffine.mat"
#> 
#> $invtransforms
#> [1] "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmppwGYy3/file6aaa53048b30GenericAffine.mat"
#> 
#> $prev_transforms
#> character(0)
#> # Applying Transformations to file
#> # Applying Transforms to other.files
#> # Writing out file
#> # Writing out other.files
#> # Removing Warping images
#> # Reading data back into R
#> # Reading in Transformed data
#> # Dropping slices not in mask
#> # Swapping Dimensions Back
#> FSLDIR='/usr/local/fsl'; PATH=${FSLDIR}/bin:${PATH};export PATH FSLDIR; sh "${FSLDIR}/etc/fslconf/fsl.sh"; FSLOUTPUTTYPE=NIFTI_GZ; export FSLOUTPUTTYPE; ${FSLDIR}/bin/fslorient -swaporient "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmppwGYy3/file6aaa14448914.nii.gz"
#> FSLDIR='/usr/local/fsl'; PATH=${FSLDIR}/bin:${PATH};export PATH FSLDIR; sh "${FSLDIR}/etc/fslconf/fsl.sh"; FSLOUTPUTTYPE=NIFTI_GZ; export FSLOUTPUTTYPE; ${FSLDIR}/bin/fslswapdim "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmppwGYy3/file6aaa14448914.nii.gz"  LR PA IS "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmppwGYy3/file6aaa1f838a71";

reg = registration(
  filename = timg,
  template.file = noneck,
  typeofTransform = "Rigid")
#> # Running Registration of file to template
#> # Applying Registration output is
#> $fwdtransforms
#> [1] "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmppwGYy3/file6aaa46ce1e1d0GenericAffine.mat"
#> 
#> $invtransforms
#> [1] "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmppwGYy3/file6aaa46ce1e1d0GenericAffine.mat"
#> 
#> $prev_transforms
#> character(0)
#> # Applying Transformations to file
#>  [1] "-d"                                                                                             
#>  [2] "3"                                                                                              
#>  [3] "-i"                                                                                             
#>  [4] "<pointer: 0x7f80bf75cd20>"                                                                      
#>  [5] "-o"                                                                                             
#>  [6] "<pointer: 0x7f80c6168d60>"                                                                      
#>  [7] "-r"                                                                                             
#>  [8] "<pointer: 0x7f80bf79e8f0>"                                                                      
#>  [9] "-n"                                                                                             
#> [10] "linear"                                                                                         
#> [11] "-t"                                                                                             
#> [12] "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmppwGYy3/file6aaa46ce1e1d0GenericAffine.mat"
#> # Writing out file
#> [1] "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmppwGYy3/file6aaa74f610e0.nii.gz"
#> # Reading data back into R
out_img = ants_apply_transforms(moving = 1 - face, 
                            fixed = noneck,
                            transformlist = reg$fwdtransforms,
                            interpolator = "nearestNeighbor")

ortho2(img, out_img)


xx = which(out_img == 1, arr.ind = TRUE)
dim3 = 1:(min(xx[,3]) + floor(  (max(xx[,3]) - min(xx[,3]))/2))
dim2 = min(xx[,2]):dim(out_img)[2]
# eg = expand.grid(dim1 = unique(xx[,1]), 
# dim2 = unique(xx[,2]), dim3 = dim3)
eg = expand.grid(dim1 = unique(xx[,1]),
                 dim2 = dim2, dim3 = dim3)
eg = as.matrix(eg)
new_mask = out_img
new_mask[eg] = 1
ortho2(img, new_mask)

ortho2(img, 1 - new_mask)

Created on 2020-02-12 by the reprex package (v0.3.0.9001)

Session info ``` r sessioninfo::session_info() #> ─ Session info ─────────────────────────────────────────────────────────────── #> setting value #> version R version 3.6.2 (2019-12-12) #> os macOS Mojave 10.14.6 #> system x86_64, darwin15.6.0 #> ui X11 #> language (EN) #> collate en_US.UTF-8 #> ctype en_US.UTF-8 #> tz America/New_York #> date 2020-02-12 #> #> ─ Packages ─────────────────────────────────────────────────────────────────── #> package * version date lib source #> abind 1.4-5 2016-07-21 [1] CRAN (R 3.6.0) #> ANTsR 0.5.4.2 2019-12-06 [1] Neuroconductor-Releases (R 3.6.1) #> ANTsRCore 0.7.3 2019-12-05 [1] Neuroconductor-Releases (R 3.6.1) #> assertthat 0.2.1 2019-03-21 [1] CRAN (R 3.6.0) #> backports 1.1.5 2019-10-02 [1] CRAN (R 3.6.0) #> bitops 1.0-6 2013-08-17 [1] CRAN (R 3.6.0) #> cli 2.0.1 2020-01-08 [1] CRAN (R 3.6.0) #> crayon 1.3.4 2017-09-16 [1] CRAN (R 3.6.0) #> curl 4.3 2019-12-02 [1] CRAN (R 3.6.0) #> digest 0.6.23 2019-11-23 [1] CRAN (R 3.6.0) #> evaluate 0.14 2019-05-28 [1] CRAN (R 3.6.0) #> extrantsr * 3.9.12.1 2020-01-03 [1] local #> fansi 0.4.1 2020-01-08 [1] CRAN (R 3.6.0) #> fs 1.3.1 2019-05-06 [1] CRAN (R 3.6.0) #> fslr * 2.24.4 2020-02-12 [1] local #> glue 1.3.1 2019-03-12 [1] CRAN (R 3.6.0) #> highr 0.8 2019-03-20 [1] CRAN (R 3.6.0) #> htmltools 0.4.0 2019-10-04 [1] CRAN (R 3.6.0) #> httr 1.4.1 2019-08-05 [1] CRAN (R 3.6.0) #> ITKR 0.5.2 2019-11-28 [1] Neuroconductor-Releases (R 3.6.1) #> knitr 1.26.1 2020-01-05 [1] Github (muschellij2/knitr@f5af631) #> lattice 0.20-38 2018-11-04 [1] CRAN (R 3.6.2) #> magrittr 1.5 2014-11-22 [1] CRAN (R 3.6.0) #> Matrix 1.2-18 2019-11-27 [1] CRAN (R 3.6.2) #> matrixStats 0.55.0 2019-09-07 [1] CRAN (R 3.6.0) #> mgcv 1.8-31 2019-11-09 [1] CRAN (R 3.6.2) #> mime 0.8 2019-12-19 [1] CRAN (R 3.6.0) #> neurobase * 1.30.0 2020-02-07 [1] local #> nlme 3.1-142 2019-11-07 [1] CRAN (R 3.6.2) #> oro.nifti * 0.10.2 2020-02-05 [1] local #> pillar 1.4.3 2019-12-20 [1] CRAN (R 3.6.0) #> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 3.6.0) #> plyr 1.8.5 2019-12-10 [1] CRAN (R 3.6.0) #> purrr 0.3.3 2019-10-18 [1] CRAN (R 3.6.0) #> R.matlab 3.6.2 2018-09-27 [1] CRAN (R 3.6.0) #> R.methodsS3 1.7.1 2016-02-16 [1] CRAN (R 3.6.0) #> R.oo 1.23.0 2019-11-03 [1] CRAN (R 3.6.0) #> R.utils 2.9.0 2019-06-13 [1] CRAN (R 3.6.0) #> R6 2.4.1 2019-11-12 [1] CRAN (R 3.6.0) #> Rcpp 1.0.3 2019-11-08 [1] CRAN (R 3.6.0) #> RcppEigen 0.3.3.7.0 2019-11-16 [1] CRAN (R 3.6.0) #> reprex 0.3.0.9001 2020-01-05 [1] Github (tidyverse/reprex@5ae0b29) #> rlang 0.4.2 2019-11-23 [1] CRAN (R 3.6.0) #> rmarkdown 2.0.7 2020-01-17 [1] Github (rstudio/rmarkdown@2faf16a) #> RNifti * 1.0.1 2019-11-27 [1] CRAN (R 3.6.0) #> rstudioapi 0.10.0-9003 2020-01-05 [1] Github (rstudio/rstudioapi@abe596d) #> sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.6.0) #> stapler 0.7.1 2020-02-07 [1] local #> stringi 1.4.5 2020-01-11 [1] CRAN (R 3.6.0) #> stringr 1.4.0 2019-02-10 [1] CRAN (R 3.6.0) #> styler 1.1.1 2019-05-06 [1] CRAN (R 3.6.0) #> tibble 2.1.3 2019-06-06 [1] CRAN (R 3.6.0) #> WhiteStripe 2.3.2 2019-10-01 [1] local #> withr 2.1.2 2018-03-15 [1] CRAN (R 3.6.0) #> xfun 0.11 2019-11-12 [1] CRAN (R 3.6.0) #> xml2 1.2.2 2019-08-09 [1] CRAN (R 3.6.0) #> yaml 2.2.0 2018-07-25 [1] CRAN (R 3.6.0) #> #> [1] /Library/Frameworks/R.framework/Versions/3.6/Resources/library ```
alexfoias commented 4 years ago

@muschellij2 The proposed solution works for me. Thank you for your help. I managed to face my dataset. We can close this issue.