#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_GENERATORS_TRIANGLE2GENERATOR_H
#define NUCLEX_GEOMETRY_AREAS_GENERATORS_TRIANGLE2GENERATOR_H

#include "Nuclex/Geometry/Config.h"
#include "Nuclex/Geometry/Areas/Triangle2.h"
#include "Nuclex/Geometry/Areas/Rectangle2.h"
#include "Nuclex/Geometry/Areas/Disc2.h"

namespace Nuclex { namespace Geometry { namespace Areas { namespace Generators {

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

  /// <summary>Point and area generation methods for 2D triangles</summary>
  template <typename TScalar>
  class Triangle2Generator {

    /// <summary>Returns a random point within the triangle</summary>
    /// <param name="triangle">Triangle in which a random point will be generated</param>
    /// <param name="random">Random number generator that will be used</param>
    /// <returns>A random point inside the triangle</returns>
    public: template <typename TRandomNumberEngine>
    static Point2<TScalar> GetRandomPointInside(
      const Triangle2<TScalar> &triangle, TRandomNumberEngine &random
    ) {
      // Treat triangle as one half of a parallelogram. Then calculate a
      // random point along the width and the height.
      //
      //         Width
      //      -----------
      //     a___________
      //     /          /  /
      //    /          /  / Height
      //   /__________/  /
      //   b         c
      //
      // Now split the parallelogram along a-c and mirror any points that
      // end up on the other side.
      //
      //   xx = rand(1.0)
      //   yy = rand(1.0)
      //
      //   if((xx + yy) > 1.0) {
      //     xx = 1.0 - xx;
      //     yy = 1.0 - yy;
      //   }
      //
      // Then calculate the absolute coordinates of the point
      //    
      //     a. 
      //     /  . 
      //    /     .   
      //   /________. 
      //   b         c
      //
      // x = b.x + xx * (c.x - b.x) + (b.x - a.x) * yy
      // y = yy * (a.y - b.y)

      // This might be using the same approach, not sure though :)
      // http://vcg.isti.cnr.it/activities/geometryegraphics/pointintetraedro.html

      // Calculate random barycentric coordinates inside the unit
      // triangle (0,0)-(1,0)-(0,1). First, we take random x and y coordinates in
      // a box with side lengths ab, ac
      TScalar x = Math<TScalar>::Random(random, 1);
      TScalar y = Math<TScalar>::Random(random, 1);

      // The triangle covers exactly half of the box our random points are distributed
      // in. Instead of rejecting these coordinates, we mirror the other half of the box
      // back into the triangle (on the bc edge of the triangle)
      if(x + y > 1) {
        x = 1 - x;
        y = 1 - y;
      }

      // Convert the barycentric coordinates back into cartesian coordinates
      TScalar z = 1 - x - y;
      return (triangle.A * x) + (triangle.B * y) + (triangle.C * z);
    }

    /// <summary>Returns a random point on the outline of the triangle</summary>
    /// <param name="triangle">
    ///   Triangle on whose outline a random point will be generated
    /// </param>
    /// <param name="random">Random number generator that will be used</param>
    /// <returns>A random point on the triangle's outline</returns>
    public: template <typename TRandomNumberEngine>
    static Point2<TScalar> GetRandomPointOnOutline(
      const Triangle2<TScalar> &triangle, TRandomNumberEngine &random
    ) {
      Vector2<TScalar> ab = (triangle.B - triangle.A);
      Vector2<TScalar> bc = (triangle.C - triangle.B);
      Vector2<TScalar> ca = (triangle.A - triangle.C);

      float lengthAB = ab.GetSquaredLength();
      float lengthBC = bc.GetSquaredLength();
      float lengthCA = ca.GetSquaredLength();
      float totalLength = lengthAB + lengthBC + lengthCA;

      float position = Math<TScalar>::Random(random, totalLength);

      // Is the position on the A->B edge?
      if(position < lengthAB) {
        float t = (position / lengthAB);
        return triangle.A + (ab * t);
      }

      // Is the position on the B->C edge?
      float lengthAC = lengthAB + lengthBC;
      if(position < lengthAC) {
        float t = ((position - lengthAB) / lengthBC);
        return triangle.B + (bc * t);
      }

      // Position is on the C->A edge!
      {
        float t = ((position - lengthAC) / lengthCA);
        return triangle.C + (ca * t);
      }
    }

