Convex hull (tools/convex_hull.hpp)

• View this file on GitHub
• View document part on GitHub
• Last update: 2024-03-23 22:25:26+09:00
• Include: #include "tools/convex_hull.hpp"

template <typename InputIterator, typename OutputIterator>
void convex_hull(InputIterator begin, InputIterator end, bool minimum, OutputIterator result);

It stores the indices of vertices or the vertices themselves which are on the edge of the convex hull of a given polygon to result. If minimum is true, only the minimal vertices required for the convex hull will be stored. On the other hand, if minimum is false, all the vertices on the edge of the convex hull will be stored.

The vertices will be stored counterclockwisely. The first sroted vertex will be the leftmost vertex. If the number of the leftmost vertices is more than one, the first stored vertex will be the lowermost vertex in them.

Constraints

• begin $\leq$ end
• The type of *begin is tools::vector2<T>.
• Either std::size_t or tools::vector2<T> is assignable to *result.

Time Complexity

• $O(n \log n)$ where $n$ is end - begin

License

• CC0

Author

• anqooqie

Depends on

• std::abs(x) extended for my library (tools/abs.hpp)
• Counter clockwise function (tools/ccw.hpp)
• Combine hash values (tools/hash_combine.hpp)
• std::less by key (tools/less_by.hpp)
• The number of nanoseconds that have elapsed since epoch (tools/now.hpp)
• Hash of std::tuple (tools/tuple_hash.hpp)
• Vector (tools/vector.hpp)
• Two dimensional vector (tools/vector2.hpp)

Verified with

• tests/convex_hull.test.cpp
• tests/polygon_2d/minimum_bounding_circle.test.cpp

Code

#ifndef TOOLS_CONVEX_HULL_HPP
#define TOOLS_CONVEX_HULL_HPP

#include <type_traits>
#include <vector>
#include <cstddef>
#include <numeric>
#include <algorithm>
#include <utility>
#include <iterator>
#include <stack>
#include "tools/vector2.hpp"
#include "tools/less_by.hpp"
#include "tools/ccw.hpp"

namespace tools {
  template <typename InputIterator, typename OutputIterator>
  void convex_hull(const InputIterator begin, const InputIterator end, bool minimum, OutputIterator result) {
    using T = ::std::decay_t<decltype(begin->x)>;

    const ::std::vector<::tools::vector2<T>> v(begin, end);
    ::std::vector<::std::size_t> a(v.size());
    ::std::iota(a.begin(), a.end(), 0);
    ::std::sort(a.begin(), a.end(), ::tools::less_by([&](const ::std::size_t i) {
      return ::std::make_pair(v[i].x, v[i].y);
    }));
    ::std::vector<::std::vector<::std::size_t>> duplicates;

    if (minimum) {
      ::std::size_t vl = 0;
      for (::std::size_t vr = 0, al = 0, ar = 0; al < a.size(); vl = vr, al = ar) {
        for (; ar < a.size() && v[a[al]].x == v[a[ar]].x; ++vr, ++ar);
        if (vl < al) ::std::move(::std::next(a.begin(), al), ::std::next(a.begin(), ar), ::std::next(a.begin(), vl));
        if (v[a[vl]].y == v[a[vr - 1]].y) {
          vr -= vr - vl - 1;
          duplicates.emplace_back();
          duplicates.back().push_back(a[vl]);
        } else {
          ::std::swap(a[vl + 1], a[vr - 1]);
          vr -= vr - vl - 2;
          duplicates.emplace_back();
          duplicates.back().push_back(a[vl]);
          duplicates.emplace_back();
          duplicates.back().push_back(a[vl + 1]);
        }
      }
      a.erase(::std::next(a.begin(), vl), a.end());
    } else {
      ::std::size_t vl = 0;
      for (::std::size_t vr = 0, al = 0, ar = 0; al < a.size(); vl = vr, al = ar) {
        for (; ar < a.size() && v[a[al]] == v[a[ar]]; ++vr, ++ar);
        if (vl < al) ::std::move(::std::next(a.begin(), al), ::std::next(a.begin(), ar), ::std::next(a.begin(), vl));
        duplicates.emplace_back();
        for (::std::size_t i = vl; i < vr; ++i) {
          duplicates.back().push_back(a[i]);
        }
        vr -= vr - vl - 1;
      }
      a.erase(::std::next(a.begin(), vl), a.end());
    }

    ::std::vector<::std::size_t> convex_hull;
    if (a.size() >= 3) {

      convex_hull.push_back(0);
      convex_hull.push_back(1);
      for (::std::size_t p3 = 2; p3 < a.size(); ++p3) {
        while (convex_hull.size() >= 2) {
          const int ccw = ::tools::ccw(v[a[convex_hull.rbegin()[1]]], v[a[convex_hull.back()]], v[a[p3]]);
          if (ccw == 1 || (!minimum && ccw == -2)) {
            break;
          }
          convex_hull.pop_back();
        }
        convex_hull.push_back(p3);
      }

      const ::std::size_t threshold = convex_hull.size() + 1;
      for (::std::size_t p3 = convex_hull.back(); p3 --> 0;) {
        while (convex_hull.size() >= threshold) {
          const int ccw = ::tools::ccw(v[a[convex_hull.rbegin()[1]]], v[a[convex_hull.back()]], v[a[p3]]);
          if (ccw == 1 || (!minimum && ccw == -2)) {
            break;
          }
          convex_hull.pop_back();
        }
        convex_hull.push_back(p3);
      }
      convex_hull.pop_back();

    } else {
      for (::std::size_t i = 0; i < a.size(); ++i) {
        convex_hull.push_back(i);
      }
    }

    for (const ::std::size_t& c : convex_hull) {
      for (const ::std::size_t& i : duplicates[c]) {
        if constexpr (::std::is_assignable_v<OutputIterator, ::std::size_t>) {
          *result = i;
        } else {
          *result = v[i];
        }
        ++result;
      }
    }
  }
}

