package de.lmu.ifi.dbs.elki.math;
/*
This file is part of ELKI:
Environment for Developing KDD-Applications Supported by Index-Structures

Copyright (C) 2011
Ludwig-Maximilians-Universität München
Lehr- und Forschungseinheit für Datenbanksysteme
ELKI Development Team

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU Affero General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program 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
GNU Affero General Public License for more details.

You should have received a copy of the GNU Affero General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

import java.util.Collections;
import java.util.Comparator;
import java.util.Iterator;
import java.util.LinkedList;
import java.util.List;
import java.util.Stack;

import de.lmu.ifi.dbs.elki.data.spatial.Polygon;
import de.lmu.ifi.dbs.elki.math.linearalgebra.Vector;
import de.lmu.ifi.dbs.elki.utilities.documentation.Reference;

/**
 * Classes to compute the convex hull of a set of points in 2D, using the
 * classic Grahams scan. Also computes a bounding box.
 * 
 * @author Erich Schubert
 */
@Reference(authors = "Paul Graham", title = "An Efficient Algorithm for Determining the Convex Hull of a Finite Planar Set", booktitle = "Information Processing Letters 1")
public class ConvexHull2D {
  /**
   * The current set of points
   */
  private List<Vector> points;

  /**
   * Min/Max in X
   */
  private DoubleMinMax minmaxX = new DoubleMinMax();

  /**
   * Min/Max in Y
   */
  private DoubleMinMax minmaxY = new DoubleMinMax();

  /**
   * Flag to indicate that the hull has been computed.
   */
  private boolean ok = false;

  /**
   * Scaling factor if we have very small polygons.
   * 
   * TODO: needed? Does this actually improve things?
   */
  private double factor = 1.0;

  /**
   * Constructor.
   */
  public ConvexHull2D() {
    this.points = new LinkedList<Vector>();
  }

  /**
   * Add a single point to the list (this does not compute the hull!)
   * 
   * @param point Point to add
   */
  public void add(Vector point) {
    if (this.ok) {
      this.points = new LinkedList<Vector>(this.points);
      this.ok = false;
    }
    this.points.add(point);
    // Update data set extends
    minmaxX.put(point.get(0));
    minmaxY.put(point.get(1));
  }

  /**
   * Compute the convex hull.
   */
  private void computeConvexHull() {
    // Trivial cases
    if(points.size() < 3) {
      this.ok = true;
      return;
    }
    // Avoid numerical instabilities by rescaling
    double maxX = Math.max(Math.abs(minmaxX.getMin()), Math.abs(minmaxX.getMax()));
    double maxY = Math.max(Math.abs(minmaxY.getMin()), Math.abs(minmaxY.getMax()));
    if(maxX < 10.0 || maxY < 10.0) {
      factor = 10 / maxX;
      if(10 / maxY > factor) {
        factor = 10 / maxY;
      }
    }
    // Find the new origin point
    findStartingPoint();
    // Sort points for the scan
    final Vector origin = points.get(0);
    Collections.sort(this.points, new Comparator<Vector>() {
      @Override
      public int compare(Vector o1, Vector o2) {
        return isLeft(o1, o2, origin) ? +1 : -1;
      }
    });

    grahamScan();
    this.ok = true;
  }

  /**
   * Find the starting point, and sort it to the beginning of the list. The
   * starting point must be on the outer hull. Any "most extreme" point will do, e.g.
   * the one with the lowest Y coordinate and for ties with the lowest X.
   */
  private void findStartingPoint() {
    // Well, we already know the best Y value...
    final double bestY = minmaxY.getMin();
    double bestX = Double.POSITIVE_INFINITY;
    int bestI = -1;
    Iterator<Vector> iter = this.points.iterator();
    for(int i = 0; iter.hasNext(); i++) {
      Vector vec = iter.next();
      if(vec.get(1) == bestY && vec.get(0) < bestX) {
        bestX = vec.get(0);
        bestI = i;
      }
    }
    assert (bestI >= 0);
    // Bring the reference object to the head.
    if(bestI > 0) {
      points.add(0, points.remove(bestI));
    }
  }

  /**
   * Get the relative X coordinate to the origin.
   * 
   * @param a
   * @param origin origin vector
   * @return relative X coordinate
   */
  private final double getRX(Vector a, Vector origin) {
    return (a.get(0) - origin.get(0)) * factor;
  }

  /**
   * Get the relative Y coordinate to the origin.
   * 
   * @param a
   * @param origin origin vector
   * @return relative Y coordinate
   */
  private final double getRY(Vector a, Vector origin) {
    return (a.get(1) - origin.get(1)) * factor;
  }

  /**
   * Test whether a point is left of the other wrt. the origin.
   * 
   * @param a Vector A
   * @param b Vector B
   * @param o Origin vector
   * @return true when left
   */
  protected final boolean isLeft(Vector a, Vector b, Vector o) {
    final double cross = getRX(a, o) * getRY(b, o) - getRY(a, o) * getRX(b, o);
    if(cross == 0) {
      // Compare manhattan distances - same angle!
      final double dista = Math.abs(getRX(a, o)) + Math.abs(getRY(a, o));
      final double distb = Math.abs(getRX(b, o)) + Math.abs(getRY(b, o));
      return dista > distb;
    }
    return cross > 0;
  }

  /**
   * Manhattan distance.
   * 
   * @param a Vector A
   * @param b Vector B
   * @return Manhattan distance
   */
  private double mdist(Vector a, Vector b) {
    return Math.abs(a.get(0) * factor - b.get(0) * factor) + Math.abs(a.get(1) * factor - b.get(1) * factor);
  }

  /**
   * Simple convexity test.
   * 
   * @param a Vector A
   * @param b Vector B
   * @param c Vector C
   * @return convexity
   */
  private final boolean isConvex(Vector a, Vector b, Vector c) {
    // We're using factor to improve numerical contrast for small polygons.
    double area = (b.get(0) * factor - a.get(0) * factor) * (c.get(1) * factor - a.get(1) * factor) - (c.get(0) * factor - a.get(0) * factor) * (b.get(1) * factor - a.get(1) * factor);
    if(area == 0) {
      return (mdist(b, c) >= mdist(a, b) + mdist(a, c));
    }
    return (area < 0);
  }

  /**
   * The actual graham scan main loop.
   */
  private void grahamScan() {
    if(points.size() < 3) {
      return;
    }
    Iterator<Vector> iter = points.iterator();
    Stack<Vector> stack = new Stack<Vector>();
    // Start with the first two points on the stack
    stack.add(iter.next());
    stack.add(iter.next());
    while(iter.hasNext()) {
      Vector next = iter.next();
      Vector curr = stack.pop();
      Vector prev = stack.peek();
      while((stack.size() > 1) && !isConvex(prev, curr, next)) {
        curr = stack.pop();
        prev = stack.peek();
      }
      stack.add(curr);
      stack.add(next);
    }
    points = stack;
  }

  /**
   * Compute the convex hull, and return the resulting polygon.
   * 
   * @return Polygon of the hull
   */
  public Polygon getHull() {
    if(!ok) {
      computeConvexHull();
    }
    return new Polygon(points, minmaxX.getMin(), minmaxX.getMax(), minmaxY.getMin(), minmaxY.getMax());
  }
}