Today's blog post was inspired by an example written by my friend and Image Processing Toolbox developer, Alex Taylor. A few years back, Alex tinkered with using toolbox algorithms to achieve a pseudo-artistic "posterization" effect, like this: ... read more >>

]]>```
imshow('eddins-horn-1000.png')
```

I thought the technique was cool and worth showing. I also noticed that the example uses several Image Processing Toolbox functions that haven't been in the product very long, including:

So I realized this would be a good chance to introduce you to some recent toolbox additions that you might not know about.

The two functions *sRGB* and *Adobe RGB (1998)*.

I have mentioned these functions in a couple of previous blog posts, including "Displaying a color gamut surface" and "Out-of-gamut colors."

The function

I = imread('rice.png'); bw = imbinarize(I,'adaptive'); mask = boundarymask(bw); subplot(1,2,1) imshow(I) title('Original image') subplot(1,2,2) imshow(mask) title('Boundary mask')

The function

```
I_overlay = imoverlay(I,mask,'yellow');
imshow(I_overlay)
```

In my first year of writing this blog (2006), I wrote my own function with this name and submitted it to the File Exchange. Now that the Image Processing Toolbox has its own version, I'll probably remove mine from the File Exchange soon.

A "superpixel" is simply a group of connected pixels that have similar colors. Computing superpixels has found a regular place in a variety of image analysis and computer vision tasks. The Image Processing Toolbox function

A = imread('kobi.png'); L = superpixels(A,1500); mask = boundarymask(L); B = imoverlay(A,mask,'cyan'); clf imshow(B)

The superpixel posterization method, as implemented by Alex, starts by using

Here's how it works.

A = imread('eddins-horn.png'); imshow(A) title('Original image')

Next, compute and visualize the superpixel clusters.

```
[L,N] = superpixels(A,1000);
BW = boundarymask(L);
imshow(imoverlay(A,BW,'cyan'))
```

In the next step, I want to replace the pixels in each superpixel cluster with the mean of the cluster's colors. But I want to compute the mean in L*a*b* space, so I start by converting from RGB to L*a*b*.

Alab = rgb2lab(A);

I'll use the function

pixel_idx = label2idx(L);

For each of the

Aplab = Alab; Ln = numel(L); for k = 1:N idx = pixel_idx{k}; Aplab(idx) = mean(Alab(idx)); Aplab(idx+Ln) = mean(Alab(idx+Ln)); Aplab(idx+2*Ln) = mean(Alab(idx+2*Ln)); end

Finally, convert back to RGB space.

Ap = lab2rgb(Aplab); imshow(Ap)]]>

Today I want to finish up my long-running discussion of Feret diameters. (Previous posts: 29-Sep-2017, 24-Oct-2017, 20-Feb-2018, and 16-Mar-2018.) ... read more >>

]]>Recall that the Feret diameter measures the extent of an object along a particular direction. The diagram below illustrates the concept. Place the object to be measured inside the jaws of a caliper, with the caliper oriented at a specified angle. Close the jaws tightly on the object while maintaining that angle. The distance between the jaws is the Feret diameter at that angle.

