convex_hull.hpp

//convex_hull.hpp

#pragma once

#include <iostream>
#include <fstream>
#include <vector>
#include <algorithm>
#include <math.h>

// Определения типов, с которыми работают функции:
typedef int TypeOfPoint; // Тип точки
typedef std::pair<TypeOfPoint, TypeOfPoint> Point; // Точка
typedef std::vector<Point> Points; // Вектор из точек

// Создание выпуклой оболочки из файла с точками:
void create_convex_hull_from_file(Points & hull, const std::string & filename);

// Создание выпуклой оболочки из вектора точек:
void create_convex_hull_from_points(Points & hull, const Points & points);

// Добавление точки в выпуклую оболочку:
void add_to_convex_hull(Points & hull, const Point & p); // Добавление точки в произвольную выпуклую оболочку
void add_to_zero_convex_hull(Points & hull, const Point & p); // В пустую выпуклую оболочку
void add_to_one_convex_hull(Points & hull, const Point & p); // В выпуклую оболочку, состоящую из одной точки
void add_to_two_convex_hull(Points & hull, const Point & p); // В выпуклую оболочку, состоящую из двух точек
void add_to_other_convex_hull(Points & hull, const Point & p); // В выпуклую оболочку, состоящую из более, чем двух точек

// Модифицирование выпуклой оболочки:
void shift_points_left(Points & hull, int k); // Циклический сдвиг всех вершин влево на k элементов
void shift_points_right(Points & hull, int k); // Циклический сдвиг всех вершин вправо на k элементов
void erase_from_points(Points & hull, int a, int b); // Вырезать из выпуклой оболочки точки начиная с индекса a и заканчивая индексом b

// Вывод вектора из точек:
void print_points(const Points & points); // Вывод на экран
void printf_points(const Points & points, const std::string & filename); // Вывод в файл

// Работа с числами в кольце по модулю n:
void to_n(int & val, int n); // Приводит значение val в диапазон значений [0, .., n-1]
void inc_n(int & val, int n, int k = 1); // Увеличивает значение val на k (по умолчанию: k=1) и вычисляет новое значение в кольце по модулю n
void dec_n(int & val, int n, int k = 1); // Уменьшает значение val на k (по умолчанию: k=1) и вычисляет новое значение в кольце по модулю n

// Работа с отдельными точками:
double dist(const Point & A, const Point & B); // Вычисление расстояния между точками A и B
double det(const Point & A, const Point & B); // Вычисление определителя, первая строка которого состоит из координат точки A, вторая из координат точки B
double square_triangle(const Point & A, const Point & B, const Point & C); // Вычисление ориентированной площади треугольника с вершинами A, B и C
bool in_interval(const Point & A, const Point & B, const Point & C); // Определение, лежит ли точка C на отрезке AB


void create_convex_hull_from_file(Points & hull, const std::string & filename){
//  Создание выпуклой оболочки из файла. Точки поступают последовательно
    if (!hull.empty()) hull.clear();

    std::ifstream fin(filename); // Открываем файл для чления
    if (!fin) return; // Если файл не открыт, выходим

    TypeOfPoint x, y;
    while(fin>>x>>y)
        add_to_convex_hull(hull, Point(x, y));
}

void create_convex_hull_from_points(Points & hull, const Points & points){
//  Создание выпуклой оболочки из вектора точек на плоскости

//  Очистим оболочку, если она не пуста:
    if (!hull.empty()) hull.clear();

//  По одной точке добавляем в выпуклую оболочку:
    for(int i = 0, n = points.size(); i < n; i++)
        add_to_convex_hull(hull, points[i]);
}

void add_to_convex_hull(Points & hull, const Point & p){
//  Добавление точки к уже существующей выпуклой оболочке

//  Выбираем случай, исходя из размера оболочки:
    switch(hull.size()){
        case 0:  add_to_zero_convex_hull(hull, p);  break; // Случай 1: оболочка пуста
        case 1:  add_to_one_convex_hull(hull, p);   break; // Случай 2: оболочка состоит из одной точки
        case 2:  add_to_two_convex_hull(hull, p);   break; // Случай 3: оболочка состоит из двух точек
        default: add_to_other_convex_hull(hull, p); break; // Случай 4: оболочка состоит более чем из двух точек
    }
}

