ANTsX / ANTsRCore

Rcpp bindings for the C++ ANTs library used by the ANTsR package
9 stars 9 forks source link

Tests show issues with some maths functions #87

Closed muschellij2 closed 5 years ago

muschellij2 commented 5 years ago

I have added the round function to the maths functions. I have added a test tests-Maths.R to check to make sure the same result is given by the ANTsRCore function and that of base R. There are some issues with gamma function (Inf vs. NaN), but these probably can be fixed or are an implementation decision. The real issues are from sinpi and tanpi. I have to look a bit more in these, but we should definitely have tests for overloaded functions.

I had to adapt round to match R behavior. You may have a better implementation than I did.

muschellij2 commented 5 years ago

These failures have been traced to https://github.com/ANTsX/ANTsRCore/blob/master/src/antsImageMath.cpp#L922 and https://github.com/ANTsX/ANTsRCore/blob/master/src/antsImageMath.cpp#L886

muschellij2 commented 5 years ago

Also from https://stackoverflow.com/questions/31502120/sin-and-cos-give-unexpected-results-for-well-known-angles, should we use tan( pi * x) or these functions from math.h in antsImageMath.cpp?

muschellij2 commented 5 years ago

I think tanpi is still an error per the help file for tanpi in R:

tanpi(0.5) is NaN. Similarly for other inputs with fractional part 0.5.

library(ANTsRCore)
#> 
#> Attaching package: 'ANTsRCore'
#> The following object is masked from 'package:stats':
#> 
#>     var
#> The following objects are masked from 'package:base':
#> 
#>     all, any, apply, max, min, prod, range, sum
outimg <- makeImage( c(2,10) , rbinom(n = 2*10, p = 0.5, size = 2))/2
arr = as.array(outimg)

The help file shows us this is true, but only for tanpi, not tan(pi*x)

tanpi(0.5)
#> Warning in tanpi(0.5): NaNs produced
#> [1] NaN
tan(pi * 0.5)
#> [1] 1.633124e+16

We see the output are not the same

ants_output = as.array(tanpi(outimg))
r_output = tanpi(arr)
#> Warning in tanpi(arr): NaNs produced
ants_output
#>              [,1]         [,2]          [,3] [,4]          [,5]
#> [1,] 1.633124e+16 0.000000e+00 -1.224647e-16    0 -1.224647e-16
#> [2,] 1.633124e+16 1.633124e+16  1.633124e+16    0  1.633124e+16
#>               [,6]          [,7]          [,8]         [,9]         [,10]
#> [1,]  1.633124e+16 -1.224647e-16 -1.224647e-16 1.633124e+16 -1.224647e-16
#> [2,] -1.224647e-16  1.633124e+16  0.000000e+00 0.000000e+00 -1.224647e-16
r_output
#>      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,]  NaN    0    0    0    0  NaN    0    0  NaN     0
#> [2,]  NaN  NaN  NaN    0  NaN    0  NaN    0    0     0

This confirms

all(ants_output == r_output)
#> [1] FALSE
arr[ants_output != r_output | is.na(ants_output)]
#>  [1] NA NA NA  1 NA  1 NA NA  1  1 NA  1 NA  1  1

# tanpi(0.5) is NaN. Similarly for other inputs with fractional part 0.5.

Created on 2019-06-17 by the reprex package (v0.2.1)