template<class T>
struct ConicProjectionInputs
{
    Vector3<T>  E, Ce, Ne, Ue, Ve;
    T           A, B;
    Vector3<T>  Cp, Np, Up, Vp;
    T           TestTheta; 

    // ...
}

template<class T>
struct ConicProjectionOutputs
{
    ConicProjection<T> projection;
    ProjectionType projectionType;
    T ThetaLocation;
};

/*
 * Implementation of:
 * https://www.geometrictools.com/Documentation/PerspectiveProjectionEllipse.pdf
 * (Originally published algorithm - the reference page has undergone a
 * significant rewrite since original publication.)
 */
template<class T>
void ProjectConicToPlane(
    const ConicProjectionInputs<T>& inputs,
    ConicProjectionOutputs<T>& outputs
)
{
  // ...
}

// 1b
Matrix3x3<T> Re33;
Re33.SetCol(0, Ue);
Re33.SetCol(1, Ve);
Re33.SetCol(2, Ne);
Matrix3x3<T> Re33_T = Transpose(Re33);

Matrix<3, 2, T> Jp32;
Jp32.SetCol(0, Up);
Jp32.SetCol(1, Vp);
Matrix<2, 3, T> Jp32_T = Transpose(Jp32);

// 2a
Vector3<T> E_bar;
Vector3<T> Ce_bar;
T e0, e1, e2;

E_bar = Re33_T * (E - Ce);
Ce_bar = Re33_T * Ce;
e0 = E_bar[0]; e1 = E_bar[1]; e2 = E_bar[2];

// 2b
Matrix3x3<T> A_bar33;
Vector3<T> B_bar({ (T)0., (T)0., (T)2. / e2 });
T c_bar = -1.;

T a00 = (T)1. / sqr(d0);
T a11 = (T)1. / sqr(d1);
T a01 = (T)0.;
T a02 = -e0 / e2 * a00;
T a12 = -e1 / e2 * a11;
T a22 = (sqr(e0) * a00 + sqr(e1) * a11 - (T)1) / sqr(e2);

A_bar33.SetCol(0, Vector3<T>{a00, a01, a02});
A_bar33.SetCol(1, Vector3<T>{a01, a11, a12});
A_bar33.SetCol(2, Vector3<T>{a02, a12, a22});

// 2c
Matrix3x3<T> A33;
Vector3<T> B;
T c;

A33 = Re33 * A_bar33 * Re33_T;
B = Re33 * (B_bar - (T)2. * (A_bar33 * Ce_bar));
c = Dot(Ce_bar, A_bar33 * Ce_bar) - Dot(B_bar, Ce_bar) + c_bar;

// 2d
Matrix<2, 2, T> A_hat22;
Vector<2, T> B_hat;
T c_hat;

A_hat22 = Jp32_T * A33 * Jp32;
B_hat = Jp32_T * (B + (T)2. * (A33 * Cp));
c_hat = Dot(Cp, A33 * Cp) + Dot(B, Cp) + c;

// 3b
Matrix<2, 2, T> A_hat22_I = Inverse(A_hat22);
Vector<2, T> k = -(A_hat22_I * ((T)0.5 * B_hat));

T d = (T)0.25 * Dot(B_hat, A_hat22_I * B_hat) - c_hat;
Matrix<2, 2, T> M = ((T)1. / d) * A_hat22;

// 3c
// Factor M2 = R*D*R^T.
SymmetricEigensolver2x2<T> eigensolver;
std::array<T, 2> S;
std::array<std::array<T, 2>, 2> evec;
eigensolver(M(0, 0), M(0, 1), M(1, 1), -1, S, evec);

T a1 = 1. / sqrt(S[0]);

bool isHyperbolic = (S[1] < 0);
T a2 = 1. / sqrt(abs(S[1]));

const auto cosT = cos(TestTheta);
const auto sinT = sin(TestTheta);
auto point = Ce + d0 * Ue * cosT + d1 * Ve * sinT;
auto tangent = -d0 * Ue * sinT + d1 * Ve * cosT;

auto toEyePoint = E - point;
Normalize(toEyePoint);

bool forwardOfCamera = Dot(toEyePoint, Np) > 0;

Line3<T> line(point, toEyePoint);
Plane3<T> plane(Np, Cp);

FIQuery<T, Line3<T>, Plane3<T>> intersectionSolver;
auto result = intersectionSolver(line, plane);

auto pointOnProjectionPlane = result.point;
pointOnProjectionPlane -= Cp;

