Playing with conic sections in homogeneous coordinates

Contents

Preparation

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

figure(1)
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,100,'filled');

click 5 points, then enter

x and y are now column vectors.

x
y

x =

  196.0000
  432.0000
  341.0000
  144.0000
   94.0000


y =

  146.0000
  253.0000
  399.0000
  282.0000
  182.0000

Estimation of conic parameters

using the technique shown in "Multiple View Geometry", page 9. questo equivale ad impostare il sistema lineare in cui la conica passa per 5 punti

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

% 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.
N = null(A);

A =

   1.0e+05 *

    0.3842    0.2862    0.2132    0.0020    0.0015    0.0000
    1.8662    1.0930    0.6401    0.0043    0.0025    0.0000
    1.1628    1.3606    1.5920    0.0034    0.0040    0.0000
    0.2074    0.4061    0.7952    0.0014    0.0028    0.0000
    0.0884    0.1711    0.3312    0.0009    0.0018    0.0000

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

cc = N(:, 1);
[a b c d e f]=deal(cc(1),cc(2),cc(3),cc(4),cc(5),cc(6));
% here is the matrix of the conic
C=[a b/2 d/2; b/2 c e/2; d/2 e/2 f]

% Remark: since the right null space has dimension one,
% the system admits an infinite number of solutions
% however, these can be expressed as lambda * n, where n
% lambda \in R and n = null(A).
% thus, they all corresponds to the same conic.

C =

    0.0000   -0.0000   -0.0003
   -0.0000    0.0000   -0.0044
   -0.0003   -0.0044    1.0000

Draw the conic

We pick every pixel, and compute the incidence relation.

im=zeros(500,500);
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.

figure(1)
imshow(im > 0);
scatter(x,y, 100, 'filled');

N.B: if we provide only 4 points, N is a matrix

which contains in each column a generator of the null space

N = null(A(1 : 4, :));

% then, any vector living in the null space give rise to a conic
% passing throguh those 4 points

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

cc = N(:, 2);
% cc = N(:, 1) + N(:, 2); % another example, any linear combination is fine
[a b c d e f]=deal(cc(1),cc(2),cc(3),cc(4),cc(5),cc(6));
% here is the matrix of the conic
C_four =[a b/2 d/2; b/2 c e/2; d/2 e/2 f]

C_four =

    0.0000    0.0000   -0.0030
    0.0000   -0.0000   -0.0011
   -0.0030   -0.0011    1.0000

Draw the conic

We pick every pixel, and compute the incidence relation.

im_four=zeros(500,500);
for i=1:500
    for j=1:500
        im_four(i,j)=[j i 1]*C_four*[j i 1]';
    end
end

figure(2),
imshow(im_four > 0);
hold on
scatter(x(1 : 4), y(1 : 4), 100, 'filled');
hold off

Tangent lines

Let's draw the tangents to our points: the tangent line in homogeneous coordinates is l the line tangent to C at the line x is defined as l = C*x

% this is the line tangent to the first point
l1 = C * [x(1); y(1); 1];
% here I compute all the tangent lines simultaneously
% stacking all the vectors in the columns of the matrix l
l=C*[x y ones(size(x))]'

l =

    0.0003    0.0019   -0.0010   -0.0019   -0.0013
   -0.0025   -0.0023    0.0029    0.0019   -0.0003
    0.3081   -0.2258   -0.8345   -0.2685    0.1810

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,:));
theta = theta';
figure(1), hold on
plot([x-cos(theta)*100 x+cos(theta)*100]',...
    [y-sin(theta)*100 y+sin(theta)*100]','b-','LineWidth',5);

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 =

  -2.7756e-16


ans =

  -2.8866e-15


ans =

  -1.1102e-16


ans =

   7.7716e-16


ans =

   6.3838e-16

Published with MATLAB® R2017b