smooth image outline in matlab by interpolation

Question

I'm plotting image outline points using this threshold method, but my outline has straight line segments. I want to plot the angle to the vertical at each point, so I really need curves.

I can get smooth curves using a convex hull.

The image was generated as follows:

B = bwboundaries(BW3);
outline = B{1,1};
plot(outline(:,2),outline(:,1),'r.','LineWidth',1) 
K = convhull(outline(:,2),outline(:,1)); 
plot(outline(K,2),outline(K,1),'b+--','LineWidth',1)

But how can I "fill in the gaps" between the convex hull points? I want a point on the blue curve for every red point.

I tried to achieve this using interp1:

outline2 = outline;
outline2(:,2)=interp1(outline(K,1),outline(K,2),outline(:,1),'spline');

but got the following error: "Error using griddedInterpolant The grid vectors must contain unique points."

I assume it's because the outline forms a loop, not a unique x point for each y. Is there a different way to fill in those missing points using a spline?

I'm also open to other ideas for finding a smooth edge.

Thanks for any help!

Answer 1

Find the initial outline with Marching Squares (https://en.wikipedia.org/wiki/Marching_squares#Isoline).

Then if you want a good estimate of the derivatives, fit an interpolating Cubic Spline (https://en.wikipedia.org/wiki/Spline_interpolation#Algorithm_to_find_the_interpolating_cubic_spline).

Here there is a little technicality: it seems that you want the slope at the centers of the pixels. But the spline obtained from marching cubes will pass through known point the edges and not through the centers. You can

• project the center point to the nearest point on the spline (unfortunately requiring the solution of a higher degree polynomial);

• implicitize the cubic arc and compute the gradient of the implicit function.

If your accuracy requirements are not stringent, you may probably do by just using the segment direction in the tagreted pixel.

Answer 2

Since your image seems smooth and well sampled, I would suggest that you find, for each edge pixel, the sub-pixel location of the true edge. With this we remove the need for the convex hull, which might be useful for your particular image, but does't generalize to arbitrary shapes.

Here is some code to accomplish what I suggest.

% A test image in the range 0-1, the true edge is assumed to be at 0.5
img = double(gaussianedgeclip(60-rr));

% Get rough outline
p = bwboundaries(img>0.5);
p = p{1,1};

% Refine outline
n = size(p,1);
q = p;          % output outline
for ii=1:n
   % Find the normal at point p(ii,:)
   if ii==1
      p1 = p(end,:);
   else
      p1 = p(ii-1,:);
   end
   if ii==n
      p2 = p(1,:);
   else
      p2 = p(ii+1,:);
   end
   g = p2-p1;
   g = (g([2,1]).*[-1,1])/norm(g);
   % Find a set of points along a line perpendicular to the outline
   s = p(ii,:) + g.*linspace(-2,2,9)';
         % NOTE: The line above requires newer versions of MATLAB. If it
         % fails, use bsxfun or repmat to compute s.
   v = interp2(img,s(:,2),s(:,1));
   % Find where this 1D sample intersects the 0.5 point,
   % using linear interpolation
   if v(1)<0.5
      j = find(v>0.5,1,'first');
   else
      j = find(v<0.5,1,'first');
   end
   x = (v(j-1)-0.5) / (v(j-1)-v(j));
   q(ii,:) = s(j-1,:) + (s(j,:)-s(j-1,:))*x;
end

% Plot
clf
imshow(img,[])
hold on
plot(p(:,2),p(:,1),'r.','LineWidth',1) 
plot(q(:,2),q(:,1),'b.-','LineWidth',1) 
set(gca,'xlim',[68,132],'ylim',[63,113])

The first line, that generates the test image, requires DIPimage, but the rest of the code uses only standard MATLAB functions, except bwboundaries, which you were using also and is from the Image Processing Toolbox.

The output set of points, q, is not sampled at integer x or y. That is a whole lot more complicated to accomplish.

Also, sorry for the one-letter variables... :)