{"id":4839,"date":"2021-09-17T15:18:39","date_gmt":"2021-09-17T19:18:39","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/?p=4839"},"modified":"2021-09-29T09:48:01","modified_gmt":"2021-09-29T13:48:01","slug":"small-peak-suppression-with-the-h-maxima-transform","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2021\/09\/17\/small-peak-suppression-with-the-h-maxima-transform\/","title":{"rendered":"Small-Peak Suppression with the H-Maxima Transform"},"content":{"rendered":"
Last time<\/span><\/a>, I introduced the idea of a <\/span>regional maximum<\/span>. Today, I want to add a concept that makes the regional maximum more useful: suppressing very small local maxima, possibly present only because of noise, that are are unimportant, before identifying the regional maxima. This \"small peak suppression\" can be accomplished using something called the <\/span>h-maxima transform.<\/span> The functions of interest today: <\/span>imregionalmax<\/span>, <\/span>imhmax<\/span>, <\/span>imextendedmax<\/span>, <\/span>imregionalmin<\/span>, <\/span>imhmin<\/span>, and <\/span>imextendedmin<\/span>. <\/span><\/div>
Let's get there by starting once again with the definition of <\/span>regional maximum<\/span>: a connected component of pixels with a constant value <\/span>h<\/span>, where every pixel that is neighbor to that connected component has a value that is lower than <\/span>h<\/span>. I'd like to elaborate on that using a one-dimensional example.<\/span><\/div>
y = peaksAndPlateaus;<\/span><\/span><\/div><\/div>
x = 1:length(y);<\/span><\/span><\/div><\/div>
plot(x,y)<\/span><\/span><\/div><\/div>
axis([-5 105 0 120])<\/span><\/span><\/div><\/div>
grid <\/span>on<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
<\/span><\/div>
I'll use dilation and erosion to identify curve samples that are on plateaus.<\/span><\/div>
plateau_mask = (y == imdilate(y,[1 1 1])) & (y == imerode(y,[1 1 1]));<\/span><\/span><\/div><\/div>
plateau_mask = imdilate(plateau_mask,[1 1 1]);<\/span><\/span><\/div><\/div>
hold <\/span>on<\/span><\/span><\/div><\/div>
plot(x(plateau_mask),y(plateau_mask),<\/span>\"*\"<\/span>)<\/span><\/span><\/div><\/div>
hold <\/span>off<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
Here are some zoomed-in views that illustrate three kinds of plateaus that are of interest.<\/span><\/div>
tiledlayout(2,2)<\/span><\/span><\/div><\/div>
nexttile<\/span><\/span><\/div><\/div>
plot(x,y)<\/span><\/span><\/div><\/div>
hold <\/span>on<\/span><\/span><\/div><\/div>
plot(x(plateau_mask),y(plateau_mask),<\/span>\"*\"<\/span>)<\/span><\/span><\/div><\/div>
hold <\/span>off<\/span><\/span><\/div><\/div>
axis([0 40 40 60])<\/span><\/span><\/div><\/div>
grid <\/span>on<\/span><\/span><\/div><\/div>
text(20,52,<\/span>\"regional maximum\"<\/span>,HorizontalAlignment=<\/span>\"center\"<\/span>)<\/span><\/span><\/div><\/div>
<\/div><\/div>
nexttile<\/span><\/span><\/div><\/div>
plot(x,y)<\/span><\/span><\/div><\/div>
hold <\/span>on<\/span><\/span><\/div><\/div>
plot(x(plateau_mask),y(plateau_mask),<\/span>\"*\"<\/span>)<\/span><\/span><\/div><\/div>
hold <\/span>off<\/span><\/span><\/div><\/div>
axis([40 60 10 30])<\/span><\/span><\/div><\/div>
grid <\/span>on<\/span><\/span><\/div><\/div>
text(47,14,[<\/span>\"regional\" \"minimum\"<\/span>],HorizontalAlignment = <\/span>\"center\"<\/span>,<\/span>...<\/span><\/span><\/div><\/div>
VerticalAlignment = <\/span>\"top\"<\/span>)<\/span><\/span><\/div><\/div>
text(55,14,[<\/span>\"regional\" \"minimum\"<\/span>],HorizontalAlignment = <\/span>\"center\"<\/span>,<\/span>...<\/span><\/span><\/div><\/div>
VerticalAlignment = <\/span>\"top\"<\/span>)<\/span><\/span><\/div><\/div>
text(51,19,[<\/span>\"regional\" \"maximum\"<\/span>],HorizontalAlignment = <\/span>\"center\"<\/span>,<\/span>...<\/span><\/span><\/div><\/div>
VerticalAlignment = <\/span>\"bottom\"<\/span>)<\/span><\/span><\/div><\/div>
<\/div><\/div>
nexttile<\/span><\/span><\/div><\/div>
plot(x,y)<\/span><\/span><\/div><\/div>
hold <\/span>on<\/span><\/span><\/div><\/div>
plot(x(plateau_mask),y(plateau_mask),<\/span>\"*\"<\/span>)<\/span><\/span><\/div><\/div>
hold <\/span>off<\/span><\/span><\/div><\/div>
axis([58 69 40 95])<\/span><\/span><\/div><\/div>
grid <\/span>on<\/span><\/span><\/div><\/div>
text(66,83,<\/span>\"plateau\"<\/span>,HorizontalAlignment = <\/span>\"right\"<\/span>)<\/span><\/span><\/div><\/div>
text(64,79,[<\/span>\"not a regional max\" \"not a regional min\"<\/span>],<\/span>...<\/span><\/span><\/div><\/div>
HorizontalAlignment = <\/span>\"left\"<\/span>, VerticalAlignment = <\/span>\"top\"<\/span>)<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
As I showed last time, the function <\/span>imregionalmax<\/span> identifies the location of regional maxima. These can be individual samples or pixels, or they can be plateaus. Take a look:<\/span><\/div>
reg_max_mask = imregionalmax(y);<\/span><\/span><\/div><\/div>
tiledlayout(3,2)<\/span><\/span><\/div><\/div>
nexttile<\/span><\/span><\/div><\/div>
plot(x,y)<\/span><\/span><\/div><\/div>
axis([-5 105 0 120])<\/span><\/span><\/div><\/div>
grid <\/span>on<\/span><\/span><\/div><\/div>
title(<\/span>\"Input curve\"<\/span>)<\/span><\/span><\/div><\/div>
<\/div><\/div>
nexttile<\/span><\/span><\/div><\/div>
plot(x,reg_max_mask)<\/span><\/span><\/div><\/div>
axis([-5 105 -0.05 1.05])<\/span><\/span><\/div><\/div>
grid <\/span>on<\/span><\/span><\/div><\/div>
title(<\/span>\"Output of imregionalmax\"<\/span>)<\/span><\/span><\/div><\/div>
<\/div><\/div>
nexttile([2 2])<\/span><\/span><\/div><\/div>
plot(x,y)<\/span><\/span><\/div><\/div>
hold <\/span>on<\/span><\/span><\/div><\/div>
plot(x(reg_max_mask),y(reg_max_mask),<\/span>\"*\"<\/span>)<\/span><\/span><\/div><\/div>
hold <\/span>off<\/span><\/span><\/div><\/div>
axis([-5 105 0 120])<\/span><\/span><\/div><\/div>
grid <\/span>on<\/span><\/span><\/div><\/div>
title(<\/span>\"Regional maxima plateaus\"<\/span>)<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
The curve I've been experimenting with is relatively smooth, with plateaus that are completely flat. Unfortunately, the <\/span>imregionalmax<\/span> is significantly less useful when the input is even a little noisy. That's because noise introduces small local maxima in many places, and <\/span>imregionalmax<\/span> picks up even the tiniest local maxima. Let's try it. If I add a very small amount of noise, you can't even see it in the plot:<\/span><\/div>
rng(17)<\/span><\/span><\/div><\/div>
yn = y + (0.01 * randn(size(y)));<\/span><\/span><\/div><\/div>
clf<\/span><\/span><\/div><\/div>
plot(x,yn)<\/span><\/span><\/div><\/div>
axis([-5 105 0 120])<\/span><\/span><\/div><\/div>
grid <\/span>on<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
But now the output of <\/span>imregionalmax<\/span> looks completely different.<\/span><\/div>
reg_max_mask_n = imregionalmax(yn);<\/span><\/span><\/div><\/div>
plot(x,yn)<\/span><\/span><\/div><\/div>
hold <\/span>on<\/span><\/span><\/div><\/div>
plot(x(reg_max_mask_n),yn(reg_max_mask_n),<\/span>\"*\"<\/span>)<\/span><\/span><\/div><\/div>
hold <\/span>off<\/span><\/span><\/div><\/div>
axis([-5 105 0 120])<\/span><\/span><\/div><\/div>
grid <\/span>on<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
Now, only a scattering of samples on the broad plateaus are identified as regional maxima, and a sample on the small \"shoulder\" plateau is marked as a regional maximum where it wasn't before. We have to zoom in pretty far to see why this is happening.<\/span><\/div>
axis([5 20 49.8 50.1])<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
axis([62 67 79.9 80.2])<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
You can see that the added noise is introducing tiny local maxima in various places, and that is throwing <\/span>imregionalmax<\/span> off the scent from what we'd like for it to do.<\/span><\/div>

imhmax<\/span><\/h2>
And here's where the function <\/span>imhmax<\/span> comes into play. It computes something called the <\/span>h-maxima<\/span> transform. This operation essentially suppresses any local maxima that are less than a certain height (<\/span>h<\/span>) above their immediate surroundings. Let's say, for example, that we are only interested in local or regional maxima that have a height of more than 5 units above their immediate surroundings. We accomplish this by calling <\/span>imhmax<\/span> like this:<\/span><\/div>
yn_h = imhmax(yn,5);<\/span><\/span><\/div><\/div>
plot(x,yn)<\/span><\/span><\/div><\/div>
hold <\/span>on<\/span><\/span><\/div><\/div>
plot(x,yn_h)<\/span><\/span><\/div><\/div>
hold <\/span>off<\/span><\/span><\/div><\/div>
axis([-5 105 0 120])<\/span><\/span><\/div><\/div>
grid <\/span>on<\/span><\/span><\/div><\/div>
legend([<\/span>\"original\" \"h-maxima transform\"<\/span>])<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
In the output of <\/span>imhmax<\/span>, the broad peaks have been flattened out, and the small peak in the middle, which has a height of only 3 above its surroundings, has been flattened out and eliminated. Now let's try <\/span>imregionalmax<\/span> again.<\/span><\/div>
reg_max_mask_n_h = imregionalmax(yn_h);<\/span><\/span><\/div><\/div>
plot(x,yn)<\/span><\/span><\/div><\/div>
hold <\/span>on<\/span><\/span><\/div><\/div>
plot(x(reg_max_mask_n_h),yn(reg_max_mask_n_h),<\/span>\"*\"<\/span>)<\/span><\/span><\/div><\/div>
hold <\/span>off<\/span><\/span><\/div><\/div>
axis([-5 105 0 120])<\/span><\/span><\/div><\/div>
grid <\/span>on<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
With that computation, we have identified only the peaks that are significantly higher than what is around them, and that's often what we are interested in. Notice that the <\/span>h-maxima<\/span> transform spreads out the peak a bit, and so the result shown above is picking out a bit of the \"shoulders\" to the left and right of each broad plateau. We could do some post-processing to correct that, if necessary.<\/span><\/div>
There is another function, <\/span>imextendedmax<\/span>, that simply combines the h-maxima transform and regional maxima steps into a single function. If you're interested in minima instead, then you can use the upside-down versions of these functions: <\/span>imhmin<\/span>, <\/span>imregionalmin<\/span>, and <\/span>imextendedmin<\/span>. <\/span><\/div>
If you have found interesting uses for these functions in your own work, I'd love to hear about it, so please leave a comment.<\/span><\/div>
<\/span><\/div>

Utility function<\/span><\/h3>
function <\/span>y = peaksAndPlateaus<\/span><\/span><\/div><\/div>
f = @(x) min(max(1-abs(x\/5),0),0.5);<\/span><\/span><\/div><\/div>
<\/div><\/div>
g1 = @(x) 100*f((x-20)\/5);<\/span><\/span><\/div><\/div>
g2 = @(x) 200*f((x-80)\/5);<\/span><\/span><\/div><\/div>
g3 = @(x) 30*f((x-50)\/3);<\/span><\/span><\/div><\/div>
g4 = @(x) 5*f(2*(x-50));<\/span><\/span><\/div><\/div>
g5 = @(x) 30*f((x-60));<\/span><\/span><\/div><\/div>
<\/div><\/div>
g = @(x) g1(x) + g2(x) + g3(x) + g4(x) + g5(x);<\/span><\/span><\/div><\/div>
<\/div><\/div>
y = round(g(0:100));<\/span><\/span><\/div><\/div>
end<\/span><\/span><\/div><\/div><\/div>\r\n<\/div>