The *maximum Feret diameter* and *minimum Feret diameter* measure the maximum and minimum width of an object. In previous posts, I talked about how to identify *antipodal vertex pairs* from the set the convex hull vertices to speed up the process of finding these maximum and minimum measures. I also considered the various assumptions one can make about the shape of an individual pixel and how those assumptions can affect the results. (In the rest of this post, I'll assume a diamond pixel shape, as suggested by Cris.)

Let's look at these several measurements for a particular object.

```
bw = imread('shape.png');
imshow(bw)
```

visualizeFeretProperties(bw)

The first diagram above shows the maximum Feret diameter and its orientation. It also shows the Feret diameter at the angle that is orthogonal to the maximum diameter.

The second diagram is similar, except that it shows the minimum Feret diameter instead of the orthogonal diameter. You can see that they are not the same, and, in general, they won't be.

The third diagram shows the minimum-area bounding box, which can be found using a search procedure similar to the minimum Feret diameter. In this diagram, notice that the bounding box is not exactly aligned with the direction of the maximum Feret dimension. In general, they won't necessarily be aligned, and then the length of the minimum-area bounding box will be less than the maximum Feret dimension.

Below is a function I wrote that adds Feret diameter properties to the table returned by

T = regionprops('table',bw,'PixelList')

T=1×1 tablePixelList _________________ [163875×2 double]

T = feretProperties(T)

T=1×12 tablePixelList MaxFeretDiameter MaxFeretDiameterEndpoints MaxFeretDiameterOrientation MinFeretDiameter MinFeretDiameterTrianglePoints MinFeretDiameterOrientation OrthogonalDiameter OrthogonalDiameterLowerPoints OrthogonalDiameterUpperPoints MinAreaBoundingBox MinAreaBoundingBoxArea _________________ ________________ _________________________ ___________________________ ________________ ______________________________ ___________________________ __________________ _____________________________ _____________________________ __________________ ______________________ [163875×2 double] 781.13 [2×2 double] -41.835 317.08 [3×2 double] -143.8 331.18 [1×2 double] [1×2 double] [5×2 double] 2.4321e+05

I hope that at least a few of you enjoyed this diversion into algorithms associated with object shape measurement.

function T = feretProperties(T) % Copyright 2017-2018 The MathWorks, Inc. maxfd = zeros(height(T),1); maxfd_endpoints = cell(height(T),1); maxfd_orientation = zeros(height(T),1); minfd = zeros(height(T),1); minfd_triangle_points = cell(height(T),1); minfd_orientation = zeros(height(T),1); minod = zeros(height(T),1); minod_lower_points = cell(height(T),1); minod_upper_points = cell(height(T),1); minbb = cell(height(T),1); minbb_a = zeros(height(T),1); for k = 1:height(T) pixels = T.PixelList{k}; V = pixelHull(pixels,'diamond'); pairs = antipodalPairs(V); [maxfd(k),maxfd_endpoints{k}] = maxFeretDiameter(V,pairs); points = maxfd_endpoints{k}; e = points(2,:) - points(1,:); maxfd_orientation(k) = atan2d(e(2),e(1)); [minfd(k),minfd_triangle_points{k}] = minFeretDiameter(V,pairs); points = minfd_triangle_points{k}; e = points(2,:) - points(1,:); thetad = atan2d(e(2),e(1)); minfd_orientation(k) = mod(thetad + 180 + 90,360) - 180; [minod(k),minod_lower_points{k},minod_upper_points{k}] = ... feretDiameter(V,maxfd_orientation(k)+90); [minbb{k},minbb_a(k)] = minAreaBoundingBox(V,pairs); end T.MaxFeretDiameter = maxfd; T.MaxFeretDiameterEndpoints = maxfd_endpoints; T.MaxFeretDiameterOrientation = maxfd_orientation; T.MinFeretDiameter = minfd; T.MinFeretDiameterTrianglePoints = minfd_triangle_points; T.MinFeretDiameterOrientation = minfd_orientation; T.OrthogonalDiameter = minod; T.OrthogonalDiameterLowerPoints = minod_lower_points; T.OrthogonalDiameterUpperPoints = minod_upper_points; T.MinAreaBoundingBox = minbb; T.MinAreaBoundingBoxArea = minbb_a; end function [bb,A] = minAreaBoundingBox(V,antipodal_pairs) % Copyright 2017-2018 The MathWorks, Inc. if nargin < 2 antipodal_pairs = antipodalPairs(V); end n = size(antipodal_pairs,1); p = antipodal_pairs(:,1); q = antipodal_pairs(:,2); A = Inf; thetad = []; for k = 1:n if k == n k1 = 1; else k1 = k+1; end pt1 = []; pt2 = []; pt3 = []; if (p(k) ~= p(k1)) && (q(k) == q(k1)) pt1 = V(p(k),:); pt2 = V(p(k1),:); pt3 = V(q(k),:); elseif (p(k) == p(k1)) && (q(k) ~= q(k1)) pt1 = V(q(k),:); pt2 = V(q(k1),:); pt3 = V(p(k),:); end if ~isempty(pt1) % Points pt1, pt2, and pt3 are possibly on the minimum-area % bounding box, with points pt1 and pt2 forming an edge coincident with % the bounding box. Call the height of the triangle with base % pt1-pt2 the height of the bounding box, h. h = triangleHeight(pt1,pt2,pt3); pt1pt2_direction = atan2d(pt2(2)-pt1(2),pt2(1)-pt1(1)); w = feretDiameter(V,pt1pt2_direction); A_k = h*w; if (A_k < A) A = A_k; thetad = pt1pt2_direction; end end end % Rotate all the points so that pt1-pt2 for the minimum bounding box points % straight up. r = 90 - thetad; cr = cosd(r); sr = sind(r); R = [cr -sr; sr cr]; Vr = V * R'; xr = Vr(:,1); yr = Vr(:,2); xmin = min(xr); xmax = max(xr); ymin = min(yr); ymax = max(yr); bb = [ ... xmin ymin xmax ymin xmax ymax xmin ymax xmin ymin ]; % Rotate the bounding box points back. bb = bb * R; end function h = triangleHeight(P1,P2,P3) % Copyright 2017-2018 The MathWorks, Inc. h = 2 * abs(signedTriangleArea(P1,P2,P3)) / norm(P1 - P2); end function area = signedTriangleArea(A,B,C) % Copyright 2017-2018 The MathWorks, Inc. area = ( (B(1) - A(1)) * (C(2) - A(2)) - ... (B(2) - A(2)) * (C(1) - A(1)) ) / 2; end function [d,V1,V2] = feretDiameter(V,theta) % Copyright 2017-2018 The MathWorks, Inc. % Rotate points so that the direction of interest is vertical. alpha = 90 - theta; ca = cosd(alpha); sa = sind(alpha); R = [ca -sa; sa ca]; % Vr = (R * V')'; Vr = V * R'; y = Vr(:,2); ymin = min(y,[],1); ymax = max(y,[],1); d = ymax - ymin; if nargout > 1 V1 = V(y == ymin,:); V2 = V(y == ymax,:); end function [d,end_points] = maxFeretDiameter(P,antipodal_pairs) % Copyright 2017-2018 The MathWorks, Inc. if nargin < 2 antipodal_pairs = antipodalPairs(P); end v = P(antipodal_pairs(:,1),:) - P(antipodal_pairs(:,2),:); D = hypot(v(:,1),v(:,2)); [d,idx] = max(D,[],1); P1 = P(antipodal_pairs(idx,1),:); P2 = P(antipodal_pairs(idx,2),:); end_points = [P1 ; P2]; end function [d,triangle_points] = minFeretDiameter(V,antipodal_pairs) % Copyright 2017-2018 The MathWorks, Inc. if nargin < 2 antipodal_pairs = antipodalPairs(V); end n = size(antipodal_pairs,1); p = antipodal_pairs(:,1); q = antipodal_pairs(:,2); d = Inf; triangle_points = []; for k = 1:n if k == n k1 = 1; else k1 = k+1; end pt1 = []; pt2 = []; pt3 = []; if (p(k) ~= p(k1)) && (q(k) == q(k1)) pt1 = V(p(k),:); pt2 = V(p(k1),:); pt3 = V(q(k),:); elseif (p(k) == p(k1)) && (q(k) ~= q(k1)) pt1 = V(q(k),:); pt2 = V(q(k1),:); pt3 = V(p(k),:); end if ~isempty(pt1) % Points pt1, pt2, and pt3 form a possible minimum Feret diameter. % Points pt1 and pt2 form an edge parallel to caliper direction. % The Feret diameter orthogonal to the pt1-pt2 edge is the height % of the triangle with base pt1-pt2. d_k = triangleHeight(pt1,pt2,pt3); if d_k < d d = d_k; triangle_points = [pt1; pt2; pt3]; end end end end]]>

What is the shape of a pixel? At various times, I have a pixel as a square (often), a point (sometimes), or a rectangle (occasionally). I recall back in grad school doing some homework where we were treating pixels as hexagons.... read more >>

]]>What is the shape of a pixel? At various times, I have a pixel as a square (often), a point (sometimes), or a rectangle (occasionally). I recall back in grad school doing some homework where we were treating pixels as hexagons.

As I haved worked through the last few posts on computing Feret diameters, though, I have started to entertain the possible usefulness of considering pixels to be circles. (See 29-Sep-2017, 24-Oct-2017, and 20-Feb-2018.) Let me try to explain why.

Here's a binary image with a single foreground blob (or "object," or "connected component.")

```
bw = imread('Martha''s Vineyard (30x20).png');
imshow(bw)
```

Most of the time, we think of image pixels as being squares with unit area.

pixelgrid

We can use `find` to get the $x$- and $y$-coordinates of the pixel centers, and then we can use `convhull` to find their convex hull. As an optimization that I think will often reduce execution time and memory, I'm going to preprocess the input binary image here by calling `bwperim`. I'm not going to show that step everywhere in this example, though.

[y,x] = find(bwperim(bw)); hold on plot(x,y,'.') hold off title('Pixel centers') h = convhull(x,y); x_hull = x(h); y_hull = y(h); hold on hull_line = plot(x_hull,y_hull,'r*','MarkerSize',12); hold off title('Pixel centers and convex hull vertices')

Notice that there are some chains of three or more colinear convex hull vertices.

```
xlim([21.5 32.5])
ylim([9.5 15.5])
title('Colinear convex hull vertices')
```

In some of the other processing steps related to Feret diameter measurements, colinear convex hull vertices can cause problems. We can eliminate these vertices directly in the call to `convhull` using the `'Simplify'` parameter.

h = convhull(x,y,'Simplify',true); x_hull = x(h); y_hull = y(h); delete(hull_line); hold on plot(x_hull,y_hull,'r*','MarkerSize',12) hold off title('Colinear hull vertices removed')

imshow(bw) hold on plot(x_hull,y_hull,'r-*','LineWidth',2,'MarkerSize',12) hold off title('A Blob''s Convex Hull and Its Vertices')

Notice, though, that there are white bits showing outside the red convex hull polygon. That's because we are only using the **pixel centers**.

Consider a simpler binary object, one that has only one row.

bw2 = false(5,15); bw2(3,5:10) = true; imshow(bw2) pixelgrid [y,x] = find(bw2);

The function `convhull` doesn't even work on colinear points.

try hull = convhull(x,y,'Simplify',true); catch e fprintf('Error message from convhull: "%s"\n', e.message); end

Error message from convhull: "Error computing the convex hull. The points may be collinear."