void add_to_zero_convex_hull(Points & hull, const Point & p){
//  Добавление указателя на точку к пустой выпуклой оболочке
    hull.push_back(p);
}

void add_to_one_convex_hull(Points & hull, const Point & p){
//  Добавление указателя на точку к выпуклой оболочке, состоящей из одной точки

//  Сравниваем уже имеющуюся точку с новой точкой. Если не совпадают, то добавляем
    if (hull[0] != p) hull.push_back(p);
}

void add_to_two_convex_hull(Points & hull, const Point & p){
//  Добавление указателя на точку к выпуклой оболочке, состоящей из двух точек

//  Считаем ориентированную площадь треугольника, состоящего из точек в оболочке и новой точки
    double sq = square_triangle(hull[0], hull[1], p);

//  Если площадь равна нулю, то проверяем, лежит ли новая точка на отрезке, соединяющем две точки выпуклой оболочки
    if (sq == 0 ){
        if (!in_interval(hull[0], hull[1], p)) // Если новая точка не лежит на отрезке, но ориентированная площадь равна нулю
            if (dist(hull[0], p) > dist(hull[1], p)) // Если расстояние от новой точки до нулевой больше, чем от новой до первой
                hull[1] = p; // То заменяем первую точку на новую
            else // Иначе
                hull[0] = p; // Заменяем вторую точку на новую
        return; // Выходим из функции
    }

//  Добавляем точку в конец
    hull.push_back(p);

//  Если значение ориентированной площади отрицательно, то меняем местами точку вначале и вконце.
//  Это позволяет сохранить направление обхода против часовой стрелки.
    if (sq < 0) swap(hull[0], hull[2]);
}

void add_to_other_convex_hull(Points & hull, const Point & p){
//  Добавляет новую точку к выпуклой оболочке, состоящей больше, чем из двух точек

//  Создаем два отрезка для хранения цепей в оболочке, которые необходимо модифицировать
    std::pair<int, int> interval[] = {std::pair<int, int>(-1, -1), std::pair<int, int>(-1, -1)};
    int n = hull.size(); // Размер оболочки
    int cur = 0; // Индекс текущего отрезка
    for(int i = 0; i < n; i++){
//      Вычислим индекс j точки, следующей сразу за текущей точкой:
        int j = i + 1;
        to_n(j, n);

//      Посчитаем ориентированную площадь треугольника, образованного точками hull[i], p, hull[j]
        double sq = square_triangle(hull[i], p, hull[j]);
        if (sq == 0 && in_interval(hull[i], hull[j], p)) // Если она равна нулю и точка лежит на отрезке между точками hull[i] и hull[j]
            return; // То выходим из функции

//      С этого момента мы считаем, что отрезок [i, j] выпуклой оболочки необходимо модифицировать

        if (sq >= 0){ // Если значение площади не отрицательно
//          Если еще не задан текущий отрезок, который нам нужно модифицировать в выпуклой оболочке
            if (interval[cur].first == interval[cur].second)
                interval[cur] = std::pair<int, int>(i, j); // То задаем его
            else if (interval[cur].second == i) // Иначе, если конец текущего отрезка совпадает с началом полученного отрезка [i, j]
                interval[cur].second = j; // Концом текущего отрезка становится конец полученного (то есть к концу текущего отрезка прибавляется полученный)
            else if (interval[cur].first == j) // Иначе, если начало текущего отрезка совпадает с концом полученного отрезка [i, j]
                interval[cur].first = i; // Началом текущего отрезка становится начало полученнго (то есть к началу текущего отрезка прибавляется полученный)
            else // Иначе, полученный отрезок [i, j] становится вторым отрезком, который необходимо модифицировать, и сразу же текущим
                interval[++cur] = std::pair<int, int>(i, j);
        }
    }

//  На выходе мы имеем два отрезка, которые необходимо соединить в один, если совпадают конец второго и начало первого:
    if (interval[1].second == interval[0].first)
        interval[0].first = interval[1].first;

//  Теперь мы будем рассматривать один целый отрезок [begin, end]
    int begin = interval[0].first;
    int end = interval[0].second;

//  Вычисляем его длину по модулю n:
    int len = (end - begin);
    to_n(len, n);

    if (begin > end && end != 0) { // Если начало отрезка больше конца
        shift_points_left(hull, end); // То сдвигаем все элементы оболочки влево на величину, равную концу отрезка
        begin -= end;
        end -= end;
    }

//  Крайние точки мы должны оставить в покое, поэтому:
    inc_n(begin, n); // Увеличиваем на единицу начало отрезка
    dec_n(end, n); // Уменьшаем на единицу конец отрезка

    switch (len){ // В зависимости от длины отрезка
        case 0: return; // Ничего не делаем, если длина отрезка равна 0
        case 1: hull.emplace(hull.begin()+begin, p); break; // Вставляем новую точку, если длина отрезка равна 1
        case 2: hull[begin] = p; break; // Заменяем точку, расположенную вначале отрезка, на новую точку
        default:
            hull[begin] = p; // Заменяем точку вначале отрезка на новую точку
            inc_n(begin, n); // Увеличиваем начало отрезка
            erase_from_points(hull, begin, end); // Вырезаем все остальные точки из отрезка [begin, end]
            break;
    }
}

