diff --git a/changelogs/diplib_next.md b/changelogs/diplib_next.md
index 45bd175c..9aa22911 100644
--- a/changelogs/diplib_next.md
+++ b/changelogs/diplib_next.md
@@ -30,6 +30,10 @@ date: 2020-00-00
 - `dip::GetOptimalDFTSize()` could loop indefinitely for inputs close to the maximum `dip::uint` value, because
   of an integer overflow in arithmetic.
 
+- The classes `dip::StatisticsAccumulator` and `dip::VarianceAccumulator` divided by 0 when adding two empty
+  accumulators together, producing NaN values as output. This caused the function `dip::SampleStatistics()`
+  to sometimes produce NaN statistics when called with a mask image.
+
 ### Updated dependencies
 
 ### Build changes
diff --git a/include/diplib/accumulators.h b/include/diplib/accumulators.h
index db8ace70..eb032815 100644
--- a/include/diplib/accumulators.h
+++ b/include/diplib/accumulators.h
@@ -87,6 +87,14 @@ class DIP_NO_EXPORT StatisticsAccumulator {
 
       /// Combine two accumulators
       StatisticsAccumulator& operator+=( StatisticsAccumulator const& b ) {
+         if( b.n_ == 0 ) {
+            return *this;
+         }
+         if( n_ == 0 ) {
+            *this = b;
+            return *this;
+         }
+         // The code below assumes n_ + b.n_ > 0
          dfloat an = static_cast< dfloat >( n_ );
          dfloat an2 = an * an;
          dfloat bn = static_cast< dfloat >( b.n_ );
@@ -208,6 +216,14 @@ class DIP_NO_EXPORT VarianceAccumulator {
 
       /// Combine two accumulators
       VarianceAccumulator& operator+=( VarianceAccumulator const& b ) {
+         if( b.n_ == 0 ) {
+            return *this;
+         }
+         if( n_ == 0 ) {
+            *this = b;
+            return *this;
+         }
+         // The code below assumes n_ + b.n_ > 0
          dfloat oldn = static_cast< dfloat >( n_ );
          n_ += b.n_;
          dfloat n = static_cast< dfloat >( n_ );
@@ -384,21 +400,25 @@ class DIP_NO_EXPORT CovarianceAccumulator {
 
       /// Combine two accumulators
       CovarianceAccumulator& operator+=( const CovarianceAccumulator& other ) {
+         if( other.n_ == 0 ) {
+            return *this;
+         }
          if( n_ == 0 ) {
             *this = other; // copy over the data
-         } else if( other.n_ > 0 ) {
-            dip::uint intN = n_ + other.n_;
-            dfloat N = static_cast< dfloat >( intN );
-            dfloat dx = other.meanx_ - meanx_;
-            dfloat dy = other.meany_ - meany_;
-            meanx_ = ( static_cast< dfloat >( n_ ) * meanx_ + static_cast< dfloat >( other.n_ ) * other.meanx_ ) / N;
-            meany_ = ( static_cast< dfloat >( n_ ) * meany_ + static_cast< dfloat >( other.n_ ) * other.meany_ ) / N;
-            dfloat fN = static_cast< dfloat >( n_ * other.n_ ) / N;
-            m2x_ += other.m2x_ + dx * dx * fN;
-            m2y_ += other.m2y_ + dy * dy * fN;
-            C_ += other.C_ + dx * dy * fN;
-            n_ = intN;
+            return *this;
          }
+         // The code below assumes n_ + b.n_ > 0
+         dip::uint intN = n_ + other.n_;
+         dfloat N = static_cast< dfloat >( intN );
+         dfloat dx = other.meanx_ - meanx_;
+         dfloat dy = other.meany_ - meany_;
+         meanx_ = ( static_cast< dfloat >( n_ ) * meanx_ + static_cast< dfloat >( other.n_ ) * other.meanx_ ) / N;
+         meany_ = ( static_cast< dfloat >( n_ ) * meany_ + static_cast< dfloat >( other.n_ ) * other.meany_ ) / N;
+         dfloat fN = static_cast< dfloat >( n_ * other.n_ ) / N;
+         m2x_ += other.m2x_ + dx * dx * fN;
+         m2y_ += other.m2y_ + dy * dy * fN;
+         C_ += other.C_ + dx * dy * fN;
+         n_ = intN;
          return *this;
       }
 
diff --git a/src/statistics/statistics.cpp b/src/statistics/statistics.cpp
index e2ac558d..43641079 100644
--- a/src/statistics/statistics.cpp
+++ b/src/statistics/statistics.cpp
@@ -374,7 +374,7 @@ class MaximumAndMinimumLineFilter : public MaximumAndMinimumLineFilterBase {
       dip::uint GetNumberOfOperations( dip::uint /**/, dip::uint /**/, dip::uint /**/ ) override { return 3; }
       void Filter( Framework::ScanLineFilterParameters const& params ) override {
          TPI const* in = static_cast< TPI const* >( params.inBuffer[ 0 ].buffer );
-         MinMaxAccumulator vars;
+         MinMaxAccumulator& vars = accArray_[ params.thread ];
          auto bufferLength = params.bufferLength;
          auto inStride = params.inBuffer[ 0 ].stride;
          if( params.inBuffer.size() > 1 ) {
@@ -401,7 +401,6 @@ class MaximumAndMinimumLineFilter : public MaximumAndMinimumLineFilterBase {
                vars.Push( static_cast< dfloat >( *in ));
             }
          }
-         accArray_[ params.thread ] += vars;
       }
       void SetNumberOfThreads( dip::uint threads ) override {
          accArray_.resize( threads );
@@ -507,7 +506,7 @@ class SampleStatisticsLineFilter : public SampleStatisticsLineFilterBase {
       dip::uint GetNumberOfOperations( dip::uint /**/, dip::uint /**/, dip::uint /**/ ) override { return 23; }
       void Filter( Framework::ScanLineFilterParameters const& params ) override {
          TPI const* in = static_cast< TPI const* >( params.inBuffer[ 0 ].buffer );
-         StatisticsAccumulator vars;
+         StatisticsAccumulator& vars = accArray_[ params.thread ];
          auto bufferLength = params.bufferLength;
          auto inStride = params.inBuffer[ 0 ].stride;
          if( params.inBuffer.size() > 1 ) {
@@ -528,7 +527,6 @@ class SampleStatisticsLineFilter : public SampleStatisticsLineFilterBase {
                in += inStride;
             }
          }
-         accArray_[ params.thread ] += vars;
       }
       void SetNumberOfThreads( dip::uint threads ) override {
          accArray_.resize( threads );
@@ -572,7 +570,7 @@ class CovarianceLineFilter : public CovarianceLineFilterBase {
       void Filter( Framework::ScanLineFilterParameters const& params ) override {
          TPI const* in1 = static_cast< TPI const* >( params.inBuffer[ 0 ].buffer );
          TPI const* in2 = static_cast< TPI const* >( params.inBuffer[ 1 ].buffer );
-         CovarianceAccumulator vars;
+         CovarianceAccumulator& vars = accArray_[ params.thread ];
          auto bufferLength = params.bufferLength;
          auto in1Stride = params.inBuffer[ 0 ].stride;
          auto in2Stride = params.inBuffer[ 1 ].stride;
@@ -596,7 +594,6 @@ class CovarianceLineFilter : public CovarianceLineFilterBase {
                in2 += in2Stride;
             }
          }
-         accArray_[ params.thread ] += vars;
       }
       void SetNumberOfThreads( dip::uint threads ) override {
          accArray_.resize( threads );
@@ -731,7 +728,7 @@ class CenterOfMassLineFilter : public CenterOfMassLineFilterBase {
       dip::uint GetNumberOfOperations( dip::uint /**/, dip::uint /**/, dip::uint /**/ ) override { return nD_ + 1; }
       void Filter( Framework::ScanLineFilterParameters const& params ) override {
          TPI const* in = static_cast< TPI const* >( params.inBuffer[ 0 ].buffer );
-         FloatArray vars( nD_ + 1, 0.0 );
+         FloatArray& vars = accArray_[ params.thread ];
          auto bufferLength = params.bufferLength;
          auto inStride = params.inBuffer[ 0 ].stride;
          UnsignedArray pos = params.position;
@@ -762,7 +759,6 @@ class CenterOfMassLineFilter : public CenterOfMassLineFilterBase {
                ++( pos[ procDim ] );
             }
          }
-         accArray_[ params.thread ] += vars;
       }
       CenterOfMassLineFilter( dip::uint nD ) : nD_( nD ) {}
       void SetNumberOfThreads( dip::uint threads ) override {
@@ -820,7 +816,7 @@ class MomentsLineFilter : public MomentsLineFilterBase {
       }
       void Filter( Framework::ScanLineFilterParameters const& params ) override {
          TPI const* in = static_cast< TPI const* >( params.inBuffer[ 0 ].buffer );
-         MomentAccumulator vars( nD_ );
+         MomentAccumulator& vars = accArray_[ params.thread ];
          auto bufferLength = params.bufferLength;
          auto inStride = params.inBuffer[ 0 ].stride;
          FloatArray pos{ params.position };
@@ -845,7 +841,6 @@ class MomentsLineFilter : public MomentsLineFilterBase {
                ++( pos[ procDim ] );
             }
          }
-         accArray_[ params.thread ] += vars;
       }
       MomentsLineFilter( dip::uint nD ) : nD_( nD ) {}
       void SetNumberOfThreads( dip::uint threads ) override {