function [L1,L2] = GetLinesFromDegenerateConic(C) detmax = abs(C(1,1)*C(2,2)-C(1,2)^2); jlast = 3; det = abs(C(2,2)*C(3,3)-C(2,3)^2); if detmax < det detmax = det; jlast = 1; end det = abs(C(1,1)*C(3,3)-C(1,3)^2); if detmax < det jlast = 2; end switch jlast case 1 a = C(2,2); b = C(3,3); c = C(2,3); p = C(2,1); q = C(3,1); r = C(1,1); case 2 a = C(1,1); b = C(3,3); c = C(1,3); p = C(1,2); q = C(3,2); r = C(2,2); case 3 a = C(1,1); b = C(2,2); c = C(1,2); p = C(1,3); q = C(2,3); r = C(3,3); end disc = (a-b)^2 + 4*c^2; if a+b > 0 d1 = (a + b + sqrt(disc))/2; d2 = (a*b - c^2)/d1; else d1 = (a + b - sqrt(disc))/2; d2 = (a*b - c^2)/d1; end if abs(a-d1)>abs(b-d1) V1x = c/sqrt(c^2 + (d1-a)^2); V1y = (d1-a)/sqrt(c^2 + (d1-a)^2); else V1x = (d1-b)/sqrt(c^2 + (d1-b)^2); V1y = c/sqrt(c^2 + (d1-b)^2); end if abs(a-d2)>abs(b-d2) V2x = c/sqrt(c^2 + (d2-a)^2); V2y = (d2-a)/sqrt(c^2 + (d2-a)^2); else V2x = (d2-b)/sqrt(c^2 + (d2-b)^2); V2y = c/sqrt(c^2 + (d2-b)^2); end if d1 >= 0 a1 = sqrt(d1)*V1x + sqrt(-d2)*V2x; b1 = sqrt(d1)*V1y + sqrt(-d2)*V2y; a2 = sqrt(d1)*V1x - sqrt(-d2)*V2x; b2 = sqrt(d1)*V1y - sqrt(-d2)*V2y; else a1 = sqrt(d2)*V2x + sqrt(-d1)*V1x; b1 = sqrt(d2)*V2y + sqrt(-d1)*V1y; a2 = sqrt(d2)*V2x - sqrt(-d1)*V1x; b2 = sqrt(d2)*V2y - sqrt(-d1)*V1y; end det = a1*b2 - a2*b1; c1 = 2*(a1*q - b1*p)/det; c2 = 2*(b2*p - a2*q)/det; switch jlast case 1 L1=[c1;a1;b1]; L2=[c2;a2;b2]; case 2 L1=[a1;c1;b1]; L2=[a2;c2;b2]; case 3 L1=[a1;b1;c1]; L2=[a2;b2;c2]; end