#endif

#line 1 "tools/convex_hull.hpp"



#include <type_traits>
#include <vector>
#include <cstddef>
#include <numeric>
#include <algorithm>
#include <utility>
#include <iterator>
#include <stack>
#line 1 "tools/vector2.hpp"



#line 1 "tools/vector.hpp"



#line 5 "tools/vector.hpp"
#include <array>
#include <initializer_list>
#include <cassert>
#include <limits>
#line 14 "tools/vector.hpp"
#include <cmath>
#include <iostream>
#include <string>
#include <functional>
#include <tuple>
#line 1 "tools/abs.hpp"



namespace tools {
  constexpr float abs(const float x) {
    return x < 0 ? -x : x;
  }
  constexpr double abs(const double x) {
    return x < 0 ? -x : x;
  }
  constexpr long double abs(const long double x) {
    return x < 0 ? -x : x;
  }
  constexpr int abs(const int x) {
    return x < 0 ? -x : x;
  }
  constexpr long abs(const long x) {
    return x < 0 ? -x : x;
  }
  constexpr long long abs(const long long x) {
    return x < 0 ? -x : x;
  }
  constexpr unsigned int abs(const unsigned int x) {
    return x;
  }
  constexpr unsigned long abs(const unsigned long x) {
    return x;
  }
  constexpr unsigned long long abs(const unsigned long long x) {
    return x;
  }
}


#line 1 "tools/tuple_hash.hpp"



#line 1 "tools/now.hpp"



#include <chrono>

namespace tools {
  inline long long now() {
    return ::std::chrono::duration_cast<::std::chrono::nanoseconds>(::std::chrono::high_resolution_clock::now().time_since_epoch()).count();
  }
}


#line 1 "tools/hash_combine.hpp"



#line 6 "tools/hash_combine.hpp"

// Source: https://github.com/google/cityhash/blob/f5dc54147fcce12cefd16548c8e760d68ac04226/src/city.h
// License: MIT
// Author: Google Inc.

// Copyright (c) 2011 Google, Inc.
//
// 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.