void shift_points_left(Points & points, int k){
//  Сдвигает все элементы выпуклой оболочки на k позиций влево
    reverse(points.begin(), points.begin()+k-1);
    reverse(points.begin()+k, points.end());
    reverse(points.begin(), points.end());
}

void shift_points_right(Points & points, int k){
//  Сдвигает все элементы выпуклой оболочки на k позиций вправо
    shift_points_left(points, points.size()-k+1);
}

void erase_from_points(Points & points, int a, int b){
//  Вырезает из выпуклой оболочки точки с позиции a до позиции b включительно
    int n = points.size();

    to_n(a, n); // Вычисляем a по модулю n
    to_n(b, n); // Вычисляем b по модулю n

    if (a > b) return; // Если левый конец отрезка больше правого, то выходим из функции

    if (a == b) // Если начало совпадает с концом, то нам нужно вырезать только один элемент
        points.erase(points.begin()+a);
    else // Иначе вырезаем весь отрезок
        points.erase(points.begin()+a, points.begin()+b+1);
}

void printf_points(const Points & points, const std::string & filename){
//  Вывод в файл выпуклой оболочки, состоящей из указателей на точки
    std::ofstream fout(filename);
    for(int i = 0, n = points.size(); i < n; i++)
        fout<<points[i].first<<" "<<points[i].second<<std::endl;
}

void print_points(const Points & points){
//  Вывод на экран вектора из точек
    for(int i = 0, n = points.size(); i < n; i++)
        std::cout<<"("<<points[i].first<<";"<<points[i].second<<") ";
    std::cout<<std::endl;
}

void to_n(int & val, int n){
//  Функция вычисляет значение val по модулю n
    if (val >= n)
        val %= n;
    else if (val < 0)
        val = n - abs(val) % n;
}

void inc_n(int & val, int n, int k){
//  Функция увеличивает значение val, а затем вычисляет его по модулю n
    val += k;
    to_n(val, n);
}

void dec_n(int & val, int n, int k){
//  Функция уменьшает значение val, а затем вычисляет его по модулю n
    val -= k;
    to_n(val, n);
}

double det(const Point & A, const Point & B){
//  Подсчет определителя, первой строкой которого идут координаты точки a, второй строкой - координаты точки b
    return (A.first * B.second - B.first * A.second);
}

double square_triangle(const Point & A, const Point & B, const Point & C){
//  Подсчет ориентированной площади треугольника, состоящего из вершин a, b и c
    return ((det(A,B)+det(B,C)+det(C,A))/2.0);
}

double dist(const Point & A, const Point & B){
//  Вычисляет расстояние от точки A до точки B

//  Вычтем точку A из точки B:
    Point C(B.first - A.first, B.second - A.second);

//  Возводим каждую компоненту в квадрат, суммируем и извлекаем корень:
    return (sqrt(C.first*C.first+C.second*C.second));
}

bool in_interval(const Point & A, const Point & B, const Point & C){
//  Проверяет, принадлежит ли точка C отрезку AB
//  Точка C точно не принадлежит этому отрезку, если сумма расстояний от его концов до точки C больше длины отрезка

//  Вычисляем длины отрезков AC, BC и AB:
    double AC = dist(A, C);
    double BC = dist(B, C);
    double AB = dist(A, B);

//  Проверяем условие принадлежности:
    return (AC + BC == AB);
}