main/src/georegression/metric/Intersection2D_F64.java
public static int intersection( LineGeneral2D_F64 line , EllipseRotated_F64 ellipse ,
									Point2D_F64 intersection0 , Point2D_F64 intersection1 , double EPS ) {

		if( EPS < 0 ) {
			EPS = GrlConstants.EPS;
		}

		// First translate the line so that coordinate origin is the same as the ellipse
		double C = line.C + (line.A*ellipse.center.x + line.B*ellipse.center.y);

		// Now rotate the line
		double cphi = Math.cos(ellipse.phi);
		double sphi = Math.sin(ellipse.phi);
		double A =  line.A*cphi + line.B*sphi;
		double B = -line.A*sphi + line.B*cphi;

		// Now solve for the intersections with the coordinate system centered and aligned to the ellipse
		// There are two different ways to solve for this.  Pick the axis with the largest slope
		// to avoid the pathological case
		double a2 = ellipse.a*ellipse.a;
		double b2 = ellipse.b*ellipse.b;

		double x0,y0;
		double x1,y1;

		int totalIntersections;

		if( Math.abs(A) > Math.abs(B) ) {
			double alpha = -C/A;
			double beta = -B/A;

			double aa = beta*beta/a2 + 1.0/b2;
			double bb = 2.0*alpha*beta/a2;
			double cc = alpha*alpha/a2 - 1.0;

			double inner = bb*bb -4.0*aa*cc;
			if( Math.abs(inner/aa) < EPS ) { // divide by aa for scale invariance
				totalIntersections = 1;
				inner = inner < 0 ? 0 : inner;
			} else if( inner < 0 ) {
				return 0;
			} else {
				totalIntersections = 2;
			}
			double right = Math.sqrt(inner);
			y0 = (-bb + right)/(2.0*aa);
			y1 = (-bb - right)/(2.0*aa);

			x0 =  -(C + B*y0)/A;
			x1 =  -(C + B*y1)/A;
		} else {
			double alpha = -C/B;
			double beta = -A/B;

			double aa = beta*beta/b2 + 1.0/a2;
			double bb = 2.0*alpha*beta/b2;
			double cc = alpha*alpha/b2-1.0;

			double inner = bb*bb -4.0*aa*cc;
			if( Math.abs(inner/aa) < EPS ) { // divide by aa for scale invariance
				totalIntersections = 1;
				inner = inner < 0 ? 0 : inner;
			} else if( inner < 0 ) {
				return 0;
			} else {
				totalIntersections = 2;
			}
			double right = Math.sqrt(inner);
			x0 = (-bb + right)/(2.0*aa);
			x1 = (-bb - right)/(2.0*aa);

			y0 = -(A*x0 + C)/B;
			y1 = -(A*x1 + C)/B;
		}

		// go back into world coordinate system
		intersection0.x = x0*cphi - y0*sphi + ellipse.center.x;
		intersection0.y = x0*sphi + y0*cphi + ellipse.center.y;

		intersection1.x = x1*cphi - y1*sphi + ellipse.center.x;
		intersection1.y = x1*sphi + y1*cphi + ellipse.center.y;

		return totalIntersections;
	}
public static int intersection( LineGeneral2D_F64 line , EllipseRotated_F64 ellipse ,
									Point2D_F64 intersection0 , Point2D_F64 intersection1 , double EPS ) {

		if( EPS < 0 ) {
			EPS = GrlConstants.EPS;
		}

		// First translate the line so that coordinate origin is the same as the ellipse
		double C = line.C + (line.A*ellipse.center.x + line.B*ellipse.center.y);

		// Now rotate the line
		double cphi = Math.cos(ellipse.phi);
		double sphi = Math.sin(ellipse.phi);
		double A =  line.A*cphi + line.B*sphi;
		double B = -line.A*sphi + line.B*cphi;

		// Now solve for the intersections with the coordinate system centered and aligned to the ellipse
		// There are two different ways to solve for this.  Pick the axis with the largest slope
		// to avoid the pathological case
		double a2 = ellipse.a*ellipse.a;
		double b2 = ellipse.b*ellipse.b;

		double x0,y0;
		double x1,y1;

		int totalIntersections;

		if( Math.abs(A) > Math.abs(B) ) {
			double alpha = -C/A;
			double beta = -B/A;

			double aa = beta*beta/a2 + 1.0/b2;
			double bb = 2.0*alpha*beta/a2;
			double cc = alpha*alpha/a2 - 1.0;

			double inner = bb*bb -4.0*aa*cc;
			if( Math.abs(inner/aa) < EPS ) { // divide by aa for scale invariance
				totalIntersections = 1;
				inner = inner < 0 ? 0 : inner;
			} else if( inner < 0 ) {
				return 0;
			} else {
				totalIntersections = 2;
			}
			double right = Math.sqrt(inner);
			y0 = (-bb + right)/(2.0*aa);
			y1 = (-bb - right)/(2.0*aa);

			x0 =  -(C + B*y0)/A;
			x1 =  -(C + B*y1)/A;
		} else {
			double alpha = -C/B;
			double beta = -A/B;

			double aa = beta*beta/b2 + 1.0/a2;
			double bb = 2.0*alpha*beta/b2;
			double cc = alpha*alpha/b2-1.0;

			double inner = bb*bb -4.0*aa*cc;
			if( Math.abs(inner/aa) < EPS ) { // divide by aa for scale invariance
				totalIntersections = 1;
				inner = inner < 0 ? 0 : inner;
			} else if( inner < 0 ) {
				return 0;
			} else {
				totalIntersections = 2;
			}
			double right = Math.sqrt(inner);
			x0 = (-bb + right)/(2.0*aa);
			x1 = (-bb - right)/(2.0*aa);

			y0 = -(A*x0 + C)/B;
			y1 = -(A*x1 + C)/B;
		}

		// go back into world coordinate system
		intersection0.x = x0*cphi - y0*sphi + ellipse.center.x;
		intersection0.y = x0*sphi + y0*cphi + ellipse.center.y;

		intersection1.x = x1*cphi - y1*sphi + ellipse.center.x;
		intersection1.y = x1*sphi + y1*cphi + ellipse.center.y;

		return totalIntersections;
	}
Determines the location(s) that a line and ellipse intersect.  Returns the number of intersections found.
NOTE: Due to floating point errors, it's possible for a single solution to returned as two points.

@param line Line
@param ellipse Ellipse
@param intersection0 Storage for first point of intersection.
@param intersection1 Storage for second point of intersection.
@param EPS Numerical precision.  Set to a negative value to use default
@return Number of intersections.  Possible values are 0, 1, or 2.