path: root/src/quadtree.hpp

/******************************************************************************
 * Project:  PROJ
 * Purpose:  Implementation of quadtree building and searching functions.
 *           Derived from shapelib, mapserver and GDAL implementations
 * Author:   Even Rouault, <even.rouault at spatialys.com>
 *
 ******************************************************************************
 * Copyright (c) 1999-2008, Frank Warmerdam
 * Copyright (c) 2008-2020, Even Rouault <even dot rouault at spatialys.com>
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 *****************************************************************************/

#ifndef QUADTREE_HPP
#define QUADTREE_HPP

#include "proj/util.hpp"

#include <functional>
#include <vector>

//! @cond Doxygen_Suppress

NS_PROJ_START

namespace QuadTree {

/* -------------------------------------------------------------------- */
/*      If the following is 0.5, psNodes will be split in half.  If it  */
/*      is 0.6 then each apSubNode will contain 60% of the parent       */
/*      psNode, with 20% representing overlap.  This can be help to     */
/*      prevent small objects on a boundary from shifting too high      */
/*      up the hQuadTree.                                               */
/* -------------------------------------------------------------------- */
constexpr double DEFAULT_SPLIT_RATIO = 0.55;

/** Describe a rectangle */
struct RectObj {
    double minx = 0; /**< Minimum x */
    double miny = 0; /**< Minimum y */
    double maxx = 0; /**< Maximum x */
    double maxy = 0; /**< Maximum y */

    /* Returns whether this rectangle is contained by other */
    inline bool isContainedBy(const RectObj &other) const {
        return minx >= other.minx && maxx <= other.maxx && miny >= other.miny &&
               maxy <= other.maxy;
    }

    /* Returns whether this rectangles overlaps other */
    inline bool overlaps(const RectObj &other) const {
        return minx <= other.maxx && maxx >= other.minx && miny <= other.maxy &&
               maxy >= other.miny;
    }

    /* Returns whether this rectangles contains the specified point */
    inline bool contains(double x, double y) const {
        return minx <= x && maxx >= x && miny <= y && maxy >= y;
    }

    /* Return whether this rectangles is different from other */
    inline bool operator!=(const RectObj &other) const {
        return minx != other.minx || miny != other.miny || maxx != other.maxx ||
               maxy != other.maxy;
    }
};

/** Quadtree */
template <class Feature> class QuadTree {

    struct Node {
        RectObj rect{}; /* area covered by this psNode */

        /* list of shapes stored at this node. */
        std::vector<std::pair<Feature, RectObj>> features{};

        std::vector<Node> subnodes{};

        explicit Node(const RectObj &rectIn) : rect(rectIn) {}
    };

    Node root{};
    unsigned nBucketCapacity = 8;
    double dfSplitRatio = DEFAULT_SPLIT_RATIO;

  public:
    /** Construct a new quadtree with the global bounds of all objects to be
    * inserted */
    explicit QuadTree(const RectObj &globalBounds) : root(globalBounds) {}

    /** Add a new feature, with its bounds specified in featureBounds */
    void insert(const Feature &feature, const RectObj &featureBounds) {
        insert(root, feature, featureBounds);
    }

#ifdef UNUSED
    /** Retrieve all features whose bounds intersects aoiRect */
    void
    search(const RectObj &aoiRect,
           std::vector<std::reference_wrapper<const Feature>> &features) const {
        search(root, aoiRect, features);
    }
#endif

    /** Retrieve all features whose bounds contains (x,y) */
    void search(double x, double y, std::vector<Feature> &features) const {
        search(root, x, y, features);
    }

  private:
    void splitBounds(const RectObj &in, RectObj &out1, RectObj &out2) {
        // The output bounds will be very similar to the input bounds,
        // so just copy over to start.
        out1 = in;
        out2 = in;

        // Split in X direction.
        if ((in.maxx - in.minx) > (in.maxy - in.miny)) {
            const double range = in.maxx - in.minx;

            out1.maxx = in.minx + range * dfSplitRatio;
            out2.minx = in.maxx - range * dfSplitRatio;
        }

        // Otherwise split in Y direction.
        else {
            const double range = in.maxy - in.miny;

            out1.maxy = in.miny + range * dfSplitRatio;
            out2.miny = in.maxy - range * dfSplitRatio;
        }
    }

    void insert(Node &node, const Feature &feature,
                const RectObj &featureBounds) {
        if (node.subnodes.empty()) {
            // If we have reached the max bucket capacity, try to insert
            // in a subnode if possible.
            if (node.features.size() >= nBucketCapacity) {
                RectObj half1;
                RectObj half2;
                RectObj quad1;
                RectObj quad2;
                RectObj quad3;
                RectObj quad4;

                splitBounds(node.rect, half1, half2);
                splitBounds(half1, quad1, quad2);
                splitBounds(half2, quad3, quad4);

                if (node.rect != quad1 && node.rect != quad2 &&
                    node.rect != quad3 && node.rect != quad4 &&
                    (featureBounds.isContainedBy(quad1) ||
                     featureBounds.isContainedBy(quad2) ||
                     featureBounds.isContainedBy(quad3) ||
                     featureBounds.isContainedBy(quad4))) {
                    node.subnodes.reserve(4);
                    node.subnodes.emplace_back(Node(quad1));
                    node.subnodes.emplace_back(Node(quad2));
                    node.subnodes.emplace_back(Node(quad3));
                    node.subnodes.emplace_back(Node(quad4));

                    auto features = std::move(node.features);
                    node.features.clear();
                    for (auto &pair : features) {
                        insert(node, pair.first, pair.second);
                    }

                    /* recurse back on this psNode now that it has apSubNodes */
                    insert(node, feature, featureBounds);
                    return;
                }
            }
        } else {
            // If we have sub nodes, then consider whether this object will
            // fit in them.
            for (auto &subnode : node.subnodes) {
                if (featureBounds.isContainedBy(subnode.rect)) {
                    insert(subnode, feature, featureBounds);
                    return;
                }
            }
        }

        // If none of that worked, just add it to this nodes list.
        node.features.push_back(
            std::pair<Feature, RectObj>(feature, featureBounds));
    }

#ifdef UNUSED
    void
    search(const Node &node, const RectObj &aoiRect,
           std::vector<std::reference_wrapper<const Feature>> &features) const {
        // Does this node overlap the area of interest at all?  If not,
        // return without adding to the list at all.
        if (!node.rect.overlaps(aoiRect))
            return;

        // Add the local features to the list.
        for (const auto &pair : node.features) {
            if (pair.second.overlaps(aoiRect)) {
                features.push_back(
                    std::reference_wrapper<const Feature>(pair.first));
            }
        }

        // Recurse to subnodes if they exist.
        for (const auto &subnode : node.subnodes) {
            search(subnode, aoiRect, features);
        }
    }
#endif

    void search(const Node &node, double x, double y,
                std::vector<Feature> &features) const {
        // Does this node overlap the area of interest at all?  If not,
        // return without adding to the list at all.
        if (!node.rect.contains(x, y))
            return;

        // Add the local features to the list.
        for (const auto &pair : node.features) {
            if (pair.second.contains(x, y)) {
                features.push_back(pair.first);
            }
        }

        // Recurse to subnodes if they exist.
        for (const auto &subnode : node.subnodes) {
            search(subnode, x, y, features);
        }
    }
};

} // namespace QuadTree

NS_PROJ_END

//! @endcond

#endif //  QUADTREE_HPP