How to detach overlapping contours but not to remove small ones. Alternative to erosion-dilation cycle

Question

Background:
I'm trying to detect trees on a sequence of aerial images by calculating a stereo disparity for consequent pairs of images. This is how disparity looks like. White spots are trees:

Next, I am locating contours on the disparity map and getting the following image:

From here I would like to detect the trees.

What I'm doing right now:
First, I'm doing an erosion:

kernel = np.ones(shape=(5, 5),
                 dtype=np.uint8)

image= cv2.erode(image,
                 kernel=kernel,
                 iterations=14)  

and getting some overlapping, e.g.:

In order to separate those overlapping contours, I do dilation this time:

image= cv2.dilate(image,
                  kernel=kernel,
                  iterations=12)

Ups, some contours became so small and almost disappeared.
Reversing damage:

image= cv2.erode(image,
                 kernel=kernel,
                 iterations=7)

Ok! Better! But.. all this process seems like a dirty workaround.. Moreover, this won't work for some other images so well (ignore the noise on the left, I should add cropping):

Questions:
- How do I make sure not to overlap the contours of the trees and not to remove them while dilating/eroding?
- Is there a better alternative to what I am doing right now?

Edit:
Probably for the full picture I should add that after obtaining a good image with separate contours, I use:

disparity, contours, hierarchy = cv2.findContours(image,
                                                  mode=cv2.RETR_LIST,
                                                  method=cv2.CHAIN_APPROX_SIMPLE)

and plot the contours on the original image:

Edit Nº2:
Adding cross-section of disparity map. Raw data available from here

Plot:

---

Update on results:

I tried a median filter as it was suggested by @A_A, and got a very good separation of trees from the ground:

ground = medfilt(disparity, kernel_size=125)
trees = disparity - ground
boolean_mask = (trees > 2.) & (trees < 15.)
plt.imshow(boolean_mask, cmap='gray')

The problem with this filter is the speed.
scipy.signal.medfilt took me more than 20 mins to calculate everything. But scipy has also a scipy.ndimage.filters.median_filter. I don't know what's the difference between them but results that I got from it were almost the same apart from some areas close to image borders, and it took somewhere from 2 to 5 mins.

Answer 1

I don't know if this is a corner case or the norm in your dataset but it is a relatively easy situation to deal with. It would be much more difficult to detect trees in an urban environment, for example, where you would probably be looking at more pattern-recognition type of solutions.

The basic idea behind this solution is surface fitting, whether to directly detect the trees or detect the "floating" level of the ground and through that information detect the trees.

The following is based on Octave and using the 1D data but generalisations to other platforms (and 2 dimensions) should be straightforward.

The first thing one can do here is kernel density estimation (KDE). In this case, you would be approximating the total available surface (both the ground and the trees) via a set of Gaussian curves. In the two dimensional case, this will give you either one Gaussian per tree or higher densities of Gaussians around the locations of trees. In one shot. No further processing required. In reality, you might find that this technique approximates your disparity image via kernel density and returns the data of a set of Gaussians that you can process further to determine where there are trees.

But, having had a quick look at the data, it contains sharp discontinuities and that might not be good for KDE convergence. While you could "blur" the disparity with a Gaussian low pass and then try to approximate it with KDE, this means that your original resolution is decreased and it might reduce your ability to discriminate between small and large trees. Anyway, this is an option too (first blur with a suitable FWHM Gaussian, then approximate with KDE).

Again, this seems a relatively easy case where you could even use simple thresholding to detect the trees. The problem there of course is that the ground height now modulates your threshold. A tree on flat ground could be anything above 25 meters but certainly not a 3 meter tall shed on top of a 22 meter hill.

So, we want to threshold but we want to take into account the ground contour. Ideally, you could fit a curve that is closer to what the ground really is (an ellipsoid) but maybe you can get away with simpler fitting schemes.

The main idea here is that we want to "undo" the modulation of the ground to then apply a simple threshold. If we flatten the ground, the 22 meter hill will disappear and then anything above 22 meters is more likely to really be a tree.

If you ignore the zeros (which I am assuming are "missing values" in your case (either high uncertainty or actually missing values)), the signal looks like this:

You can threshold this at $\approx 62$ but you can see that the second tree is lower than the first but further up the slope of what looks like a hill.

Let's flatten the hill. Flattening the hill here is basically equivalent to subtracting the slowly varying component (the hill) from the signal itself and possibly even controlling the amount of removal with a coefficient. Here, I will use simple polynomial fitting but anything that could approximate the slow varying component of the height field would do.

xrange=0:(length(b)-1);
sv = polyval(polyfit(xrange,b,2),xrange); # Fit and evaluate a second order polynomial. You might want to go higher if you have ridges or depressions but not too high because we just want to approximate the slow varying ground changes.
plot(b);hold on;plot(sv,'r');grid on; # Plot everything

This now looks like this:

Obviously, we can simply subtract the red curve from the blue and get back the "flattened" version of the signal, but as you can see, the trees have "pulled the ground up" slightly and this is why if you do this blindingly, you will get this:

In terms of comparing between the trees, we are now definitely in a better condition because we have removed the ground and we can also tell that the second tree is lower than the first, but, we now have a dip from trying to remove the hill. To adjust for situations like this, you can introduce a coefficient by which to control how much subtraction you want to perform, for example:

m=0.43;plot(m.*b-(1-m).*sv);

Here, $m \in [0 . . 1.0]$ controls how much of the real signal and how much of its approximation we are taking into account. I determined 0.43 after trial and error and this results to:

Where basically, the trees start at almost the same height.

So, now you can threshold anything above -8.2 and be more certain that you are selecting trees.

There is a way by which you can flatten the ground almost completely but it requires some knowledge about the size of the tree canopy.

This works on the same premise but approximates the ground using a different method.

Basically, you can use a median filter to "chop the peaks" that the trees are causing to the height field. Anything that remains after that is basically ground variation. The median filter is applied in a "rolling window" way and substitutes a sample with the median value of $N$ values surrounding it.

Now, to apply this, you need to determine that $N$ which should be almost as wide as your widest tree canopy. I have determined this to be approximately 125 (pixels) in this example, so by doing:

plot(b-medfilt1(b,125)); grid on;

You get something like this:

Which is almost too good. At least the trees start at exactly the same point. Again, it is now the point of applying a simple threshold (anything above 1 (?)) to detect where the trees are.

This process will give you a binary mask which you can then shrink (or even simply cluster, with simple methods) to determine a centroid that you would then associate with the location of a tree. Obviously, processing the size of each "cluster" could probably also tell you how wide is the tree canopy.

In terms of generalising this, you can either fit a 2D polynomial or a 2D NURB or even simply apply a 2D median with medfilt.

Hope this helps.