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

A MATLAB user recently contacted tech support with a question about the output of `fft2`. The user had a function, q(x,y), evaluated on an (x,y) grid. The grid evaluation produced a matrix that the user passed to `fft2`. The user's question: what are the spatial frequencies associated with the elements of the output from `fft2`?

The question was forwarded to the MATLAB Math team. Chris Turnes, a fellow Georgia Tech grad, answered it. Here is Chris' answer, lightly edited.

*To map discrete Fourier coefficients into samples of the corresponding Discrete Time Fourier Transform, we need to know about the sampling pattern of the spatial signal. Since the customer presumably has the vectors x and y that she used to query the continuous function q(x,y), this ought to be easy to determine.*

*To determine the sampling frequencies, we use the spacing of the sampling grid. The sampling frequencies in each dimension are fsx = 1/(x(2)-x(1)) and fsy = 1/(y(2)-y(1)) if the grid is uniform in each dimension.*

*Next, the normalized DFT frequencies of a (shifted) N-point DFT (where N here is even for convenience) are fa = ((-N/2):(N/2-1))/N. Therefore, the corresponding DTFT frequency grid would be composed of the frequencies fcx = fa*fsx and fcy = fa*fsy. That would be the correct way to do the mapping between FFT coefficients and the corresponding continuous spatial frequency.*

Thanks, Chris.

Get
the MATLAB code

Published with MATLAB® R2017a

A few weeks ago, I wrote about a solver algorithm for the Eight Queens Problem. The post included diagrams like this one.... read more >>

]]>A few weeks ago, I wrote about a solver algorithm for the Eight Queens Problem. The post included diagrams like this one.

Today I want to show you how I made that diagram in MATLAB. First, let's talk about the board. Back in 2011, I wrote about a variety of ways to make a checkerboard (or chessboard) pattern. In that post, I played games with `repmat`, powers of -1, `floor`, and `round`. It got a little crazy.

I'm still fond of using integer powers of -1, like this.

$(-1)^{(n_1 + n_2)}$

In 2011, I used `ndgrid` to make a matrix of integer powers, $n_1 + n_2$.

[n1,n2] = ndgrid(0:4); (-1).^(n1 + n2)

ans = 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1

Since the introduction of implicit expansion in R2016b, though, I no longer need to use `ndgrid` to explicitly form the matrix of powers.

n = 0:4; (-1).^(n + n')

ans = 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1

In my eight queens code, here's how I got the shades of gray for a chessboard. (There are lots of ways to do this.)

N = 8; dark_square_color = .7; light_square_color = .9; color_range = light_square_color - dark_square_color; f = (-1).^((1:N)' + (1:N)); f = (f + 1) * color_range / 2; f = f + dark_square_color; f(1:4,1:4)

ans = 0.9000 0.7000 0.9000 0.7000 0.7000 0.9000 0.7000 0.9000 0.9000 0.7000 0.9000 0.7000 0.7000 0.9000 0.7000 0.9000

Next, I needed to replicate each element of `f` to make an image with larger squares. Do you know how to replicate elements of a matrix? Some experienced MATLAB users would do it using the `kron` function. Now, however, you can just use the new `repelem` function. The reference page tells you when this function was introduced.

I'm going to make my squares 60 pixels wide, and then I'll turn it into a truecolor image so that the pixel colors are independent of the figure's colormap.

```
pixel_size = 60;
f = repelem(f,pixel_size,pixel_size);
f = repmat(f,1,1,3);
h = imshow(f,'InitialMagnification',100);
```

Now I want to show you a coordinate system trick. If you turn on the axes display, you can see the image pixel coordinates. (In Image Processing Toolbox terminology, these are called *intrinsic coordinates*.)

```
axis on
```

When I get to the part about displaying the queens, however, I would like to be able to place them based on the coordinates of each chessboard square, independent of the number of pixels per square. I can do that by manipulating the `XData` and `YData` properties of the image. These properties assign spatial coordinates to the left/right and top/bottom edges of the image.

h.XData = [0.5 N+0.5]; h.YData = [0.5 N+0.5]; axis([0.5 N+0.5 0.5 N+0.5])

Now let's put some queens on the board. This turns out to be pretty easy because you can just do it with a `text` object. Since R2014b, you've been able to draw text in MATLAB graphics using Unicode characters. And the Unicode character set includes chess symbols! The black queen symbol is at code point 9819 (decimal).

queen = char(9819)

queen = '♛'

Let's put a queen on the square in the second row, third column. (That's $x=3$, $y=2$.)

hq = text(3,2,queen);

Oops, our queen is pretty small, and she also appears to be off-center. Let's fix that.

```
hq.FontSize = 40;
hq.HorizontalAlignment = 'center';
```

That's better.

I have one more little coding trick to show you. When I wrote my eight queens animation code, I wanted to have an easy to show and remove a queen at any square on the board. So, I made an matrix of text objects, one at each square. Then I could just index into the matrix of text objects and turn the `Visible` property on and off. Here's how. (First, let me delete the text object I just made.)

delete(hq) for r = 1:N for c = 1:N queens(r,c) = text(c,r,char(9819),... 'HorizontalAlignment','center',... 'Visible','off',... 'FontSize',40); end end

OK, let's turn on the queen displays on a couple of squares.

queens(5,3).Visible = 'on'; queens(2,4).Visible = 'on';

For a final bit of fun, the function `set` lets you modify properties of a whole array of graphics objects at once.

Turn them all on!

set(queens,'Visible','on')

(The board above, by the way, is NOT a solution to the Eight Queens Problem. In case you were wondering.)

Get
the MATLAB code

Published with MATLAB® R2017a

On April 20, I wrote about an algorithm for solving the Eight Queens Problem. In this problem, the task is to place eight queens on a chessboard so that none of the queens is attacking any other queen. Here is one solution.... read more >>

]]>On April 20, I wrote about an algorithm for solving the Eight Queens Problem. In this problem, the task is to place eight queens on a chessboard so that none of the queens is attacking any other queen. Here is one solution.

