search

muschellij2 / WhiteStripe

WhiteStripe Package
9 stars 6 forks source link

Error while running WhiteStripe #5

Closed chihyanghsu0805 closed 7 years ago

chihyanghsu0805 commented 7 years ago

Hi, thank you for the awesome package.

As I was using WhiteStripe to process my data with the following line, wsIndex = whitestripe(img = ws, type = ImageType, stripped = FALSE)

one of my cases raised the following error. Error in model.frame.default(formula = y ~ X - 1, drop.unused.levels = TRUE)

All my images are without skull stripping, but if I used TRUE for the argument stripped, the error-causing case would actually finish.

Therefore, I was wondering how significant does skull-stripping affect the normalization. Thank you for your time.

muschellij2 commented 7 years ago

I think I need some more details. Do you have an example image you can share? What is ws in your example?

chihyanghsu0805 commented 7 years ago

ws is an NIFTI object read using the following line.

ws = readNIfTI(filename)

Sorry, but I can't share my images due to confidentiality.

For some of my cases, I also run into the following error

Error in gam(y ~ X - 1, paraPen = list(X = D), method = method, ...) : Model has more coefficients than data

I ma wondering whether it is because of my image registration or bias correction.

chihyanghsu0805 commented 7 years ago

After carefully reading the instructions, I found that it was noted that

This method performs white matter mean and standard deviation estimates on data that has been rigidly-registered to the MNI template and uses histogram-based methods.

I have not register to MNI template, I will try it.

muschellij2 commented 7 years ago

Hope this worked. Reopen if not.

MLTOL commented 3 years ago

Hi,

I have a similar issue. I would like to apply the normalisation to the DWI images from the ISLES 2015 challange (SISS). In addition to their preprocess I have applied rigid registration to register the images to MNI space. I also have strippid the background. I get the following error:

Smoothing Histogram

Error in model.frame.default(formula = y ~ X - 1, drop.unused.levels = TRUE) : invalid type (list) for variable 'X' Timing stopped at: 0 0 0

First, I thought it was a type problem, since the images were stored as single. However, converting to double did not help. No NaNs are present in the image array either.

Applying the normalisation to images from our own center that went through the same pre-process did not raise any issues.

Do you have any suggestion on how to solve this issue?

Thank you.

muschellij2 commented 3 years ago

Can you make a reproducible example? It’ll be great if you can share the data from the challenge. I’ve done the challenge as well so should have access to the data but it’s much easier if you make a reproducible example. Also indicate which identifiers are failing please

On Thu, Mar 11, 2021 at 6:56 AM MLTOL @.***> wrote:

Hi,

I have a similar issue. I would like to apply the normalisation to the DWI images from the ISLES 2015 challange (SISS). In addition to their preprocess I have applied rigid registration to register the images to MNI space. I also have strippid the background. I get the following error:

Smoothing Histogram

Error in model.frame.default(formula = y ~ X - 1, drop.unused.levels = TRUE) : invalid type (list) for variable 'X' Timing stopped at: 0 0 0

First, I thought it was a type problem, since the images were stored as single. However, converting to double did not help. No NaNs are present in the image array either.

Applying the normalisation to images from our own center that went through the same pre-process did not raise any issues.

Do you have any suggestion on how to solve this issue?

Thank you.

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/muschellij2/WhiteStripe/issues/5#issuecomment-796683903, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAIGPLSIJ337PL4WYDY5KO3TDCVVPANCNFSM4DHKDTQQ .

-- Best, John

MLTOL commented 3 years ago

Thank you for your quick response.

I have shared one of the images.

as formality: This SICAS Medical Image Repository Data is made available under the Open Database License: http://opendatacommons.org/licenses/odbl/1.0/. Any rights in individual contents of the database are licensed under the Database Contents License: http://opendatacommons.org/licenses/dbcl/1.0/.

rb1000_backgroundstrip.nii.gz

Rigid registration and background strip has been applied to the original image using the spm8 toolbox, which has placed the image in MNI space.

For normalisation, I have used the following code:

img=readNIfTI(fname = "rb1000_backgroundstrip.nii.gz" , reorient = FALSE)

apply whitestripe algorithm

