dbscan.h

#ifndef DBSCAN_H 
#define DBSCAN_H 
 
#ifdef __cplusplus
#include <iostream>
#include <vector>
#include <map>
#include <math.h>
#include <utility>
#include <cstdint>
#include <stdexcept>
#endif
 
class DBSCAN{
private:
  std::vector<std::pair<double, double>> setOfPoints;
  double Eps;
  int64_t MinPts;
  int64_t cluster_id{0};
  std::map<std::pair<double, double>, int64_t> points;
 
public:
  
  explicit DBSCAN(std::vector<std::pair<double, double>> setOfPoints, double Eps, int64_t MinPts) noexcept : setOfPoints(setOfPoints), Eps(Eps), MinPts(MinPts) {
    // cluster_id is by default noise
    
    for(size_t i = 0; i<setOfPoints.size(); ++i){
      if(points.find(setOfPoints[i]) == points.end()){
        if(ExpandCluster(setOfPoints, setOfPoints[i], cluster_id, Eps, MinPts)){
          cluster_id = nextId(cluster_id);
        }
      }
    }
  }
  
  int64_t nextId(int64_t cluster_id);
 
  bool ExpandCluster(std::vector<std::pair<double, double> > setOfPoints, std::pair<double, double> point, int64_t cluster_id, double Eps, int64_t MinPts);
  
  std::vector<std::pair<double, double> > get_query(std::vector<std::pair<double, double> > setOfPoints, std::pair<double, double> point, double Eps);
  
  double dist(std::pair<double, double> a, std::pair<double, double> b);
  
  std::map<std::pair<double, double>, int64_t> get_clusters();
  
  std::vector<std::pair<double, double> > get_noise();
};
 
int64_t DBSCAN::nextId(int64_t cluster_id){
  try{
    if(cluster_id < -1){
      throw std::logic_error("cluster_id can't be less than -1");
    }
  }
  catch(std::logic_error &e){
    std::cerr << e.what() << '\n';
    return -2;
  }
  cluster_id++;
  return cluster_id;
}
 
bool DBSCAN::ExpandCluster(std::vector<std::pair<double, double> > setOfPoints, std::pair<double, double> point, int64_t cluster_id, double Eps, int64_t MinPts){
  std::vector<std::pair<double, double> > seeds = get_query(setOfPoints, point, Eps);
  if(seeds.size() < MinPts){
    // no core point 
    points[point] = -1;
    return false;
  }
  else{
    // we have a core point
    // all the points in seeds are density-reachable from point
    for(auto & x : seeds){
      points[x] = cluster_id;
    }
    for(auto it = seeds.begin(); it != seeds.end(); it++){
      if((*it).first == point.first && (*it).second == point.second){
        seeds.erase(it);
        break;
      }
    }
 
    while(!seeds.empty()){
      auto current = seeds[0];
      std::vector<std::pair<double, double> > result = get_query(setOfPoints, current, Eps);
 
      if(result.size() >= MinPts){
        for(size_t i = 0; i<result.size(); i++){
          std::pair<double, double> result_p = result[i];
          if(points.find(result_p) == points.end() || points[result_p] == -1){
            if(points.find(result_p) == points.end()){
              seeds.push_back(result_p);
            }
            points[result_p] = cluster_id;
          } // unclassified or noise
        }      
      }
      for(auto it = seeds.begin(); it != seeds.end(); it++){
        if((*it).first == current.first && (*it).second == current.second){
          seeds.erase(it);
          break;
        }
      }
    }
    return true;
  }
  return false;
}
 
std::vector<std::pair<double, double> > DBSCAN::get_query(std::vector<std::pair<double, double> > setOfPoints, std::pair<double, double> point, double Eps){
  std::vector<std::pair<double, double> > ans;
  for(size_t i = 0; i<setOfPoints.size(); i++){
    std::pair<double, double> curr = setOfPoints[i];
    if(dist(point, curr) <= Eps){
      ans.push_back(curr);
    }
  }
  return ans;
}
 
double DBSCAN::dist(std::pair<double, double> a, std::pair<double, double> b){
  return sqrt( pow(b.second - a.second, 2) + pow(b.first - a.first, 2) );
}
 
std::map< std::pair<double, double>, int64_t> DBSCAN::get_clusters(){
  std::map< std::pair<double, double>, int64_t> ans;
  for(size_t i = 0; i<setOfPoints.size(); i++){
    if(points[setOfPoints[i]] != -1){
      ans[setOfPoints[i]] = points[setOfPoints[i]];
    }
  }
  return ans;
}
 
std::vector< std::pair<double, double>> DBSCAN::get_noise(){
  std::vector< std::pair<double, double> > ans;
  for(size_t i = 0; i<setOfPoints.size(); i++){
    if(points[setOfPoints[i]] == -1){
      ans.push_back(setOfPoints[i]);
    }
  }
  return ans;
}
 
 
#endif