Implicit functions and fun with contours

Stefan Güttel, 30th July 2012

(Chebfun example fun/ContourFun.m)
[Tags: #implicitfunction, #fun, #contour]

close all
clear all
LW = 'linewidth'; lw = 2.5;
MS = 'markersize'; ms = 10;

Matlab provides a function called contours, which returns a vector of points on level lines formed by the entries of a matrix. This function can easily be used to define chebfuns satisfying an implicit equation, at least approximately. To demonstrate this, suppose we want to create a chebfun that corresponds to a circle of radius 0.8 centered at the origin. The points on this circle satisfy z(x,y) = 0 with

z = @(x,y) x.^2 + y.^2 - 0.8.^2;

We now sample z at 20 x 20 grid points on the square [-1,1] x [1,1], and use Matlab's contours to detect the level lines of the circle:

x = linspace(-1,1,20);
[X,Y] = meshgrid(x,x);
Z = z(X,Y);
C = contours(X,Y,Z,[0,0]);

The following plot shows that the algorithm does quite a good job by actually interpolating between the grid points:

plot(.8*exp(2i*pi*linspace(0,1,200)),'k-',LW,lw)
hold on, axis square
plot(X(:),Y(:),'b+',MS,ms)
plot(C(1,2:end),C(2,2:end),'ro',MS,ms)

The following lines will convert the coordinates of the contours to a quasimatrix f which each column corresponding to a contour (in this case, there is only one contour):

j = 1; f = chebfun();
while j < length(C),
    k = j + C(2,j); D = C(:,j+1:k); j = k + 1;
    f = [ f , chebfun(D(1,:)+1i*D(2,:)) ];
end
figure
plot(f,'r',LW,lw)
axis square, grid on

Of course, we cannot expect this to be a perfect circle as we have only sampled at 20 grid points in each coordinate direction, but for some practical purposes this accuracy may be sufficient:

norm(z(real(f),imag(f)),inf)

ans =
   0.005386578687440

One example where high precision is certainly not necessary is when a function has do be defined from an image. To demonstrate this, let's read an image showing a whiteboard in our visitors office. I'm grateful to Yuji Nakatsukasa for taking this snapshot.

A = imread('stefun.jpg');
imshow(A)
axis on, hold on

In the following we will actually operate on a grayscale version of this image, which is computed by the following two lines:

A = 0.299*double(A(:,:,1)) + 0.587*double(A(:,:,2)) + 0.114*double(A(:,:,3));
A = A/max(max(A));

We crop the image to the region of the graph on the upper right, and compute contour lines using the contours command:

x = 185:463;
y = 60:210;
B = A(y,x);
[X,Y] = meshgrid(x,y);
C = contours(X,Y,B,[0.8,0.8]);

As above, the following lines will convert the coordinates of the contours to a quasimatrix f:

j = 1; f = chebfun();
while j < length(C),
    k = j + C(2,j); D = C(:,j+1:k); j = k + 1;
    f = [ f , chebfun(D(1,:)+1i*D(2,:)) ];
end

It turns out that f has four columns, and there are pairs of chebfuns with function values very close to each other. These pairs correspond to the upper and lower edge of a single line on the whiteboard. By taking the arithmetic mean of two neighbouring chebfuns we obtain a pretty good agreement with the curves on the board:

xaxis = (f(:,1) + f(:,2))/2;
curve = (f(:,3) + f(:,4))/2;
plot([xaxis,curve],LW,lw)

Since these are chebfuns, we can do all kinds of things with them. For example, find the area enclosed by the x-axis and the curve. Note that the real part of each chebfun will correspond to the x-coordinate, and the imaginary part is the y-coordinate measured in pixels. The unit of area therefore is square-pixels.

meanaxis = 1i*mean(imag(xaxis));
curve2 = curve - meanaxis;
Area = abs(sum(imag(curve2).*real(curve2)/2));
sprintf('The enclosed area is %0.1f square-pixels.',Area)

ans =
The enclosed area is 34484.6 square-pixels.

Or, we could compute a degree 3 best polynomial approximation to the curve, overlay it with the plot, and show the intersections of both functions:

p = 1i*remez(imag(curve),3) + real(curve);
plot(p,'r--',LW,lw)
rts = roots(curve - p);
plot(p(rts),'ro',MS,ms)

One may wonder, how close is the lower curve on the whiteboard to a circle? In order to answer this question we again crop to the region of interest, and compute contours, which in this case gives a quasimatrix with two columns corresponding to the edges of the painted curve. Let's take the first column only as a prediction to the curve:

x = 310:470;
y = 240:410;
B = A(y,x);
[X,Y] = meshgrid(x,y);
C = contours(X,Y,B,[0.8,0.8]);
j = 1; f = chebfun();
while j < length(C),
    k = j + C(2,j); D = C(:,j+1:k); j = k + 1;
    f = [ f , chebfun(D(1,:)+1i*D(2,:)) ];
end
curve = f(:,1);
plot(curve,'m',LW,lw)

We can compute the area and centroid of this curve (see also [1]), and draw a perfect circle with the same area for comparison:

Area = sum(real(curve).*diff(imag(curve)));
centroid = sum(diff(curve).*curve.*conj(curve))/(2i*Area);
radius = sqrt(Area/pi);
plot(centroid,'m+',MS,ms)
circle = chebfun('exp(1i*pi*x)');
circle = radius*circle + centroid;
plot(circle,'k--',LW,lw)

Last and least, we go crazy and approximate ourselfs via a collection of quasimatrices plotted with three shades of gray:

x = 10:170;
y = 90:310;
B = A(y,x);
[X,Y] = meshgrid(x,y);
levels = [.2 .4 .8];
C = contours(X,Y,B,levels);
j = 1; f = chebfun(); level = [];
while j < length(C),
    k = j + C(2,j);  D = C(:,j+1:k); level = [ level , C(1,j) ];
    j = k + 1; f = [ f , chebfun(D(1,:)+1i*D(2,:)) ];
end
figure
for j = 1:length(levels),
    ind = find(level == levels(j));
    plot(f(:,ind),'Color',[levels(j) levels(j) levels(j)],LW,lw)
    hold on
end
axis ij, axis image

I refrain from doing any calculations on this!

References:

[1] http://www2.maths.ox.ac.uk/chebfun/examples/geom/html/Area.shtml