namespace tools {
  template <typename T>
  void hash_combine(::std::size_t& seed, const T& v) {
    static const ::std::hash<T> hasher;
    static constexpr ::std::size_t k_mul = 0x9ddfea08eb382d69ULL;
    ::std::size_t a = (hasher(v) ^ seed) * k_mul;
    a ^= (a >> 47);
    ::std::size_t b = (seed ^ a) * k_mul;
    b ^= (b >> 47);
    seed = b * k_mul;
  }
}


#line 11 "tools/tuple_hash.hpp"

namespace tools {
  template <typename... Ts>
  struct tuple_hash {
    template <::std::size_t I = sizeof...(Ts) - 1>
    ::std::size_t operator()(const ::std::tuple<Ts...>& key) const {
      if constexpr (I == ::std::numeric_limits<::std::size_t>::max()) {
        static const ::std::size_t seed = ::tools::now();
        return seed;
      } else {
        ::std::size_t seed = this->operator()<I - 1>(key);
        ::tools::hash_combine(seed, ::std::get<I>(key));
        return seed;
      }
    }
  };
}


#line 21 "tools/vector.hpp"

namespace tools {
  namespace detail {
    namespace vector {
      template <typename T, ::std::size_t N>
      class members {
      protected:
        constexpr static bool variable_sized = false;
        constexpr static bool has_aliases = false;
        ::std::array<T, N> m_values;
        members() : m_values() {}
        members(const ::std::initializer_list<T> il) : m_values(il) {
          assert(il.size() == N);
        }
      };

      template <typename T>
      class members<T, 2> {
      protected:
        constexpr static bool variable_sized = false;
        constexpr static bool has_aliases = true;
        ::std::array<T*, 2> m_values;
        members() : m_values{&this->x, &this->y} {}
        members(const T& x, const T& y) : m_values{&this->x, &this->y}, x(x), y(y) {}
        members(const ::std::initializer_list<T> il) : m_values{&this->x, &this->y}, x(il.begin()[0]), y(il.begin()[1]) {
          assert(il.size() == 2);
        }

      public:
        T x;
        T y;
      };

      template <typename T>
      class members<T, 3> {
      protected:
        constexpr static bool variable_sized = false;
        constexpr static bool has_aliases = true;
        ::std::array<T*, 3> m_values;
        members() : m_values{&this->x, &this->y, &this->z} {}
        members(const T& x, const T& y, const T& z) : m_values{&this->x, &this->y, &this->z}, x(x), y(y), z(z) {}
        members(const ::std::initializer_list<T> il) : m_values{&this->x, &this->y, &this->z}, x(il.begin()[0]), y(il.begin()[1]), z(il.begin()[2]) {
          assert(il.size() == 3);
        }

      public:
        T x;
        T y;
        T z;
      };

      template <typename T>
      class members<T, 4> {
      protected:
        constexpr static bool variable_sized = false;
        constexpr static bool has_aliases = true;
        ::std::array<T*, 4> m_values;
        members() : m_values{&this->x, &this->y, &this->z, &this->w} {}
        members(const T& x, const T& y, const T& z, const T& w) : m_values{&this->x, &this->y, &this->z, &this->w}, x(x), y(y), z(z), w(w) {}
        members(const ::std::initializer_list<T> il) : m_values{&this->x, &this->y, &this->z, &this->w}, x(il.begin()[0]), y(il.begin()[1]), z(il.begin()[2]), w(il.begin()[3]) {
          assert(il.size() == 4);
        }

      public:
        T x;
        T y;
        T z;
        T w;
      };

      template <typename T>
      class members<T, ::std::numeric_limits<::std::size_t>::max()> {
      protected:
        constexpr static bool variable_sized = true;
        constexpr static bool has_aliases = false;
        ::std::vector<T> m_values;
        members() = default;
        members(const ::std::size_t n) : m_values(n) {}
        members(const ::std::size_t n, const T& value) : m_values(n, value) {}
        template <typename InputIter>
        members(const InputIter first, const InputIter last) : m_values(first, last) {}
        members(const ::std::initializer_list<T> il) : m_values(il) {}
      };
    }
  }

