PCL
RobustChauvenetRejection.h
Go to the documentation of this file.
1 // ____ ______ __
2 // / __ \ / ____// /
3 // / /_/ // / / /
4 // / ____// /___ / /___ PixInsight Class Library
5 // /_/ \____//_____/ PCL 2.7.0
6 // ----------------------------------------------------------------------------
7 // pcl/RobustChauvenetRejection.h - Released 2024-06-18T15:48:54Z
8 // ----------------------------------------------------------------------------
9 // This file is part of the PixInsight Class Library (PCL).
10 // PCL is a multiplatform C++ framework for development of PixInsight modules.
11 //
12 // Copyright (c) 2003-2024 Pleiades Astrophoto S.L. All Rights Reserved.
13 //
14 // Redistribution and use in both source and binary forms, with or without
15 // modification, is permitted provided that the following conditions are met:
16 //
17 // 1. All redistributions of source code must retain the above copyright
18 // notice, this list of conditions and the following disclaimer.
19 //
20 // 2. All redistributions in binary form must reproduce the above copyright
21 // notice, this list of conditions and the following disclaimer in the
22 // documentation and/or other materials provided with the distribution.
23 //
24 // 3. Neither the names "PixInsight" and "Pleiades Astrophoto", nor the names
25 // of their contributors, may be used to endorse or promote products derived
26 // from this software without specific prior written permission. For written
27 // permission, please contact info@pixinsight.com.
28 //
29 // 4. All products derived from this software, in any form whatsoever, must
30 // reproduce the following acknowledgment in the end-user documentation
31 // and/or other materials provided with the product:
32 //
33 // "This product is based on software from the PixInsight project, developed
34 // by Pleiades Astrophoto and its contributors (https://pixinsight.com/)."
35 //
36 // Alternatively, if that is where third-party acknowledgments normally
37 // appear, this acknowledgment must be reproduced in the product itself.
38 //
39 // THIS SOFTWARE IS PROVIDED BY PLEIADES ASTROPHOTO AND ITS CONTRIBUTORS
40 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
41 // TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
42 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL PLEIADES ASTROPHOTO OR ITS
43 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
44 // EXEMPLARY OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, BUSINESS
45 // INTERRUPTION; PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; AND LOSS OF USE,
46 // DATA OR PROFITS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
47 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
48 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
49 // POSSIBILITY OF SUCH DAMAGE.
50 // ----------------------------------------------------------------------------
51 
52 #ifndef __PCL_RobustChauvenetRejection_h
53 #define __PCL_RobustChauvenetRejection_h
54 
56 
57 #include <pcl/Defs.h>
58 #include <pcl/Diagnostics.h>
59 
60 #include <pcl/Constants.h>
61 #include <pcl/LinearFit.h>
62 
63 namespace pcl
64 {
65 
66 // ----------------------------------------------------------------------------
67 
82 class PCL_CLASS RobustChauvenetRejection
83 {
84 public:
85 
90 
95 
100  {
101  }
102 
106  RobustChauvenetRejection& operator =( const RobustChauvenetRejection& ) = default;
107 
122  float RejectionLimit() const
123  {
124  return m_rejectionLimit;
125  }
126 
131  void SetRejectionLimit( float r )
132  {
133  PCL_PRECONDITION( r > 0 && r < 1 )
134  m_rejectionLimit = Range( r, 0.0F, 1.0F );
135  }
136 
152  int LargeSampleSize() const
153  {
154  return m_largeSampleSize;
155  }
156 
161  void SetLargeSampleSize( int n )
162  {
163  m_largeSampleSize = Max( 1, n );
164  }
165 
216  template <class C, typename I>
217  void PerformRejection( I& i, I& j, double& mean, double& sigma, C& data ) const
218  {
219  data.Sort();
220  i = I( 0 );
221  j = I( data.Length() );
222 
223  for ( int phase = 0; phase < 3; ++phase )
224  for ( ;; )
225  {
226  switch ( phase )
227  {
228  case 0: // + robustness / - precision
229  mean = Median( data.At( i ), data.At( j ) );
230  sigma = LineFitDeviation( data, i, j, mean );
231  break;
232  case 1: // = robustness / = precision
233  mean = Median( data.At( i ), data.At( j ) );
234  sigma = SampleDeviation( data, i, j, mean );
235  break;
236  case 2: // - robustness / + precision
237  mean = Mean( data.At( i ), data.At( j ) );
238  sigma = StdDev( data.At( i ), data.At( j ), mean );
239  break;
240  }
241 
242  if ( 1 + sigma == 1 )
243  return;
244 
245  I n = j - i;
246  if ( n < I( 3 ) )
247  return;
248 
249  if ( n <= I( m_largeSampleSize ) + I( m_largeSampleSize >> 1 ) )
250  {
251  /*
252  * Optimal single-item rejection iteration for 'small' samples.
253  */
254  double d0 = n*QF( (mean - data[i])/sigma );
255  double d1 = n*QF( (data[j-1] - mean)/sigma );
256  if ( d0 >= m_rejectionLimit )
257  if ( d1 >= m_rejectionLimit )
258  break;
259  if ( d1 < d0 )
260  --j;
261  else
262  ++i;
263  }
264  else
265  {
266  /*
267  * Accelerated bulk rejection for large samples.
268  */
269  int nc = RoundInt( double( n )/m_largeSampleSize );
270  int c = 0;
271  for ( int it = 0; it < nc && n*QF( (mean - data[i])/sigma ) < m_rejectionLimit; ++i, ++c, ++it ) {}
272  for ( int it = 0; it < nc && n*QF( (data[j-1] - mean)/sigma ) < m_rejectionLimit; --j, ++c, ++it ) {}
273  if ( c == 0 )
274  break;
275  }
276  }
277  }
278 
287  template <class C, typename I>
288  void operator()( I& i, I& j, double& mean, double& sigma, C& data ) const
289  {
290  return PerformRejection( i, j, mean, sigma, data );
291  }
292 
293 private:
294 
295  float m_rejectionLimit = 0.5; // 0.5 = Chauvenet rejection criterion
296  int m_largeSampleSize = 20000;
297 
298  /*
299  * The Q-function: the probability for a value from a normal distribution of
300  * being more than x standard deviations from the mean.
301  */
302  static double QF( double x )
303  {
304  return 0.5*(1 - Erf( x/Const<double>::sqrt2() ));
305  }
306 
307  /*
308  * Correction factor for 68.3-percentile deviations to avoid overaggressive
309  * rejection for small samples.
310  */
311  template <typename I>
312  static double FN( I N )
313  {
314  return 1/(1 - 2.9442*Pow( double( N ), -1.073 ));
315  }
316 
317  /*
318  * The 68.3-percentile value from the distribution of absolute deviations.
319  */
320  template <class C, typename I>
321  static double SampleDeviation( const C& R, I i, I j, double m )
322  {
323  int N = int( j - i );
324  if ( N < 2 )
325  return 0;
326  Vector D( N );
327  for ( int k = 0; i < j; ++i, ++k )
328  D[k] = Abs( double( R[i] ) - m );
329  return FN( N ) * *pcl::Select( D.Begin(), D.End(), TruncInt( 0.683*D.Length() ) );
330  }
331 
332  /*
333  * 68.3-percentile deviation by fitting a zero-intercept line to the vector
334  * of absolute differences.
335  */
336  template <class C, typename I>
337  static double LineFitDeviation( const C& R, I i, I j, double m )
338  {
339  int N = int( j - i );
340  int n = int( 0.683*N + 0.317 );
341  if ( n < 8 )
342  return SampleDeviation( R, i, j, m );
343 
344  Vector y( N );
345  int k = 0;
346  for ( I ii = i; ii < j; ++ii, ++k )
347  y[k] = Abs( double( R[ii] ) - m );
348  y.Sort();
349  y = Vector( y.Begin(), n );
350 
351  Vector x( n );
352  for ( int i = 0; i < n; ++i )
353  x[i] = Const<double>::sqrt2() * ErfInv( (i + 1 - 0.317)/N );
354 
355  double s;
356  LinearFit f( x, y );
357  if ( f.IsValid() )
358  s = FN( N ) * f( 1.0 );
359  else
360  s = SampleDeviation( R, i, j, m );
361 
362  return s;
363  }
364 };
365 
366 // ----------------------------------------------------------------------------
367 
368 } // pcl
369 
370 #endif // __PCL_RobustChauvenetRejection_h
371 
372 // ----------------------------------------------------------------------------
373 // EOF pcl/RobustChauvenetRejection.h - Released 2024-06-18T15:48:54Z
Fundamental numeric constants.
Definition: Constants.h:73
static constexpr T sqrt2()
Definition: Constants.h:179
Robust Chauvenet outlier rejection algorithm.
RobustChauvenetRejection(const RobustChauvenetRejection &)=default
void PerformRejection(I &i, I &j, double &mean, double &sigma, C &data) const
void operator()(I &i, I &j, double &mean, double &sigma, C &data) const
64-bit floating point real vector.
Complex< T1 > Pow(const Complex< T1 > &c, T2 x) noexcept
Definition: Complex.h:747
T Abs(const Complex< T > &c) noexcept
Definition: Complex.h:429
constexpr T ErfInv(T x) noexcept
Definition: Math.h:4616
constexpr T Erf(T x) noexcept
Definition: Math.h:4599
int RoundInt(T x) noexcept
Definition: Math.h:1503
int TruncInt(T x) noexcept
Definition: Math.h:1132
RI Select(RI i, RI j, distance_type k)
Definition: Selection.h:165
double StdDev(const T *__restrict__ i, const T *__restrict__ j, double center) noexcept
Definition: Math.h:2813
double Mean(const T *__restrict__ i, const T *__restrict__ j) noexcept
Definition: Math.h:2709
double Median(const T *__restrict__ i, const T *__restrict__ j)
Definition: Math.h:2878
constexpr const T & Range(const T &x, const T &a, const T &b) noexcept
Definition: Utility.h:190
constexpr const T & Max(const T &a, const T &b) noexcept
Definition: Utility.h:119
PCL root namespace.
Definition: AbstractImage.h:77