Как вы можете практически решить выпуклый корпус, когда есть проблемы с плавающей точностью?

Предположим, что у вас есть 100000 точек на кривой y = x^2, Вы хотите найти выпуклую оболочку этих точек. Все координаты являются плавающими числами.

В моей реализации сканирования Грэма единственное место, где я оперирую плавающими числами, — это когда я сначала сортирую все точки по их координатам, а затем у меня есть одна функция, которая определяет, поворачивают ли три точки влево или вправо.

Точки:

struct point {
double x;
double y;
};

Сортировочный компаратор:

inline bool operator() (const point &p1, const point &p2) {
return (p1.x < p2.x) || (p1.x == p2.x && p1.y > p2.y);
}

Левый / правый поворот:

inline int ccw(point *p1, point *p2, point *p3) {
double left  = (p1->x - p3->x)*(p2->y - p3->y);
double right = (p1->y - p3->y)*(p2->x - p3->x);
double res = left - right;
return res > 0;
}

Моя программа говорит, что из 100 000 точек только 68894 являются частью выпуклой оболочки. Но поскольку они находятся на кривой, все они должны быть частью выпуклой оболочки.

Для твоего глаза это не имеет никакого значения. Смотрите рисунок ниже. Красные точки являются частью выпуклой оболочки.

изображение http://oi57.tinypic.com/t8lsvs.jpg

Но если вы посмотрите достаточно близко и приблизитесь к точкам, вы увидите, что некоторые из них синего цвета, поэтому они не входят в выпуклый корпус.

изображение http://oi61.tinypic.com/2eol37a.jpg

Теперь мое первоначальное предположение состоит в том, что ошибки с плавающей точкой вызывают эту проблему.

Я думаю, я мог бы использовать внешнюю библиотеку, которая имеет произвольную точность для чисел с плавающей запятой, но меня больше интересуют простые типы данных, которые мы имеем, например, в C ++.

Как я мог увеличить точность? Я читал о эпсилонах, но как использование эпсилон здесь поможет? Я бы по-прежнему предполагал, что некоторые точки, близкие друг к другу, одинаковы, поэтому точность не приблизится к 100%.

Как лучше всего подойти к этой проблеме?

2

Решение

Вы правы, что все точки должны быть на выпуклой оболочке, если вы действительно используете точки вида (x, x^2), Однако три точки могут быть коллинеарными. Если вы перемещаете их или делаете что-то еще странное, это уходит в окно.

Если вы выберете 100000 баллов, я бы предложил использовать целые числа в [-50000,49999]. Ваш ccw функция будет вычислять left а также right быть целыми числами меньше по абсолютной величине, чем 2.5e14 < 2 ^ 53, поэтому округление не произойдет.

Ваша сортировка на основе координат будет работать правильно независимо от ввода.

Для общих входных данных, следующее ccw предикат глючит:

inline int ccw(point *p1, point *p2, point *p3) {
double left  = (p1->x - p3->x)*(p2->y - p3->y);
double right = (p1->y - p3->y)*(p2->x - p3->x);
double res = left - right;
return res > 0;
}

Может быть округление как в вычитаниях, так и в умножениях. Если все ваши точки лежат в ограничительной рамке H * W, различия по координате x будут вычислены с абсолютной ошибкой около H * eps / 2, а различия по координате y будут вычислены с абсолютной ошибкой вокруг W * eps / 2. Поэтому продукты будут вычисляться с абсолютной погрешностью около H * W * eps / 2. Если fabs(left - right) < 3*H*W*eps/2нужно оценить left а также right точнее. eps здесь 2-52.

Я, вероятно, рекомендую использовать MPFR, если double Сравнение ничего вам не говорит. Вы можете сделать это без, однако. Трюк с суммированием Кахана поможет вам получить небольшие разницы от различий, а 227+1 трюк может помочь вам точно вычислить продукты.

0

Другие решения

Очень часто в математике с плавающей точкой нужно вводить понятие «толерантности», иногда обозначаемое как эпсилон. В вашем случае вы могли бы сделать свой ccw() функция трехзначная: истина / ложь / неопределенная. Затем, когда вы пытаетесь выяснить, может ли новая точка быть частью выпуклой оболочки, вы спрашиваете «является ли она ccw = true или неопределенным», и в любом случае вы принимаете точку. Неопределенный результат может произойти, когда наклон слишком близок к прямой линии, чтобы принять решение.

0