  template <typename T, ::std::size_t N = ::std::numeric_limits<::std::size_t>::max()>
  class vector : public ::tools::detail::vector::members<T, N> {
  private:
    using Base = ::tools::detail::vector::members<T, N>;
    using F = ::std::conditional_t<::std::is_floating_point_v<T>, T, double>;
    using V = ::tools::vector<T, N>;
    constexpr static bool variable_sized = Base::variable_sized;
    constexpr static bool has_aliases = Base::has_aliases;

  public:
    using reference = T&;
    using const_reference = const T&;
    using size_type = ::std::size_t;
    using difference_type = ::std::ptrdiff_t;
    using pointer = T*;
    using const_pointer = const T*;
    using value_type = T;
    class iterator {
    private:
      V* m_parent;
      size_type m_i;

    public:
      using difference_type = ::std::ptrdiff_t;
      using value_type = T;
      using reference = T&;
      using pointer = T*;
      using iterator_category = ::std::random_access_iterator_tag;

      iterator() = default;
      iterator(const iterator&) = default;
      iterator(iterator&&) = default;
      ~iterator() = default;
      iterator& operator=(const iterator&) = default;
      iterator& operator=(iterator&&) = default;

      iterator(V * const parent, const size_type i) : m_parent(parent), m_i(i) {}

      reference operator*() const {
        return (*this->m_parent)[this->m_i];
      }
      pointer operator->() const {
        return &(*(*this));
      }

      iterator& operator++() {
        ++this->m_i;
        return *this;
      }
      iterator operator++(int) {
        const auto self = *this;
        ++*this;
        return self;
      }
      iterator& operator--() {
        --this->m_i;
        return *this;
      }
      iterator operator--(int) {
        const auto self = *this;
        --*this;
        return self;
      }
      iterator& operator+=(const difference_type n) {
        this->m_i += n;
        return *this;
      }
      iterator& operator-=(const difference_type n) {
        this->m_i -= n;
        return *this;
      }
      friend iterator operator+(const iterator& self, const difference_type n) {
        return iterator(self) += n;
      }
      friend iterator operator+(const difference_type n, const iterator& self) {
        return iterator(self) += n;
      }
      friend iterator operator-(const iterator& self, const difference_type n) {
        return iterator(self) -= n;
      }
      friend difference_type operator-(const iterator& lhs, const iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return static_cast<difference_type>(lhs.m_i) - static_cast<difference_type>(rhs.m_i);
      }
      reference operator[](const difference_type n) const {
        return *(*this + n);
      }

      friend bool operator==(const iterator& lhs, const iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return lhs.m_i == rhs.m_i;
      }
      friend bool operator!=(const iterator& lhs, const iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return lhs.m_i != rhs.m_i;
      }
      friend bool operator<(const iterator& lhs, const iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return lhs.m_i < rhs.m_i;
      }
      friend bool operator<=(const iterator& lhs, const iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return lhs.m_i <= rhs.m_i;
      }
      friend bool operator>(const iterator& lhs, const iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return lhs.m_i > rhs.m_i;
      }
      friend bool operator>=(const iterator& lhs, const iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return lhs.m_i >= rhs.m_i;
      }
    };
    class const_iterator {
    private:
      const V *m_parent;
      size_type m_i;

    public:
      using difference_type = ::std::ptrdiff_t;
      using value_type = T;
      using reference = const T&;
      using pointer = const T*;
      using iterator_category = ::std::random_access_iterator_tag;

      const_iterator() = default;
      const_iterator(const const_iterator&) = default;
      const_iterator(const_iterator&&) = default;
      ~const_iterator() = default;
      const_iterator& operator=(const const_iterator&) = default;
      const_iterator& operator=(const_iterator&&) = default;

      const_iterator(const V * const parent, const size_type i) : m_parent(parent), m_i(i) {}

      reference operator*() const {
        return (*this->m_parent)[this->m_i];
      }
      pointer operator->() const {
        return &(*(*this));
      }

