sugr_plot_conic_explicit

% explicit plot of conic (ellipse or hyperbola)

 h = plot_conic_explicit( C, ... ) plots the conic defined by the 3x3 matrix
 C explicitly. 

 * C = conic
 * center_type = string specifying plot type of center
 * bound_type = string specifying plot type of conic
 * linewidth = width of line of conic
 * factor = magification of conic compared to position
 * se = [start,end] of hyperbola for line segments
        se=[0,0]'  then line segment (+-sqrt(2))
        se=[1,1]'  then line until xlim,ylim (default)

 $Log: plot_conic_explicit.m,v $
 Revision 1.1  2009/08/03 11:10:33  Meidow
 *** empty log message ***

 adapted
 Wolfgang Förstner 1/2011
 wfoerstn@uni-bonn.de

0001 %% explicit plot of conic (ellipse or hyperbola)
0002 %
0003 % h = plot_conic_explicit( C, ... ) plots the conic defined by the 3x3 matrix
0004 % C explicitly.
0005 %
0006 % * C = conic
0007 % * center_type = string specifying plot type of center
0008 % * bound_type = string specifying plot type of conic
0009 % * linewidth = width of line of conic
0010 % * factor = magification of conic compared to position
0011 % * se = [start,end] of hyperbola for line segments
0012 %        se=[0,0]'  then line segment (+-sqrt(2))
0013 %        se=[1,1]'  then line until xlim,ylim (default)
0014 %
0015 % $Log: plot_conic_explicit.m,v $
0016 % Revision 1.1  2009/08/03 11:10:33  Meidow
0017 % *** empty log message ***
0018 %
0019 % adapted
0020 % Wolfgang Förstner 1/2011
0021 % wfoerstn@uni-bonn.de
0022 
0023 function sugr_plot_conic_explicit...
0024     (C,center_type,bound_type,linewidth,factor,se)
0025 
0026 warning off
0027 
0028 if nargin < 4
0029     linewidth=2;
0030 end
0031 if nargin < 5
0032     factor=1;
0033 end
0034 if nargin < 6
0035     se = [1,1];
0036 end
0037 
0038 C = C/norm(C+eps);
0039 N0=50;
0040 N = 2*N0;  % number of supporting points
0041 
0042 
0043 Chh = C(1:2,1:2);
0044 ch0 = C(1:2,  3);
0045 c00 = C(  3,  3);
0046 
0047 d = det( Chh +eps*eye(2));
0048 x0 = -inv(Chh+eps*eye(2)) *ch0;            % centre point
0049 plot( x0(1), x0(2), center_type);
0050 
0051 c00q = c00 -ch0'*inv(Chh+eps*eye(2))*ch0;                                  %#ok<MINV>
0052 [Rho,Lambda] = eig( -Chh/ (c00q+eps) );
0053 
0054 lambda = diag(Lambda);
0055 
0056 hold on
0057 switch sign(d)
0058     case +1
0059         %% Ellipse
0060         plot( x0(1), x0(2), center_type);
0061         %  (1) plot ellipse
0062         if lambda(1)<0
0063             lambda = -lambda;
0064         end
0065         
0066         lambda=lambda/factor^2;
0067 
0068         t = linspace(0,2*pi, N);
0069 
0070         x = sqrt(1/lambda(1)) * sin(t);
0071         y = sqrt(1/lambda(2)) * cos(t);
0072         x = bsxfun(@plus, Rho * [x; y], x0);
0073         plot( x(1,:), x(2,:), bound_type,'LineWidth',linewidth);
0074 
0075     case -1
0076         %% Hyperbola
0077         % (2) plot hyperbola
0078         if lambda(2)<0
0079             lambda = flipud(lambda);
0080             Rho = fliplr(Rho);
0081         end
0082                 
0083         lambda(2)=lambda(2)/factor^2;
0084 
0085         %         t = linspace(-asinh(a), +asinh(a), N);
0086         %         % t = linspace(-pi,pi, N);
0087         %
0088         %         x = sqrt(1/-lambda(1)) *sinh(t);
0089         %         y = sqrt(1/+lambda(2))* cosh(t);
0090          
0091         
0092         if se(1) ~= 0 && se(2) ~= 0 || se(1) ~= 1 && se(2) ~= 1     % bounds given
0093             t = linspace(atan(se(1)*sqrt(-lambda(1))),atan(se(2)*sqrt(-lambda(1))),N);
0094         end
0095         if se(1) == 1 && se(2) == 1            % bounds by xlim and ylim
0096             dir = Rho(:,1);
0097             xl  = xlim;
0098             yl  = ylim;
0099             s1  = dir'*[xl(1)-x0(1);yl(1)-x0(2)];
0100             s2  = dir'*[xl(1)-x0(1);yl(2)-x0(2)];
0101             s3  = dir'*[xl(2)-x0(1);yl(2)-x0(2)];
0102             s4  = dir'*[xl(2)-x0(1);yl(1)-x0(2)];
0103             se(1) = min([s1,s2,s3,s4])-1;
0104             se(2) = max([s1,s2,s3,s4])+1;
0105             t = linspace(atan(se(1)*sqrt(-lambda(1))),atan(se(2)*sqrt(-lambda(1))),N);
0106         end
0107         if  se(1) == 0 && se(2) == 0                      % line segment
0108             t = linspace(-pi/4,pi/4,N);
0109         end
0110         x = sqrt(1/-lambda(1)) * tan(t);
0111         y = sqrt(1/+lambda(2)) ./ cos(t);
0112 
0113         x1 = repmat(x0, 1, N) + Rho*[x; +y];
0114         x2 = repmat(x0, 1, N) + Rho*[x; -y];
0115         plot( x1(1,:), x1(2,:), bound_type,'LineWidth',linewidth);
0116         plot( x2(1,:), x2(2,:), bound_type,'LineWidth',linewidth);
0117         plot( [(x1(1,1)+x2(1,1))/2,(x1(1,N)+x2(1,N))/2],...
0118             [(x1(2,1)+x2(2,1))/2,(x1(2,N)+x2(2,N))/2] ,center_type);
0119 
0120 
0121     otherwise
0122         error('degenrated conic')
0123 end