auto planeSpaceTangentCoords = Jp32_T * tangent;
auto planeSpaceCoords = Jp32_T * pointOnProjectionPlane;
planeSpaceCoords -= k;

Matrix<2, 2, T> R_T;
R_T.SetRow(0, evec[0]);
R_T.SetRow(1, evec[1]);

// And, now we know the location of this point in the conic's coord system!
auto conicSpaceCoords = R_T * planeSpaceCoords;

if ((conicSpaceCoords[0] > 0) ^ (!forwardOfCamera))
{
    projectionType = ProjectionType::PositiveHyperbola;
}
else
{
    projectionType = ProjectionType::NegativeHyperbola;
}

// Ensure the hyperbola is defined such that Cross(v1,v2) > 0
auto conicSpaceTangent = R_T * planeSpaceTangentCoords;
bool velocityPointsUpwards = conicSpaceTangent[1] > 0;
if ((projectionType == ProjectionType::PositiveHyperbola) != velocityPointsUpwards)
{
    a2 *= -1;
}

// Avoid dividing by zero, ... which realistically isn't going to happen
// because then we're not hyperbolic... Which, is why there's a singularity here
// that blows up the math anyways.
if (forwardOfCamera)
{
    ThetaLocation = atanh(conicSpaceCoords[1] * a1 / conicSpaceCoords[0] / a2);
}
else
{
    // The current object is off the projected hyperbola.
    // We can at least find out which side it's closer to...
    ThetaLocation = -std::numeric_limits<T>::infinity();
    if (atanh(conicSpaceCoords[1] * a1 / conicSpaceCoords[0] / a2) > 0)
    {
        ThetaLocation *= -1;
    }
}

if (forwardOfCamera)
{
    projectionType = ProjectionType::Ellipse;
}
else
{
    projectionType = ProjectionType::NotVisible;
}

// Ensure the ellipse is defined such that Cross(v1,v2) > 0
if (planeSpaceCoords[0] * planeSpaceTangentCoords[1] - planeSpaceCoords[1] * planeSpaceTangentCoords[0] < 0)
{
    a2 *= -1.;

    ThetaLocation = atan2(-conicSpaceCoords[1] * a1, -conicSpaceCoords[0] * a2);
}
else
{
    ThetaLocation = atan2(+conicSpaceCoords[1] * a1, +conicSpaceCoords[0] * a2);
}
if (ThetaLocation < 0) ThetaLocation += twopi<T>;

// SegmentList - X < Y for non-wrap around segments
//               X > Y if the segment wraps around 0   
template<class T>
bool ClipEllipseToFrustum(
    const FrustumParameters<T>& frustum,
    const ConicProjection<T>& projection,
    vector<ConicSegment<T>>& SegmentList
)
{
  // ...
}

struct ClipChange
{
    int     Whom;
    bool    HiddenChange;
    T       When;
    int     HiddenState;
    ClipChange(int whom, bool hidden, T when)
    {
        Whom = whom;
        HiddenChange = hidden;
        When = normalizeRadians0toTwoPi(when);
        HiddenState = 0;
    }
};

vector<ClipChange> clipPoints;
SegmentList.clear();

T maxX = z * tan(fov / 2);
T maxY = maxX / aspectRatio;
T a, b, c;
T theta1, dt1;
T theta2, dt2;
bool interceptFound;

// Right
a = a1 * v1[0]; b = a2 * v2[0]; c = k[0];
if (ClipEllipse(a, b, c, (T)+maxX, interceptFound, theta1, dt1, theta2, dt2))
{
    if (interceptFound)
    {
        //        UE_LOG(LogTemp, Warning, TEXT("X+ Clip at th = %f (%f)/ %f (%f)"), theta1, dt1, theta2, dt2);
        bool hidden = dt1 >= 0;
        clipPoints.push_back(ClipChange(0, hidden, theta1));
        clipPoints.push_back(ClipChange(0, !hidden, theta2));
    }
}
else
{
    // Outside of frustum... early out!
    return false;
}

// Bottom
a = a1 * v1[1]; b = a2 * v2[1]; c = k[1];
if (ClipEllipse(a, b, c, (T)-maxY, interceptFound, theta2, dt2, theta1, dt1))
{
    if (interceptFound)
    {
        bool hidden = dt2 >= 0;

        clipPoints.push_back(ClipChange(1, hidden, theta1));
        clipPoints.push_back(ClipChange(1, !hidden, theta2));
    }
}
else
{
    // Outside of frustum... early out!
    return false;
}