      const_iterator& operator++() {
        ++this->m_i;
        return *this;
      }
      const_iterator operator++(int) {
        const auto self = *this;
        ++*this;
        return self;
      }
      const_iterator& operator--() {
        --this->m_i;
        return *this;
      }
      const_iterator operator--(int) {
        const auto self = *this;
        --*this;
        return self;
      }
      const_iterator& operator+=(const difference_type n) {
        this->m_i += n;
        return *this;
      }
      const_iterator& operator-=(const difference_type n) {
        this->m_i -= n;
        return *this;
      }
      friend const_iterator operator+(const const_iterator& self, const difference_type n) {
        return const_iterator(self) += n;
      }
      friend const_iterator operator+(const difference_type n, const const_iterator& self) {
        return const_iterator(self) += n;
      }
      friend const_iterator operator-(const const_iterator& self, const difference_type n) {
        return const_iterator(self) -= n;
      }
      friend difference_type operator-(const const_iterator& lhs, const const_iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return static_cast<difference_type>(lhs.m_i) - static_cast<difference_type>(rhs.m_i);
      }
      reference operator[](const difference_type n) const {
        return *(*this + n);
      }

      friend bool operator==(const const_iterator& lhs, const const_iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return lhs.m_i == rhs.m_i;
      }
      friend bool operator!=(const const_iterator& lhs, const const_iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return lhs.m_i != rhs.m_i;
      }
      friend bool operator<(const const_iterator& lhs, const const_iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return lhs.m_i < rhs.m_i;
      }
      friend bool operator<=(const const_iterator& lhs, const const_iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return lhs.m_i <= rhs.m_i;
      }
      friend bool operator>(const const_iterator& lhs, const const_iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return lhs.m_i > rhs.m_i;
      }
      friend bool operator>=(const const_iterator& lhs, const const_iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return lhs.m_i >= rhs.m_i;
      }
    };
    using reverse_iterator = ::std::reverse_iterator<iterator>;
    using const_reverse_iterator = ::std::reverse_iterator<const_iterator>;

    vector() = default;
    vector(const V& other) : Base() {
      if constexpr (has_aliases) {
        ::std::copy(other.begin(), other.end(), this->begin());
      } else {
        this->m_values = other.m_values;
      }
    }
    vector(V&& other) noexcept {
      if constexpr (has_aliases) {
        ::std::copy(other.begin(), other.end(), this->begin());
      } else {
        this->m_values = ::std::move(other.m_values);
      }
    }
    ~vector() = default;
    V& operator=(const V& other) {
      if constexpr (has_aliases) {
        ::std::copy(other.begin(), other.end(), this->begin());
      } else {
        this->m_values = other.m_values;
      }
      return *this;
    }
    V& operator=(V&& other) noexcept {
      if constexpr (has_aliases) {
        ::std::copy(other.begin(), other.end(), this->begin());
      } else {
        this->m_values = ::std::move(other.m_values);
      }
      return *this;
    }

    template <bool SFINAE = variable_sized, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    explicit vector(size_type n) : Base(n) {}
    template <bool SFINAE = variable_sized, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    vector(size_type n, const_reference value) : Base(n, value) {}
    template <typename InputIter, bool SFINAE = variable_sized, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    vector(const InputIter first, const InputIter last) : Base(first, last) {}
    template <bool SFINAE = N == 2, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    vector(const T& x, const T& y) : Base(x, y) {}
    template <bool SFINAE = N == 3, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    vector(const T& x, const T& y, const T& z) : Base(x, y, z) {}
    template <bool SFINAE = N == 4, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    vector(const T& x, const T& y, const T& z, const T& w) : Base(x, y, z, w) {}
    vector(const ::std::initializer_list<T> il) : Base(il) {}

