Ellipse.h

Go to the documentation of this file.

// Jolt Physics Library (https://github.com/jrouwe/JoltPhysics)
// SPDX-FileCopyrightText: 2021 Jorrit Rouwe
// SPDX-License-Identifier: MIT
 
#pragma once
 
#include <Jolt/Math/Float2.h>
 
JPH_NAMESPACE_BEGIN
 
class Ellipse
{
public:
    JPH_OVERRIDE_NEW_DELETE
 
                    Ellipse(float inA, float inB) : mA(inA), mB(inB) { JPH_ASSERT(inA > 0.0f); JPH_ASSERT(inB > 0.0f); }
 
    bool            IsInside(const Float2 &inPoint) const
    {
        return Square(inPoint.x / mA) + Square(inPoint.y / mB) <= 1.0f;
    }
 
    Float2          GetClosestPoint(const Float2 &inPoint) const
    {
        float a_sq = Square(mA);
        float b_sq = Square(mB);
 
        // Equation of ellipse: f(x, y) = (x/a)^2 + (y/b)^2 - 1 = 0                                         [1]
        // Normal on surface: (df/dx, df/dy) = (2 x / a^2, 2 y / b^2)
        // Closest point (x', y') on ellipse to point (x, y): (x', y') + t (x / a^2, y / b^2) = (x, y)
        // <=> (x', y') = (a^2 x / (t + a^2), b^2 y / (t + b^2))
        // Requiring point to be on ellipse (substituting into [1]): g(t) = (a x / (t + a^2))^2 + (b y / (t + b^2))^2 - 1 = 0
 
        // Newton raphson iteration, starting at t = 0
        float t = 0.0f;
        for (;;)
        {
            // Calculate g(t)
            float t_plus_a_sq = t + a_sq;
            float t_plus_b_sq = t + b_sq;
            float gt = Square(mA * inPoint.x / t_plus_a_sq) + Square(mB * inPoint.y / t_plus_b_sq) - 1.0f;
 
            // Check if g(t) it is close enough to zero
            if (abs(gt) < 1.0e-6f)
                return Float2(a_sq * inPoint.x / t_plus_a_sq, b_sq * inPoint.y / t_plus_b_sq);
 
            // Get derivative dg/dt = g'(t) = -2 (b^2 y^2 / (t + b^2)^3 + a^2 x^2 / (t + a^2)^3)
            float gt_accent = -2.0f *
                (a_sq * Square(inPoint.x) / Cubed(t_plus_a_sq)
                + b_sq * Square(inPoint.y) / Cubed(t_plus_b_sq));
 
            // Calculate t for next iteration: tn+1 = tn - g(t) / g'(t)
            float tn = t - gt / gt_accent;
            t = tn;
        }
    }
 
    Float2          GetNormal(const Float2 &inPoint) const
    {
        // Calculated by [d/dx f(x, y), d/dy f(x, y)], where f(x, y) is the ellipse equation from above
        return Float2(inPoint.x / Square(mA), inPoint.y / Square(mB));
    }
 
private:
    float           mA;             
    float           mB;             
};
 
JPH_NAMESPACE_END