/*
 *    GeoTools - The Open Source Java GIS Toolkit
 *    http://geotools.org
 *    
 *    (C) 2001-2006  Vivid Solutions
 *    (C) 2001-2008, Open Source Geospatial Foundation (OSGeo)
 *    
 *    This library is free software; you can redistribute it and/or
 *    modify it under the terms of the GNU Lesser General Public
 *    License as published by the Free Software Foundation;
 *    version 2.1 of the License.
 *
 *    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 GNU
 *    Lesser General Public License for more details.
 */
package org.geotools.geometry.iso.util.algorithm2D;

import java.util.Iterator;
import java.util.List;

import org.geotools.geometry.iso.aggregate.MultiSurfaceImpl;
import org.geotools.geometry.iso.coordinate.DirectPositionImpl;
import org.geotools.geometry.iso.primitive.CurveImpl;
import org.geotools.geometry.iso.primitive.RingImpl;
import org.geotools.geometry.iso.primitive.RingImplUnsafe;
import org.geotools.geometry.iso.primitive.SurfaceBoundaryImpl;
import org.geotools.geometry.iso.primitive.SurfaceImpl;
import org.geotools.geometry.iso.root.GeometryImpl;
import org.opengis.geometry.DirectPosition;
import org.opengis.geometry.primitive.OrientableSurface;
import org.opengis.referencing.crs.CoordinateReferenceSystem;



/**
 * Computes the centroid of an area geometry.
 * <h2>Algorithm</h2>
 * Based on the usual algorithm for calculating the centroid as a weighted sum
 * of the centroids of a decomposition of the area into (possibly overlapping)
 * triangles. The algorithm has been extended to handle holes and
 * multi-polygons. See
 * <code>http://www.faqs.org/faqs/graphics/algorithms-faq/</code> for further
 * details of the basic approach.
 *
 *
 * @source $URL$
 */
public class CentroidArea2D {

    static final int X = 0;
    static final int Y = 1;
    static final int Z = 2;
	//private FeatGeomFactoryImpl factory = null;
	private CoordinateReferenceSystem crs = null;
	
	// the point all triangles are based at
	private DirectPosition basePt = null;

	// partial area sum
	private double areasum2 = 0;

	// partial centroid sum
	//private DirectPositionImpl cg3 = new DirectPositionImpl();
	double centSumX = 0.0;
	double centSumY = 0.0;
	double centSumZ = 0.0;
	

	/**
	 * Creates a new Centroid operation
	 * 
	 * @param crs
	 */
	public CentroidArea2D(CoordinateReferenceSystem crs) {
		this.crs = crs;
		this.basePt = null;
		
	}

	/**
	 * Adds the area defined by a Geometry to the centroid total. If the
	 * geometry has no area it does not contribute to the centroid.
	 * 
	 * @param geom
	 *            the geometry to add
	 */
	public void add(GeometryImpl geom) {
		if (geom instanceof SurfaceImpl) {
			SurfaceBoundaryImpl sb = ((SurfaceImpl) geom).getBoundary();
			this.setBasePoint( ((CurveImpl)sb.getExterior().getGenerators().iterator().next()).getStartPoint());
			this.addSurface(sb);
		} else if (geom instanceof MultiSurfaceImpl) {
			Iterator<OrientableSurface> surfaces = ((MultiSurfaceImpl) geom).getElements().iterator();
			while (surfaces.hasNext()) {
				this.add((GeometryImpl) surfaces.next());
			}
		}
	}

	public DirectPositionImpl getCentroid() {
		DirectPositionImpl centroid =  new DirectPositionImpl(crs); //this.factory.getGeometryFactoryImpl().createDirectPosition();
		centroid.setX(this.centSumX / 3 / this.areasum2);
		centroid.setY(this.centSumY / 3 / this.areasum2);
		return centroid;
	}

	private void setBasePoint(DirectPosition basePt) {
		if (this.basePt == null)
			this.basePt = basePt;
	}

	private void addSurface(SurfaceBoundaryImpl sb) {
		this.addShell(((RingImplUnsafe)sb.getExterior()).asDirectPositions());
		for (int i = 0; i < sb.getInteriors().size(); i++) {
			this.addHole(((RingImplUnsafe)sb.getInteriors().get(i)).asDirectPositions());
		}
	}

	private void addShell(List<DirectPosition> pts) {
		boolean isPositiveArea = !CGAlgorithms.isCCW(pts);
		for (int i = 0; i < pts.size() - 1; i++) {
			addTriangle(basePt, pts.get(i), pts.get(i+1), isPositiveArea);
		}
	}

	private void addHole(List<DirectPosition> pts) {
		boolean isPositiveArea = CGAlgorithms.isCCW(pts);
		for (int i = 0; i < pts.size() - 1; i++) {
			addTriangle(basePt, pts.get(i), pts.get(i+1), isPositiveArea);
		}
	}

	private void addTriangle(DirectPosition p0, DirectPosition p1, DirectPosition p2,
			boolean isPositiveArea) {
		double sign = (isPositiveArea) ? 1.0 : -1.0;

		//this.centroid3(p0, p1, p2);
		double tempSumX = p0.getOrdinate(X) + p1.getOrdinate(X) + p2.getOrdinate(X);
		double tempSumY = p0.getOrdinate(Y) + p1.getOrdinate(Y)+ p2.getOrdinate(Y);
		//double tempSumZ = 0.0;

		double area2 = area2(p0, p1, p2);
		this.centSumX += sign * area2 * tempSumX;
		this.centSumY += sign * area2 * tempSumY;
		this.areasum2 += sign * area2;
	}

	/**
	 * Returns twice the signed area of the triangle p1-p2-p3, positive if a,b,c
	 * are oriented ccw, and negative if cw.
	 */	
	private static double area2(DirectPosition p1, DirectPosition p2, DirectPosition  p3) {
		return (p2.getOrdinate(X) - p1.getOrdinate(X)) * (p3.getOrdinate(Y) - p1.getOrdinate(Y))
		    - (p3.getOrdinate(X) - p1.getOrdinate(X)) * (p2.getOrdinate(Y) - p1.getOrdinate(Y));
	}

}