    iterator begin() noexcept { return iterator(this, 0); }
    const_iterator begin() const noexcept { return const_iterator(this, 0); }
    const_iterator cbegin() const noexcept { return const_iterator(this, 0); }
    iterator end() noexcept { return iterator(this, this->size()); }
    const_iterator end() const noexcept { return const_iterator(this, this->size()); }
    const_iterator cend() const noexcept { return const_iterator(this, this->size()); }
    reverse_iterator rbegin() noexcept { return ::std::make_reverse_iterator(this->end()); }
    const_reverse_iterator rbegin() const noexcept { return ::std::make_reverse_iterator(this->end()); }
    const_reverse_iterator crbegin() const noexcept { return ::std::make_reverse_iterator(this->cend()); }
    reverse_iterator rend() noexcept { return ::std::make_reverse_iterator(this->begin()); }
    const_reverse_iterator rend() const noexcept { return ::std::make_reverse_iterator(this->begin()); }
    const_reverse_iterator crend() const noexcept { return ::std::make_reverse_iterator(this->cbegin()); }

    size_type size() const noexcept { return this->m_values.size(); }
    bool empty() const noexcept { return this->m_values.empty(); }

    reference operator[](const size_type n) {
      if constexpr (has_aliases) {
        return *this->m_values[n];
      } else {
        return this->m_values[n];
      }
    }
    const_reference operator[](const size_type n) const {
      if constexpr (has_aliases) {
        return *this->m_values[n];
      } else {
        return this->m_values[n];
      }
    }
    reference front() { return *this->begin(); }
    const_reference front() const { return *this->begin(); }
    reference back() { return *this->rbegin(); }
    const_reference back() const { return *this->rbegin(); }

    V operator+() const {
      return *this;
    }
    V operator-() const {
      V res = *this;
      for (auto& v : res) v = -v;
      return res;
    }
    V& operator+=(const V& other) {
      assert(this->size() == other.size());
      for (::std::size_t i = 0; i < this->size(); ++i) {
        (*this)[i] += other[i];
      }
      return *this;
    }
    friend V operator+(const V& lhs, const V& rhs) {
      return V(lhs) += rhs;
    }
    V& operator-=(const V& other) {
      assert(this->size() == other.size());
      for (::std::size_t i = 0; i < this->size(); ++i) {
        (*this)[i] -= other[i];
      }
      return *this;
    }
    friend V operator-(const V& lhs, const V& rhs) {
      return V(lhs) -= rhs;
    }
    V& operator*=(const T& c) {
      for (auto& v : *this) v *= c;
      return *this;
    }
    friend V operator*(const T& lhs, const V& rhs) {
      return V(rhs) *= lhs;
    }
    friend V operator*(const V& lhs, const T& rhs) {
      return V(lhs) *= rhs;
    }
    V& operator/=(const T& c) {
      for (auto& v : *this) v /= c;
      return *this;
    }
    friend V operator/(const V& lhs, const T& rhs) {
      return V(lhs) /= rhs;
    }

    friend bool operator==(const V& lhs, const V& rhs) {
      return ::std::equal(lhs.begin(), lhs.end(), rhs.begin(), rhs.end());
    }
    friend bool operator!=(const V& lhs, const V& rhs) {
      return !(lhs == rhs);
    }

    T inner_product(const V& other) const {
      assert(this->size() == other.size());
      T res{};
      for (::std::size_t i = 0; i < this->size(); ++i) {
        res += (*this)[i] * other[i];
      }
      return res;
    }
    T l1_norm() const {
      T res{};
      for (const auto& v : *this) {
        res += ::tools::abs(v);
      }
      return res;
    }
    T squared_l2_norm() const {
      return this->inner_product(*this);
    }
    F l2_norm() const {
      return ::std::sqrt(static_cast<F>(this->squared_l2_norm()));
    }
    template <bool SFINAE = ::std::is_floating_point_v<T>, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    V normalized() const {
      return *this / this->l2_norm();
    }

    friend ::std::ostream& operator<<(::std::ostream& os, const V& self) {
      os << '(';
      ::std::string delimiter = "";
      for (const auto& v : self) {
        os << delimiter << v;
        delimiter = ", ";
      }
      return os << ')';
    }
    friend ::std::istream& operator>>(::std::istream& is, V& self) {
      for (auto& v : self) {
        is >> v;
      }
      return is;
    }

