meanshift.cpp

#include <iostream>
#include <memory>
#include <vector>
#include <cmath>
#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>


namespace py = pybind11;
typedef double cdt; // coordinate data type, float or double
typedef std::vector<std::pair<cdt,cdt>> ctype; // type for coordinates

class Grouper{

public:
Grouper() : _kernelSize(1.0), _numberOfIterations(10), _lookDistance(1.0), _groupingDistance(1.0){

};
void setKernelSize(double kernelSize){_kernelSize = kernelSize;}
void setNumberOfIterations (int numberOfIterations){_numberOfIterations = numberOfIterations;}
void setLookDistance (double lookDistance){_lookDistance = lookDistance;}
void setGroupingDistance (double groupingDistance){_groupingDistance = groupingDistance;}
void readCoordinates(py::array_t<cdt> array, py::array_t<cdt> farray); 
void run();
void group(py::list l);

private:
double gaussian_kernel(double distance);
 ctype _coordinates;
 std::vector<cdt> _fluxes;
 double _kernelSize;
 int _numberOfIterations;
 double _lookDistance;
 double _groupingDistance;

};


void Grouper::readCoordinates(py::array_t<cdt> array, py::array_t<cdt> farray){  
  auto arraybuf    = array.request();
  auto farraybuf = farray.request();
  cdt* coordinates = (cdt*) arraybuf.ptr;
  cdt* fluxes = (cdt*) farraybuf.ptr;
  if (arraybuf.ndim != 2 || arraybuf.shape[1] != 2){
    py::print("Grouper::readCoordinates: expected 2D array with two columns!");
  }
  if (farraybuf.ndim != 1){
    py::print("Grouper::readCoordinates: expected 1D array!");
  }
  unsigned N = arraybuf.shape[0];
  if (farraybuf.shape[0] != N){
    py::print("Grouper::readCoordinates: array of fluxes does not have the same length as array of coordinates!");
  }
  for (unsigned n=0; n<N; ++n){
    _coordinates.push_back(std::pair<cdt, cdt>(coordinates[2*n],coordinates[2*n+1]));
    _fluxes.push_back(fluxes[n]);  
  }

}

static double euclid_distance(std::pair<cdt, cdt> coordinate1, std::pair<cdt,cdt> coordinate2){
// Euclidian distance between two points
  double result = 0.0;
  double distx = coordinate1.first - coordinate2.first;
  double disty = coordinate1.second - coordinate2.second;
  return sqrt(distx * distx + disty * disty);
}

static std::vector<unsigned> neighbourhood_points(std::pair<cdt,cdt> centroid, ctype coordinates, double max_distance){
std::vector<unsigned> result;