ws = whitestripe(img = img, type = "T1", stripped = TRUE) norm = whitestripe_norm(img = img, indices=ws$whitestipe.ind) std = ws$mu.whitestripe-ws$whitestripe[1]

apply normalisation

norm_img = (img - ws$mu.whitestripe)/std

datatype(norm_img) <- 64 bitpix(norm_img) <-64

save normalised image

writeNIfTI(norm_img, "normalised_image", gzipped=FALSE)

muschellij2 commented 3 years ago

The image is all zero:

library(WhiteStripe)
library(neurobase)
#> Loading required package: oro.nifti
#> oro.nifti 0.11.0
img=readNIfTI(fname = "~/Downloads/rb1000_backgroundstrip.nii.gz" , reorient = FALSE)
range(img)
#> [1] 0 0
ortho2(img)

Created on 2021-03-11 by the reprex package (v1.0.0)

Session info ``` r sessioninfo::session_info() #> ─ Session info ─────────────────────────────────────────────────────────────── #> setting value #> version R version 4.0.2 (2020-06-22) #> os macOS Catalina 10.15.7 #> system x86_64, darwin17.0 #> ui X11 #> language (EN) #> collate en_US.UTF-8 #> ctype en_US.UTF-8 #> tz America/New_York #> date 2021-03-11 #> #> ─ Packages ─────────────────────────────────────────────────────────────────── #> package * version date lib source #> abind 1.4-5 2016-07-21 [2] CRAN (R 4.0.0) #> assertthat 0.2.1 2019-03-21 [2] CRAN (R 4.0.0) #> backports 1.2.1 2020-12-09 [1] CRAN (R 4.0.2) #> bitops 1.0-6 2013-08-17 [2] CRAN (R 4.0.0) #> cli 2.3.0 2021-01-31 [1] CRAN (R 4.0.2) #> crayon 1.4.1 2021-02-08 [1] CRAN (R 4.0.2) #> curl 4.3 2019-12-02 [2] CRAN (R 4.0.0) #> debugme 1.1.0 2017-10-22 [2] CRAN (R 4.0.0) #> digest 0.6.27 2020-10-24 [1] CRAN (R 4.0.2) #> ellipsis 0.3.1 2020-05-15 [2] CRAN (R 4.0.0) #> evaluate 0.14 2019-05-28 [2] CRAN (R 4.0.0) #> fansi 0.4.2 2021-01-15 [1] CRAN (R 4.0.2) #> fs 1.5.0 2020-07-31 [2] CRAN (R 4.0.2) #> glue 1.4.2 2020-08-27 [1] CRAN (R 4.0.2) #> highr 0.8 2019-03-20 [2] CRAN (R 4.0.0) #> htmltools 0.5.1.1 2021-01-22 [1] CRAN (R 4.0.2) #> httr 1.4.2 2020-07-20 [2] CRAN (R 4.0.2) #> knitr 1.31 2021-01-27 [1] CRAN (R 4.0.2) #> lattice 0.20-41 2020-04-02 [2] CRAN (R 4.0.2) #> lifecycle 1.0.0 2021-02-15 [1] CRAN (R 4.0.2) #> magrittr 2.0.1 2020-11-17 [1] CRAN (R 4.0.2) #> Matrix 1.3-2 2021-01-06 [1] CRAN (R 4.0.2) #> matrixStats 0.58.0 2021-01-29 [1] CRAN (R 4.0.2) #> mgcv 1.8-33 2020-08-27 [1] CRAN (R 4.0.2) #> mime 0.10 2021-02-13 [1] CRAN (R 4.0.2) #> neurobase * 1.31.0 2020-10-07 [1] local #> nlme 3.1-152 2021-02-04 [1] CRAN (R 4.0.2) #> oro.nifti * 0.11.0 2020-09-04 [2] local #> pillar 1.5.1 2021-03-05 [1] CRAN (R 4.0.2) #> pkgconfig 2.0.3 2019-09-22 [2] CRAN (R 4.0.0) #> purrr 0.3.4 2020-04-17 [2] CRAN (R 4.0.0) #> R.methodsS3 1.8.1 2020-08-26 [1] CRAN (R 4.0.2) #> R.oo 1.24.0 2020-08-26 [1] CRAN (R 4.0.2) #> R.utils 2.10.1 2020-08-26 [1] CRAN (R 4.0.2) #> R6 2.5.0 2020-10-28 [1] CRAN (R 4.0.2) #> Rcpp 1.0.6 2021-01-15 [1] CRAN (R 4.0.2) #> reprex 1.0.0 2021-01-27 [1] CRAN (R 4.0.2) #> rlang 0.4.10 2020-12-30 [1] CRAN (R 4.0.2) #> rmarkdown 2.6 2020-12-14 [1] CRAN (R 4.0.2) #> RNifti 1.3.0 2020-12-04 [1] CRAN (R 4.0.2) #> sessioninfo 1.1.1 2018-11-05 [2] CRAN (R 4.0.0) #> stringi 1.5.3 2020-09-09 [1] CRAN (R 4.0.2) #> stringr 1.4.0 2019-02-10 [2] CRAN (R 4.0.0) #> styler 1.3.2 2020-02-23 [2] CRAN (R 4.0.0) #> tibble 3.1.0 2021-02-25 [1] CRAN (R 4.0.2) #> utf8 1.1.4 2018-05-24 [2] CRAN (R 4.0.0) #> vctrs 0.3.6 2020-12-17 [1] CRAN (R 4.0.2) #> WhiteStripe * 2.3.2 2019-10-01 [2] CRAN (R 4.0.0) #> withr 2.4.1 2021-01-26 [1] CRAN (R 4.0.2) #> xfun 0.21 2021-02-10 [1] CRAN (R 4.0.2) #> xml2 1.3.2 2020-04-23 [2] CRAN (R 4.0.0) #> yaml 2.2.1 2020-02-01 [2] CRAN (R 4.0.0) #> #> [1] /Users/johnmuschelli/Library/R/4.0/library #> [2] /Library/Frameworks/R.framework/Versions/4.0/Resources/library ```
MLTOL commented 3 years ago

