diff --git a/README.md b/README.md
index 5c9ae705e61ab28b2b962f3a99c4775c7cc5f74c..829c38f84ce4fad9f10d4273a4f4d6fa569bd6b6 100644
--- a/README.md
+++ b/README.md
@@ -64,6 +64,12 @@ Then install with:
     cd LSMTool
     python setup.py install
 
+If you have a C++11-compliant compiler, you can build a faster
+version of the mean shift grouping algorithm with:
+
+    cd LSMTool
+    python setup.py install --build_c_extentions
+
 ### Testing
 
 You can test that the installation worked with:
diff --git a/meanshift.cpp b/lsmtool/operations/_grouper.cpp
similarity index 95%
rename from meanshift.cpp
rename to lsmtool/operations/_grouper.cpp
index 10ef88c672064d943867a590bd854dea3dbdc507..dbac7b5c05ed4c543411245a619bfb4ed01afbb6 100644
--- a/meanshift.cpp
+++ b/lsmtool/operations/_grouper.cpp
@@ -20,7 +20,7 @@ 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 readCoordinates(py::array_t<cdt> array, py::array_t<cdt> farray);
 void run();
 void group(py::list l);
 
@@ -36,7 +36,7 @@ double gaussian_kernel(double distance);
 };
 
 
-void Grouper::readCoordinates(py::array_t<cdt> array, py::array_t<cdt> farray){  
+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;
@@ -53,7 +53,7 @@ void Grouper::readCoordinates(py::array_t<cdt> array, py::array_t<cdt> farray){
   }
   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]);  
+    _fluxes.push_back(fluxes[n]);
   }
 
 }
@@ -96,7 +96,7 @@ void Grouper::run(){
 
   for (unsigned it = 0; it < _numberOfIterations; ++it){
    py::print("starting iteration ", it);
-#pragma omp parallel for 
+#pragma omp parallel for
     for (unsigned n=0; n<_coordinates.size(); ++n){
       std::vector<unsigned> idx_neighbours = neighbourhood_points(_coordinates[n], _coordinates, _lookDistance);
       double denominator = 0.0;
@@ -119,14 +119,14 @@ void Grouper::run(){
       if (diff > maxdiff){
         maxdiff = diff;
       }
-      ++n; 
+      ++n;
     }
-    _coordinates = newcoords;    
+    _coordinates = newcoords;
     if ((it > 1) && (maxdiff < _groupingDistance / 2.0)){
       break;
     }
 
-  } 
+  }
 
 
 }
@@ -158,7 +158,7 @@ void Grouper::group(py::list l){
   }
 }
 
-PYBIND11_MODULE(grouper, m)    
+PYBIND11_MODULE(_grouper, m)
 {
   py::class_<Grouper>(m, "Grouper")
     .def(py::init<>())
@@ -167,7 +167,7 @@ PYBIND11_MODULE(grouper, m)
     .def("setNumberOfIterations", &Grouper::setNumberOfIterations)
     .def("setLookDistance", &Grouper::setLookDistance)
     .def("setGroupingDistance", &Grouper::setGroupingDistance)
-    .def("readCoordinates", &Grouper::readCoordinates)  
+    .def("readCoordinates", &Grouper::readCoordinates)
     .def("group", &Grouper::group)
 ;
 }
diff --git a/lsmtool/operations/_meanshiftc.py b/lsmtool/operations/_meanshiftc.py
index 901579e01a4b641619e7672be8aa8e778dc9b6c9..56b85843370cc6d84814ae5d585245ee3ae57280 100644
--- a/lsmtool/operations/_meanshiftc.py
+++ b/lsmtool/operations/_meanshiftc.py
@@ -18,7 +18,7 @@
 # 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
+from . import _grouper
 import logging
 log = logging.getLogger('LSMTool.Group')
 
@@ -52,7 +52,7 @@ class Grouper(object):
         self.grouping_distance = grouping_distance
         self.past_coords = [np.copy(self.coords)]
         self.clusters = []
