1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
|
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* https://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#ifndef itkCompensatedSummation_h
#define itkCompensatedSummation_h
#include "itkNumericTraits.h"
#include "itkConceptChecking.h"
namespace itk
{
/** \class CompensatedSummation
* \brief Perform more precise accumulation of floating point numbers.
*
* The \c float and \c double datatypes only have finite precision. When
* performing a running sum of floats, the accumulated errors get progressively
* worse as the magnitude of the sum gets large relative to new elements.
*
* From Wikipedia, https://en.wikipedia.org/wiki/Kahan_summation_algorithm
*
* "In numerical analysis, the Kahan summation algorithm (also known as
* compensated summation) significantly reduces the numerical error in the total
* obtained by adding a sequence of finite precision floating point numbers,
* compared to the obvious approach. This is done by keeping a separate running
* compensation (a variable to accumulate small errors)."
*
* For example, instead of
\code
double sum = 0.0;
for( unsigned int i = 0; i < array.Size(); ++i )
{
sum += array.GetElement(i);
}
\endcode
*
* do
*
\code
using CompensatedSummationType = CompensatedSummation<double>;
CompensatedSummationType compensatedSummation;
for( unsigned int i = 0; i < array.Size(); ++i )
{
compensatedSummation += array.GetElement(i);
}
double sum = compensatedSummation.GetSum();
\endcode
*
* \ingroup ITKCommon
*/
template <typename TFloat>
class ITK_TEMPLATE_EXPORT CompensatedSummation
{
public:
/** Type of the input elements. */
using FloatType = TFloat;
/** Type used for the sum and compensation. */
using AccumulateType = typename NumericTraits<FloatType>::AccumulateType;
/** Standard class type aliases. */
using Self = CompensatedSummation;
/** Constructors. */
CompensatedSummation() = default;
CompensatedSummation(FloatType value);
/** Copy constructor. */
CompensatedSummation(const Self & rhs);
/** Assignment operator. */
Self &
operator=(const Self & rhs);
/** Add an element to the sum. */
void
AddElement(const FloatType & element);
Self &
operator+=(const FloatType & rhs);
Self &
operator+=(const Self & rhs);
/** Subtract an element from the sum. */
Self &
operator-=(const FloatType & rhs);
/** Division and multiplication. These do not provide any numerical advantages
* relative to vanilla division and multiplication. */
Self &
operator*=(const FloatType & rhs);
Self &
operator/=(const FloatType & rhs);
/** Reset the sum and compensation to zero. */
void
ResetToZero();
/** Reset the sum to the given value and the compensation to zero. */
Self &
operator=(const FloatType & rhs);
/** Get the sum. */
const AccumulateType &
GetSum() const;
/** explicit conversion */
explicit operator FloatType() const;
private:
AccumulateType m_Sum{};
AccumulateType m_Compensation{};
// Maybe support more types in the future with template specialization.
#ifdef ITK_USE_CONCEPT_CHECKING
itkConceptMacro(OnlyDefinedForFloatingPointTypes, (itk::Concept::IsFloatingPoint<TFloat>));
#endif // ITK_USE_CONCEPT_CHECKING
};
void ITKCommon_EXPORT
CompensatedSummationAddElement(float & compensation, float & sum, const float element);
void ITKCommon_EXPORT
CompensatedSummationAddElement(double & compensation, double & sum, const double element);
} // end namespace itk
#ifndef ITK_MANUAL_INSTANTIATION
# include "itkCompensatedSummation.hxx"
#endif
#endif
|