// Left
a = a1 * v1[0]; b = a2 * v2[0]; c = k[0];
if (ClipEllipse(a, b, c, (T)-maxX, interceptFound, theta2, dt2, theta1, dt1))
{
    if (interceptFound)
    {
        bool hidden = dt2 >= 0;

        clipPoints.push_back(ClipChange(2, hidden, theta1));
        clipPoints.push_back(ClipChange(2, !hidden, theta2));
    }
}
else
{
    // Outside of frustum... early out!
    return false;
}

// Top
a = a1 * v1[1]; b = a2 * v2[1]; c = k[1];
if (ClipEllipse(a, b, c, (T)+maxY, interceptFound, theta1, dt1, theta2, dt2))
{
    if (interceptFound)
    {
        bool hidden = dt1 >= 0;

        clipPoints.push_back(ClipChange(3, hidden, theta1));
        clipPoints.push_back(ClipChange(3, !hidden, theta2));
    }
}
else
{
    // Outside of frustum... early out!
    return false;
}

// If we survived here, the ellipse is visible in at least one place...
// No clip points means it never crosses a frustum line/plane, so it must be fully visible
bool allVisible = clipPoints.empty();

if (allVisible)
{
    return true;
}

// The ellipse is hidden in at least one place...
// It will hidden anytime it's crossed 1 or more clip line.
// So, for each clip line noted, remove all states (via removing associated state
// changes) where that clip segment becomes hidden.

// We're effectively ref counting the segment 'hides'.
// Segments are visible anywhere there's zero hidden segments on the ref count.

// Start, by sorting the state changes by time...
std::sort(
    clipPoints.begin(), clipPoints.end(),
    [](const ClipChange& a, const ClipChange& b)
    {
        return a.When < b.When;
    });

// Loop twice through the list to wrap around the states...
int state = 0;
for (int _rawindex = 0; _rawindex < 2 * clipPoints.size(); ++_rawindex)
{
    int index = _rawindex % clipPoints.size();

    // We can't early out if the state is non-zero, because there may be
    // outstanding ref's to tally...

    // If this is the start of a segment, clear the hidden flag
    bool segmentIsHidden = clipPoints[index].HiddenChange;
    int whom = clipPoints[index].Whom;
    int bitFlag = 1 << whom;

    if (segmentIsHidden)
    {
        state |= bitFlag;
    }
    else
    {
        state &= (-1) ^ bitFlag;
    }

    clipPoints[index].HiddenState = state;
}

// loop through twice again
// It would be moar elegant, logically, if these loops were combined.
// ... but its still only O(n) and this isn't an inner loop.
bool wasVisible = false;
T segmentStart = 0;
for (int _rawindex = 0; _rawindex < 2 * clipPoints.size(); ++_rawindex)
{
    int index = _rawindex % clipPoints.size();
    int hiddenState = clipPoints[index].HiddenState;

    // We can early out if this one's already been done...
    if (hiddenState == -1) break;

    bool visible = hiddenState == 0;
    if (!wasVisible && visible)
    {
        segmentStart = clipPoints[index].When;
    }
    else if (wasVisible && !visible)
    {
        T segmentEnd = clipPoints[index].When;
        T segmentLength = segmentEnd - segmentStart;
        if (segmentLength < 0) segmentLength += twopi<T>;
        if (segmentStart < 0) segmentStart += twopi<T>;
        SegmentList.push_back(ConicSegment<T>(normalizeRadians0toTwoPi<T>(segmentStart), segmentLength));
    }

    // Mark the state so we don't reconsider it... If, it was segmentized
    if (visible)
    {
        clipPoints[index].HiddenState = -1;
    }

    wasVisible = visible;
}

/*
*   Clip a screen-space projected ellipse to the viewable area
*/

// Intercept:
//  a*cos(t) + b*sin(t) + c == d
template<class T>
bool ClipEllipse(T a, T b, T c, T d, bool& interceptsFustumLine, T& theta1, T& dt1, T& theta2, T& dt2)
{
    // https://www.mathcentre.ac.uk/resources/uploaded/mc-ty-rcostheta-alpha-2009-1.pdf
    T R = (T)(sqrt(sqr(a) + sqr(b)));
    T act = (d - c) / R;

    // Is there an intercept?  Welp, can we take the arcos of act?
    if (act <= (T)1. && act >= (T)-1.)
    {
        // Yes, there's an intercept!  Two, in fact!

        // Compute the phase shift
        T alpha = (T)atan2((double)b, (double)a);

        T arcos = (T)acos(act);
        theta1 = alpha + arcos;
        theta2 = alpha - arcos;

        dt1 = -R * sin(arcos);
        dt2 = +R * sin(arcos);
        interceptsFustumLine = true;
        return true;
    }
    else
    {
        interceptsFustumLine = false;
        return d * act > 0;
    }

    return false;
}