-        self.g = grouper.Grouper()
+        self.g = _grouper.Grouper()
         self.g.readCoordinates(self.coords, self.fluxes)
         self.g.setKernelSize(self.kernel_size)
         self.g.setNumberOfIterations(self.n_iterations)
diff --git a/lsmtool/operations/group.py b/lsmtool/operations/group.py
index 91fa0bd5e908ebd466b8496c7d40e05c2cdae693..fcc9977f6ca3ab4eec4e7603d766b133b073c125 100644
--- a/lsmtool/operations/group.py
+++ b/lsmtool/operations/group.py
@@ -150,7 +150,10 @@ def group(LSM, algorithm, targetFlux=None, weightBySize=False, numClusters=100,
     from . import _tessellate
     from . import _cluster
     from . import _threshold
-    from . import _meanshift
+    try:
+        from . import _meanshiftc as _meanshift
+    except ImportError:
+        from . import _meanshift
     import numpy as np
     import os
     from itertools import groupby
diff --git a/setup.py b/setup.py
index 92b5368477d0901afd7af8643f5eb4f34538a8cb..61335ebbfd79668a7cc11a7a533c3d02c938c0cd 100644
--- a/setup.py
+++ b/setup.py
@@ -1,12 +1,35 @@
 from __future__ import print_function
-from setuptools import setup, Command
-import os
+from setuptools import setup, Command, Extension, Distribution
+from setuptools.command.build_ext import build_ext
 import sys
 import lsmtool._version
 
 
+# Flag that determines whether to build the optional (but faster) C++
+# extensions. Set to False to install only the pure Python versions.
+if "--build_c_extentions" in sys.argv:
+    build_c_extentions = True
+    sys.argv.remove("--build_c_extentions")
+else:
+    build_c_extentions = False
+
+# Handle Python 3-only dependencies
+if sys.version_info < (3, 0):
+    reqlist = ['numpy', 'astropy >= 0.4, <3.0']
+else:
+    reqlist = ['numpy', 'astropy >= 0.4']
+if build_c_extentions:
+    reqlist.append('pybind11>=2.2.0')
+    ext_modules = [Extension('lsmtool.operations._grouper',
+                             ['lsmtool/operations/_grouper.cpp'],
+                             language='c++')]
+else:
+    ext_modules = []
+
+
 class PyTest(Command):
     user_options = []
+
     def initialize_options(self):
         pass
 
@@ -14,15 +37,36 @@ class PyTest(Command):
         pass
 
     def run(self):
-        import sys,subprocess
+        import sys
+        import subprocess
         errno = subprocess.call([sys.executable, 'runtests.py'])
         raise SystemExit(errno)
 
-# Handle Python 3-only dependencies
-if sys.version_info < (3, 0):
-    reqlist = ['numpy','astropy >= 0.4, <3.0']
-else:
-    reqlist = ['numpy','astropy >= 0.4']
+
+class LSMToolDistribution(Distribution):
+
+    def is_pure(self):
+        if self.pure:
+            return True
+
+    def has_ext_modules(self):
+        return not self.pure
+
+    global_options = Distribution.global_options + [
+        ('pure', None, "use pure Python code instead of C++ extensions")]
+    pure = False
+
+
+class BuildExt(build_ext):
+
+    def build_extensions(self):
+        opts = ['-std=c++11']
+        if sys.platform == 'darwin':
+            opts += ['-stdlib=libc++']
+        for ext in self.extensions:
+            ext.extra_compile_args = opts
+        build_ext.build_extensions(self)
+
 
 setup(
     name='lsmtool',
@@ -41,7 +85,10 @@ setup(
         ],
     install_requires=reqlist,
     scripts=['bin/lsmtool'],
-    packages=['lsmtool','lsmtool.operations'],
+    ext_modules=ext_modules,
+    cmdclass={'build_ext': BuildExt},
+    distclass=LSMToolDistribution,
+    packages=['lsmtool', 'lsmtool.operations'],
     setup_requires=['pytest-runner'],
     tests_require=['pytest']
     )