    /// <summary>Determines the closest point to another point inside the triangle</summary>
    /// <param name="triangle">Triangle in which the closest point will be found</param>
    /// <param name="point">Point to which the closest point will be determined</param>
    /// <returns>The closest point to the specified point inside the triangle</returns>
    public: static Point2<TScalar> GetClosestPointInside(
      const Triangle2<TScalar> &triangle, const Point2<TScalar> &point
    ) {
      Vector2<TScalar> offsetAB = triangle.B - triangle.A;
      Vector2<TScalar> offsetAPoint = point - triangle.A;
      TScalar crossAPoint = (offsetAB.X * offsetAPoint.Y) - (offsetAB.Y * offsetAPoint.X);

      Vector2<TScalar> offsetBC = triangle.C - triangle.B;
      Vector2<TScalar> offsetBPoint = point - triangle.B;
      TScalar crossBPoint = (offsetBC.X * offsetBPoint.Y) - (offsetBC.Y * offsetBPoint.X);

      Vector2<TScalar> offsetCA = triangle.A - triangle.C;
      Vector2<TScalar> offsetCPoint = point - triangle.C;
      TScalar crossCPoint = (offsetCA.X * offsetCPoint.Y) - (offsetCA.Y * offsetCPoint.X);

      //      A
      //      O
      //     / \
      // AB /   \ CA
      //   /     \
      //  O-------O
      // B   BC    C

      if(crossCPoint < 0) { // Point lies outside on the CA edge

        // Note: the point can't be outside all three edges, so we ignore this possibility!

        if(crossAPoint < 0) { // Point lies outside the AB and CA edges -> vertex A
          return triangle.A;
        } else if(crossBPoint < 0) { // Point lies outside the BC and CA edges -> vertex C
          return triangle.C;
        } else { // Point lies inside the AB and BC edges -> along edge CA
          TScalar wd = Math<TScalar>::Clamp(
            ((offsetCA.X * offsetCPoint.X) + (offsetCA.Y * offsetCPoint.Y)) /
            ((offsetCA.X * offsetCA.X) + (offsetCA.Y * offsetCA.Y)),
            0, 1
          );
          return Point2<TScalar>(
            triangle.C.X + (triangle.A.X - triangle.C.X) * wd,
            triangle.C.Y + (triangle.A.Y - triangle.C.Y) * wd
          );
        }

      } else { // Point lies inside the CA edge

        if(crossAPoint < 0) {
          if(crossBPoint < 0) { // Point lies outside the AB and BC edges -> vertex B
            return triangle.B;
          } else { // Point lies inside the BC and CA edges -> along edge AB
            TScalar wd = Math<TScalar>::Clamp(
              ((offsetAB.X * offsetAPoint.X) + (offsetAB.Y * offsetAPoint.Y)) /
              ((offsetAB.X * offsetAB.X) + (offsetAB.Y * offsetAB.Y)),
              0, 1
            );
            return Point2<TScalar>(
              triangle.A.X + (triangle.B.X - triangle.A.X) * wd,
              triangle.A.Y + (triangle.B.Y - triangle.A.Y) * wd
            );
          }
        } else {
          if(crossBPoint < 0) { // Point lies inside the AB and CA edges -> along edge BC
            TScalar wd = Math<TScalar>::Clamp(
              ((offsetBC.X * offsetBPoint.X) + (offsetBC.Y * offsetBPoint.Y)) /
              ((offsetBC.X * offsetBC.X) + (offsetBC.Y * offsetBC.Y)),
              0, 1
            );
            return Point2<TScalar>(
              triangle.B.X + (triangle.C.X - triangle.B.X) * wd,
              triangle.B.Y + (triangle.C.Y - triangle.B.Y) * wd
            );
          } else { // Point lies inside all edges -> contained

            // If the point is parallel to one edge and yet inside the triangle,
            // the triangle is degenerate and the point hit it
            if(crossAPoint == 0) {
              std::pair<TScalar, TScalar> rangeX = Math<TScalar>::Extremes(
                triangle.A.X, triangle.B.X, triangle.C.X
              );
              std::pair<TScalar, TScalar> rangeY = Math<TScalar>::Extremes(
                triangle.A.Y, triangle.B.Y, triangle.C.Y
              );
              return Point2<TScalar>(
                Math<TScalar>::Clamp(point.X, rangeX.first, rangeX.second),
                Math<TScalar>::Clamp(point.Y, rangeY.first, rangeY.second)
              );
            } else { // Point is inside a non-degenerate triangle
              return point;
            }

          }
        }
      }
    }

    /// <summary>Determines the closest point to another point on the triangle's outline</summary>
    /// <param name="triangle">Triangle on whose outline the closest point will be found</param>
    /// <param name="point">Point to which the closest point will be determined</param>
    /// <returns>The closest point to the specified point on the triangle's outline</returns>
    public: static Point2<TScalar> GetClosestPointOnOutline(
      const Triangle2<TScalar> &triangle, const Point2<TScalar> &point
    ) {
      Point2<TScalar> closestAB = getClosestPointOnSegment(
        triangle.A, triangle.B, point - triangle.A
      );
      Point2<TScalar> closestBC = getClosestPointOnSegment(
        triangle.B, triangle.C, point - triangle.B
      );
      Point2<TScalar> closestCA = getClosestPointOnSegment(
        triangle.C, triangle.A, point - triangle.C
      );

      TScalar squaredDistanceAB = point.SquaredDistanceTo(closestAB);
      TScalar squaredDistanceBC = point.SquaredDistanceTo(closestBC);
      TScalar squaredDistanceCA = point.SquaredDistanceTo(closestCA);

      if(squaredDistanceAB < squaredDistanceBC) {
        if(squaredDistanceAB < squaredDistanceCA) {
          return closestAB;
        } else {
          return closestCA;
        }
      } else {
        if(squaredDistanceBC < squaredDistanceCA) {
          return closestBC;
        } else {
          return closestCA;
        }
      }
    }

