quad2d.m

this a SVM toolbox,it is very useful for someone who just learn SVM.In order to be undestood easily,

function [X1,Y1,X2,Y2]=quad2d(A,B,C,win,interp)
% QUAD2D computes contour of quadratic function in 2D.
%  [X1,Y1,X2,Y2]=quad2d(A,B,C,win,interp)
%
% QUAD2D computes points on countour of 2D quadratic function
%         x'*A*x + B'*x + C = 0,       for x in win.
%
%  The computed curve(s) are returned in the format useful 
%  for Matlab 'plot' function.
%
% Input:
%   A [2x2], B[2x1], C[1x1] parameters of 2D quadratic function.
%   win [1x4] defines function domains.
%   interp [1x1] number of lines used for interpolation.
%
% Output:
%   X1 [1x(interp+1)], Y1[1x(interp+1)] points on the 1st countour.
%   X2 [1x(interp+1)], Y2[1x(interp+1)] points on the 2nd countour.
%
% Example:
%
%   A=-eye(2,2); B=[2 2]; C=[1];
%   [X1,Y1,X2,Y2]=quad2d(A,B,C,[-1 1 -1 1]);
%   figure; hold on;
%   plot(X1,Y1);
%   plot(X2,Y2);
%
% See also L2Q2D, PQUAD2D, QTRANSF.
%

% Statistical Pattern Recognition Toolbox, Vojtech Franc, Vaclav Hlavac
% (c) Czech Technical University Prague, http://cmp.felk.cvut.cz
% Written Vojtech Franc (diploma thesis) 19.11.1999, 23.12.1999, 5.4.2000
% Modifications
% 26-June-2001, V.Franc, comments repared.
% 24. 6.00 V. Hlavac, comments polished.

% default setting
if nargin < 4,
   win=[0 1 0 1];
end
if nargin < 5,
   interp=75;
end

% x-limits
xmax=win(2);
xmin=win(1);

% y-limits
ymax=win(4);
ymin=win(3);


%%%%%%%%%%%%%%%%%%%%%%%%5
a=A(1,1);
b=B(1);
c=A(1,2)+A(2,1);
d=B(2);
e=A(2,2);
f=C;

alpha=c^2-4*e*a;
beta1=2*c*d-4*e*b;
beta2=2*b*c-4*a*d;
gama1=d^2-4*e*f;
gama2=b^2-4*a*f;

Ddx=beta1^2-4*alpha*gama1;
Ddy=beta2^2-4*alpha*gama2;

if Ddx >0 & Ddy > 0,
   yy1=(-beta2+sqrt(Ddy))/(2*alpha);
   yy2=(-beta2-sqrt(Ddy))/(2*alpha);
   if alpha < 0,
      ymin=max(ymin,yy1);
      ymax=min(ymax,yy2);
   else
      ymin2=max(ymin,min(yy2,ymax));
      ymax2=min(ymax,max(yy1,ymin));
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%5

% curvature
X1=[]; X2=[];
Y1=[]; Y2=[];

% go trough y value
if Ddx < 0,  % Ddx > 0 -> x=(-inf,inf)
   step=(xmax-xmin)/interp;
   for x=xmin:step:xmax,
      ay=e;
      by=c*x+d;
      cy=a*x^2+b*x+f;

      D=by^2-4*ay*cy;
      if D > 0,
         X1=[X1,x];
         X2=[X2,x];
         Y1=[Y1,(-by-sqrt(D))/(2*ay)];
         Y2=[Y2,(-by+sqrt(D))/(2*ay)];
      elseif D==0,
         X1=[X1,x];
         X2=[X2,x];
         if ay~=0,
            Y1=[Y1,-by/(2*ay)];
            Y2=[Y2,-by/(2*ay)];
         else
            Y1=[Y1,-cy/by];
            Y2=[Y2,-cy/by];
         end
      end
   end
elseif Ddy < 0 | alpha < 0, % Ddy < 0 -> y=(-inf,inf)
                            % or alpha < 0 & Ddy > 0 & Ddx > 0 -> ellipse
   step=(ymax-ymin)/interp;
   for y=ymin:step:ymax,
      ax=a;
      bx=b+c*y;
      cx=e*y^2+d*y+f;

      D=bx^2-4*ax*cx;
      if D > 0,
         Y1=[Y1,y];
         Y2=[Y2,y];
         X1=[X1,(-bx-sqrt(D))/(2*ax)];
         X2=[X2,(-bx+sqrt(D))/(2*ax)];
      elseif D==0,
         Y1=[Y1,y];
         Y2=[Y2,y];
         if ax~=0,
            X1=[X1,-bx/(2*ax)];
            X2=[X2,-bx/(2*ax)];
         else
            X1=[X1,-cx/bx];
            X2=[X2,-cx/bx];
         end
      end
   end
   % connect curvature
   if Ddx > 0 & Ddy > 0 & size(X2,2) > 0,
      if ymin > win(3),
         X1=[X2(1),X1];
         Y1=[Y2(1),Y1];
      end
      if ymax < win(4),
         X1=[X1,X2(size(X2,2))];
         Y1=[Y1,Y2(size(Y2,2))];
      end
   end

elseif alpha > 0,
   step=(ymin2-ymin)/interp;
   for y=ymin:step:ymin2,
      upc=1;
      ax=a;
      bx=b+c*y;
      cx=e*y^2+d*y+f;

      D=bx^2-4*ax*cx;
      if D > 0,
         Y1=[Y1,y];
         Y2=[Y2,y];
         X1=[X1,(-bx-sqrt(D))/(2*ax)];
         X2=[X2,(-bx+sqrt(D))/(2*ax)];
      elseif D==0,
         Y1=[Y1,y];
         Y2=[Y2,y];
         if ax~=0,
            X1=[X1,-bx/(2*ax)];
            X2=[X2,-bx/(2*ax)];
         else
            X1=[X1,-cx/bx];
            X2=[X2,-cx/bx];
         end
      end
   end
   % connect curvature
   uindex=size(X2,2);
   if step > 0 & ymin2 < ymax,
      X1=[X1,X2(uindex)];
      Y1=[Y1,Y2(uindex)];
   end
   step=(ymax-ymax2)/interp;
   first=1;
   for y=ymax2:step:ymax,
      ax=a;
      bx=b+c*y;
      cx=e*y^2+d*y+f;

      D=bx^2-4*ax*cx;
      if D > 0,
         if first==1 & ymin < ymax2,
            first=0;
            Y1=[Y1,y];
            X1=[X1,(-bx+sqrt(D))/(2*ax)];
         end
         Y1=[Y1,y];
         Y2=[Y2,y];
         X1=[X1,(-bx-sqrt(D))/(2*ax)];
         X2=[X2,(-bx+sqrt(D))/(2*ax)];
      elseif D==0,
         if first==1 & ymin < ymax2,
            first=0;
            Y1=[Y1,y];
            X1=[X1,-bx/(2*ax)];
         end
         Y1=[Y1,y];
         Y2=[Y2,y];
         if ax~=0,
            X1=[X1,-bx/(2*ax)];
            X2=[X2,-bx/(2*ax)];
         else
            X1=[X1,-cx/bx];
            X2=[X2,-cx/bx];
         end
      end

   end
end


return;