#pragma region CPL License
/*
Nuclex Native Framework
Copyright (C) 2002-2015 Nuclex Development Labs

This library is free software; you can redistribute it and/or
modify it under the terms of the IBM Common Public License as
published by the IBM Corporation; either version 1.0 of the
License, or (at your option) any later version.

This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
IBM Common Public License for more details.

You should have received a copy of the IBM Common Public
License along with this library
*/
#pragma endregion // CPL License

#ifndef NUCLEX_GEOMETRY_AREAS_TRIANGLE2_H
#define NUCLEX_GEOMETRY_AREAS_TRIANGLE2_H

#include "Nuclex/Geometry/Config.h"
#include "Nuclex/Geometry/Point2.h"
#include "Nuclex/Geometry/Vector2.h"

namespace Nuclex { namespace Geometry { namespace Areas {

  // ------------------------------------------------------------------------------------------- //

  /// <summary>Triangle in 2D space constructed from its 3 corner vertices</summary>
  template <typename TScalar>
  struct Triangle2 {
  
    /// <summary>Initializes a triangle using the specified corner points</summary>
    /// <param name="a">First corner of the triangle</param>
    /// <param name="b">Second corner of the triangle</param>
    /// <param name="c">Third corner of the triangle</param>
    public: Triangle2(
      const Point2<TScalar> &a,
      const Point2<TScalar> &b,
      const Point2<TScalar> &c
    ) : A(a), B(b), C(c) {}

    /// <summary>Returns the center point of the triangle</summary>
    /// <returns>The center of the triangle</returns>
    /// <remarks>
    ///   This is also called the centroid - a cut-out triangle could be balanced on this
    ///   center point, as it is the center of gravity.
    /// </remarks>
    public: Point2<TScalar> GetCenter() const {
      return Point2<TScalar>(
        (this->A.X + this->B.X + this->C.X) / 3,
        (this->A.Y + this->B.Y + this->C.Y) / 3
      );
    }

    /// <summary>Returns the circumcenter of the triangle</summarx>
    /// <returns>The triangle's circumcenter</returns>
    /// <remarks>
    ///   The circumcenter is a point inside the triangle to which all three vertices
    ///   have the same distance. Therefore, it also is the center of the smallest
    ///   enclosing circle for triangles with a height of more than half their width.
    /// </remarks>
    public: Point2<TScalar> GetCircumcenter() const {
      TScalar magnitude = (
        this->A.X * (this->C.Y - this->B.Y) +
        this->B.X * (this->A.Y - this->C.Y) +
        this->C.X * (this->B.Y - this->A.Y)
      );
 
      Point2<TScalar> circumcenter;
      if(magnitude == 0) { // Triangle is a flat degenerate, center within bounds

        std::pair<TScalar, TScalar> extremeX = Math<TScalar>::Extremes(
          this->A.X, this->B.X, this->C.X
        );
        circumcenter.X = (extremeX.first + extremeX.second) / 2;

        std::pair<TScalar, TScalar> extremeY = Math<TScalar>::Extremes(
          this->A.Y, this->B.Y, this->C.Y
        );
        circumcenter.Y = (extremeY.first + extremeY.second) / 2;

      } else { // Otherwise calculate the circumcenter

        TScalar a = (this->A.X * this->A.X) + (this->A.Y * this->A.Y);
        TScalar b = (this->B.X * this->B.X) + (this->B.Y * this->B.Y);
        TScalar c = (this->C.X * this->C.X) + (this->C.Y * this->C.Y);
        magnitude *= 2;
      
        circumcenter.X = (
          a * (this->C.Y - this->B.Y) +
          b * (this->A.Y - this->C.Y) +
          c * (this->B.Y - this->A.Y)
        ) / magnitude;

        circumcenter.Y = -(
          a * (this->C.X - this->B.X) +
          b * (this->A.X - this->C.X) +
          c * (this->B.X - this->A.X)
        ) / magnitude;

      }

      return circumcenter;
    }

    /// <summary>Returns the incenter of the triangle</summary>
    /// <returns>The triangle's incenter</returns>
    /// <remarks>
    ///   The incenter is a point inside the triangle to which all three edges have
    ///   the same distance. Therefore, it also serves as the center of the largest
    ///   contained circle.
    /// </remarks>
    public: Point2<TScalar> GetIncenter() const {
      TScalar distanceAB = Point2<TScalar>::DistanceBetween(this->A, this->B);
      TScalar distanceBC = Point2<TScalar>::DistanceBetween(this->B, this->C);
      TScalar distanceCA = Point2<TScalar>::DistanceBetween(this->C, this->A);

      TScalar perimeter = distanceAB + distanceBC + distanceCA;
      if(perimeter == 0) { // Are all of the triangle's points are in the same place?
        return this->A; 
      } else { // Triangle area > 0
        TScalar x = (distanceBC * this->A.X + distanceCA * this->B.X + distanceAB * this->C.X);
        TScalar y = (distanceBC * this->A.Y + distanceCA * this->B.Y + distanceAB * this->C.Y);
        return Point2<TScalar>(x / perimeter, y / perimeter);
      }
    }

