How to detect light spots in image?

Question

In order to detect inselbergs from a geographic region, I've downloaded topographic imagery (relief and declivity) from that region. The declivity images seems best suited for the task.

After applying a gaussian blur (or a much faster common blur, with ImageMagick), the image seems ready for automatic detection.

Now I'm wondering the best/fastest way to detect these white stains on the black background. My first idea is to use a simple function (no external library) that works like the "bucket paint" function of drawing programs, calculating the area of an object lighter than a predefined threshold. Problem is: since the images are quite large, 5400x3600 pixels, the usual recursion function will probably cause a stack overflow, specially on massive mountain ranges like below.

So, any suggestions? I think the ideal language may be C++ (perhaps JavaScript). I'm not used to Python. Maybe it gets easier using a library like OpenCV, but maybe the problem is too trivial to ask for an external library (including the implied learning curve).

---

The TIFF images came from here (I can convert them to other format for processing). The example image is at the 17S42 quadricule ("Declividade" option), close to the coordinates 18°S 41°W. In the middle image (above), each "white ball" is an inselberg. The precision will depend on the chosen threshold of gray.

Answer 1

In short, you'll get a very poor estimate of size using your method. The edges of the blurred regions, no matter what threshold you pick, do not correspond to the edges of the inselbergs you want to measure. Instead, I suggest you follow the recipe below. I'm using DIPlib in Python for this (disclaimer: I'm an author). The Python bindings are a thin layer on the C++ library, it's fairly simple to translate the Python code below to C++ (it's just easier for me to develop it interactively in Python). Install with pip install diplib.

I downloaded the original height data from the link you provided (rather than the declivity). DIPlib can directly read floating-point valued TIFF files, so there is no need for any special conversion. I cropped a region similar to the one used by OP for this demo, but there is no reason to not apply the method to the whole tile.

import diplib as dip

height = dip.ImageRead('17S42_ZN.tif')
height.SetPixelSize(0.000278, 'rad')  # not really radian, but we don't have degrees
height = height[3049:3684, 2895:3513];

The code also sets the pixel size according to data in the TIFF file (using units of radian, because DIPlib doesn't do degrees).

Next, I apply the top-hat filter with a specific diameter (25 pixels). This isolates all peaks with a diameter of 25 pixels or less. Adjust this size according to what you think the maximum width an inselberg should be.

local_height = dip.Tophat(height, 25)

In effect, the result is a local height, the height above some baseline determined by the size of the filter.

Next, I apply a hysteresis threshold (double threshold). This yields a binary image, thresholded at 100m above the baseline, where the terrain goes above 200m above that baseline. That is, I decided that an inselberg should be at least 200m above the baseline, but cut each of them off at 100m. At this height we'll be measuring the size (area). Again, adjust the thresholds as you see fit.

inselbergs = dip.HysteresisThreshold(local_height, 100, 200)

Now all that is left is measuring the regions we found:

labels = dip.Label(inselbergs)
result = dip.MeasurementTool.Measure(labels, features=['Size', 'Center'])
print(result)

This outputs:

   |       Size |                  Center | 
-- | ---------- | ----------------------- | 
   |            |       dim0 |       dim1 | 
   |     (rad²) |      (rad) |      (rad) | 
-- | ---------- | ---------- | ---------- | 
 1 |  1.863e-05 |     0.1514 |    0.01798 | 
 2 |  4.220e-05 |     0.1376 |    0.02080 | 
 3 |  6.214e-05 |    0.09849 |    0.04429 | 
 4 |  6.492e-06 |     0.1282 |    0.04710 | 
 5 |  3.022e-05 |     0.1354 |    0.04925 | 
 6 |  4.274e-05 |     0.1510 |    0.05420 | 
 7 |  2.218e-05 |     0.1228 |    0.05802 | 
 8 |  1.932e-05 |     0.1420 |    0.05689 | 
 9 |  7.690e-05 |     0.1493 |    0.06960 | 
10 |  3.285e-05 |     0.1120 |    0.07089 | 
11 |  5.248e-05 |     0.1389 |    0.07851 | 
12 |  4.637e-05 |     0.1096 |    0.09016 | 
13 |  3.787e-05 |    0.07146 |     0.1012 | 
14 |  2.133e-05 |    0.09046 |    0.09908 | 
15 |  3.895e-05 |    0.08553 |     0.1064 | 
16 |  3.308e-05 |    0.09972 |     0.1143 | 
17 |  3.277e-05 |    0.05312 |     0.1174 | 
18 |  2.581e-05 |    0.07298 |     0.1167 | 
19 |  1.955e-05 |    0.04038 |     0.1304 | 
20 |  4.846e-05 |    0.03657 |     0.1448 | 

(Remember, where it says 'rad' it is really degrees.) An area in square degrees is a bit weird, but you can convert this to square meters since you know the location on the globe. It might in fact be easier to translate the pixel sizes to meters before the computations.

The values given here for 'Center' are with respect to the top-left pixel, if we hadn't cropped the tile to begin with, we could have added the coordinates for the tile (as can be obtained from the corresponding tag in the TIFF file): (-42.0, -17.0).

---

In C++ the code should look like this:

#include <diplib/simple_file_io.h>
#include <diplib/morphology.h>
#include <diplib/segmentation.h>
#include <diplib/regions.h>
#include <diplib/measurement.h>

//...

dip::Image height = dip::ImageRead("17S42_ZN.tif");
height.SetPixelSize(0.000278 * dip::Units::Radian());
height = height.At(dip::Range(3049, 3684), dip::Range(2895, 3513));

dip::Image local_height = dip::Tophat(height, 25);

dip::Image inselbergs = dip::HysteresisThreshold(local_height, 100, 200);

dip::Image labels = dip::Label(inselbergs);
dip::MeasurementTool measurementTool;
dip::Measurement result = measurementTool.Measure(labels, {}, {"Size", "Center"});
std::cout << result;