    template <bool SFINAE = N == 2, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    T outer_product(const V& other) const {
      return this->x * other.y - this->y * other.x;
    }
    template <bool SFINAE = N == 2, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    V turned90() const {
      return V{-this->y, this->x};
    }
    template <bool SFINAE = N == 2, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    V turned270() const {
      return V{this->y, -this->x};
    }

    template <bool SFINAE = N == 2, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    static const ::std::array<V, 4>& four_directions() {
      static const ::std::array<V, 4> res = {
        V{T(1), T(0)},
        V{T(0), T(1)},
        V{T(-1), T(0)},
        V{T(0), T(-1)}
      };
      return res;
    }
    template <bool SFINAE = N == 2, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    static const ::std::array<V, 8>& eight_directions() {
      static const ::std::array<V, 8> res = {
        V{T(1), T(0)},
        V{T(1), T(1)},
        V{T(0), T(1)},
        V{T(-1), T(1)},
        V{T(-1), T(0)},
        V{T(-1), T(-1)},
        V{T(0), T(-1)},
        V{T(1), T(-1)}
      };
      return res;
    }

    template <bool SFINAE = N == 3, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    V outer_product(const V& other) const {
      return V{this->y * other.z - this->z * other.y, this->z * other.x - this->x * other.z, this->x * other.y - this->y * other.x};
    }
    template <bool SFINAE = N == 3 && ::std::is_floating_point_v<T>, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    ::std::array<V, 3> orthonormal_basis() const {
      assert((*this != V{0, 0, 0}));

      ::std::array<V, 3> v;
      v[0] = *this;
      v[1] = V{0, this->z, -this->y};
      if (v[1] == V{0, 0, 0}) {
        v[1] = V{-this->z, 0, this->x};
      }
      v[1] -= v[0].inner_product(v[1]) / v[0].inner_product(v[0]) * v[0];

      v[0] = v[0].normalized();
      v[1] = v[1].normalized();
      v[2] = v[0].outer_product(v[1]);

      return v;
    }
  };
}

namespace std {
  template <typename T>
  struct hash<::tools::vector<T, 2>> {
    using result_type = ::std::size_t;
    using argument_type = ::tools::vector<T, 2>;
    result_type operator()(const argument_type& key) const {
      static const ::tools::tuple_hash<T, T> hasher;
      return hasher(::std::make_tuple(key.x, key.y));
    }
  };
  template <typename T>
  struct hash<::tools::vector<T, 3>> {
    using result_type = ::std::size_t;
    using argument_type = ::tools::vector<T, 3>;
    result_type operator()(const argument_type& key) const {
      static const ::tools::tuple_hash<T, T, T> hasher;
      return hasher(::std::make_tuple(key.x, key.y, key.z));
    }
  };
  template <typename T>
  struct hash<::tools::vector<T, 4>> {
    using result_type = ::std::size_t;
    using argument_type = ::tools::vector<T, 4>;
    result_type operator()(const argument_type& key) const {
      static const ::tools::tuple_hash<T, T, T, T> hasher;
      return hasher(::std::make_tuple(key.x, key.y, key.z, key.w));
    }
  };
}


#line 5 "tools/vector2.hpp"

namespace tools {
  template <typename T>
  using vector2 = ::tools::vector<T, 2>;
}


#line 1 "tools/less_by.hpp"



namespace tools {

  template <class F>
  class less_by {
  private:
    F selector;

  public:
    less_by(const F& selector) : selector(selector) {
    }

    template <class T>
    bool operator()(const T& x, const T& y) const {
      return selector(x) < selector(y);
    }
  };
}


#line 1 "tools/ccw.hpp"



#line 5 "tools/ccw.hpp"

namespace tools {
  template <typename T>
  int ccw(const ::tools::vector2<T>& a, ::tools::vector2<T> b, ::tools::vector2<T> c) {
    b -= a;
    c -= a;
    if (b.outer_product(c) > T(0)) return +1;
    if (b.outer_product(c) < T(0)) return -1;
    if (b.inner_product(c) < T(0)) return +2;
    if (b.squared_l2_norm() < c.squared_l2_norm()) return -2;
    return 0;
  }
}