    /// <summary>Calculates the area of the triangle</summary>
    /// <returns>The area of the triangle</returns>
    ///   <para>
    ///     Heron's triangle area formular states that, given s = (a + b + c) / 2
    ///     (semiperimeter), the area of a triangle can be calculated as
    ///   </para>
    ///   <code>
    ///              _________________________________
    ///     area = \/ s * (s - a) * (s - b) * (s - c)   
    ///
    ///   </code>
    ///   <para>
    ///     In a paper by W. Kahan this method is proven to be numerically unstable
    ///     for floating point numbers. He recommends to use the following formula
    ///     instead, where the lengths a, b and c have to be sorted in ascending order.
    ///   </para>
    ///   <code>
    ///                     ______________________________________________________________
    ///     area = 0.25 * \/ (a + (b + c)) * (c - (a - b)) * (c + (a - b)) * (a + b - c))
    ///
    ///   </code>
    /// </remarks>
    public: TScalar GetArea() const {
      TScalar a = Point2<TScalar>::DistanceBetween(A, B);
      TScalar b = Point2<TScalar>::DistanceBetween(B, C);
      TScalar c = Point2<TScalar>::DistanceBetween(C, A);

      TScalar s = (a + b + c) / 2;
      return Math<TScalar>::SquareRoot(s * (s - a) * (s - b) * (s - c));
    }

    /// <summary>Calculates the perimeter of the triangle</summary>
    /// <returns>The triangle's perimeter</returns>
    public: TScalar GetPerimeter() const {
      return (
        Point2<TScalar>::DistanceBetween(A, B) +
        Point2<TScalar>::DistanceBetween(B, C) +
        Point2<TScalar>::DistanceBetween(C, A)
      );
    }

    /// <summary>Offsets the entire triangle by the specified amount</summary>
    /// <param name="amountX">Amount the triangle will be offset on the X axis</param>
    /// <param name="amountY">Amount the triangle will be offset on the Y axis</param>
    public: void ShiftBy(TScalar amountX, TScalar amountY) {
      this->A.X += amountX;
      this->B.X += amountX;
      this->C.X += amountX;
      this->A.Y += amountY;
      this->B.Y += amountY;
      this->C.Y += amountY;
    }

    /// <summary>Offsets the entire triangle by the specified amount</summary>
    /// <param name="amount">Amount the triangle will be offset</param>
    public: void ShiftBy(const Vector2<TScalar> &amount) {
      this->A += amount;
      this->B += amount;
      this->C += amount;
    }

    /// <summary>Checks whether the triangle is winding in clockwise order</summary>
    /// <returns>True if the triangle is winding in clockwise order</returns>
    public: bool IsWindingClockwise() const {
      float determinant = (this->B.X - this->A.X) * (this->C.Y - this->A.Y);
      determinant -= (this->B.Y - this->A.Y) * (this->C.X - this->A.X);
      return (determinant >= 0);
    }

    /// <summary>Checks whether the triangle is winding in counter-clockwise order</summary>
    /// <returns>True if the triangle is winding in counter-clockwise order</returns>
    public: bool IsWindingCounterClockwise() const {
      float determinant = (this->B.X - this->A.X) * (this->C.Y - this->A.Y);
      determinant -= (this->B.Y - this->A.Y) * (this->C.X - this->A.X);
      return (determinant < 0);
    }

    /// <summary>Fails an assertion if the triangle's points are not counter-clockwise</summary>
    public: void EnforceCounterClockwiseWinding() const {
      using namespace std;
      assert(
        IsWindingCounterClockwise() && "Triangle2 points have to be in counter-clockwise order"
      );
    }

    /// <summary>Fails an assertion if the triangle's points are not clockwise</summary>
    public: void EnforceClockwiseWinding() const {
      using namespace std;
      assert(
        IsWindingClockwise() && "Triangle2 points have to be in clockwise order"
      );
    }

    /// <summary>Coordinates of the triangle's first corner</summary>
    public: Point2<TScalar> A;
    /// <summary>Coordinates of the triangle's second corner</summary>
    public: Point2<TScalar> B;
    /// <summary>Coordinates of the triangle's third corner</summary>
    public: Point2<TScalar> C;

  };

  // ------------------------------------------------------------------------------------------- //

}}} // namespace Nuclex::Geometry::Areas

#endif // NUCLEX_GEOMETRY_AREAS_TRIANGLE2_H