But even if it did return an answer, the answer would be a degenerate polygon with length 5 (even though the number of foreground pixels is 6) and zero area.

hold on plot(x,y,'r-*','MarkerSize',12,'LineWidth',2) hold off title('Degenerate convex hull polygon')

We can solve this degeneracy problem by using **square pixels**.

In the computation of the convex hull above, we treated each pixel as a **point**. We can, instead, treat each pixel as a **square** by computing the convex hull of all the **corners** of every pixel. Here's one way to perform that computation.

offsets = [ ... 0.5 -0.5 0.5 0.5 -0.5 -0.5 -0.5 0.5 ]'; offsets = reshape(offsets,1,2,[]); P = [x y]; Q = P + offsets; R = permute(Q,[1 3 2]); S = reshape(R,[],2); h = convhull(S,'Simplify',true); x_hull = S(h,1); y_hull = S(h,2);

imshow(bw2) pixelgrid hold on plot(x_hull,y_hull,'r-*','MarkerSize',12,'LineWidth',2) hold off title('Convex hull of square pixels')

This result looks good at first glance. However, it loses some of its appeal when you consider the implications for computing the maximum Feret diameter.

points = [x_hull y_hull]; [d,end_points] = maxFeretDiameter(points,antipodalPairs(points)) hold on plot(end_points(:,1),end_points(:,2),'k','LineWidth',3) hold off title('The maximum Feret diameter is not horizontal')

d = 6.0828 end_points = 10.5000 2.5000 4.5000 3.5000

The maximum Feret distance of this horizontal segment of is 6.0828 ($\sqrt{37}$) instead of 6, and the corresponding orientation in degrees is:

atan2d(1,6)

ans = 9.4623

instead of 0.

Another worthy attempt is to use **diamond** pixels.

Instead of using the four corners of each pixel, let's try using the middle of each pixel edge. Once we define the offsets, the code is exactly the same as for square pixels.

offsets = [ ... 0.5 0.0 0.0 0.5 -0.5 0.0 0.0 -0.5 ]'; offsets = reshape(offsets,1,2,[]); P = [x y]; Q = P + offsets; R = permute(Q,[1 3 2]); S = reshape(R,[],2); h = convhull(S,'Simplify',true); x_hull = S(h,1); y_hull = S(h,2);

imshow(bw2) pixelgrid hold on plot(x_hull,y_hull,'r-*','MarkerSize',12,'LineWidth',2) hold off title('Convex hull of diamond pixels')

Now the max Feret diameter result looks better for the horizontal row of pixels.

points = [x_hull y_hull]; [d,end_points] = maxFeretDiameter(points,antipodalPairs(points)) hold on plot(end_points(:,1),end_points(:,2),'k','LineWidth',3) hold off

d = 6 end_points = 10.5000 3.0000 4.5000 3.0000

Hold on, though. Consider a square blob.

```
bw3 = false(9,9);
bw3(3:7,3:7) = true;
imshow(bw3)
pixelgrid
[y,x] = find(bw3);
P = [x y];
Q = P + offsets;
R = permute(Q,[1 3 2]);
S = reshape(R,[],2);
h = convhull(S,'Simplify',true);
x_hull = S(h,1);
y_hull = S(h,2);
```

hold on plot(x_hull,y_hull,'r-*','MarkerSize',12,'LineWidth',2) points = [x_hull y_hull]; [d,end_points] = maxFeretDiameter(points,antipodalPairs(points)) plot(end_points(:,1),end_points(:,2),'k','LineWidth',3) hold off title('The max Feret diameter is not at 45 degrees')

d = 6.4031 end_points = 7.5000 3.0000 2.5000 7.0000

We'd like to see the max Feret diameter oriented at 45 degrees, and clearly we don't.

OK, I'm going to make one more attempt. I'm going to treat each pixel as **approximately** a **circle**. I'm going to approximate a circle using 24 points that are spaced at 15-degree intervals along the circumference.

thetad = 0:15:345; offsets = 0.5*[cosd(thetad) ; sind(thetad)]; offsets = reshape(offsets,1,2,[]); Q = P + offsets; R = permute(Q,[1 3 2]); S = reshape(R,[],2); h = convhull(S,'Simplify',true); x_hull = S(h,1); y_hull = S(h,2); imshow(bw3) pixelgrid hold on plot(x_hull,y_hull,'r-*','MarkerSize',12,'LineWidth',2) points = [x_hull y_hull]; [d,end_points] = maxFeretDiameter(points,antipodalPairs(points)) plot(end_points(:,1),end_points(:,2),'k','LineWidth',3) axis on hold off

d = 6.6569 end_points = 7.3536 7.3536 2.6464 2.6464

Now the max Feret diameter orientation is what we would naturally expect, which is $\pm 45^{\circ}$. The orientation would also be as expected for a horizontal or vertical segment of pixels.

Still, a circular approximation might not always give exactly what a user might expect. Let's go back to the Martha's Vinyard blob that I started with. I wrote a function called `pixelHull` that can compute the convex hull of binary image pixels in a variety of different ways. The call `pixelHull(bw,24)` computes the pixel hull using a 24-point circle approximation.

Here's the maximum Feret diameter using that approximation.

imshow(bw) V = pixelHull(bw,24); hold on plot(V(:,1),V(:,2),'r-','LineWidth',2,'MarkerSize',12) [d,end_points] = maxFeretDiameter(V,antipodalPairs(V)); plot(end_points(:,1),end_points(:,2),'m','LineWidth',3) axis on pixelgrid hold off

I think many people might expect the maximum Feret diameter to go corner-to-corner in this case, but it doesn't exactly do that.

xlim([22.07 31.92]) ylim([8.63 15.20])

You have to use square pixels to get corner-to-corner.

imshow(bw) V = pixelHull(bw,'square'); hold on plot(V(:,1),V(:,2),'r-','LineWidth',2,'MarkerSize',12) [d,end_points] = maxFeretDiameter(V,antipodalPairs(V)); plot(end_points(:,1),end_points(:,2),'m','LineWidth',3) axis on pixelgrid hold off

xlim([22.07 31.92]) ylim([8.63 15.20])

After all this, I'm still not completely certain which shape assumption will generally work best. My only firm conclusion is that the point approximation is the worst choice. The degeneracies associated with point pixels are just too troublesome.

If you have an opinion, please share it in the comments. (Note: A comment that says, "Steve, you're totally overthinking this" would be totally legit.)

*The rest of the post contains functions used by the code above.*

function V = pixelHull(P,type) if nargin < 2 type = 24; end if islogical(P) P = bwperim(P); [i,j] = find(P); P = [j i]; end if strcmp(type,'square') offsets = [ ... 0.5 -0.5 0.5 0.5 -0.5 0.5 -0.5 -0.5 ]; elseif strcmp(type,'diamond') offsets = [ ... 0.5 0 0 0.5 -0.5 0 0 -0.5 ]; else % type is number of angles for sampling a circle of diameter 1. thetad = linspace(0,360,type+1)'; thetad(end) = []; offsets = 0.5*[cosd(thetad) sind(thetad)]; end offsets = offsets'; offsets = reshape(offsets,1,2,[]); Q = P + offsets; R = permute(Q,[1 3 2]); S = reshape(R,[],2); k = convhull(S,'Simplify',true); V = S(k,:); end

Get
the MATLAB code

Published with MATLAB® R2017b

