RANSAC.h

Go to the documentation of this file.

// Copyright (c) 2002-2023, The OpenMS Team -- EKU Tuebingen, ETH Zurich, and FU Berlin
// SPDX-License-Identifier: BSD-3-Clause
//
// --------------------------------------------------------------------------
// $Maintainer: George Rosenberger $
// $Authors: George Rosenberger, Hannes Roest, Chris Bielow $
// --------------------------------------------------------------------------
 
#pragma once
 
#include <OpenMS/config.h>
 
#include <OpenMS/MATH/MISC/RANSACModel.h>
 
#include <OpenMS/CONCEPT/Exception.h>
#include <OpenMS/DATASTRUCTURES/String.h>
#include <OpenMS/MATH/MISC/RANSACModelLinear.h>
#include <OpenMS/MATH/MISC/MathFunctions.h>
 
#include <limits>       // std::numeric_limits
#include <vector>       // std::vector
#include <sstream>      // stringstream
 
namespace OpenMS
{
 
  namespace Math
  {
    struct RANSACParam
    {
      RANSACParam()
        : n(0), k(0), t(0), d(0), relative_d(false)
        {
        }
      RANSACParam(size_t p_n, size_t p_k, double p_t, size_t p_d, bool p_relative_d = false)
        : n(p_n), k(p_k), t(p_t), d(p_d), relative_d(p_relative_d)
      {
        if (relative_d)
        {
          if (d >= 100) throw Exception::Precondition(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, String("RANSAC: Relative 'd' >= 100% given. Use a lower value; the more outliers you expect, the lower it should be."));
        }
      }
 
      [[nodiscard]] std::string toString() const
      {
        std::stringstream r;
        r << "RANSAC param:\n  n: " << n << "\n  k: " << k << " iterations\n  t: " << t << " threshold\n  d: " << d << " inliers\n\n";
        return r.str();
      }
 
      size_t n; 
      size_t k; 
      double t; 
      size_t d; 
      bool relative_d; 
    };
 
    template<typename TModelType = RansacModelLinear>
    class RANSAC
    {
public:
 
      explicit RANSAC(uint64_t seed = time(nullptr)):
      shuffler_(seed)
      {}
 
      ~RANSAC() = default;
 
 
      void setSeed(uint64_t seed)
      {
        shuffler_.seed(seed);
      }
 
      std::vector<std::pair<double, double> > ransac(
        const std::vector<std::pair<double, double> >& pairs,
        const RANSACParam& p)
      {
        return ransac(pairs, p.n, p.k, p.t, p.d, p.relative_d);
      }
 
      std::vector<std::pair<double, double> > ransac(
          const std::vector<std::pair<double, double> >& pairs, 
          size_t n, 
          size_t k, 
          double t, 
          size_t d, 
          bool relative_d = false)
      {
        // translate relative percentages into actual numbers
        if (relative_d)
        {
          if (d >= 100) throw Exception::Precondition(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, String("RANSAC: Relative 'd' >= 100% given. Use a lower value; the more outliers you expect, the lower it should be."));
          d = pairs.size() * d / 100;
        }
 
        // implementation of the RANSAC algorithm according to http://wiki.scipy.org/Cookbook/RANSAC.
 
        if (pairs.size() <= n)
        {
          throw Exception::Precondition(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION,
                                        String("RANSAC: Number of total data points (") + String(pairs.size()) + ") must be larger than number of initial points (n=" + String(n) + ").");
        }
 
        TModelType model;
 
        std::vector< std::pair<double, double> > alsoinliers, betterdata, bestdata;
        std::vector<std::pair<double, double> > pairs_shuffled = pairs;  // mutable data. will be shuffled in every iteration
        double besterror = std::numeric_limits<double>::max();
        typename TModelType::ModelParameters coeff;
    #ifdef DEBUG_RANSAC
        std::pair<double, double > bestcoeff;
        double betterrsq = 0;
        double bestrsq = 0;
    #endif
 
        for (size_t ransac_int=0; ransac_int<k; ransac_int++)
        {
          // check if the model already includes all points
          if (bestdata.size() == pairs.size()) break;
 
          // use portable RNG in test mode
          shuffler_.portable_random_shuffle(pairs_shuffled.begin(), pairs_shuffled.end());
 
          // test 'maybeinliers'
          try
          { // fitting might throw UnableToFit if points are 'unfortunate'
            coeff = model.rm_fit(pairs_shuffled.begin(), pairs_shuffled.begin()+n);
          }
          catch (...)
          {
            continue;
          }
          // apply model to remaining data; pick inliers
          alsoinliers = model.rm_inliers(pairs_shuffled.begin()+n, pairs_shuffled.end(), coeff, t);
          // ... and add data
          if (alsoinliers.size() > d 
              || alsoinliers.size() >= (pairs_shuffled.size()-n)) // maximum number of inliers we can possibly have (i.e. remaining data)
          {
            betterdata.clear();
            std::copy( pairs_shuffled.begin(), pairs_shuffled.begin()+n, back_inserter(betterdata) );
            betterdata.insert( betterdata.end(), alsoinliers.begin(), alsoinliers.end() );
            typename TModelType::ModelParameters bettercoeff = model.rm_fit(betterdata.begin(), betterdata.end());
            double bettererror = model.rm_rss(betterdata.begin(), betterdata.end(), bettercoeff);
    #ifdef DEBUG_RANSAC
            betterrsq = model.rm_rsq(betterdata);
    #endif
 
            // If the current model explains more points, we assume its better (these points pass the error threshold 't', so they should be ok);
            // If the number of points is equal, we trust rss.
            // E.g. imagine gaining a zillion more points (which pass the threshold!) -- then rss will automatically be worse, no matter how good
            //      these points fit, since its a simple absolute SUM() of residual error over all points.
            if (betterdata.size() > bestdata.size() || (betterdata.size() == bestdata.size() && (bettererror < besterror)))
            {
              besterror = bettererror;
              bestdata = betterdata;
    #ifdef DEBUG_RANSAC
              bestcoeff = bettercoeff;
              bestrsq = betterrsq;
              std::cout << "RANSAC " << ransac_int << ": Points: " << betterdata.size() << " RSQ: " << bestrsq << " Error: " << besterror << " c0: " << bestcoeff.first << " c1: " << bestcoeff.second << std::endl;
    #endif
            }
          }
        }
 
    #ifdef DEBUG_RANSAC
        std::cout << "=======STARTPOINTS=======" << std::endl;
        for (std::vector<std::pair<double, double> >::iterator it = bestdata.begin(); it != bestdata.end(); ++it)
        {
          std::cout << it->first << "\t" << it->second << std::endl;
        }
        std::cout << "=======ENDPOINTS=======" << std::endl;
    #endif
 
        return(bestdata);
      } // ransac()
 
    private:
      Math::RandomShuffler shuffler_{};
    }; // class
  
  } // namespace Math
 
 
} // namespace OpenMS