c++ version of mean shift algorithm

CMakeLists.txt

+cmake_minimum_required(VERSION 3.16)
+project(grouper)      
+
+set(CMAKE_CXX_STANDARD 14)
+
+find_package(Python 3.8 COMPONENTS Interpreter Development REQUIRED)
+find_package(OpenMP)
+
+if(Python_FOUND)
+    message(STATUS "Python Found: ${Python_EXECUTABLE}")
+    message(STATUS "Python Found: ${Python_INCLUDE_DIRS}")
+    message(STATUS "Python Found: ${Python_LIBRARIES}")
+    message(STATUS "Python Found: ${Python_LIBRARY_DIRS}")
+    include_directories(${Python_INCLUDE_DIRS})
+endif()
+
+add_library(grouper   SHARED meanshift.cpp)  
+SET_TARGET_PROPERTIES(grouper   PROPERTIES PREFIX "")
+if(Python_FOUND)
+    target_link_libraries(grouper ${Python_LIBRARIES} OpenMP::OpenMP_CXX)                     
+endif()
+

lsmtool/operations/_meanshiftc.py

+# -*- coding: utf-8 -*-
+#
+# Defines meanshift grouper functions used by the group operation. Based on
+# implementation of Francesco de Gasperin (see
+# https://github.com/revoltek/LiLF/blob/master/lib_dd.py)
+#
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+import numpy as np
+import grouper
+import logging
+log = logging.getLogger('LSMTool.Group')
+
+
+class Grouper(object):
+    """
+    Class to group a list of coordinates and fluxes by the meanshift algorighm.
+    Based on: http://www.chioka.in/meanshift-algorithm-for-the-rest-of-us-python/
+    Parameters
+    ----------
+    coords : list
+        List of (x, y) pixel coordinates for source positions
+    fluxes: list
+        total flux density for each source
+    kernel_size : float
+        Kernel size in pixels
+    n_iterations : int
+        Number of iterations
+    look_distance : float
+        Look distance in pixels
+    grouping_distance : float
+        Grouping distance in pixels
+    """
+    def __init__(self, coords, fluxes, kernel_size, n_iterations, look_distance,
+                 grouping_distance):
+        self.coords = np.array(coords)
+        self.fluxes = fluxes
+        self.kernel_size = kernel_size
+        self.n_iterations = n_iterations
+        self.look_distance = look_distance
+        self.grouping_distance = grouping_distance
+        self.past_coords = [np.copy(self.coords)]
+        self.clusters = []
+        self.g = grouper.Grouper()
+        self.g.readCoordinates(self.coords, self.fluxes)
+        self.g.setKernelSize(self.kernel_size)
+        self.g.setNumberOfIterations(self.n_iterations)
+        self.g.setLookDistance(self.look_distance)
+        self.g.setGroupingDistance(self.grouping_distance)
+
+
+    def run(self):
+        self.g.run()
+
+    def grouping(self):
+        """
+        Take the last coords set and group sources nearby, then return a list of lists.
+        Each list has the index of one cluster.
+        """
+
+        self.g.group(self.clusters)
+
+        log.info('Creating %i groups.' % len(self.clusters))
+        return self.clusters

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;
+  for (unsigned n=0; n<coordinates.size(); ++n){
+    double distance = euclid_distance(centroid, coordinates[n]);
+    if (distance < max_distance){
+      result.push_back(n);
+    }
+  }
+return result;
+}
+
+double Grouper::gaussian_kernel(double distance){
+// Gaussian kernel
+  double result = 1.0 / (_kernelSize * std::sqrt(2.0*3.14159265359));
+  result *= std::exp(-0.5 * distance * distance / (_kernelSize * _kernelSize));
+  return result;
+}
+
+void Grouper::run(){
+  // run the algorithm
+
+  if (_coordinates.size() == 0){
+    py::print("Grouper::run: coordinates have not been read!");
+  }
+
+
+  ctype newcoords = _coordinates;
+
+  for (unsigned it = 0; it < _numberOfIterations; ++it){
+   py::print("starting iteration ", it); 
+#pragma omp parallel
+    for (unsigned n=0; n<_coordinates.size(); ++n){
+      std::vector<unsigned> idx_neighbours = neighbourhood_points(_coordinates[n], _coordinates, _lookDistance);
+      double denominator = 0.0;
+      double numx = 0.0;
+      double numy = 0.0;
+      for (auto i : idx_neighbours){
+         double distance = euclid_distance(_coordinates[i], _coordinates[n]);
+         double weight = gaussian_kernel(distance);
+         weight *= _fluxes[i];
+         numx += weight * _coordinates[i].first;
+         numy += weight * _coordinates[i].second;
+         denominator +=  weight;
+      }
+       newcoords[n] = std::pair<cdt,cdt>(numx/denominator, numy/denominator);
+    }
+    double maxdiff = 0.0;
+    unsigned n=0;
+    for (auto coordinate : _coordinates){
+      double diff = euclid_distance(coordinate, newcoords[n]);
+      if (diff > maxdiff){
+        maxdiff = diff;
+      }
+      ++n; 
+    }
+    _coordinates = newcoords;    
+    if ((it > 1) && (maxdiff < _groupingDistance / 2.0)){
+      break;
+    }
+
+  } 
+
+
+}
+
+void Grouper::group(py::list l){
+
+  ctype coords_to_check = _coordinates;
+
+  while (coords_to_check.size() > 0){
+    auto idx_cluster = neighbourhood_points(coords_to_check[0], _coordinates, _groupingDistance);
+    auto idx_cluster_to_remove = neighbourhood_points(coords_to_check[0], coords_to_check, _groupingDistance);
+    unsigned n = 0;
+    for (auto idx : idx_cluster_to_remove){
+      coords_to_check.erase(coords_to_check.begin() + idx - n);
+      ++n;
+    }
+
+    py::list newcluster;
+    for (auto idx : idx_cluster){
+      newcluster.append(py::int_(idx));
+    }
+    l.append(newcluster);
+
+  }
+}
+
+PYBIND11_MODULE(grouper, m)    
+{
+  py::class_<Grouper>(m, "Grouper")
+    .def(py::init<>())
+    .def("run", &Grouper::run)
+    .def("setKernelSize", &Grouper::setKernelSize)
+    .def("setNumberOfIterations", &Grouper::setNumberOfIterations)
+    .def("setLookDistance", &Grouper::setLookDistance)
+    .def("setGroupingDistance", &Grouper::setGroupingDistance)
+    .def("readCoordinates", &Grouper::readCoordinates)  
+    .def("group", &Grouper::group)
+;
+}
+
+
+int main(){
+
+  Py_Initialize();
+
+  return 0;
+}
+

test.py

+import grouper
+import numpy as np
+
+a = np.array([[1.0, 2.0],[3.0,4.0],[5.0,6.0]])
+f = np.array([1.0,2.0,3.0])
+
+g = grouper.Grouper()
+
+g.readCoordinates(a,f)
+
+g.setLookDistance(5.0)
+g.setGroupingDistance(5.0)
+
+g.run()
+b = []
+g.group(b)
+
+print(b)

testlsm.py

+import lsmtool
+
+s = lsmtool.load('test_large.sky')
+
+s.group('meanshift', byPatch=True,lookDistance=0.075,groupingDistance=0.01)
+