After I published that post, I became curious to see what others might have submitted about this problem on the File Exchange. A little searching brought me to this submission by uu tii. When I looked at the submission, I was surprised to how little code it contained. This is it:

function M = eightqueens() N=8; T=N-2; rows=1:N; cols=perms(rows); S=size(cols,1); M=zeros(N,N,S); linearInd = sub2ind(size(M), repmat(rows',1,S), cols', repmat(1:S,N,1)); M(linearInd) = 1; dv=arrayfun(@(k)max([arrayfun(@(x)sum(diag(M(:,:,k),x)),-T:T),arrayfun(@(x)sum(diag(rot90(M(:,:,k)),x)),-T:T)]),1:S); M(:,:,dv>1)=[];

When I ran it, it returned an 8x8x92 array. The size caught my eye because there are exactly 92 solutions to the Eight Queens Problem (including rotations and flips).

```
M = eightqueens;
whos M
```

Name Size Bytes Class Attributes M 8x8x92 47104 double

Let's peek at one of the planes of `M`.

M(:,:,45)

ans = 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0

Indeed, that is a solution. So, what's going on in this small chunk of code? Well, as I see it, this code has three conceptual chunks. The first chunk, consisting of eight lines, is the longest. Let's run it.

```
N=8;
T=N-2;
rows=1:N;
cols=perms(rows);
S=size(cols,1);
M=zeros(N,N,S);
linearInd = sub2ind(size(M), repmat(rows',1,S), cols', repmat(1:S,N,1));
M(linearInd) = 1;
whos M
```

Name Size Bytes Class Attributes M 8x8x40320 20643840 double

This seems to be a collection of 40320 boards. There aren't that many solutions, so the collection must contain many boards that aren't valid solutions. For example:

M(:,:,45)

ans = 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0

There are four queens lined up on the main antidiagonal, so this board is not a solution.

This set of boards has been constructed by making the observation that any solution must have exactly one queen on each row and exactly one queen on each column. The array `M` contains every board that has this property. All the computational magic here is in the line `cols=perms(rows)`. The rest of this chunk is some indexing computations to initialize and fill in the array `M` with queens at the locations given by `rows` and `cols`.

The second chunk is this line:

dv=arrayfun(@(k)max([arrayfun(@(x)sum(diag(M(:,:,k),x)),-T:T),arrayfun(@(x)sum(diag(rot90(M(:,:,k)),x)),-T:T)]),1:S);

Wow, there's a lot going on there. I would summarize it this way: For every board `k` in M, sum every diagonal and every antidiagonal and take the maximum of those sums. The idea is that `M` was constructed so that every board has exactly one queen on each row and exactly one queen on each column, so you really only have to check the diagonals of a board to determine whether it is a solution. Let's look at `dv` for the 45th board, which I displayed above.

dv(45)

ans = 4

Remember my comment above that board 45 has 4 queens lined up on the main antidiagonal? That's why `dv(45)` is 4.

The last conceptual chunk is eliminate all the boards in the `M` that are not solutions. We do that by eliminating all boards for which `dv > 1`.

```
M(:,:,dv>1)=[];
whos M
```

Name Size Bytes Class Attributes M 8x8x92 47104 double

And there it is. That's a very clever use of `perms` and indexing and `arrayfun` to generate all possible solutions. Thanks, uu tii.

Get
the MATLAB code

Published with MATLAB® R2017a