diff --git a/dlib/image_processing/scan_image.h b/dlib/image_processing/scan_image.h
index 1a9c46eda..6cc74cfd3 100644
--- a/dlib/image_processing/scan_image.h
+++ b/dlib/image_processing/scan_image.h
@@ -11,6 +11,7 @@
 #include "../rand.h"
 #include "../array2d.h"
 #include "../image_transforms/spatial_filtering.h"
+#include "../image_transforms/thresholding.h"
 
 namespace dlib
 {
@@ -235,6 +236,143 @@ namespace dlib
         }
     }
 
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename image_type
+        >
+    std::vector<point> find_peaks (
+        const image_type& img_,
+        const double non_max_suppression_radius,
+        const typename pixel_traits<typename image_traits<image_type>::pixel_type>::basic_pixel_type& thresh
+    )
+    {
+        DLIB_CASSERT(non_max_suppression_radius >= 0);
+        const_image_view<image_type> img(img_);
+
+        using basic_pixel_type = typename pixel_traits<typename image_traits<image_type>::pixel_type>::basic_pixel_type;
+
+        std::vector<std::pair<basic_pixel_type,point>> peaks;
+
+        for (long r = 1; r+1 < img.nr(); ++r)
+        {
+            for (long c = 1; c+1 < img.nc(); ++c)
+            {
+                auto val = img[r][c];
+                if (val < thresh)
+                    continue;
+
+                if (
+                    val <= img[r-1][c] ||
+                    val <= img[r+1][c] ||
+                    val <= img[r][c+1] ||
+                    val <= img[r][c-1] ||
+                    val <= img[r-1][c-1] ||
+                    val <= img[r+1][c+1] ||
+                    val <= img[r-1][c+1] ||
+                    val <= img[r+1][c-1]
+                )
+                {
+                    continue;
+                }
+
+                peaks.emplace_back(val,point(c,r));
+            }
+        }
+
+
+        // now do non-max suppression of the peaks according to the supplied radius.
+        using pt = std::pair<basic_pixel_type,point>;
+        // First sort the peaks so the strongest peaks come first.  We will greedily accept
+        // them and then do the normal peak sorting/non-max suppression thing.
+        std::sort(peaks.rbegin(), peaks.rend(), [](pt& a, pt&b ){ return a.first < b.first; });
+        std::vector<point> final_peaks;
+        const double radius_sqr = non_max_suppression_radius*non_max_suppression_radius;
+
+        // If there are a lot of peaks then we will make a mask image and use that to do
+        // the non-max suppression since this is fast when peaks.size() is large.  Otherwise we
+        // will do the simpler thing in the else block that doesn't require us to allocate a
+        // temporary mask image.
+        if (peaks.size() > 500 && radius_sqr != 0)
+        {
+            // hit will record which areas of the image have already been accounted for by some
+            // peak.  So it is our mask image.
+            matrix<unsigned char> hit(img.nr(), img.nc());
+            // initially nothing has been hit.
+            hit = 0;
+            const unsigned long win_size = std::round(2*non_max_suppression_radius);
+            const rectangle area = get_rect(img);
+            for (auto& pp : peaks)
+            {
+                auto& p = pp.second;
+                if (!hit(p.y(),p.x()))
+                {
+                    final_peaks.emplace_back(p);
+
+                    // mask out a circle around this new peak
+                    rectangle win = centered_rect(p, win_size, win_size).intersect(area); 
+                    for (long r = win.top(); r <= win.bottom(); ++r)
+                    {
+                        for (long c = win.left(); c <= win.right(); ++c)
+                        {
+                            if (length_squared(point(c,r)-p) <= radius_sqr)
+                                hit(r,c) = 1;
+                        }
+                    }
+                }
+            }
+        }
+        else
+        {
+            // if peaks.size() is relatively small then this is a faster way to do the non-max
+            // suppression.
+            for (auto& p : peaks)
+            {
+                bool hits_any_existing_peak = false;
+                // If the user set the radius to 0 then just copy the peaks to the output without
+                // checking anything.
+                if (radius_sqr != 0)
+                {
+                    for (auto& v : final_peaks)
+                    {
+                        if (length_squared(p.second-v) <= radius_sqr)
+                        {
+                            hits_any_existing_peak = true;
+                            break;
+                        }
+                    }
+                }
+                if (!hits_any_existing_peak)
+                {
+                    final_peaks.emplace_back(p.second);
+                }
+            }
+        }
+
+        return final_peaks;
+    }
+
+    template <
+        typename image_type
+        >
+    std::vector<point> find_peaks (
+        const image_type& img
+    )
+    {
+        return find_peaks(img, 0, partition_pixels(img));
+    }
+
+    template <
+        typename image_type
+        >
+    std::vector<point> find_peaks (
+        const image_type& img,
+        const double non_max_suppression_radius
+    )
+    {
+        return find_peaks(img, non_max_suppression_radius, partition_pixels(img));
+    }
+
 // ----------------------------------------------------------------------------------------
 
     template <
diff --git a/dlib/image_processing/scan_image_abstract.h b/dlib/image_processing/scan_image_abstract.h
index fe2fc51ac..6ff377fa3 100644
--- a/dlib/image_processing/scan_image_abstract.h
+++ b/dlib/image_processing/scan_image_abstract.h
@@ -132,6 +132,62 @@ namespace dlib
                 - #dets == all the points which passed the threshold test.
     !*/
 
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename image_type
+        >
+    std::vector<point> find_peaks (
+        const image_type& img,
+        const double non_max_suppression_radius,
+        const typename pixel_traits<typename image_traits<image_type>::pixel_type>::basic_pixel_type& thresh
+    );
+    /*!
+        requires
+            - image_type == an image object that implements the interface defined in
+              dlib/image_processing/generic_image.h.  Moreover, these it must contain a
+              scalar pixel type (e.g. int rather than rgb_pixel)
+            - non_max_suppression_radius >= 0
+        ensures
+            - Scans the given image and finds all pixels with values >= thresh that are
+              also local maximums within their 8-connected neighborhood of the image.  Such
+              pixels are collected, sorted in decreasing order of their pixel values, and
+              then non-maximum suppression is applied to this list of points using the
+              given non_max_suppression_radius.  The final list of peaks is then returned.
+
+              Therefore, the returned list, V, will have these properties:
+                - V.size() == the number of peaks found in the image.
+                - When measured in image coordinates, no elements of V are within
+                  non_max_suppression_radius distance of each other.  That is, for all valid i!=j
+                  it is true that length(V[i]-V[j]) > non_max_suppression_radius.
+                - For each element of V, that element has the maximum pixel value of all
+                  pixels in the ball centered on that pixel with radius
+                  non_max_suppression_radius.
+    !*/
+
+    template <
+        typename image_type
+        >
+    std::vector<point> find_peaks (
+        const image_type& img
+    );
+    /*!
+        ensures
+            - performs: return find_peaks(img, 0, partition_pixels(img))
+    !*/
+
+    template <
+        typename image_type
+        >
+    std::vector<point> find_peaks (
+        const image_type& img,
+        const double non_max_suppression_radius
+    );
+    /*!
+        ensures
+            - performs: return find_peaks(img, non_max_suppression_radius, partition_pixels(img))
+    !*/
+
 // ----------------------------------------------------------------------------------------
 
     template <