My apologies, something went wrong during gzipping. This image should work.

rb1000_backgroundstrip.nii.gz

muschellij2 commented 3 years ago

So the issue is this https://github.com/muschellij2/WhiteStripe/blob/fd993fad86e3e3b30a7771f36bcecea94375d5f8/R/get.largest.mode.R#L78: 

  if (remove.tail) {
    which.rare <- y < (rare.prop * max(y))
    y = y[!which.rare]
    x = x[!which.rare]
  }

And your data has a HUGE spike around zero. The data you have is not integer, which is interesting given it should likely be a raw-ish image. I would definitely look at rounding issues or rounding the data to some value. Here I show you exactly what I mean and have put a warning in there just in case. There are 2 solutions presented:

Load in Packages

library(RNifti)
library(neurobase)
#> Loading required package: oro.nifti
#> oro.nifti 0.11.0
#> 
#> Attaching package: 'oro.nifti'
#> The following object is masked from 'package:RNifti':
#> 
#>     origin
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
library(WhiteStripe)
fname = "~/Downloads/rb1000_backgroundstrip.nii.gz"

Read image and explore

img = readNifti(fname)
range(img)
#> [1] -193.8515 2141.6748
img=readNIfTI(fname = fname , reorient = FALSE)
range(img)
#> [1] -193.8515 2141.6748
ortho2(img)

LARGEE SPike at zero

This is likely due to really small (but not zero values)

hist(c(img), breaks = 100)

VOI That’s being created

This is from the whitestripe code and defaults

breaks = 2000
img.voi <- img[img > 0]

See small values

range(img.voi)
#> [1] 2.283842e-05 2.141675e+03
img.hist = hist(img.voi, breaks = breaks, plot = FALSE)
y.in = img.hist$counts
x.in = img.hist$mids
x.in = x.in[!is.na(y.in)]
y.in = y.in[!is.na(y.in)]

The largest spike is x of 0.5

x.in[which.max(y.in)]
#> [1] 0.5

range(y.in)
#> [1]     0 94659

This causes HUGE SPIKE

plot(img.hist)

see get.last.mode

The get.last.mode uses remove.tail and the proportion below the max times rare.prop is removed

range(x.in)
#> [1]    0.5 2141.5
max(y.in) 
#> [1] 94659
max(y.in) * 1/5
#> [1] 18931.8
quantile(y.in, probs = 0.999)
#>    99.9% 
#> 6335.542
max(y.in) * 1/5
#> [1] 18931.8