Last time (if you can remember that long ago), I talked about how to find the maximum *Feret diameter* of a shape. The Feret diameter, sometimes called the *caliper diameter*, is illustrated by the diagram below. In a virtual sense, place the object to be measured inside the jaws of a caliper, with the caliper oriented at a specified angle. Close the jaws tightly on the object while maintaining that angle. The distance between the jaws is the Feret diameter at that angle.... read more >>

Last time (if you can remember that long ago), I talked about how to find the maximum *Feret diameter* of a shape. The Feret diameter, sometimes called the *caliper diameter*, is illustrated by the diagram below. In a virtual sense, place the object to be measured inside the jaws of a caliper, with the caliper oriented at a specified angle. Close the jaws tightly on the object while maintaining that angle. The distance between the jaws is the Feret diameter at that angle.

In the last post, I demonstrated how to find all the *antipodal* vertex pairs of a shape, which is a useful step in finding the maximum Feret diameters of a shape.

Here is a convex shape with all of the antipodal vertex pairs shown.

It turns out that the **minimum** Feret diameter can also be found by looking at these antipodal pairs. Furthermore, it happens that the minimum-distance Feret calipers touch the shape at **three** of these vertices. I'm not going to try to prove that here, but let me illustrate it with a couple of diagrams.

Here is a shape with a couple of "caliper lines" drawn at -30 degrees. (Note that the y-axis is reversed; that's why the angle is negative.)

hull = [ 2.5000 5.5000 3.5000 4.5000 6.5000 2.5000 9.5000 1.5000 10.5000 1.5000 10.5000 3.5000 9.5000 5.5000 5.5000 7.5000 2.5000 7.5000 2.5000 5.5000 ]; plot(hull(:,1),hull(:,2),'r','LineWidth',2) axis equal axis ij axis([0 15 0 10]) hold on plot(hull(:,1),hull(:,2),'r*') [x1,y1] = fullLine(gca,[9.5 5.5],-30); [x2,y2] = fullLine(gca,[6.5 2.5],-30); caliper_lines(1) = plot(x1,y1,'k'); caliper_lines(2) = plot(x2,y2,'k'); hold off

You can always rotate these two caliper lines, by the same angle, until at least one of them touches the **next** antipodal vertex. When you do that rotation, the distance between the lines shrinks.

set(caliper_lines,'Color',[0.8 0.8 0.8]); [x3,y3] = fullLine(gca,[9.5 5.5],thetad); [x4,y4] = fullLine(gca,[6.5 2.5],thetad); hold on plot(x3,y3,'k') plot(x4,y4,'k') hold off

The Feret diameter at that rotated angle, then, is the height of the triangle formed by the three vertices, with the base of the triangle defined by the two vertices touched by the same caliper line.

delete(caliper_lines); hold on plot([5.5 9.5 6.5 5.5],[7.5 5.5 2.5 7.5],'Color','b','LineWidth',2) hold off

The function `minFeretDiameter`, which appears at the end of this post, simply walks around the set of adjacent-vertex triangles formed from antipodal vertex pairs, looking for the triangle with the minimum height.

Let's look for the minimum Feret diameter in a slightly more interesting shape.

load shape plot(bx,by,'LineWidth',1) axis equal axis ij axis([0 650 0 650])

Find the convex hull, making sure to simplify it, and then find the antipodal pairs. (I showed the function `antipodalPairs` in my previous post.)

```
hull = convhull(bx,by,'Simplify',true);
S = [bx(hull) by(hull)];
pairs = antipodalPairs(S);
```

Now we can call `minFeretDiameter`.

[d,tri_points] = minFeretDiameter(S,pairs)

d = 318.7890 tri_points = 613 211 459 443 239 198

Superimpose the minimum-height triangle that `minFeretDiameter` found.

tri_points = [tri_points; tri_points(1,:)]; hold on plot(tri_points(:,1),tri_points(:,2),'LineWidth',2) hold off

Now calculate the caliper angle so that we can visualize the parallel lines of support corresponding to the minimum diameter.

dx = tri_points(2,1) - tri_points(1,1); dy = tri_points(2,2) - tri_points(1,2); angle = atan2d(dy,dx); [x5,y5] = fullLine(gca,tri_points(1,:),angle); [x6,y6] = fullLine(gca,tri_points(3,:),angle);

hold on plot(x5,y5,'k') plot(x6,y6,'k') hold off

For next time, I'm tentatively planning to put together some of these concepts and use them to compute the minimum bounding box for an object.

function [x,y] = fullLine(ax,point,angle_degrees) % Steve Eddins limits = axis(ax); width = abs(limits(2) - limits(1)); height = abs(limits(4) - limits(3)); d = 2*hypot(width,height); x1 = point(1) - d*cosd(angle_degrees); x2 = point(1) + d*cosd(angle_degrees); y1 = point(2) - d*sind(angle_degrees); y2 = point(2) + d*sind(angle_degrees); x = [x1 x2]; y = [y1 y2]; end function [d,triangle_points] = minFeretDiameter(V,antipodal_pairs) % Steve Eddins if nargin < 2 antipodal_pairs = antipodalPairs(V); end n = size(antipodal_pairs,1); p = antipodal_pairs(:,1); q = antipodal_pairs(:,2); d = Inf; triangle_points = []; for k = 1:n if k == n k1 = 1; else k1 = k+1; end pt1 = []; pt2 = []; pt3 = []; if (p(k) ~= p(k1)) && (q(k) == q(k1)) pt1 = V(p(k),:); pt2 = V(p(k1),:); pt3 = V(q(k),:); elseif (p(k) == p(k1)) && (q(k) ~= q(k1)) pt1 = V(q(k),:); pt2 = V(q(k1),:); pt3 = V(p(k),:); end if ~isempty(pt1) % Points pt1, pt2, and pt3 form a possible minimum Feret diameter. % Points pt1 and pt2 form an edge parallel to caliper direction. % The Feret diameter orthogonal to the pt1-pt2 edge is the height % of the triangle with base pt1-pt2. d_k = triangleHeight(pt1,pt2,pt3); if d_k < d d = d_k; triangle_points = [pt1; pt2; pt3]; end end end end

Copyright 2017 The MathWorks, Inc.

Get
the MATLAB code

Published with MATLAB® R2017b

Last time, I wrote about finding the maximum Feret diameter for an object in a binary image, ending up with this figure:... read more >>

]]>Last time, I wrote about finding the maximum Feret diameter for an object in a binary image, ending up with this figure:

I had computed the convex hull of all the pixel corners, and then I computed the pairwise distance between every pair of convex hull vertices to find the maximum distance.

The procedure would work fine in many cases, but the time required to find the maximum distance this way grows with the square of the number of convex hull vertices. With modern digital image resolutions, it's not hard to imagine having thousands of vertices and therefore millions of pairwise distances to compute.

There is a procedure for reducing the number of vertex pairs we can consider. It is based on this theorem:

*The diameter of a convex figure is the greatest distance between parallel lines of support.* [Theorem 4,18, Preparata and Shamos, *Computational Geometry*, 1985]

A *line of support* for a polygon is a line that contains a vertex of the polygon, with the polygon lying entirely on one side of the line.

Let me show you a picture using the convex hull points from last time.

