public class GridGenerator {

  // constants
  private static final int dimensions = 2;
  private static final RectDomain<1> Ddims = [1:dimensions], Ddirs = [-1:1:2];

  // parameters for splitting
  private static double MinFillRatio = 0.7;
  private static int MinInflection = 3;

  // domain that must contain all grids returned
  private static Domain<2> properDomain;

  // list of grids to be returned
  private static BoxedList_RectDomain2 recsGood;

  /*

    Returns an array of RectDomains that cover the given Domain.

    Uses algorithm of Berger & Rigoutsos,
    ``An Algorithm for Point Clustering and Grid Generation'',
    IEEE Trans. Syst., Man, Cybern., vol. 21, pp. 1278--1286 (1991).

  */

  public static final RectDomain<2> [1d] MakeGrids(Domain<2> tagged,
						   Domain<2> pnd) {

    // Return a set of grids covering the tagged domain,
    // all contained in pnd.

    properDomain = pnd;
    // Initialize recsGood with an empty list.
    recsGood = new BoxedList_RectDomain2();

    MakeGridsRecur(tagged);

    RectDomain<2> [1d] recsArray = recsGood.toArray();
    return recsArray;
  }


  private static final void MakeGridsRecur(Domain<2> tagged) {
    RectDomain<2> rec = tagged.boundingBox();
    Point<2> recMin = rec.min(), recMax = rec.max();

    // Check if grid is full enough.
    double fillratio = (float) tagged.size() / (float) rec.size();

    // Need to check also that grid does not contain points outside
    // properDomain.
    if (fillratio >= MinFillRatio && (rec - properDomain).isNull()) {
      // Full enough.  We'll use it.
      recsGood.push(rec);
      Point<2> ext = recMax - recMin + Point<2>.all(1);

//      System.out.println("Accepted " + rec.toString() +
//			 " dimensions " + ext.toString() +
//			 " fill ratio " + fillratio);

    } else {
      // Compute signatures and then cut somewhere.
      int [1d][1d] sig = signature(tagged);

      if (!cutAtHole(tagged, sig, recMin, recMax))
	if (!cutAtInflection(tagged, sig, recMin, recMax))
	  cutAtHalf(tagged, sig, recMin, recMax);
    }
  }

    
  private static final boolean
    cutAtHole(Domain<2> tagged, int [1d][1d] sig,
	      Point<2> recMin, Point<2> recMax) {
    // Look for a hole.
    // Split at the first hole in the first dimension where you find one.
    foreach (pdim in Ddims) {
      int dim = pdim[1];
      for (int hole = recMin[dim]; hole <= recMax[dim]; hole++) {
	if (sig[dim][hole] == 0) {
//	  System.out.println("splitting in dimension " + dim +
//			     " at hole at " + hole);
	  cut(tagged, dim, hole-1, hole+1);
	  return true;
	}
      }
    }
    return false;
  }


  private static final boolean
    cutAtInflection(Domain<2> tagged, int [1d][1d] sig,
		    Point<2> recMin, Point<2> recMax) {
    // Look for a strong enough inflection point (set MinInflection).
    // Split at the point of greatest inflection.
    boolean [1d] existInfl = new boolean[Ddims];
    boolean existInflAny = false;
    int [1d] bestLoc = new int[Ddims];
    int [1d] bestVal = new int[Ddims];
    foreach (pdim in Ddims) {
      int dim = pdim[1];
      RectDomain<1> Ddiff = [[recMin[dim] + 1]: [recMax[dim] - 2]];
      existInfl[pdim] = false;
      if (!Ddiff.isNull()) {
	bestVal[pdim] = 0;
	int [1d] sigdim = sig[pdim];
	foreach (p in Ddiff) {
	  int diffLaplacian =
	    Math.abs(sigdim[p + [2]] - 3 * sigdim[p + [1]]
		     + 3 * sigdim[p] - sigdim[p - [1]]);
	  if (diffLaplacian >= MinInflection) {
	    existInfl[pdim] = true;
	    existInflAny = true;
	    if (diffLaplacian > bestVal[pdim]) {
	      bestVal[pdim] = diffLaplacian;
	      bestLoc[pdim] = p[1];
	    }
	  }
	}
      }
    }
    if (existInflAny) {
      // Find best inflection point, and split there.
      int dimSplit = 1;
      for (int dim = 2; dim <= dimensions; dim++)
	if (existInfl[dim] && bestVal[dim] > bestVal[dimSplit])
	  dimSplit = dim;

//      System.out.println("splitting in dimension " + dimSplit +
//			     " at inflection at " + bestLoc[dimSplit] +
//			 " with value " + bestVal[dimSplit]);
      cut(tagged, dimSplit, bestLoc[dimSplit], bestLoc[dimSplit] + 1);
      return true;
    }
    return false;
  }


  private static final void
    cutAtHalf(Domain<2> tagged, int [1d][1d] sig,
	      Point<2> recMin, Point<2> recMax) {
    // Find longest side, and split there.
    Point<2> lengths = recMax - recMin;
    int dimSplit = 1;
    for (int dim = 2; dim <= dimensions; dim++)
      if (lengths[dim] > lengths[dimSplit]) dimSplit = dim;
    int splitpoint = (recMin[dimSplit] + recMax[dimSplit]) / 2;

//    System.out.println("splitting in dimension " +dimSplit+ " at half");
    cut(tagged, dimSplit, splitpoint, splitpoint + 1);
  }


  private static final void
    cut(Domain<2> tagged, int dim, int end1, int start2) {

    // Cut the domain tagged along the dimension dim,
    // with the first part ending at coordinate end1 and the second
    // part beginning at coordinate start2.

    Point<2> tagMin = tagged.min(), tagMax = tagged.max();
    Point<2> unit = Point<2>.direction(dim);
    Point<2> wtlim = Point<2>.all(1) - unit;

    Point<2> pend1 = wtlim * tagMax + unit * end1;
    Point<2> pstart2 = wtlim * tagMin + unit * start2;

    MakeGridsRecur(tagged * [tagMin : pend1]);
    MakeGridsRecur(tagged * [pstart2 : tagMax]);
  }


  private static final int [1d][1d] signature(Domain<2> tagged) {
    int [1d][1d] sig = new int[Ddims][1d];
    // sig[1:dimensions] returns an array of int counting the number of
    // filled-in cells in each cross-section.

    Point<2> tagMin = tagged.min();
    Point<2> tagMax = tagged.max();
    foreach (pdim in Ddims) {
      int dim = pdim[1];
      RectDomain<1> tagProj = [tagMin[dim] : tagMax[dim]];
      sig[dim] = new int[tagProj];
      foreach (p1 in tagProj) sig[dim][p1] = 0;
    }

    foreach (p in tagged) foreach (pdim in Ddims) sig[pdim][p[pdim[1]]]++;

    return sig;
  }
}