REMEMBER - y.in represents counts from the histogram, not observed values # Fail

ws = whitestripe(img = img, type = "T1", stripped = TRUE)
#> Making Image VOI
#> Making T1 Histogram
#> Warning in whitestripe(img = img, type = "T1", stripped = TRUE): Stripped data has very small, but > 0 values, 
#>                    probably rounding needed, such as img[abs(img) <= 0.1] = 0
#> Getting T1 Modes
#> Smoothing Histogram
#> Error in model.frame.default(formula = y ~ X - 1, drop.unused.levels = TRUE): invalid type (list) for variable 'X'
#> Timing stopped at: 0.008 0 0.008

Using a different rare proportion

max(y.in) * 1/100
#> [1] 946.59

Could also use a quantile on the counts

quantile(y.in, probs = 0.999) * 1/5
#>    99.9% 
#> 1267.108

See where the tail removal works

plot(img.hist)
abline(h = max(y.in) * 1/5, col = "red")
abline(h = max(y.in) * 1/100, col = "blue")
abline(h = quantile(y.in, probs = 0.999) * 1/5, col = "green")

You can change rare.prop, but not a great solution necessarily

ws = whitestripe(img = img, type = "T1", stripped = TRUE, rare.prop = 1/100)
#> Making Image VOI
#> Making T1 Histogram
#> Warning in whitestripe(img = img, type = "T1", stripped = TRUE, rare.prop = 1/100): Stripped data has very small, but > 0 values, 
#>                    probably rounding needed, such as img[abs(img) <= 0.1] = 0
#> Getting T1 Modes
#> Smoothing Histogram
#> Smoothing Derivative
#> Quantile T1 VOI

Or you should likely round your data at zero (recommended)

img[img < 0] = 0
img[ abs(img) < 1]= 0
ws = whitestripe(img = img, type = "T1", stripped = TRUE)
#> Making Image VOI
#> Making T1 Histogram
#> Getting T1 Modes
#> Smoothing Histogram
#> Smoothing Derivative
#> Quantile T1 VOI

Created on 2021-03-11 by the reprex package (v1.0.0)