hull = [ 2.5000 5.5000 3.5000 4.5000 6.5000 2.5000 9.5000 1.5000 10.5000 1.5000 10.5000 3.5000 9.5000 5.5000 5.5000 7.5000 2.5000 7.5000 2.5000 5.5000 ]; plot(hull(:,1),hull(:,2),'r','LineWidth',2) hold on plot(hull(:,1),hull(:,2),'r*') hold off axis equal axis ij axis([0 15 0 10]) for theta = -55:5:-35 drawFullLine(gca,[9.5 5.5],theta,'Color',[.6 .6 .6]); end

The plot above shows five different lines of support drawn through the (9.5,5.5) vertex. Now I'll add some of the lines of support through the (2.5,7.5) vertex.

for theta = 10:10:80 drawFullLine(gca,[2.5 7.5],theta,'Color',[.6 .6 .6]); end

You can see that there is no line of support for the (2.5,7.5) vertex that is parallel to a line of support for the (9.5,5.5) vertex. That rules out this pair of vertices for computing the maximum Feret diameter.

Now I'll draw lines of support for a different pair of vertices.

plot(hull(:,1),hull(:,2),'r','LineWidth',2) hold on plot(hull(:,1),hull(:,2),'r*') hold off axis equal axis ij axis([0 15 0 10]) drawFullLine(gca,[6.5 2.5],-30,'Color',[.6 .6 .6]); drawFullLine(gca,[9.5 5.5],-30,'Color',[.6 .6 .6]);

Because the vertices (6.5,2.5) and (9.5,5.5) have parallel lines of support, they are called *antipodal vertices*. There is an algorithm in the Preparata and Shamos book (referenced above) that finds all the antipodal vertices for a convex polygon. There's an implementation of the algorithm in a function at the bottom of this post. I'll use it to find all the antipodal pairs of the convex hull vertices.

pq = antipodalPairs(hull);

**Important note**: the algorithm in `antipodalPairs` assumes that its input is convex. Further, it assumes that the input does not contain any vertices that are on the straight line between the two adjacent vertices. To satisfy this condition, compute the convex hull using the call: `k = convhull(P,'Simplify',true)`.

Now let's plot the line segments joining each antipodal pair.

plot(hull(:,1),hull(:,2),'r','LineWidth',2) hold on plot(hull(:,1),hull(:,2),'r*') axis equal axis ij axis([0 15 0 10]) for k = 1:size(pq,1) x = [hull(pq(k,1),1) hull(pq(k,2),1)]; y = [hull(pq(k,1),2) hull(pq(k,2),2)]; plot(x,y,'Color',[.6 .6 .6]); end hold off

To find the maximum Feret diameter, we only have to check the distances of these 10 segments, instead of 10*9 = 90 line segments as before.

p1 = hull(pq(:,1),:); p2 = hull(pq(:,2),:); v = p1 - p2; d = hypot(v(:,1),v(:,2)); [d_max,idx] = max(d); p1_max = p1(idx,:); p2_max = p2(idx,:); hold on plot([p1_max(1) p2_max(1)],[p1_max(2) p2_max(2)],'b','LineWidth',2) hold off

d_max

d_max = 10

And there is the maximum Feret distance.

I am sure that, in many cases, it would be quicker to find the maximum distance by brute force comparison of all pairs off convex hull vertices. The computation of the antipodal pairs does take some time, after all. I have not done a performance study to investigate this question further. However, the antipodal pairs computation is useful for another measurement: the minimum Feret distance.

I will look at that next time.

*End of post*

*Algorithm, visualization, and utility functions used by the code above*

function h = drawFullLine(ax,point,angle_degrees,varargin) %drawFullLine Draw a line that spans the entire plot % % drawFullLine(ax,point,angle_degrees) draws a line in the % specified axes that goes through the specified point at the % specified angle (in degrees). The line is drawn to span the % entire plot. % % drawFullLine(___,Name,Value) passes name-value parameter pairs % to the line function. % Steve Eddins limits = axis(ax); width = abs(limits(2) - limits(1)); height = abs(limits(4) - limits(3)); d = 2*hypot(width,height); x1 = point(1) - d*cosd(angle_degrees); x2 = point(1) + d*cosd(angle_degrees); y1 = point(2) - d*sind(angle_degrees); y2 = point(2) + d*sind(angle_degrees); h = line(ax,'XData',[x1 x2],'YData',[y1 y2],varargin{:}); end function pq = antipodalPairs(S) % antipodalPairs Antipodal vertex pairs of simple, convex polygon. % % pq = antipodalPairs(S) computes the antipodal vertex pairs of a simple, % convex polygon. S is a Px2 matrix of (x,y) vertex coordinates for the % polygon. S must be simple and convex without repeated vertices. It is % not checked for satisfying these conditions. S can either be closed or % not. The output, pq, is an Mx2 matrix representing pairs of vertices in % S. The coordinates of the k-th antipodal pair are S(pq(k,1),:) and % S(pq(k,2),:). % % TERMINOLOGY % % For a convex polygon, an antipodal pair of vertices is one where you % can draw distinct lines of support through each vertex such that the % lines of support are parallel. % % A line of support is a line that goes through a polygon vertex such % that the interior of the polygon lies entirely on one side of the line. % % EXAMPLE % % Compute antipodal vertices of a polygon and plot the corresponding % line segments. % % x = [0 0 1 3 5 4 0]; % y = [0 1 4 5 4 1 0]; % S = [x' y']; % pq = antipodalPairs(S); % % plot(S(:,1),S(:,2)) % hold on % for k = 1:size(pq,1) % xk = [S(pq(k,1),1) S(pq(k,2),1)]; % yk = [S(pq(k,1),2) S(pq(k,2),2)]; % plot(xk,yk,'LineStyle','--','Marker','o','Color',[0.7 0.7 0.7]) % end % hold off % axis equal % % ALGORITHM NOTES % % This function uses the "ANTIPODAL PAIRS" algorithm, Preparata and % Shamos, Computational Geometry: An Introduction, Springer-Verlag, 1985, % p. 174. % Steve Eddins n = size(S,1); if isequal(S(1,:),S(n,:)) % The input polygon is closed. Remove the duplicate vertex from the % end. S(n,:) = []; n = n - 1; end % The algorithm assumes the input vertices are in counterclockwise order. % If the vertices are in clockwise order, reverse the vertices. clockwise = simplePolygonOrientation(S) < 0; if clockwise S = flipud(S); end % The following variables, including the two anonymous functions, are set % up to follow the notation in the pseudocode on page 174 of Preparata and % Shamos. p and q are indices (1-based) that identify vertices of S. p0 and % q0 identify starting vertices for the algorithm. area(i,j,k) is the area % of the triangle with the corresponding vertices from S: S(i,:), S(j,:), % and S(k,:). next(p) returns the index of the next vertex of S. % % The initialization of p0 is missing from the Preparata and Shamos text. area = @(i,j,k) signedTriangleArea(S(i,:),S(j,:),S(k,:)); next = @(i) mod(i,n) + 1; % mod((i-1) + 1,n) + 1 p = n; p0 = next(p); q = next(p); % The list of antipodal vertices will be built up in the vectors pp and qq. pp = zeros(0,1); qq = zeros(0,1); % ANTIPODAL PAIRS step 3. while (area(p,next(p),next(q)) > area(p,next(p),q)) q = next(q); end q0 = q; % Step 4. while (q ~= p0) % Step 5. p = next(p); % Step 6. % Step 7. (p,q) is an antipodal pair. pp = [pp ; p]; qq = [qq ; q]; % Step 8. while (area(p,next(p),next(q)) > area(p,next(p),q)) q = next(q); % Step 9. if ~isequal([p q],[q0,p0]) % Step 10. pp = [pp ; p]; qq = [qq ; q]; else % This loop break is omitted from the Preparata and Shamos % text. break end end % Step 11. Check for parallel edges. if (area(p,next(p),next(q)) == area(p,next(p),q)) if ~isequal([p q],[q0 n]) % Step 12. (p,next(q)) is an antipodal pair. pp = [pp ; p]; qq = [qq ; next(q)]; else % This loop break is omitted from the Preparata and Shamos % text. break end end end if clockwise % Compensate for the flipping of the polygon vertices. pp = n + 1 - pp; qq = n + 1 - qq; end pq = [pp qq]; end function s = vertexOrientation(P0,P1,P2) % vertexOrientation Orientation of a vertex with respect to line segment. % % s = vertexOrientation(P0,P1,P2) returns a positive number if P2 is to % the left of the line through P0 to P1. It returns 0 if P2 is on the % line. It returns a negative number if P2 is to the right of the line. % % Stating it another way, a positive output corresponds to a % counterclockwise traversal from P0 to P1 to P2. % % P0, P1, and P2 are two-element vectors containing (x,y) coordinates. % % Reference: http://geomalgorithms.com/a01-_area.html, function isLeft() % Steve Eddins s = (P1(1) - P0(1)) * (P2(2) - P0(2)) - ... (P2(1) - P0(1)) * (P1(2) - P0(2)); end function s = simplePolygonOrientation(V) % simplePolygonOrientation Determine vertex order for simple polygon. % % s = simplePolygonOrientation(V) returns a positive number if the simple % polygon V is counterclockwise. It returns a negative number of the % polygon is clockwise. It returns 0 for degenerate cases. V is a Px2 % matrix of (x,y) vertex coordinates. % % Reference: http://geomalgorithms.com/a01-_area.html, function % orientation2D_Polygon() % Steve Eddins n = size(V,1); if n < 3 s = 0; return end % Find rightmost lowest vertext of the polygon. x = V(:,1); y = V(:,2); ymin = min(y,[],1); y_idx = find(y == ymin); if isscalar(y_idx) idx = y_idx; else [~,x_idx] = max(x(y_idx),[],1); idx = y_idx(x_idx(1)); end % The polygon is counterclockwise if the edge leaving V(idx,:) is left of % the entering edge. if idx == 1 s = vertexOrientation(V(n,:), V(1,:), V(2,:)); elseif idx == n s = vertexOrientation(V(n-1,:), V(n,:), V(1,:)); else s = vertexOrientation(V(idx-1,:), V(idx,:), V(idx+1,:)); end end