template<class T>
void ClipHyperbolaToFrustum(
    const FrustumParameters<T>& frustum,
    const ConicProjection<T>& projection,
    bool positiveOrientation,
    vector<ConicSegment<T>>& SegmentList
)
{
    // ...
}

T maxX = z * tan(fov / 2);
T maxY = maxX / aspectRatio;
const Vector<2, T> screenDimentions({ maxX, maxY });

// { {direction}, {center} }
static const pair<Vector<2, T>, Vector<2, T>> clipLines[4] =
{
    // Top Line
    {{ +1,  0 }, { 0, +1 }},
    // Right Light
    {{  0, -1 }, { +1, 0 }},
    // Bottom
    {{ -1,  0 }, { 0, -1 }},
    // Left
    {{  0, +1 }, { -1, 0 }}
};

// rotate line's coordinate system coordinate system so:
// y ->  v1
// x -> -v2
// (we do this for convenience in intercepting the hyperbolic equation, in its standard form)
Matrix<2, 2, T> R;
R.SetRow(0, v1);
R.SetRow(1, v2);

SegmentList.clear();

// bit field of hidden/showing states....
int states = 0;
vector<pair <int, T>> stateChangeList;

const bool hidden = true;
const bool showing = false;

for (int i = 0; i < 4; ++i)
{
    const pair<Vector<2, T>, Vector<2, T>>& line = clipLines[i];

    bool initialState;
    T stateChanges[2];
    int n;

    ClipHyperbolaToLine(line, screenDimentions, R, positiveOrientation, k, a, b, initialState, stateChanges, n);

    if (initialState == hidden && n == 0)
    {
        // The hyperbola is never on the visible side of this line...
        /// We can bail right now...
        return;
    }
    
    int stateFlag = initialState ? hidden : showing;

    // Remember if this line initially was hidden or showing
    states |= stateFlag << i;

    // Add all state changes...
    if (n > 0)
    {
        stateChangeList.push_back({ i, stateChanges[0] });
        if (n > 1)
        {
            stateChangeList.push_back({ i, stateChanges[1] });
        }
    }
}

if (stateChangeList.size() > 0)
{
    std::sort(
        stateChangeList.begin(), stateChangeList.end(),
        [](const  pair <int, T>& a, const  pair <int, T>& b)
        {
            return a.second < b.second;
        });

    // Ignore any cases where the line is visible to infinity (initially showing)
    int lastState = states;
    T segmentBegin {};

    // Loop through the state changes looking for transitions from
    // one or more hidden clip line bisections.
    for (int i = 0; i < stateChangeList.size(); ++i)
    {
        const pair <int, T>& change = stateChangeList[i];

        int who = change.first;
        T when = change.second;

        int bitMask = (1 << who);

        // Flip the state...
        states ^= bitMask;

        bool visible = !states;

        if (!states && lastState)
        {
            segmentBegin = when;
        }
        else if (states && !lastState)
        {
            SegmentList.push_back({ segmentBegin, when-segmentBegin });
        }

        lastState = states;
    }
}