Session info ``` r sessioninfo::session_info() #> ─ Session info ─────────────────────────────────────────────────────────────── #> setting value #> version R version 4.0.2 (2020-06-22) #> os macOS Catalina 10.15.7 #> system x86_64, darwin17.0 #> ui X11 #> language (EN) #> collate en_US.UTF-8 #> ctype en_US.UTF-8 #> tz America/New_York #> date 2021-03-11 #> #> ─ Packages ─────────────────────────────────────────────────────────────────── #> package * version date lib source #> abind 1.4-5 2016-07-21 [2] CRAN (R 4.0.0) #> ANTsR 0.5.7.5 2021-02-16 [1] Github (ANTsX/ANTsR@4abee8d) #> ANTsRCore 0.7.4.9 2021-02-16 [1] Github (ANTsX/ANTsRCore@ddf445b) #> assertthat 0.2.1 2019-03-21 [2] CRAN (R 4.0.0) #> backports 1.2.1 2020-12-09 [1] CRAN (R 4.0.2) #> bitops 1.0-6 2013-08-17 [2] CRAN (R 4.0.0) #> cli 2.3.0 2021-01-31 [1] CRAN (R 4.0.2) #> crayon 1.4.1 2021-02-08 [1] CRAN (R 4.0.2) #> curl 4.3 2019-12-02 [2] CRAN (R 4.0.0) #> debugme 1.1.0 2017-10-22 [2] CRAN (R 4.0.0) #> digest 0.6.27 2020-10-24 [1] CRAN (R 4.0.2) #> ellipsis 0.3.1 2020-05-15 [2] CRAN (R 4.0.0) #> evaluate 0.14 2019-05-28 [2] CRAN (R 4.0.0) #> extrantsr * 3.9.13.1 2020-09-03 [2] Github (muschellij2/extrantsr@00c75ad) #> fansi 0.4.2 2021-01-15 [1] CRAN (R 4.0.2) #> fs 1.5.0 2020-07-31 [2] CRAN (R 4.0.2) #> fslr 2.25.0 2021-02-16 [1] local #> glue 1.4.2 2020-08-27 [1] CRAN (R 4.0.2) #> highr 0.8 2019-03-20 [2] CRAN (R 4.0.0) #> htmltools 0.5.1.1 2021-01-22 [1] CRAN (R 4.0.2) #> httr 1.4.2 2020-07-20 [2] CRAN (R 4.0.2) #> ITKR 0.5.3.3.0 2021-02-15 [1] Github (stnava/ITKR@ea0ac19) #> knitr 1.31 2021-01-27 [1] CRAN (R 4.0.2) #> lattice 0.20-41 2020-04-02 [2] CRAN (R 4.0.2) #> lifecycle 1.0.0 2021-02-15 [1] CRAN (R 4.0.2) #> magrittr 2.0.1 2020-11-17 [1] CRAN (R 4.0.2) #> Matrix 1.3-2 2021-01-06 [1] CRAN (R 4.0.2) #> matrixStats 0.58.0 2021-01-29 [1] CRAN (R 4.0.2) #> mgcv 1.8-33 2020-08-27 [1] CRAN (R 4.0.2) #> mime 0.10 2021-02-13 [1] CRAN (R 4.0.2) #> neurobase * 1.31.0 2020-10-07 [1] local #> nlme 3.1-152 2021-02-04 [1] CRAN (R 4.0.2) #> oro.nifti * 0.11.0 2020-09-04 [2] local #> pillar 1.5.1 2021-03-05 [1] CRAN (R 4.0.2) #> pkgconfig 2.0.3 2019-09-22 [2] CRAN (R 4.0.0) #> plyr 1.8.6 2020-03-03 [2] CRAN (R 4.0.0) #> purrr 0.3.4 2020-04-17 [2] CRAN (R 4.0.0) #> R.matlab 3.6.2 2018-09-27 [2] CRAN (R 4.0.0) #> R.methodsS3 1.8.1 2020-08-26 [1] CRAN (R 4.0.2) #> R.oo 1.24.0 2020-08-26 [1] CRAN (R 4.0.2) #> R.utils 2.10.1 2020-08-26 [1] CRAN (R 4.0.2) #> R6 2.5.0 2020-10-28 [1] CRAN (R 4.0.2) #> Rcpp 1.0.6 2021-01-15 [1] CRAN (R 4.0.2) #> RcppEigen 0.3.3.9.1 2020-12-17 [1] CRAN (R 4.0.2) #> reprex 1.0.0 2021-01-27 [1] CRAN (R 4.0.2) #> rlang 0.4.10 2020-12-30 [1] CRAN (R 4.0.2) #> rmarkdown 2.6 2020-12-14 [1] CRAN (R 4.0.2) #> RNifti * 1.3.0 2020-12-04 [1] CRAN (R 4.0.2) #> sessioninfo 1.1.1 2018-11-05 [2] CRAN (R 4.0.0) #> stapler 0.7.2 2020-07-09 [2] Github (muschellij2/stapler@79e23d2) #> stringi 1.5.3 2020-09-09 [1] CRAN (R 4.0.2) #> stringr 1.4.0 2019-02-10 [2] CRAN (R 4.0.0) #> styler 1.3.2 2020-02-23 [2] CRAN (R 4.0.0) #> tibble 3.1.0 2021-02-25 [1] CRAN (R 4.0.2) #> utf8 1.1.4 2018-05-24 [2] CRAN (R 4.0.0) #> vctrs 0.3.6 2020-12-17 [1] CRAN (R 4.0.2) #> WhiteStripe * 2.4.0 2021-03-11 [1] local #> withr 2.4.1 2021-01-26 [1] CRAN (R 4.0.2) #> xfun 0.21 2021-02-10 [1] CRAN (R 4.0.2) #> xml2 1.3.2 2020-04-23 [2] CRAN (R 4.0.0) #> yaml 2.2.1 2020-02-01 [2] CRAN (R 4.0.0) #> #> [1] /Users/johnmuschelli/Library/R/4.0/library #> [2] /Library/Frameworks/R.framework/Versions/4.0/Resources/library ```
MLTOL commented 3 years ago

Great, thank you for pointing out that the images were not integer. I think that is a bit strange also. After rounding the images, the issue was solved. Thank you!