points_from_conic(A,B,C,D,E,F, domain,number_points, noise_power)

% Juan Manuel Perez Rua ([e-mail])

% This function generates points from a generic conic given by
% parameters A,B,C,D,E and F of its canonical equation:
% A*x^2 + 2*B*x*y + C*y^2 + 2*D*x + 2*E*y + F = 0
% The points can be noisy, and the power of the noise can be configured 
% with noise_power variable (std. deviation of gaussian noise).
% The number of output points can be configured with number_points

% xpoints are computed depending on the given domain and number_points
% ypoitns are chosen to solve the canonic equation with xpoints.
% If a combinations of parameters give a unsolvable condition on the real
% field, the complex numbers are ignored. Thus, sometimes a smaller number
% of points can be get.

% EXAMPLES:
% [x,y] = points_from_conic;
% plot(x,y,'g.'); title('CIRCLE')

% [x,y]=points_from_conic(3,4,5,3,0,1,[-10 0],200,0.0);
% plot(x,y,'r*');  title('HYPERBOLA')

% [x,y]=points_from_conic(1,0,1,0,0,-1,[-2 2],100,0.05);
% plot(x,y,'bo');  title('Noisy ELLIPSE')

function [xpoints,ypoints] = points_from_conic(A,B,C,D,E,F, domain,number_points, noise_power)
    % Parsing parameters
    switch nargin
        case 0
            A=1;
            B=0;
            C=1;
            D=0;
            E=0;
            F=-1;
            domain=[-1,1];
            number_points=50;
            noise_power=0;
        case 7
            number_points=50;
            noise_power=0;
        case 8
            noise_power=0;
        case 9
        otherwise
            error('Wrong number of parameters. Expected 0, 6, 7 or 8 parameters.')
    end
    
    % Generating domain vector
    x = domain(1):(domain(2)-domain(1))/round(number_points/2):domain(2)-1/round(number_points/2);
    % Expression to solve
    exp = strcat(num2str(A),'*x^2+2*',num2str(B),'*x*y+',num2str(C),'*y^2+2*',...
        num2str(D),'*x+2*',num2str(E),'*y+',num2str(F),'=0');
    yy = solve(exp,'y');
    y = eval(yy);
 
    % Point generation
    xpoints=0;
    ypoints=0;
    cont=0;
    for k=1:size(y,2)
        for m=1:size(y,1)
           if ( imag(y(m,k))==0 )
               cont=cont+1;
               xpoints(cont)=x(k)    + noise_power*randn;
               ypoints(cont)=y(m,k)  + noise_power*randn;
           end
        end
    end
    
    
end