polygonCurvature

POLYGONCURVATURE Estimate curvature on polygon vertices using polynomial fit.

   CURV = polygonCurvature(POLY, M)
   Estimate the curvature for each vertex of a polygon, using polynomial
   fit from the M verties located around current vertex. M is usually an
   odd value, resulting in a symmetric neighborhood.

   Polynomial fitting is of degree 2 by default.
   

   Example
     img = imread('circles.png');
     img = imfill(img, 'holes');
     imgf = imfilter(double(img), fspecial('gaussian', 7, 2));
     figure(1), imshow(imgf);
     contours = imContours(imgf, .5); poly = contours{1};
     poly2 = smoothPolygon(poly, 7);
     hold on; drawPolygon(poly2);
     curv = polygonCurvature(poly2, 11);
     figure; plot(curv);
     minima = bwlabel(imextendedmin(curv, .05));
     centroids = imCentroid(minima);
     inds = round(centroids(:,2));
     figure(1); hold on; drawPoint(poly2(inds, :), 'g*')

   See also
     polygons2d

0001 function curv = polygonCurvature(poly, M)
0002 %POLYGONCURVATURE Estimate curvature on polygon vertices using polynomial fit.
0003 %
0004 %   CURV = polygonCurvature(POLY, M)
0005 %   Estimate the curvature for each vertex of a polygon, using polynomial
0006 %   fit from the M verties located around current vertex. M is usually an
0007 %   odd value, resulting in a symmetric neighborhood.
0008 %
0009 %   Polynomial fitting is of degree 2 by default.
0010 %
0011 %
0012 %   Example
0013 %     img = imread('circles.png');
0014 %     img = imfill(img, 'holes');
0015 %     imgf = imfilter(double(img), fspecial('gaussian', 7, 2));
0016 %     figure(1), imshow(imgf);
0017 %     contours = imContours(imgf, .5); poly = contours{1};
0018 %     poly2 = smoothPolygon(poly, 7);
0019 %     hold on; drawPolygon(poly2);
0020 %     curv = polygonCurvature(poly2, 11);
0021 %     figure; plot(curv);
0022 %     minima = bwlabel(imextendedmin(curv, .05));
0023 %     centroids = imCentroid(minima);
0024 %     inds = round(centroids(:,2));
0025 %     figure(1); hold on; drawPoint(poly2(inds, :), 'g*')
0026 %
0027 %   See also
0028 %     polygons2d
0029  
0030 % ------
0031 % Author: David Legland
0032 % e-mail: david.legland@inra.fr
0033 % Created: 2018-03-02,    using Matlab 9.3.0.713579 (R2017b)
0034 % Copyright 2018 INRA - Cepia Software Platform.
0035 
0036 % number of vertices of polygon
0037 n = size(poly, 1);
0038 
0039 % allocate memory for result
0040 curv = zeros(n, 1);
0041 
0042 % number of vertices before and after current vertex
0043 s1 = floor((M - 1) / 2);
0044 s2 = ceil((M - 1) / 2);
0045 
0046 % parametrisation basis
0047 % As we recenter the points, the constant factor is omitted
0048 ti = (-s1:s2)';
0049 X = [ti ti.^2];
0050     
0051 % Iteration on vertex indices
0052 for i = 1:n
0053     % coordinate of current vertex, for recentring neighbor vertices
0054     x0 = poly(i,1);
0055     y0 = poly(i,2);
0056     
0057     % indices of neighbors
0058     inds = i-s1:i+s2;
0059     inds = mod(inds-1, n) + 1;
0060     
0061     % Least square estimation using mrdivide
0062     xc = X \ (poly(inds,1) - x0);
0063     yc = X \ (poly(inds,2) - y0);
0064     
0065     % compute curvature
0066     curv(i) = 2 * (xc(1)*yc(2) - xc(2)*yc(1) ) / power(xc(1)^2 + yc(1)^2, 3/2);
0067 end
0068 
0069