/**********************************************************************

 *

 * GEOS - Geometry Engine Open Source

 * http://geos.osgeo.org

 *

 * Copyright (C) 2020 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: operation/overlayng/ElevationModel.java 4c88fea52

 *

 **********************************************************************/

#include <geos/operation/overlayng/ElevationModel.h>

#include <geos/geom/Coordinate.h>

#include <geos/geom/Geometry.h>

#include <geos/geom/CoordinateSequence.h>

#include <geos/geom/CoordinateFilter.h>

#include <geos/geom/CoordinateSequenceFilter.h>

#include <geos/geom/Envelope.h>

#include <memory> // for unique_ptr

#ifndef GEOS_DEBUG

# define GEOS_DEBUG 0

#endif

#if GEOS_DEBUG

# include <iostream>

#endif

namespace geos {      // geos.

namespace operation { // geos.operation

namespace overlayng { // geos.operation.overlayng

namespace {

    template<typename T>

    T clamp(const T& v, const T& low, const T& high) {

        return v < low ? low : high < v ? high : v;

    }

}

using geos::geom::Coordinate;

using geos::geom::CoordinateSequence;

using geos::geom::Geometry;

using geos::geom::Envelope;

/* public static */

std::unique_ptr<ElevationModel>

ElevationModel::create(const geom::Geometry& geom1, const geom::Geometry& geom2)

{

    Envelope extent;

    if (! geom1.isEmpty() ) {

      extent.expandToInclude(geom1.getEnvelopeInternal());

    }

    if (! geom2.isEmpty() ) {

      extent.expandToInclude(geom2.getEnvelopeInternal());

    }

    std::unique_ptr<ElevationModel> model ( new ElevationModel(extent, DEFAULT_CELL_NUM, DEFAULT_CELL_NUM) );

    if (! geom1.isEmpty() ) model->add(geom1);

    if (! geom2.isEmpty() ) model->add(geom2);

    return model;

}

/* public static */

std::unique_ptr<ElevationModel>

ElevationModel::create(const geom::Geometry& geom1)

{

    Envelope extent;

    if (! geom1.isEmpty() ) {

      extent.expandToInclude(geom1.getEnvelopeInternal());

    }

    std::unique_ptr<ElevationModel> model ( new ElevationModel(extent, DEFAULT_CELL_NUM, DEFAULT_CELL_NUM) );

    if (! geom1.isEmpty() ) model->add(geom1);

    return model;

}

const int ElevationModel::DEFAULT_CELL_NUM = 3;

ElevationModel::ElevationModel(const Envelope& nExtent, int nNumCellX, int nNumCellY)

    :

    extent(nExtent),

    numCellX(nNumCellX),

    numCellY(nNumCellY)

{

    cellSizeX = extent.getWidth() / numCellX;

    cellSizeY = extent.getHeight() / numCellY;

    if(cellSizeX <= 0.0) {

        numCellX = 1;

    }

    if(cellSizeY <= 0.0) {

        numCellY = 1;

    }

    cells.resize(static_cast<std::size_t>(numCellX) * static_cast<std::size_t>(numCellY));

}

/* public */

void

ElevationModel::add(const Geometry& geom)

{

    class Filter: public geom::CoordinateSequenceFilter {

        ElevationModel& model;

        bool hasZ;

    public:

        Filter(ElevationModel& nModel) : model(nModel), hasZ(true) {}

        void filter_ro(const geom::CoordinateSequence& seq, std::size_t i) override

        {

            if (! seq.hasZ()) {

                hasZ = false;

                return;

            }

            const Coordinate& c = seq.getAt(i);

#if GEOS_DEBUG

            std::cerr << "Coordinate " << i << " of added geom is: "

<< c << std::endl;

#endif

            model.add(c.x, c.y, c.z);

        }

        bool isDone() const override {

            // no need to scan if no Z present

            return ! hasZ;

        }

        bool isGeometryChanged() const override {

            return false;

        }

    };

    Filter filter(*this);

    geom.apply_ro(filter);

}

/* protected */

void

ElevationModel::add(double x, double y, double z)

{

    if (std::isnan(z))

      return;

    hasZValue = true;

    ElevationCell& cell = getCell(x, y); //, true);

    cell.add(z);

}

/* private */

void

ElevationModel::init()

{

    isInitialized = true;

    int numCells = 0;

    double sumZ = 0.0;

#if GEOS_DEBUG

    int offset = 0;

#endif

    for (ElevationCell& cell : cells )

    {

        if (!cell.isNull()) {

          cell.compute();

#if GEOS_DEBUG

          std::cerr << "init: cell" << offset

                    << " getZ: " << cell.getZ()

                    << std::endl;

#endif

          numCells++;

          sumZ += cell.getZ();

        }

#if GEOS_DEBUG

        ++offset;

#endif

    }

    averageZ = DoubleNotANumber;

    if (numCells > 0) {

      averageZ = sumZ / numCells;

    }

#if GEOS_DEBUG

    std::cerr << "init: numCells: " << numCells

              << " averageZ: " << averageZ << std::endl;

#endif

}

/* public */

double

ElevationModel::getZ(double x, double y)

{

    if (! isInitialized)

    {

      init();

    }

    const ElevationModel::ElevationCell& cell = getCell(x, y); //, false);

    if (cell.isNull())

    {

      return averageZ;

    }

    return cell.getZ();

}

/* public */

void

ElevationModel::populateZ(Geometry& geom)

{

    // short-circuit if no Zs are present in model

    if (! hasZValue)

      return;

    if (! isInitialized)

      init();

    class Filter: public geom::CoordinateFilter {

        ElevationModel& model;

    public:

        Filter(ElevationModel& nModel) : model(nModel) {}

        void filter_rw(geom::Coordinate* c) const override

        {

#if GEOS_DEBUG

            std::cerr << "Input coord:  " << *c << std::endl;

#endif

            if (std::isnan( c->z )) {

                c->z = model.getZ(c->x, c->y);

#if GEOS_DEBUG

                std::cerr << "New coord: " << *c << std::endl;

#endif

            }

        }

    };

    Filter filter(*this);

#if GEOS_DEBUG

    std::cerr << "ElevationModel::poplateZ calling apply_rw(CoordinateSequenceFilter&) against" <<

              //std::type_name<decltype(ci)>()

              typeid(geom).name()

              << std::endl;

#endif

    geom.apply_rw(&filter);

}

/* private */

ElevationModel::ElevationCell&

ElevationModel::getCell(double x, double y) //, bool isCreateIfMissing

{

    int ix = 0;

    if (numCellX > 1) {

      ix = (int) ((x - extent.getMinX()) / cellSizeX);

      ix = clamp(ix, 0, numCellX - 1);

    }

    int iy = 0;

    if (numCellY > 1) {

      iy = (int) ((y - extent.getMinY()) / cellSizeY);

      iy = clamp(iy, 0, numCellY - 1);

    }

    int cellOffset = getCellOffset(ix, iy);

#if GEOS_DEBUG

    std::cerr << "Cell of coordinate " << x << "," << y

              << " is " << ix << "," << iy

              << " offset " << cellOffset << std::endl;

#endif

    assert ( cellOffset < numCellX * numCellY );

    ElevationModel::ElevationCell& cell = cells[static_cast<std::size_t>(cellOffset)];

    return cell;

}

} // namespace geos.operation.overlayng

} // namespace geos.operation

} // namespace geos