#line 15 "tools/convex_hull.hpp"

namespace tools {
  template <typename InputIterator, typename OutputIterator>
  void convex_hull(const InputIterator begin, const InputIterator end, bool minimum, OutputIterator result) {
    using T = ::std::decay_t<decltype(begin->x)>;

    const ::std::vector<::tools::vector2<T>> v(begin, end);
    ::std::vector<::std::size_t> a(v.size());
    ::std::iota(a.begin(), a.end(), 0);
    ::std::sort(a.begin(), a.end(), ::tools::less_by([&](const ::std::size_t i) {
      return ::std::make_pair(v[i].x, v[i].y);
    }));
    ::std::vector<::std::vector<::std::size_t>> duplicates;

    if (minimum) {
      ::std::size_t vl = 0;
      for (::std::size_t vr = 0, al = 0, ar = 0; al < a.size(); vl = vr, al = ar) {
        for (; ar < a.size() && v[a[al]].x == v[a[ar]].x; ++vr, ++ar);
        if (vl < al) ::std::move(::std::next(a.begin(), al), ::std::next(a.begin(), ar), ::std::next(a.begin(), vl));
        if (v[a[vl]].y == v[a[vr - 1]].y) {
          vr -= vr - vl - 1;
          duplicates.emplace_back();
          duplicates.back().push_back(a[vl]);
        } else {
          ::std::swap(a[vl + 1], a[vr - 1]);
          vr -= vr - vl - 2;
          duplicates.emplace_back();
          duplicates.back().push_back(a[vl]);
          duplicates.emplace_back();
          duplicates.back().push_back(a[vl + 1]);
        }
      }
      a.erase(::std::next(a.begin(), vl), a.end());
    } else {
      ::std::size_t vl = 0;
      for (::std::size_t vr = 0, al = 0, ar = 0; al < a.size(); vl = vr, al = ar) {
        for (; ar < a.size() && v[a[al]] == v[a[ar]]; ++vr, ++ar);
        if (vl < al) ::std::move(::std::next(a.begin(), al), ::std::next(a.begin(), ar), ::std::next(a.begin(), vl));
        duplicates.emplace_back();
        for (::std::size_t i = vl; i < vr; ++i) {
          duplicates.back().push_back(a[i]);
        }
        vr -= vr - vl - 1;
      }
      a.erase(::std::next(a.begin(), vl), a.end());
    }

    ::std::vector<::std::size_t> convex_hull;
    if (a.size() >= 3) {

      convex_hull.push_back(0);
      convex_hull.push_back(1);
      for (::std::size_t p3 = 2; p3 < a.size(); ++p3) {
        while (convex_hull.size() >= 2) {
          const int ccw = ::tools::ccw(v[a[convex_hull.rbegin()[1]]], v[a[convex_hull.back()]], v[a[p3]]);
          if (ccw == 1 || (!minimum && ccw == -2)) {
            break;
          }
          convex_hull.pop_back();
        }
        convex_hull.push_back(p3);
      }

      const ::std::size_t threshold = convex_hull.size() + 1;
      for (::std::size_t p3 = convex_hull.back(); p3 --> 0;) {
        while (convex_hull.size() >= threshold) {
          const int ccw = ::tools::ccw(v[a[convex_hull.rbegin()[1]]], v[a[convex_hull.back()]], v[a[p3]]);
          if (ccw == 1 || (!minimum && ccw == -2)) {
            break;
          }
          convex_hull.pop_back();
        }
        convex_hull.push_back(p3);
      }
      convex_hull.pop_back();

    } else {
      for (::std::size_t i = 0; i < a.size(); ++i) {
        convex_hull.push_back(i);
      }
    }

    for (const ::std::size_t& c : convex_hull) {
      for (const ::std::size_t& i : duplicates[c]) {
        if constexpr (::std::is_assignable_v<OutputIterator, ::std::size_t>) {
          *result = i;
        } else {
          *result = v[i];
        }
        ++result;
      }
    }
  }
}

Back to top page