// Hyperbola:
//      sqr(x)/sqr(a) - sqr(y)/sqr(b) = 1
// 
// Line:
//      A*x + B*y = C
//
// find the intersection... And discard cases where X < 0
template<class T>
int BisectHyperbola(T A, T B, T C, bool positiveOrientation, T a, T b, bool& initiallyHidden, T(&points)[2])
{
    double sign = positiveOrientation ? +1. : -1;

    // The intersection is solved by a simple quadratic equation....
    // (Easily derived by substititing the line equation into the hyperbola eq)
    T A_quad = sqr(B) - sqr(a * A / b);
    T B_quad = 2 * B * C;
    T C_quad = sqr(C) - sqr(a * A);

    T d = sqr(B_quad) - 4. * A_quad * C_quad;

    int n = 0;

    const Vector<2, T> lineNormal({ A, B });

    // First Asymptote
    const Vector<2, T> asymptote({ sign * a, -sign * b });

    // Does the first asymptote cross INTO the hidden side?  If so, it starts off hidden
    // By definition, the lines have a positive cross product in this case...
    // Which means their normal vector points towards the hidden side...
    // Which means the dot product of it and another line will be positive if it's heading off to the hidden side
    initiallyHidden = Dot(lineNormal, asymptote) > 0;

    // if the discriminant is > 0, there's two real roots
    if (d > 0)
    {
        // But only one of them is against the x>0 half of the hyperbola we care about...
        T x, y;

        // (Wash)
        // Solve the quatratic, using the '-' discriminant
        T q = (-B_quad - sqrt(d)) / (2. * A_quad);

        // what does the quadratic yield?  the point intersection point (x(y), y).
        y = q;
        // knowing y, we can resolve z...
        x = (-B * y - C) / A;

        // (Rinse)
        // We only care about the hyperbola in the same side as (x,y)...
        // Any -sign*x < 0 intersections are on the mirrored hyperbola, which is of no interest
        // to us in this use case.
        if (sign * x >= 0)
            // Intersetction at: {x, y}
            points[n++] = atanh(y / x * a / b);

        // (Repeat.)
        // using the '+' discriminant
        q = (-B_quad + sqrt(d)) / (2. * A_quad);
        y = q;
        x = (-B * y - C) / A;

        if (sign * x >= 0)
            points[n++] = atanh(y / x * a / b);
    }
    // vvv -------------------
    // Included for clarity.
    // (The compiler is smart enough to optimize the branches away, ... Except for
    // unoptimized builds, which means you may want to set a breakpoint here anyways
    // since you're likely debugging.)
    else if (d == 0)
    {
        // This case only happens if the intersection is a tangent...
        // In the case of a tanget, a hyperbola never crosses the line...
        // So, we don't care about this case.... as the hyperbola will never
        // need clipped here.
    }
    else if (d < 0)
    {
        // There are no real roots, to the quadratic, which means the line does
        // not intersect the hyperbola at all.
    }
    // ^^^ -------------------

    return n;
}

template<class T>
void ClipHyperbolaToLine(
    const pair<Vector<2, T>, Vector<2, T>>& init,
    const Vector<2, T>& screenDimentions,
    const Matrix<2, 2, T>& R,
    bool positiveOrientation,
    const Vector<2, T>& k,
    T a, T b,
    bool& initialState,
    T (&stateChanges)[2],
    int& nStateChangeCount
)
{
    // Transform the line from screen space into the hyperbola's space...
    const Vector<2, T >& direction = init.first;
    const Vector<2, T >& center = init.second * screenDimentions;

    const Vector<2, T> direction_prime = R * direction;
    const Vector<2, T> center_prime = R * (center - k);

    // Get linear coefficients of the line in the form:
    // Ax + By + C = 0
    double A, B, C;
    A = -direction_prime[1];
    B = direction_prime[0];
    C = -A * center_prime[0] - B * center_prime[1];

    // Now the line is position in they hyperbola's space, so the intersection can be solved
    // via an intersection against a "standard" form hyperbola: x^2/a^2 - y^2/b^2 == 1
    nStateChangeCount = BisectHyperbola(A, B, C, positiveOrientation, a, b, initialState, stateChanges);
}

float thetaB = thetaA + segmentLength;
FVector2D c1 = Conic(thetaB);

FVector B = Center + c1.X * Axis1 + c1.Y * Axis2;

FVector AtoB = B - A;

const FVector TangentX = AtoB.GetSafeNormal();
const FVector TangentZ = CameraDirection.GetSafeNormal();
const FVector TangentY = (TangentX ^ TangentZ).GetSafeNormal();

Vertex.SetTangents((FVector3f)TangentX, (FVector3f)TangentY, (FVector3f)TangentZ);

FVector Perp = TangentY;
Perp *= LineThickness * 0.5f;

const FVector WorldProjectionPlanePositionA1 = A + Perp;
const FVector WorldProjectionPlanePositionA2 = A - Perp;

// 'Conic' defines the 2D shape of our conic section.
// The possibilities of 'Circle' and 'Parabola' are treated as elliptical and hyperbolic.
FVector2D(*Conic)(float theta) = nullptr;

FVector2D UConicRendererComponent::Ellipse(float theta)
{
    return FVector2D(cos(theta), sin(theta));
}

FVector2D UConicRendererComponent::PositiveHyperbola(float theta)
{
    return FVector2D(cosh(theta), sinh(theta));
}

FVector2D UConicRendererComponent::NegativeHyperbola(float theta)
{
    return FVector2D(-cosh(theta), -sinh(theta));
}