У меня есть функция сравнения для сравнения двух точек в 2 (или более) измерениях на основе угла каждой точки в стандартных полярных координатах. Другими словами, точка a следует сравнивать меньше точки b если угол, проходящий через начало координат от оси + x к a меньше соответствующего угла для b. Для начала координат нормализованный угол равен 0 НО для здравомыслия сравнивается с любой точкой.
Это можно сделать довольно просто, используя atan2 функция, которая имеет встроенную обработку для разных квадрантов и просто возвращает нормализованный угол (от 0 включительно до 2 пи не считая) для любой заданной точки, кроме начала координат. К сожалению, триггер медленный, поэтому лучше сравнить b.x*a.y <=> a.x*b.y. Это также хорошо, так как с целыми числами он работает лучше, чем триггерный. Это сравнение может быть получено из того факта, что tan строго возрастает на интервале от 0 до pi / 2. Однако это работает, только если точки a а также b находятся в том же квадранте, поэтому нам нужно провести некоторую предварительную обработку, прежде чем мы сможем ее применить.
Здесь и возникает моя дилемма с моим текущим кодом. Изначально у меня был довольно длинный и трудный для чтения лес if заявления, которые обрабатывали все случаи, когда a а также b были не в том же квадранте, но я заменил его кодом, который вы видите здесь. Сначала я получаю «номер региона» для каждой точки. a а также b, а затем используйте b.x*a.y <=> a.x*b.y сравнение при необходимости. Я определяю номер региона как 0 для происхождения и 1-8 для оси + x через квадрант IV против часовой стрелки. А именно, 2 это квадрант I, 3 ось + y, 4 это квадрант II, 5 ось -x, 6 — квадрант III, и 7 ось -y.
Таким образом, я могу сравнить две точки по
- получение номера региона для каждого
- если любой из них
0, возвращаться0так как одна из точек является началом координат, поэтому они должны сравнивать равные, поскольку - в противном случае, если их номера регионов не совпадают, вернуть
-1если номер региона дляaменьше и1если номер региона дляbменьше - в противном случае номера регионов для
aа такжеbодинаковы. Если номер региона дляaстранно, возврат0потому что точки находятся на одной оси в одном направлении (+ x, + y, -x или -y) - иначе,
aа такжеbоба находятся внутри одного и того же квадранта, поэтому мы можем использоватьb.x*a.y <=> a.x*b.yсравнение
Нахождение номера региона для каждой точки в некоторых случаях делает некоторую избыточную работу, но, вероятно, это то, что компилятор может оптимизировать. Я не на том этапе, когда мне нужно профилировать и оптимизировать это, хотя я хочу избежать вычислений с плавающей запятой.
Вот код функции сравнения:
/// Return an identifier for the region of the xy plane containing a point
///
/// Returns 0 for origin or 1 to 8 for the +/- axis and quadrants in
/// counterclockwise order starting from the +x axis (ie 1 for +x axis,
/// 2 for quadrant I, etc, and 8 for quadrant IV)
static inline int find_x_y_region(const int64_t *a){
if(a[0] == 0){// a is on the y axis
if(a[1] == 0){// a is the origin
return 0;
}
return a[1] > 0 ? 3 : 7;// +y axis and -y axis respectively
}else if(a[1] == 0){// a is on the x axis (but not the origin)
return a[0] > 0 ? 1 : 5;// +x axis and -x axis respectively
}else if(a[1] > 0){// a is in the upper half plane
return a[0] > 0 ? 2 : 4;// quadrant I and quadrant II respectively
}// otherwise a is in the lower half plane
return a[0] > 0 ? 8 : 6;// quadrant IV and quadrant III respectively
}
static inline int cmp_x_y(const void *_a, const void *_b){
const int64_t *a = _a, *b = _b;
// consider the points a and b projected into the xy plane (ie ignore their z components)
// we want to find which point, a or b, has a smaller angle in its standard polar representation
// return -1 if a has a smaller angle, 1 if b does, 0 if they have the same angle
// we do not need to do trig for this, but there are a lot of cases
// for a fixed z, all xy points have a fixed radius, so comparison of points at the same z
// can be done instantly with one coordinate. However, when a and b have different z coordinates we
// need to work harder.
int a_region = find_x_y_region(a), b_region = find_x_y_region(b);
if(a_region == 0 || b_region == 0){
return 0;
}else if(a_region < b_region){
return -1;
}else if(a_region > b_region){
return 1;
}else if(a_region & 1){// both regions are the same, if they are odd both points are on the same axis (+/- x/y)
return 0;
}// otherwise both points are in the same quadrant
int64_t ord = b[0]*a[1] - a[0]*b[1];
if(ord < 0){
return -1;
}else if(ord == 0){
return 0;
}
return 1;
}
Факторинг find_x_y_region очень помогло, но могу ли я еще больше упростить эту функцию сравнения? Кроме того, следует ли мне беспокоиться о переполнении при вычислении ord? Как указано в комментариях, я собираюсь использовать это, чтобы сделать дерево KD из точек на поверхности сферы, поэтому я уже требую, чтобы квадрат любой координаты был представлен в виде int64_t и для таких точек переполнение не проблема. Я мог бы объединить пару ветвей в функции сравнения с тернарными операторами, которые сократят ее на 4 строки, но это не улучшит читаемость.
2 ответа
должен ли я беспокоиться о переполнении при вычислении ord?
Да. Разница может вылиться через край.
Вместо этого сравните:
// Compute cross product sign (see next bullet)
int64_t ba = b[0]*a[1];
int64_t ab = a[0]*b[1];
if(ba < ab) return -1;
if(ba == ab) return 0;
return 1;
могу ли я еще больше упростить эту функцию сравнения?
Сделано чрезмерное сравнение.
Чтобы расширить идею @vnp:
… Если оба вектора находятся в верхней полуплоскости, перекрестное произведение работает. Если они оба находятся в нижней полуплоскости, тоже работает (в обратном направлении). Если они находятся в разных полуплоскостях, то в верхней меньше.
int cmp_x_y(const int64_t a[2], const int64_t b[2]) {
if (a[1] < 0 == b[1] < 0) { // same upper/lower plane?
int cps = cross_product_sign(a,b);
if (a[1] < 0) cps = -cps;
return cps;
}
return a[1] < 0 ? -1 : 1;
}
Вышеупомянутое может потребовать некоторого обзора крайнего случая, но в нем представлена идея использования перекрестного произведения в качестве общего пути решения.
void * против. int64_t *
*a = _a приводит к неопределенное поведение (UB) следует *_a не соответствуют требованиям выравнивания int64_t. void * также теряет некоторую проверку типов.
Предложите простое использование int64_t * или же int64_t a[2].
//static inline int cmp_x_y(const void *_a, const void *_b){
// const int64_t *a = _a, *b = _b;
static inline int cmp_x_y(const int64_t *a, const int64_t *b) {
// or
static inline int cmp_x_y(const int64_t a[2], const int64_t b[2]) {
Незначительное: ошибка обсуждения
Это можно сделать довольно тривиально, используя функцию atan2, которая имеет встроенную обработку для разных квадрантов и просто возвращает нормализованный угол (от 0 включительно до 2 пи исключая).
Выше неверно: atan2(y, x) возвращает «arctan y = x в интервале [-π π] радианы «.
Раннее возвращение
У вас их куча. Я за них (не все), но вы непоследовательны в своих else стиль после скорейшего возвращения. Иногда вы делаете это:
if(a[1] == 0){
return a[0] > 0 ? 1 : 5;
}else // ...
и иногда вы обходитесь без else:
if(a[1] > 0){
return a[0] > 0 ? 2 : 4;
}
return a[0] > 0 ? 8 : 6;
Я предпочитаю второй стиль, так как else является избыточным.
В соответствии
Даже если компилятор уважает ваше явное inline (чего обычно не происходит), вы бы этого не захотели. Компиляторы лучше людей решают, что следует автоматически вставлять, а что нет. Просто брось inline.
Дальнейшее чтение, например: https://www.greenend.org.uk/rjk/tech/inline.html, особенно
Иногда компилятору необходимо создать автономную копию объектного кода для функции, даже если это встроенная функция.
[…]
У вас может быть отдельное (не встроенное) определение в другой единице перевода, и компилятор может выбрать либо это, либо встроенное определение.
а также https://stackoverflow.com/questions/2082551/what-does-inline-mean, в том числе
в настоящее время большинство компиляторов даже не считают это подсказкой для встраивания функции.