Copyright 2017 The MathWorks, Inc.

Get
the MATLAB code

Published with MATLAB® R2017b

This is the first of a few blog posts about object measurements based on a concept called the Feret diameter, sometimes called the caliper diameter. The diagram below illustrates the concept. Place the object to be measured inside the jaws of a caliper, with the caliper oriented at a specified... read more >>

]]>This is the first of a few blog posts about object measurements based on a concept called the *Feret diameter*, sometimes called the *caliper diameter*. The diagram below illustrates the concept. Place the object to be measured inside the jaws of a caliper, with the caliper oriented at a specified angle. Close the jaws tightly on the object while maintaining that angle. The distance between the jaws is the Feret diameter at that angle.

Eventually, we'll examine several related measurements (and algorithms for computing them), like these.

Today, though, I'm just going to play around with a very small binary image. (Note: the `pixelgrid` function is a utility that I wrote. It's at the bottom of this post.)

```
bw = [ ...
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 1 1 1 1 0
0 0 0 0 0 1 1 1 1 0 0
0 0 0 1 1 1 1 1 1 0 0
0 0 1 1 1 0 0 0 0 0 0
0 0 1 1 1 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 ];
imshow(bw)
pixelgrid
```

I want to brute-force my to finding the two corners of the white blob that are farthest apart. Let's start by finding all the pixel corners. (I'll be using implicit expansion here, which works in R2016b or later.)

[y,x] = find(bw); hold on plot(x,y,'xb','MarkerSize',5) corners = [x y] + cat(3,[.5 .5],[.5 -.5],[-.5 .5],[-.5 -.5]); corners = permute(corners,[1 3 2]); corners = reshape(corners,[],2); corners = unique(corners,'rows'); plot(corners(:,1),corners(:,2),'sr','MarkerSize',5) hold off

Well, let's get away from pure brute force just a little. We really only need to look at the points that are part of the *convex hull* of all the corner points.

hold on k = convhull(corners); hull_corners = corners(k,:); plot(hull_corners(:,1),hull_corners(:,2),'r','LineWidth',3) plot(hull_corners(:,1),hull_corners(:,2),'ro','MarkerSize',10,'MarkerFaceColor','r') hold off

Now let's find which two points on the convex hull are farthest apart.

dx = hull_corners(:,1) - hull_corners(:,1)'; dy = hull_corners(:,2) - hull_corners(:,2)'; pairwise_dist = hypot(dx,dy); [max_dist,j] = max(pairwise_dist(:)); [k1,k2] = ind2sub(size(pairwise_dist),j); point1 = hull_corners(k1,:); point2 = hull_corners(k2,:); hold on plot([point1(1) point2(1)],[point1(2) point2(2)],'-db','LineWidth',3,'MarkerFaceColor','b') hold off

And the maximum distance is ...

max_dist

max_dist = 10

(I thought I had made a mistake when I first saw this number come out to be an exact integer. But no, the answer is really 10.)

That's a pretty good start. I don't know exactly where I'll go with this topic next time. I'm making it up as I go.

