Thanks for reporting this.

Because aggregate-images uses internally 1D histograms for the samples at each voxel, my suspicion is that this is caused by histogram bins that are not aligned with the input label discretization, but bin limits that may be slightly off. In particular, the default number of -bins used for these histograms, i.e.,

I would assume for your particular input data, the variable values should be:

• min_value = 0.0f (though cast to float may result in 1e-12f or so...)
• max_value = 95.0f
• bins = 95

Image values are internally cast to double, added to the histogram, and the mode retrieved as the center value of the bin with highest sample frequency. As you obtain non-integral output values, it seems these bin centers do not coincide with integral numbers.

Another possibility, though less likely here I hope, would be errors introduced by the conversion to double. More likely the mapping Histogram1D::BinToVal used by Histogram1D::Mode:

Now, looking at BinToVal, turns out it adds half the bin width to get the bin center. This suggests that the histogram parameters are calculated incorrectly in aggregate-images, because without the + _width/2, I would expect that the output values you would get in agg.nii.gz would be more or less integral (certainly with output -dtype uchar).

I guess the correct parameters for Histogram1D constructor should be:

• min = min_value - width / 2
• max = max_value + width / 2
• width = (max_value - min_value) / max(1, bins - 1)

where bins = (max_value - min_value) + 1.

Which basically means width = 1 and thus simplifies to:

• min = min_value - 0.5
• max = max_value + 0.5
• width = 1

Positive side effect, when min_value == max_value, width == 1 and (max - min) == 1, so no degenerate histogram bins even if input images have constant value.

@soundray Could you check the linked PR #692 to see if that fixes the issue? Thanks!

Also please feel free to review the changes and comment if anything seems odd. Given the lack of testing, this would be very good to have a second pair of eyes check the validity of the change.

Thank you for the PR. Here are the updated test results:

$ info seg95.nii.gz | grep type
Voxel type is unsigned char
$ calculate-element-wise seg95.nii.gz
Mean = 9.9747
Standard deviation = 20.819
Extrema = [0, 95]
Range = 95
$ aggregate-images mode seg95.nii.gz seg95.nii.gz seg95.nii.gz -output agg.nii.gz
$ info agg.nii.gz | grep type
Voxel type is float
$ calculate-element-wise agg.nii.gz
Mean = 9.9728
Standard deviation = 20.796
Extrema = [0, 95]
Range = 95
$ convert-image agg.nii.gz agg-char.nii.gz -char
$ info agg-char.nii.gz | grep type
Voxel type is char
$ calculate-element-wise agg-char.nii.gz 
Mean = 9.9866
Standard deviation = 20.766
Extrema = [0, 95]
Range = 95

It is still strange that the output from aggregating three identical input images is not identical to these inputs. So I break out seg_stats from NiftySeg to get a count of each label (does MIRTK have a tool for this?):

$ seg_stats seg95.nii.gz -Nl temp.csv ; cat temp.csv ; printf '\n'
4.2417e+06,2411,1875,1517,1514,7144,9791,4384,4173,5485,5697,16245,16785,21066,14997,4056,4861,68092,69993,25288,3091,3480,44426,52010,6844,5954,8160,8572,55722,52882,54776,52119,17504,12999,3839,3936,312,267,5041,4812,8067,8091,1103,1104,20371,5966,8978,677,347,894,36815,37495,3785,4425,7102,7429,20247,18707,59516,61095,31407,33339,47643,49308,13118,12389,12009,12236,6312,5944,3655,4011,5617,5747,338,328,1617,1313,281,215,875,432,4108,4443,32989,28377,2895,2371,820,1058,1851,2224,2447,2556,2484,2884
rahec@cx1-103-7-1 [stp:1] { ...ild/seg95-compare }$ seg_stats agg-char.nii.gz -Nl temp.csv ; cat temp.csv ; printf '\n'
4.2417e+06,0,4286,1517,0,8658,9791,0,8557,5485,0,21942,16785,0,36063,4056,0,72953,69993,0,28379,3480,0,96436,6844,0,14114,8572,0,108604,54776,0,69623,12999,0,7775,312,0,5308,4812,0,16158,1103,0,21475,5966,0,9655,1241,0,36815,41280,0,4425,14531,0,20247,78223,0,61095,64746,0,47643,62426,0,12389,24245,0,6312,9599,0,4011,11364,0,338,1945,0,1313,496,0,875,4540,0,4443,61366,0,2895,3191,0,1058,4075,0,2447,5040,0,2884

Pardon the long lines, but the essence is that every third label is as expected (those divisible by 3), and the other two are 0 or otherwise non-matching.
seg95.nii.gz

Hm, one step closer. Thanks for checking and sharing your input file. I will have to add proper tests for these applications, starting with this one...

PS: I thought there was a command similar to IRTK's labelStats, but apparently I didn't migrate it.

There is evaluate-overlap but it gives me strange results, too. I'll investigate and report.

Sorry about that, blame it on my PhD! :-)