ChainHull.cpp

Go to the documentation of this file.

/*
 * Copyright 2011-2016 Universiteit Antwerpen
 *
 * Licensed under the EUPL, Version 1.1 or  as soon they will be approved by
 * the European Commission - subsequent versions of the EUPL (the "Licence");
 * You may not use this work except in compliance with the Licence.
 * You may obtain a copy of the Licence at: http://ec.europa.eu/idabc/eupl5
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the Licence is distributed on an "AS IS" basis,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the Licence for the specific language governing
 * permissions and limitations under the Licence.
 */
#include "ChainHull.h"
#include <algorithm>

namespace {
        using SimPT_Sim::ChainHull;

        inline float isLeft(ChainHull::Point P0, ChainHull::Point P1, ChainHull::Point P2)
        {
                return (P1.GetX() - P0.GetX())*(P2.GetY() - P0.GetY()) - (P2.GetX() - P0.GetX())*(P1.GetY() - P0.GetY());
        }
}

namespace SimPT_Sim {

std::vector<std::array<double, 2> > ChainHull::Compute(const std::vector<std::array<double, 2> >& polygon)
{
        // Step 1: Prepare data for 2D hull code - get boundary polygon
        int pc = 0;
        Point* p = new Point[polygon.size() + 1];
        for (auto const& vertex : polygon) {
                p[pc++] = Point(vertex[0], vertex[1]);
        }

        // chainHull algorithm requires sorted points
        std::sort(p, p + pc);

        // Step 2: Call 2D hull code
        int np = polygon.size();
        Point* hull = new Point[np + 1];
        int nph = chainHull_2D(p, np, hull);

        // Step 3: Convert hull array to vector
        std::vector<std::array<double, 2> > result;
        for (int i = 0; i < nph - 1; i++) {
                result.push_back({ {hull[i].GetX(), hull[i].GetY()} });
        }

        delete[] p;
        delete[] hull;

        return result;
}

// required to sort points (e.g. Qt's qSort())
bool operator<(ChainHull::Point const& p1, ChainHull::Point const& p2)
{
        bool smaller = true;
        if (p1.GetY() < p2.GetY()) {
                smaller = true;
        } else if (p1.GetY() > p2.GetY()) {
                smaller = false;
        } else if (p1.GetX() < p2.GetX()) {
                smaller = true;
        } else {
                smaller = false;
        }
        return smaller;
}

int ChainHull::chainHull_2D(Point* P, int n, Point* H)
{
        // the output array H[] will be used as the stack
        int bot = 0;
        int top = (-1);  // indices for bottom and top of the stack
        int i;           // array scan index

        // Get the indices of points with min x-coord and min|max y-coord
        int minmin = 0;
        int minmax;
        float xmin = P[0].GetX();
        for (i = 1; i < n; i++) {
                if (P[i].GetX() != xmin) {
                        break;
                }
        }
        minmax = i - 1;
        if (minmax == n - 1) {       // degenerate case: all x-coords == xmin
                H[++top] = P[minmin];
                if (P[minmax].GetY() != P[minmin].GetY()) { // a nontrivial segment
                        H[++top] = P[minmax];
                }
                H[++top] = P[minmin];           // add polygon endpoint
                return top + 1;
        }

        // Get the indices of points with max x-coord and min|max y-coord
        int maxmin;
        int maxmax = n - 1;
        float xmax = P[n - 1].GetX();
        for (i = n - 2; i >= 0; i--) {
                if (P[i].GetX() != xmax) {
                        break;
                }
        }
        maxmin = i + 1;

        // Compute the lower hull on the stack H
        H[++top] = P[minmin];      // push minmin point onto stack
        i = minmax;
        while (++i <= maxmin) {
                // the lower line joins P[minmin] with P[maxmin]
                if (isLeft(P[minmin], P[maxmin], P[i]) >= 0 && i < maxmin) {
                        continue;      // ignore P[i] above or on the lower line
                }
                while (top > 0)        // there are at least 2 points on the stack
                {
                        // test if P[i] is left of the line at the stack top
                        if (isLeft(H[top - 1], H[top], P[i]) > 0)
                                break;         // P[i] is a new hull vertex
                        else
                                top--;         // pop top point off stack
                }
                H[++top] = P[i];       // push P[i] onto stack
        }

        // Next, compute upper hull on stack H above bottom hull
        if (maxmax != maxmin) {        // if distinct xmax points
                H[++top] = P[maxmax];  // push maxmax point onto stack
        }
        bot = top;                     // the bottom point of the upper hull stack
        i = maxmin;
        while (--i >= minmax) {
                // the upper line joins P[maxmax] with P[minmax]
                if (isLeft(P[maxmax], P[minmax], P[i]) >= 0 && i > minmax) {
                        continue;      // ignore P[i] below or on the upper line
                }
                while (top > bot)      // at least 2 points on the upper stack
                {
                        // test if P[i] is left of the line at the stack top
                        if (isLeft(H[top - 1], H[top], P[i]) > 0)
                                break;         // P[i] is a new hull vertex
                        else
                                top--;         // pop top point off stack
                }
                H[++top] = P[i];       // push P[i] onto stack
        }
        if (minmax != minmin) {
                H[++top] = P[minmin];  // push joining endpoint onto stack
        }

        return top + 1;
}

} // namespace