    /// <summary>Calculates the bounding rectangle of a triangle</summary>
    /// <param name="triangle">Triangle of which the bounding rectangle will be calculated</param>
    /// <returns>The triangle's bounding rectangle</returns>
    public: static Rectangle2<TScalar> GetBoundingRectangle(
      const Triangle2<TScalar> &triangle
    ) {
      Rectangle2<TScalar> boundingRectangle;

      {
        std::pair<TScalar, TScalar> extremesX = Math<TScalar>::Extremes(
          triangle.A.X, triangle.B.X, triangle.C.X
        );
        boundingRectangle.Min.X = extremesX.first;
        boundingRectangle.Max.X = extremesX.second;
      }
      {
        std::pair<TScalar, TScalar> extremesY = Math<TScalar>::Extremes(
          triangle.A.Y, triangle.B.Y, triangle.C.Y
        );
        boundingRectangle.Min.Y = extremesY.first;
        boundingRectangle.Max.Y = extremesY.second;
      }

      return boundingRectangle;
    }

    /// <summary>Calculates the bounding disc of a triangle</summary>
    /// <param name="triangle">Triangle of which the bounding disc will be calculated</param>
    /// <returns>The triangle's bounding disc</returns>
    /// <remarks>
    ///   Based on &quot;Minimum bounding circle for a triangle&quot; article by
    ///   Christer Ericson published at http://realtimecollisiondetection.net/blog/?p=20
    /// </remarks>
    public: static Disc2<TScalar> GetBoundingDisc(const Triangle2<TScalar> &triangle) {
      Vector2<TScalar> ab = triangle.B - triangle.A;
      Vector2<TScalar> ac = triangle.C - triangle.A;

      TScalar dotABAB = Vector2<TScalar>::Dot(ab, ab);
      TScalar dotABAC = Vector2<TScalar>::Dot(ab, ac);
      TScalar dotACAC = Vector2<TScalar>::Dot(ac, ac);

      Disc2<TScalar> boundingDisc;

      TScalar d = 2 * (dotABAB * dotACAC - dotABAC * dotABAC);
      if(d == 0) { // Triangle is a flat degenerate, center within bounds

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

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

        boundingDisc.Radius = boundingDisc.Center.DistanceTo(
          extremeX.first, extremeY.first
        );

      } else { // Otherwise use either the circumcenter or the longest edge

        TScalar s = (dotABAB * dotACAC - dotACAC * dotABAC) / d;
        TScalar t = (dotACAC * dotABAB - dotABAB * dotABAC) / d;

        if(s <= 0) {
          boundingDisc.Center.X = (triangle.A.X + triangle.C.X) / 2;
          boundingDisc.Center.Y = (triangle.A.Y + triangle.C.Y) / 2;
          boundingDisc.Radius = boundingDisc.Center.DistanceTo(triangle.A);
        } else if(t <= 0) {
          boundingDisc.Center.X = (triangle.A.X + triangle.B.X) / 2;
          boundingDisc.Center.Y = (triangle.A.Y + triangle.B.Y) / 2;
          boundingDisc.Radius = boundingDisc.Center.DistanceTo(triangle.A);
        } else if(s + t >= 1) {
          boundingDisc.Center.X = (triangle.B.X + triangle.C.X) / 2;
          boundingDisc.Center.Y = (triangle.B.Y + triangle.C.Y) / 2;
          boundingDisc.Radius = boundingDisc.Center.DistanceTo(triangle.B);
        } else {
          boundingDisc.Center = triangle.A + s * ab + t * ac;
          boundingDisc.Radius = boundingDisc.Center.DistanceTo(triangle.A);
        }

      }

      return boundingDisc;
    }

    // Could add a GetEnclosedDisc() method here...

    /// <summary>Finds the closest point to an offset on a line segment</summary>
    /// <param name="start">Point at which the line segment begins</param>
    /// <param name="end">Point at which the line segment ends</param>
    /// <param name="offset">Offset of the line segment from the line's start</param>
    /// <returns>The closest point on the line segment to the offset</returns>
    private: static Point2<TScalar> getClosestPointOnSegment(
      const Point2<TScalar> &start, const Point2<TScalar> &end,
      const Vector2<TScalar> &offset
    ) {
      Vector2<TScalar> direction = (end - start);

      TScalar t = Vector2<TScalar>::Dot(offset, direction);
      if(t < 0) {
        return start;
      } else {
        TScalar squaredSegmentLength = direction.GetSquaredLength();
        if(t > squaredSegmentLength) {
          return end;
        } else {
          return start + (direction * (t / squaredSegmentLength));
        }
      }
    }

  };

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

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

#endif // NUCLEX_GEOMETRY_AREAS_GENERATORS_TRIANGLE2GENERATOR_H