Playing with conic sections in homogeneous coordinates

Preparation
Estimation of conic parameters
Draw the conic
Improvements

Preparation

First, we create an empty image and display it (will be our "canvas")

im=zeros(500,500);
imshow(im);
hold on;

Data input: enter five points (then enter)

disp('click 5 points, then enter');
[x y]=getpts;
scatter(x,y,25,'filled');

click 5 points, then enter

x and y are now column vectors.

x
y

Estimation of conic parameters

using the technique shown in "Multiple View Geometry", page 31.

A=[x.^2 x.*y y.^2 x y ones(size(x))]

A =

   1.0e+05 *

  Columns 1 through 5

    0.3725    0.6774    1.2320    0.0019    0.0035
    0.0902    0.2375    0.6250    0.0009    0.0025
    0.1588    0.1663    0.1742    0.0013    0.0013
    0.4580    0.2247    0.1102    0.0021    0.0010
    1.6160    1.0934    0.7398    0.0040    0.0027

  Column 6

    0.0000
    0.0000
    0.0000
    0.0000
    0.0000

A is 5x6: we find the parameter vector (containing [a b c d e f]) as the right null space of A. This returns a vector cc such that A*cc=0. Note that this expresses the fact that the conic passes through all the points we inserted.

cc=null(A);

Let's assign the values to variables and build the 3x3 conic matrix (which is symmetrical)

[a b c d e f]=deal(cc(1),cc(2),cc(3),cc(4),cc(5),cc(6));
C=[a b/2 d/2; b/2 c e/2; d/2 e/2 f]

C =

    0.0000   -0.0000   -0.0021
   -0.0000    0.0000   -0.0031
   -0.0021   -0.0031    1.0000

Draw the conic

We pick every pixel, and compute the incidence relation.

for i=1:500 for j=1:500
im(i,j)=[j i 1]*C*[j i 1]';
end
end

Let's draw black for pixels less than zero, white for pixels greater than 0.

imshow(im>0);
scatter(x,y,25,'filled');

Improvements

Let's draw the tangents to our points: the tangent line in homogeneous coordinates is l

l=C*[x y ones(size(x))]'

l =

   -0.0009   -0.0018   -0.0011    0.0000    0.0017
    0.0020    0.0006   -0.0014   -0.0021    0.0001
   -0.5162    0.0111    0.3153    0.2114   -0.7161

in order to draw it we find its angle, then plot a line of length 200 around it. In order to write less lines of code, I'm doing a bit of hacks with matrices and vectors here: theta will be a column vector of angles, and with a single plot command I'll plot all the lines. Check out the documentation of plot to understand what's happening.

theta=atan2(-l(1,:),l(2,:))
plot([x-cos(theta')*100 x+cos(theta')*100]',[y-sin(theta')*100 y+sin(theta')*100]','b-','LineWidth',5);

theta =

    0.4356    1.2299    2.4781   -3.1321   -1.5074

Let's also draw the pole of the line at the infinity w.r.t. the conic.

linf=[0 0 1]';
center=inv(C)*linf;
center=center/center(3);
scatter(center(1),center(2),25,'filled');

Let's compute the dual conic and check that all the tangent lines we computed previously actually satisfy the relation (up to numeric precision issues).

Cstar=inv(C);
l(:,1)'*Cstar*l(:,1)
l(:,2)'*Cstar*l(:,2)
l(:,3)'*Cstar*l(:,3)
l(:,4)'*Cstar*l(:,4)
l(:,5)'*Cstar*l(:,5)

ans =

   1.1102e-16


ans =

   4.1633e-16


ans =

  -2.7756e-16


ans =

  -8.3267e-17


ans =

   1.1449e-16

Playing with conic sections in homogeneous coordinates

Contents

Preparation

Estimation of conic parameters

Draw the conic

Improvements