{"id":5001,"date":"2021-10-04T16:51:25","date_gmt":"2021-10-04T20:51:25","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/?p=5001"},"modified":"2021-10-04T16:51:25","modified_gmt":"2021-10-04T20:51:25","slug":"peaks-peak-filtering-and-gray-scale-weighted-centroids","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2021\/10\/04\/peaks-peak-filtering-and-gray-scale-weighted-centroids\/","title":{"rendered":"Peaks, Peak Filtering, and Gray-Scale Weighted Centroids"},"content":{"rendered":"
Today, I'm continuing my recent theme of thinking about peak-finding in images. When I wrote the first one (<\/span>19-Aug-2021<\/span><\/a>), I didn't realize it was going to turn into a series. This might be the last one\u2014but no promises!<\/span><\/div>
My previous post (<\/span>17-Sep-2021<\/span><\/a>) was based on 1-D examples. Today's post focuses on an image example (in 2-D), and it connects to using <\/span>regionprops<\/span> to compute gray-weighted centroids of peaks. I'll also throw in some 3-D surface visualization, just for fun. Here is today's image:<\/span><\/div>
url = <\/span>\"https:\/\/blogs.mathworks.com\/steve\/files\/snowflakes2.png\"<\/span>;<\/span><\/span><\/div><\/div>
A = imread(url);<\/span><\/span><\/div><\/div>
imshow(A)<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
I'd like to see if we can focus algorithmically on these <\/span>white blobby things<\/span> (that's a technical term that I learned in engineering school). First, let's try computing the <\/span>regional maxima<\/span>.<\/span><\/div>
A_regmax = imregionalmax(A);<\/span><\/span><\/div><\/div>
imshow(A_regmax)<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
At first glance, that does <\/span>not<\/span> look like what I'm after. Let's find out what's going on here.<\/span><\/div>
A_regmax_overlay = imoverlay(A,A_regmax,<\/span>\"green\"<\/span>);<\/span><\/span><\/div><\/div>
imshow(A_regmax_overlay)<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
<\/span><\/div>
I want to look closely at one white blobby thing in particular.<\/span><\/div>
imshow(A_regmax_overlay(30:47,263:285,:))<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
One, or maybe two, of these locations looks like the kind of peak that we are interested in. What about the others? Let me show you two different ways to explore further. The first is with a 3-D visualization, where we treat the image pixel values as heights of a gray surface. The code below shows that surface, and then it displays the detected regional maxima as blue-and-yellow dots just above the surface.<\/span><\/div>
figure<\/span><\/span><\/div><\/div>
A_cropped = A(30:47,263:285);<\/span><\/span><\/div><\/div>
A_regmax_cropped = A_regmax(30:47,263:285);<\/span><\/span><\/div><\/div>
surf(A_cropped, EdgeColor = <\/span>\"none\"<\/span>)<\/span><\/span><\/div><\/div>
colormap(<\/span>\"gray\"<\/span>)<\/span><\/span><\/div><\/div>
lighting <\/span>gouraud<\/span><\/span><\/div><\/div>
daspect([1 1 15])<\/span><\/span><\/div><\/div>
[y,x,~] = find(A_regmax_cropped);<\/span><\/span><\/div><\/div>
z = A_cropped(A_regmax_cropped);<\/span><\/span><\/div><\/div>
hold <\/span>on<\/span><\/span><\/div><\/div>
plot3(x,y,z+1,LineStyle=<\/span>\"none\"<\/span>,Marker=<\/span>\"o\"<\/span>,<\/span>...<\/span><\/span><\/div><\/div>
MarkerEdgeColor=<\/span>\"b\"<\/span>,<\/span>...<\/span><\/span><\/div><\/div>
MarkerFaceColor=<\/span>\"y\"<\/span>)<\/span><\/span><\/div><\/div>
xlabel(<\/span>\"x\"<\/span>)<\/span><\/span><\/div><\/div>
ylabel(<\/span>\"y\"<\/span>)<\/span><\/span><\/div><\/div>
set(gca,YDir = <\/span>\"reverse\"<\/span>)<\/span><\/span><\/div><\/div>
view(-35,45)<\/span><\/span><\/div><\/div>
camlight<\/span><\/span><\/div><\/div>
hold <\/span>off<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
A second way is to go back to 2-D, but using an extremely zoomed-in view, with the actual pixel values superimposed on the individual pixels. To accomplish, I used the <\/span>Image Viewer<\/span><\/a> and the <\/span>Pixel Region Tool<\/span><\/a>. Below is a screen shot of the Pixel Region Tool. I have annotated the screen shot with the locations of the regional maxima.<\/span><\/div>
\"\"<\/img><\/div>
With either technique, you can get the idea that several of these maxima locations are just small, uninteresting fluctuations.<\/span><\/div>
As I discussed in the previous posts, you can use the <\/span>h-maxima transform<\/span> to filter out these small peaks.<\/span><\/div>
h = 75;<\/span><\/span><\/div><\/div>
B = imhmax(A,h);<\/span><\/span><\/div><\/div>
imshow(B(30:47,263:285))<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
B_regmax = imregionalmax(B);<\/span><\/span><\/div><\/div>
imshow(B_regmax(30:47,263:285))<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
B_regmax_overlay = imoverlay(B,B_regmax,<\/span>\"green\"<\/span>);<\/span><\/span><\/div><\/div>
imshow(B_regmax_overlay)<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
axis([225 315 30 60])<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
I want to finish up with a connection to the <\/span>regionprops<\/span><\/a> function. This function is used most often with a binary image input, and it computes various geometries properties of objects (connected components) in the binary image. But you can also provide a gray-scale image at the same time. With that additional info, you can compute other properties. Here, I am interested in the <\/span>weighted centroids<\/span> of the various detected peaks. The original image pixel values will be used as the weights. Here's how to do it:<\/span><\/div>
props = regionprops(<\/span>\"table\"<\/span>,B_regmax,A,<\/span>\"WeightedCentroid\"<\/span>)<\/span><\/span><\/div>
props = <\/span>99\u00d71 table <\/span><\/div>
 <\/span><\/th>WeightedCentroid<\/span><\/th><\/tr><\/thead>
1<\/span><\/th>1.6702<\/span><\/td>24.6139<\/span><\/td><\/tr>
2<\/span><\/th>2.4948<\/span><\/td>43.3596<\/span><\/td><\/tr>
3<\/span><\/th>1.4239<\/span><\/td>97.7269<\/span><\/td><\/tr>
4<\/span><\/th>14.5813<\/span><\/td>52.3127<\/span><\/td><\/tr>
5<\/span><\/th>20.5567<\/span><\/td>100.2771<\/span><\/td><\/tr>
6<\/span><\/th>28.5663<\/span><\/td>89.2373<\/span><\/td><\/tr>
7<\/span><\/th>29.9773<\/span><\/td>18.0773<\/span><\/td><\/tr>
8<\/span><\/th>31.1731<\/span><\/td>50.3474<\/span><\/td><\/tr>
9<\/span><\/th>41.2062<\/span><\/td>2.7718<\/span><\/td><\/tr>
10<\/span><\/th>45.1645<\/span><\/td>39.6400<\/span><\/td><\/tr>
11<\/span><\/th>44.1222<\/span><\/td>25.3352<\/span><\/td><\/tr>
12<\/span><\/th>43.9098<\/span><\/td>67.8786<\/span><\/td><\/tr>
13<\/span><\/th>44.8320<\/span><\/td>103.4289<\/span><\/td><\/tr>
14<\/span><\/th>50.9091<\/span><\/td>16.5648<\/span><\/td><\/tr>
\u22ee<\/th><\/tr><\/tbody><\/table><\/div><\/div><\/div><\/div><\/div><\/div>
And here's one way to visualize the results:<\/span><\/div>
imshow(A)<\/span><\/span><\/div><\/div>
hold <\/span>on<\/span><\/span><\/div><\/div>
plot(props.WeightedCentroid(:,1),props.WeightedCentroid(:,2),<\/span>...<\/span><\/span><\/div><\/div>
LineStyle=<\/span>\"none\"<\/span>, Marker=<\/span>\"o\"<\/span>, MarkerEdgeColor=<\/span>\"y\"<\/span>,<\/span>...<\/span><\/span><\/div><\/div>
MarkerFaceColor=<\/span>\"b\"<\/span>)<\/span><\/span><\/div><\/div>
hold <\/span>off<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
axis([225 315 30 60])<\/span><\/span><\/div><\/div>
title(<\/span>\"Gray-scale weighted centroids\"<\/span>)<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>\r\n<\/div>