function hout = pixelgrid(h) %pixelgrid Superimpose a grid of pixel edges on an image % pixelgrid superimposes a grid of pixel edges on the image in the % current axes. % % pixelgrid(h_ax) superimposes the grid on the image in the specified % axes. % % pixelgrid(h_im) superimposes the grid on the specified image. % % group = pixelgrid(___) returns an hggroup object that contains the two % lines that are used to draw the grid. One line is thick and darker, and % the other is thin and lighter. Using two contrasting lines in this way % guarantees that the grid will be visible regardless of pixel colors. % % EXAMPLES % % Superimpose pixel grid on color image. Zoom in. % % rgb = imread('peppers.png'); % imshow(rgb) % pixelgrid % axis([440 455 240 250]) % % Change the colors and line widths of the pixel grid. % % rgb = imread('peppers.png'); % imshow(rgb) % h = pixelgrid; % axis([440 455 240 250]) % h.Children(1).Color = [178 223 138]/255; % h.Children(1).LineWidth = 2; % h.Children(2).Color = [31 120 180]/255; % h.Children(2).LineWidth = 4; % % LIMITATIONS % % This function is intended for use when looking at a zoomed-in image % region with a relatively small number of rows and columns. If you use % this function on a typical, full-size image without zooming in, the % image will not be visible under the grid lines. % % This function attempts to determine if MATLAB is running with a high % DPI display. If so, then it displays the pixel grid using thinner % lines. However, the method for detecting a high DPI display on works on % a Mac when Java AWT is available. % Steve Eddins % Copyright The MathWorks, Inc. 2017 if nargin < 1 him = findobj(gca,'type','image'); elseif strcmp(h.Type,'axes') him = findobj(h,'type','image'); elseif strcmp(h.Type,'image') him = h; else error('Invalid graphics object.') end if isempty(him) error('Image not found.'); end hax = ancestor(him,'axes'); xdata = him.XData; ydata = him.YData; [M,N,~] = size(him.CData); if M > 1 pixel_height = diff(ydata) / (M-1); else % Special case. Assume unit height. pixel_height = 1; end if N > 1 pixel_width = diff(xdata) / (N-1); else % Special case. Assume unit width. pixel_width = 1; end y_top = ydata(1) - (pixel_height/2); y_bottom = ydata(2) + (pixel_height/2); y = linspace(y_top, y_bottom, M+1); x_left = xdata(1) - (pixel_width/2); x_right = xdata(2) + (pixel_width/2); x = linspace(x_left, x_right, N+1); % Construct xv1 and yv1 to draw all the vertical line segments. Separate % the line segments by NaN to avoid drawing diagonal line segments from the % bottom of one line to the top of the next line over. xv1 = NaN(1,3*numel(x)); xv1(1:3:end) = x; xv1(2:3:end) = x; yv1 = repmat([y(1) y(end) NaN], 1, numel(x)); % Construct xv2 and yv2 to draw all the horizontal line segments. yv2 = NaN(1,3*numel(y)); yv2(1:3:end) = y; yv2(2:3:end) = y; xv2 = repmat([x(1) x(end) NaN], 1, numel(y)); % Put all the vertices together so that they can be drawn with a single % call to line(). xv = [xv1(:) ; xv2(:)]; yv = [yv1(:) ; yv2(:)]; hh = hggroup(hax); dark_gray = [0.3 0.3 0.3]; light_gray = [0.8 0.8 0.8]; if isHighDPI bottom_line_width = 1.5; top_line_width = 0.5; else bottom_line_width = 3; top_line_width = 1; end % When creating the lines, use AlignVertexCenters to avoid antialias % effects that would cause some lines in the grid to appear brighter than % others. line(... 'Parent',hh,... 'XData',xv,... 'YData',yv,... 'LineWidth',bottom_line_width,... 'Color',dark_gray,... 'LineStyle','-',... 'AlignVertexCenters','on'); line(... 'Parent',hh,... 'XData',xv,... 'YData',yv,... 'LineWidth',top_line_width,... 'Color',light_gray,... 'LineStyle','-',... 'AlignVertexCenters','on'); % Only return an output if requested. if nargout > 0 hout = hh; end end function tf = isHighDPI % Returns true if running with a high DPI default display. % % LIMITATION: Only works on a Mac with Java AWT available. In other % situations, this function always returns false. tf = false; if ismac if usejava('awt') sf = java.awt.GraphicsEnvironment.getLocalGraphicsEnvironment().getDefaultScreenDevice.getScaleFactor(); tf = sf > 1; end end end

Get
the MATLAB code

Published with MATLAB® R2017b

I happen to know the author, and he's an OK fellow.

]]>Before I dive into today's topic, I want to invite you to come back to the MATLAB Central Blogs page in a couple of weeks. Look for something new to appear related to deep learning.... read more >>

]]>*Before I dive into today's topic, I want to invite you to come back to the MATLAB Central Blogs page in a couple of weeks. Look for something new to appear related to deep learning*.

Earlier today, I was looking at recent image processing questions on MATLAB Answers, and I came across this question:

"Using regionprops to change intensity of blobs"

I thought to myself, "I bet I've written about something like before." It turns out that I have -- but it was more than ten years ago! It is ever-so-slightly possible that a couple of you have not been reading this blog that long, so I thought I would update the post for you today.

The function `regionprops` is very useful for measuring the properties of shapes in a binary image. There are documentation examples and product demos showing how to do this, and I've shown this function in action several times in this blog.

But sometimes I get questions about how to process pixel values in the "original" gray scale image. In other words, suppose your process is something like this:

- Segment gray scale image to get a binary image of objects
- Do something with the original gray scale pixel values corresponding to each blob in the binary image

This post is about how to do that last step. Let's start by creating some simple sample images to work with.

```
I = imread('coins.png');
imshow(I)
```

bw = imbinarize(I); bw = imfill(bw, 'holes'); imshow(bw) title('Thresholded, with holes filled')

The key to processing the pixels in `I` using the blobs in `bw` is to use `regionprops` to get the `'PixelIdxList'` for each blob.

```
s = regionprops(L, 'PixelIdxList')
```

s = 10×1 struct array with fields: PixelIdxList

s(1)

ans = struct with fields: PixelIdxList: [2651×1 double]

`s` is a struct array. `s(k).PixelIdxList` is a vector of the linear indices of the k-th region. You can use these values to index directly into `I`. For example, here is the mean pixel value of the 3rd coin:

idx = s(3).PixelIdxList; mean(I(idx))

ans = 172.8769

You can also replace object pixel values in `I` by indexing on the left side of an assignment statement, like this:

```
% Turn the 4th coin white and the 5th coin black:
J = I;
J(s(4).PixelIdxList) = 255;
J(s(5).PixelIdxList) = 0;
imshow(J)
```

Here's a loop that replaces the pixel values for each coin by the mean pixel value for that coin:

F = I; for k = 1:numel(s) idx = s(k).PixelIdxList; F(idx) = mean(I(idx)); end imshow(F)

I'll finish with a question that Mat asked in a blog comment: "How can I find the location of the maximum pixel value in each labeled region?" Here's how to find the locations and then plot them on top of the original image. (In addition to `'PixelIdxList'`, I'll also ask `regionprops` for `'PixelList'`, which gives the x and y coordinates of each pixel in the region.)

s = regionprops(L, 'PixelIdxList', 'PixelList'); imshow(I) hold on for k = 1:numel(s) idx = s(k).PixelIdxList; pixels_in_region_k = I(idx); [max_value, max_idx] = max(pixels_in_region_k); max_location = s(k).PixelList(max_idx, :); plot(max_location(1),max_location(2),'o','MarkerSize',10,... 'MarkerFaceColor','w','MarkerEdgeColor','k') end hold off

Get
the MATLAB code

Published with MATLAB® R2017a

Today I want to tell you how and why I made these images:... read more >>

]]>Today I want to tell you how and why I made these images:

After the MATLAB R2014b release, I wrote several blog posts (Part 1, Part 2, Part 3, and Part 4) about the new default colormap, parula, that was introduced in that release.

Sometime later, I came across some material by Peter Kovesi about designing perceptually uniform colormaps (or colourmaps, as Peter writes it).

I was particularly intrigued by a test image that Peter designed for the purpose of visually evaluating the perceptual characteristics of a colormap. Here is the test image:

The image is constructed by superimposing a sinusoid on a linear ramp, with the amplitude of the sinusoid getting smaller as you move away from the top row. Here are three cross-sections of the image: row 1, row 64, and row 128.

url = 'http://peterkovesi.com/projects/colourmaps/cm_MATLAB_gray256.png'; I = im2double(imread(url)); subplot(3,1,1) plot(I(1,:)) axis([1 512 0 1]) title('Row 1') subplot(3,1,2) plot(I(64,:)) axis([1 512 0 1]) title('Row 64') subplot(3,1,3) plot(I(128,:)) axis([1 512 0 1]) title('Row 128')

