/**********************************************************************
 *
 * GEOS - Geometry Engine Open Source
 * http://geos.osgeo.org
 *
 * Copyright (C) 2013 Sandro Santilli <strk@kbt.io>
 *
 * This is free software; you can redistribute and/or modify it under
 * the terms of the GNU Lesser General Public Licence as published
 * by the Free Software Foundation.
 * See the COPYING file for more information.
 *
 **********************************************************************
 *
 * Last port: algorithm/Centroid.java r728 (JTS-0.13+)
 *
 **********************************************************************/

#include <geos/algorithm/Centroid.h>
#include <geos/algorithm/Orientation.h>
#include <geos/geom/CoordinateSequence.h>
#include <geos/geom/Geometry.h>
#include <geos/geom/Point.h>
#include <geos/geom/Polygon.h>
#include <geos/geom/GeometryCollection.h>
#include <geos/geom/LineString.h>


#include <cmath> // for std::abs

#include "geos/util.h"

using namespace geos::geom;

namespace geos {
namespace algorithm { // geos.algorithm

/* static public */
bool
Centroid::getCentroid(const Geometry& geom, CoordinateXY& pt)
{
    Centroid cent(geom);
    return cent.getCentroid(pt);
}

/* public */
bool
Centroid::getCentroid(CoordinateXY& cent) const
{
    if(std::abs(areasum2) > 0.0) {
        cent.x = cg3.x / 3 / areasum2;
        cent.y = cg3.y / 3 / areasum2;
    }
    else if(totalLength > 0.0) {
        // if polygon was degenerate, compute linear centroid instead
        cent.x = lineCentSum.x / totalLength;
        cent.y = lineCentSum.y / totalLength;
    }
    else if(ptCount > 0) {
        cent.x = ptCentSum.x / ptCount;
        cent.y = ptCentSum.y / ptCount;
    }
    else {
        return false;
    }
    return true;
}

/* private */
void
Centroid::add(const Geometry& geom)
{
    util::ensureNoCurvedComponents(geom);

    if(geom.isEmpty()) {
        return;
    }

    if(const Point* pt = dynamic_cast<const Point*>(&geom)) {
        addPoint(*pt->getCoordinate());
    }
    else if(const LineString* ls = dynamic_cast<const LineString*>(&geom)) {
        addLineSegments(*ls->getCoordinatesRO());
    }
    else if(const Polygon* p = dynamic_cast<const Polygon*>(&geom)) {
        add(*p);
    }
    else if(const GeometryCollection* g = dynamic_cast<const GeometryCollection*>(&geom)) {
        for(std::size_t i = 0; i < g->getNumGeometries(); i++) {
            add(*g->getGeometryN(i));
        }
    }
}

/* private */
void
Centroid::setAreaBasePoint(const CoordinateXY& basePt)
{
    areaBasePt.reset(new CoordinateXY(basePt));
}

/* private */
void
Centroid::add(const Polygon& poly)
{
    addShell(*poly.getExteriorRing()->getCoordinatesRO());
    for(std::size_t i = 0; i < poly.getNumInteriorRing(); i++) {
        addHole(*poly.getInteriorRingN(i)->getCoordinatesRO());
    }
}

/* private */
void
Centroid::addShell(const CoordinateSequence& pts)
{
    std::size_t len = pts.size();
    if(len > 0) {
        setAreaBasePoint(pts.getAt<CoordinateXY>(0));
    }
    bool isPositiveArea = ! Orientation::isCCW(&pts);
    for(std::size_t i = 0; i < len - 1; ++i) {
        addTriangle(*areaBasePt, pts.getAt<CoordinateXY>(i), pts.getAt<CoordinateXY>(i + 1), isPositiveArea);
    }
    addLineSegments(pts);
}

/* private */
void
Centroid::addHole(const CoordinateSequence& pts)
{
    if (pts.isEmpty()) {
        return;
    }

    bool isPositiveArea = Orientation::isCCW(&pts);
    for(std::size_t i = 0, e = pts.size() - 1; i < e; ++i) {
        addTriangle(*areaBasePt, pts.getAt<CoordinateXY>(i), pts.getAt<CoordinateXY>(i + 1), isPositiveArea);
    }
    addLineSegments(pts);
}

/* private */
void
Centroid::addTriangle(const CoordinateXY& p0, const CoordinateXY& p1, const CoordinateXY& p2, bool isPositiveArea)
{
    double sign = (isPositiveArea) ? 1.0 : -1.0;
    centroid3(p0, p1, p2, triangleCent3);
    double a2 = area2(p0, p1, p2);
    cg3.x += sign * a2 * triangleCent3.x;
    cg3.y += sign * a2 * triangleCent3.y;
    areasum2 += sign * a2;
}

/* static private */
void
Centroid::centroid3(const CoordinateXY& p1, const CoordinateXY& p2, const CoordinateXY& p3, CoordinateXY& c)
{
    c.x = p1.x + p2.x + p3.x;
    c.y = p1.y + p2.y + p3.y;
    return;
}

/* static private */
double
Centroid::area2(const CoordinateXY& p1, const CoordinateXY& p2, const CoordinateXY& p3)
{
    return
        (p2.x - p1.x) * (p3.y - p1.y) -
        (p3.x - p1.x) * (p2.y - p1.y);
}

/* private */
void
Centroid::addLineSegments(const CoordinateSequence& pts)
{
    std::size_t npts = pts.size();
    double lineLen = 0.0;
    for(std::size_t i = 0; i < npts - 1; i++) {
        double segmentLen = pts.getAt<CoordinateXY>(i).distance(pts.getAt<CoordinateXY>(i + 1));
        if(segmentLen == 0.0) {
            continue;
        }

        lineLen += segmentLen;

        double midx = (pts.getAt<CoordinateXY>(i).x + pts.getAt<CoordinateXY>(i + 1).x) / 2;
        lineCentSum.x += segmentLen * midx;
        double midy = (pts.getAt<CoordinateXY>(i).y + pts.getAt<CoordinateXY>(i + 1).y) / 2;
        lineCentSum.y += segmentLen * midy;
    }
    totalLength += lineLen;
    if(lineLen == 0.0 && npts > 0) {
        addPoint(pts[0]);
    }
}

/* private */
void
Centroid::addPoint(const CoordinateXY& pt)
{
    ptCount += 1;
    ptCentSum.x += pt.x;
    ptCentSum.y += pt.y;
}

} // namespace geos.algorithm
} // namespace geos