Line #170 in LabelStatisticsImageFilter implementation here:
is causing a floating point exception and crashing my application. I think the reason is a numerically unstable calculation of variance, resulting in a negative value passed to std::sqrt(). In my case, I have a signed 16-bit image with a uniform intensity of -100 for a large background region, the variance of which is computed as approximately -2e-12 by the filter code.
I’m using ITKv4, 4.10.1 to be specific. My image is 512x512x236. Most of it (60838487 pixels) is background. This means the code is computing variance as (608384870000 - -6083848700*-6083848700/60838487) / ( 60838487 - 1.0 ).
Kahan Summation will not resolve this issue, will it? I think the use of the formula above is the issue.
Please find a fix for the sort of negative value here:
Thank you for the link to the Walford method of computing variance. However, ITK uses a parallel divide and conquer method to compute the variance with multiple threads. It is not clear to me how this recurrent method can be adapted for chunking of that data.
The data set is not mine (or even my company’s) to share unfortunately. It should be very easy to reproduce though by creating a large image with constant intensity value of -100 (or some such number) and running the LabelStatisticsImageFilter on it.
Create a 512x512x250 image with mostly (at least 60 million) background pixels of constant intensity value -100, and a small foreground region. Run it through itk::LabelMap to create a label map of the image, and pass the image and the label map to the statistics filter.