Here is the basic code for making this test image. I'm going to vary Kovesi's image slightly. I'll add an extra half-cycle of the sinusoid so that it reaches a peak at the right side of the image, and I'll add a full-range linear ramp section at the bottom. (If you look closely at the cross-section curves above, you'll see that the linear ramp goes from 5% to 95% of the range.)

First, compute the ramp.

num_cycles = 64.5; pixels_per_cycle = 8; A = 0.05; width = pixels_per_cycle * num_cycles + 1; height = round((width - 1) / 4); ramp = linspace(A, 1-A, width);

Next, compute the sinusoid.

k = 0:(width-1); x = -A*cos((2*pi/pixels_per_cycle) * k);

Now, vary the amplitude of the sinusoid with the square of the distance from the bottom of the image.

q = 0:(height-1); y = ((height - q) / (height - 1)).^2; I1 = (y') .* x;

Superimpose the sinusoid on the ramp.

I = I1 + ramp;

Finally, add a full-range linear ramp section to the bottom of the image.

```
I = [I ; repmat(linspace(0,1,width), round(height/4), 1)];
clf
imshow(I)
title('Colormap test image')
```

Last week, I posted Colormap Test Image to the File Exchange. It contains the function `colormapTestImage`, which does all this for you.

I = colormapTestImage;

The function has another syntax, too. If you pass it the name of a colormap, it will display the test image using that colormap. For example, here is the test image with the old MATLAB default colormap, `jet`.

```
colormapTestImage('jet')
```

This test image illustrates why we replaced `jet` as the default MATLAB colormap. I have annotated the image below to show some of the issues.

Now compare with the new default colormap, `parula`.

```
colormapTestImage('parula')
```

I think that illustrates what we were trying to achieve with `parula`: perceptual fidelity to the data.

Since I'm talking about `parula`, I'll finish by mentioning that we need some very subtle tweaks to `parula` in the R2017b release. So you can compare, I'll show you the original version that shipped with R2014b.

```
colormapTestImage('parula_original')
```

Readers, can you tell what is different? Let us know in the comments.

Get
the MATLAB code

Published with MATLAB® R2017a

I saw a presentation last month that mentioned a user request to have the ability to customize regionprops. That is, a user wanted to be able to add their own measurement to regionprops.... read more >>

]]>I saw a presentation last month that mentioned a user request to have the ability to customize `regionprops`. That is, a user wanted to be able to add their own measurement to `regionprops`.

Today, I'll show you how to do this yourself.

First, here's a brief recap on what `regionprops` does. The function computes measurements of image regions. Some of these measurements are based purely on a region's shape, while others incorporate pixel values within the regions. Here's an example using the coins.png sample image.

```
I = imread('coins.png');
imshow(I)
```

Let's convert this image to binary, using adaptive thresholding, filling holes, and removing small "noise" pixels.

bw = imbinarize(I,'adaptive'); bw = imfill(bw,'holes'); bw = bwareafilt(bw,[100 Inf]); imshow(bw)

You can count the "blobs" (object) yourself; there are 10 of them.

The simplest `regionprops` call, `regionprops(bw)` computes the `Area`, `Centroid`, and `BoundingBox` for each object.

s = regionprops(bw)

s = 10×1 struct array with fields: Area Centroid BoundingBox

But I don't think this is the best way to call `regionprops` anymore. You can now tell `regionprops` to return the results as a table.

```
t = regionprops('table',bw)
```

t = 10×3 table Area Centroid BoundingBox ____ ________________ ____________ 2635 37.133 106.85 [1x4 double] 1846 56.131 49.693 [1x4 double] 2672 96.199 146.05 [1x4 double] 1839 109.97 84.848 [1x4 double] 2744 120.37 208.73 [1x4 double] 2520 148.57 34.404 [1x4 double] 2589 174.83 120.01 [1x4 double] 2518 216.81 70.649 [1x4 double] 1857 236.03 173.36 [1x4 double] 1829 265.96 102.64 [1x4 double]

The table form is a lot more convenient for many tasks. For today's topic, one especially nice thing thing about tables is how easy it is to add your own variables to the table.

To illustrate, let's add a measurement that I've seen called `Roundness`. One definition for roundness is:

$R = \frac{4A}{\pi L^2}$

where $A$ is the object area and $L$ is the major axis length of the best-fit ellipse for the object. Here's how to compute roundness and add it directly to the measurements returned by `regionprops`.

First, note that both `Area` and `MajorAxisLength` are supported by `regionprops`, so let's start with those.

t = regionprops('table',bw,'Area','MajorAxisLength')

t = 10×2 table Area MajorAxisLength ____ _______________ 2635 60.08 1846 50.178 2672 59.792 1839 49.674 2744 60.374 2520 58.08 2589 58.676 2518 58.162 1857 49.77 1829 49.564

You access table variables using dot notation, like `t.Area`. Similarly, you can create a new table variable using dot notation and assignment, like `t.MyVariable = ...`. So adding `Roundness` to the table returned by `regionprops` is this simple.

t.Roundness = 4 * t.Area ./ (pi * t.MajorAxisLength.^2)

t = 10×3 table Area MajorAxisLength Roundness ____ _______________ _________ 2635 60.08 0.92945 1846 50.178 0.93352 2672 59.792 0.9516 1839 49.674 0.94893 2744 60.374 0.9585 2520 58.08 0.95118 2589 58.676 0.95745 2518 58.162 0.94772 1857 49.77 0.95453 1829 49.564 0.94798

Let's try this computation with an image containing objects that are not quite as round.

```
I2 = imread('rice.png');
imshow(I2)
```

bw2 = imbinarize(I2,'adaptive'); bw2 = imfill(bw2,'holes'); bw2 = bwareafilt(bw2,[100 Inf]); imshow(bw2)

t2 = regionprops('table',bw2,'Area','MajorAxisLength'); t2.Roundness = 4 * t2.Area ./ (pi * t2.MajorAxisLength.^2); head(t2)

ans = 8×3 table Area MajorAxisLength Roundness ____ _______________ _________ 138 23.594 0.31562 120 18.152 0.4637 169 28.123 0.27207 157 23.793 0.3531 284 43.757 0.18885 200 26.259 0.36929 141 21.647 0.38311 177 29.087 0.26636

I'm a big fan of the (relatively) new `histogram` function in MATLAB, so let's use it to compare our roundness numbers. I will follow the advice given in the `histogram` reference page for normalizing multiple histograms so that they can be more easily compared. I'll set the y-axis limits to [0 1], which is appropriate for probability normalization, and I'll set the x-axis limits to [0 1], which is the range for `Roundness`.

h1 = histogram(t.Roundness); hold on h2 = histogram(t2.Roundness); hold off h1.Normalization = 'probability'; h2.Normalization = 'probability'; h1.BinWidth = 0.02; h2.BinWidth = 0.02; xlim([0 1]); ylim([0 1]); title('Histogram of roundness (probability normalization)') legend('coins','rice')

There you have it. You can add your own object measurements to the output of `regionprops`. It's especially easy if you tell `regionprops` to return a table.

I'll leave you with this question, dear reader: Are there measurements you would like us to add to `regionprops`? I am aware of an enhancement request for the Feret diameter. What else would you like to see?

Get
the MATLAB code

Published with MATLAB® R2017a