Алгоритм прыгающего пузыря для наименьшей объемлющей сферы

Я понимаю, что этому посту уже год, но я реализовал алгоритм прыгающего пузыря на C++ и подумал, что, возможно, кому-то он покажется полезным. Это основано на статье Тянь Бо, которую я купил за 18 долларов.

BoundSphere calculateBoundSphere(vector<Vertex> vertices){
    Vector3D center = vertices[0].position;
    float radius = 0.0001f;
    Vector3D pos, diff;
    float len, alpha, alphaSq;
    vector<Vertex>::iterator it;

    for (int i = 0; i < 2; i++){
        for (it = vertices.begin(); it != vertices.end(); it++){
            pos = it->position;
            diff = pos - center;
            len = diff.length();
            if (len > radius){
                alpha = len / radius;
                alphaSq = alpha * alpha;
                radius = 0.5f * (alpha + 1 / alpha) * radius;
                center = 0.5f * ((1 + 1 / alphaSq) * center + (1 - 1 / alphaSq) * pos);
            }
        }
    }

    for (it = vertices.begin(); it != vertices.end(); it++){
        pos = it->position;
        diff = pos - center;
        len = diff.length();
        if (len > radius){
            radius = (radius + len) / 2.0f;
            center = center + ((len - radius) / len * diff);
        }
    }

    return BoundSphere(center, radius);
}

Вот реализация, о которой вы просили. Что касается аналитической части, вы можете получить представление, прочитав следующее.

Этот анализ предназначен для круга (2d), но расширение до сферы (3d) должно быть легким. Я уверен, что анализ поможет.

Алгоритм времени O(n2)

• Нарисуйте окружность в центре с так, чтобы точки данного множества лежали внутри окружности. Ясно, что этот круг можно уменьшить (иначе у нас есть решение).
• Уменьшите круг, найдя точку A, наиболее удаленную от центра круга, и нарисуйте новый круг с тем же центром, проходящий через точку A. Эти два шага создают меньший окружающий круг. Причина того, что новый круг меньше, заключается в том, что этот новый круг по-прежнему содержит все точки заданного набора, за исключением того, что теперь он проходит через самую дальнюю точку x, а не охватывает ее.
• Если окружность проходит через 2 или более точек, перейдите к шагу 4. В противном случае уменьшите окружность, переместив центр в сторону точки А, пока окружность не коснется другой точки В из набора.
• At this stage, we a circle, C, that passes through two or more points of a given set. If the circle contains an interval (point-free interval) of arc greater than half the circle's circumference on which no points lie, the circle can be made smaller. Let D and E be the points at the ends of this point-free interval. While keeping D and E on the circle's boundary, reduce the diameter of the circle until we have either case (a) or case (b).
  • Case (a) The diameter is the distance DE.
    We are done!
  • Случай (b) Окружность C касается другой точки из множества F.
    Проверить, не выходят ли бесточечные дуговые интервалы длиной более половины длины окружности C.
    ЕСЛИ
    нет таких бесточечных дуговые интервалы exit THEN Готово!
    Иначе
    Перейти к шагу 4. В этом случае три точки должны лежать на дуге длиной менее половины окружности. Повторяем шаг 4 для двух внешних из трех точек дуги.

Анализ

Этот алгоритм прост, но дорог в реализации. Шаги 1, 2 и 3 требуют линейного времени по количеству точек в данном наборе. На шаге 4 выше требуется линейное время, чтобы найти каждую новую точку F. Однако нахождение точки F не гарантирует завершение алгоритма. Шаг 4 необходимо повторять до тех пор, пока в окружности не останется ни одного бесточечного интервала длиннее половины ее окружности. В худшем случае для этого требуется (n-2) итераций шага 4, что означает, что общее время, затраченное на шаг 4, может быть порядка квадрата размера набора точек.

Следовательно, этот алгоритм равен O(n2).

Алгоритм O(n)

Нимрод Мегиддо предлагает алгоритм, и он использует методы сокращения и поиска для линейного программирования, чтобы найти минимальную охватывающую окружность за время O (n).

Сокращение и поиск

Суть алгоритма Мегиддо заключается в сокращении и поиске. В алгоритме сокращения и поиска на каждом шаге выполняется линейный объем работы, чтобы уменьшить размер входных данных на постоянную долю f. Если это может быть достигнуто, то общий объем проделанной работы уменьшится до O(n)*(1 + (1-f) + (1-f)2 + ...). В этой формуле бесконечный ряд является геометрическим и в сумме дает постоянное значение, поэтому общее время работы равно O(n).

Например, предположим, что, проверив наш набор из n входных элементов, мы можем отбросить 1/4 из них как не относящиеся к решению. Повторно применяя эту проверку к оставшимся элементам, мы можем уменьшить входные данные до размера, который несложно решить, скажем, n≤3. Общее время, необходимое для достижения этого сокращения, будет пропорционально (n + 3n/4 + 9n/16 + ...). Легко показать, что этот ряд приближается к пределу 4n и никогда не превышает его. Следовательно, общее время работы пропорционально n, что и требуется.

Идея использования геометрического ряда для сведения алгоритма к линейному времени возникла еще до работы Мегиддо; в частности, ранее он использовался для разработки алгоритмов нахождения медианы O (n). Однако он был первым, кто применил его к ряду фундаментальных задач вычислительной геометрии.

Алгоритм линейного времени Мегиддо

Чтобы найти минимальную охватывающую окружность (MEC) набора точек, алгоритм Мегиддо отбрасывает не менее n/16 точек на каждой (линейной) итерации. То есть для заданного набора S из n точек алгоритм идентифицирует n/16 точек, которые можно удалить из S, не влияя на MEC S. Эту процедуру можно многократно применять до тех пор, пока не будет достигнут какой-то тривиальный базовый случай (например, n=3). ), при этом общее время работы пропорционально (n + 15n/16 + 225n/256 + ...) = 16n.

Чтобы найти n/16 точек для отбрасывания, требуется большая смекалка. Алгоритм активно использует две подпрограммы:

• median(S, >): принимает набор S элементов и метрику для их попарного сравнения и возвращает медиану элементов.
• MEC-center(S, L): берет набор S точек и линию L и определяет, на какой стороне L находится центр MEC S.

Как упоминалось выше, median() предшествует работе Мегиддо, тогда как алгоритм, описанный здесь как MEC-center(), был представлен как часть его статьи 1983 года. Подробное изучение этих процедур выходит за рамки этого плана, но каждая из них использует сокращение и поиск для выполнения за линейное время. Алгоритм, используемый MEC-center(), представляет собой упрощенную версию алгоритма в целом.

С учетом этих примитивов алгоритм отбрасывания n/16 входных точек работает следующим образом:

• Произвольно соедините n точек в S, чтобы получить n/2 пары.
• Постройте биссектрису для каждой пары точек, чтобы получить всего n/2 биссектрисы.
• Вызовите median(), чтобы найти биссектрису с медианным наклоном. Назовите этот наклон mmid.
• Соедините каждую биссектрису с уклоном ≥ mmid с другой с наклоном ‹ mmid, чтобы получить n/4 точки пересечения. Назовем множество этих точек I.
• Вызовите median(), чтобы найти точку в I с медианным значением y. Назовите это значение y ymid.
• Вызовите MEC-center(), чтобы определить, на какой стороне линии y=ymid находится MEC-центр C. (Не ограничивая общности, предположим, что он лежит выше.)
• Пусть I' будет подмножеством точек I, значения y которых меньше ymid. (I' содержит n/8 точек.)
• Найдите прямую L с наклоном mmid так, чтобы половина точек I' лежала слева от L, а половина - справа.
• Вызовите MEC-center() на L. Без ограничения общности предположим, что C лежит справа от L.
• Пусть I'' будет подмножеством I', точки которого лежат слева от L. (I'' содержит n/16 точек.)

Теперь мы можем отбросить одну точку в S для каждой из n/16 точек пересечения в I''. Рассуждение выглядит следующим образом. После двух вызовов MEC-center() мы обнаружили, что MEC-центр C должен лежать выше ymid и правее L, тогда как любая точка в I'' находится ниже ymid и левее L.

Каждая точка I'' находится на пересечении двух биссектрис. Одна из этих биссектрис должна иметь наклон ≥ mmid и, следовательно, никогда не должна проходить через квадрант, в котором, как мы знаем, лежит C. Назовем эту биссектрису B. Теперь мы знаем, какая сторона B лежит на C, поэтому из двух точек, биссектриса которых равна B, пусть PC будет той, которая лежит на той же стороне, что и C, а другая пусть будет PX.

Легко показать, что PC должен быть ближе к C, чем PX. Отсюда следует, что PC не может лежать на минимальной окружности, и поэтому мы можем безопасно отбросить точку PC для каждой из n/16 точек пересечения в I''.

Мы не обсуждали здесь, как этот алгоритм может работать с вырожденными входными данными (параллельные биссектрисы, коллинеарные точки и т. д.), но оказывается, что мы получаем такие же гарантии производительности для таких случаев. Дело в том, что при вырожденном вводе алгоритм способен отбросить более n/16 точек. Короче говоря, алгоритм Мегиддо гарантирует удаление не менее n/16 точек на каждой итерации независимо от ввода.

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

Дополнительные сведения об алгоритме O(n) см. в этой статье.

ОБНОВЛЕНИЕ. После более внимательного изучения анимации и статьи в Википедии я решил, что мой алгоритм — не Прыгающий Пузырь. Из анимации видно, что Bouncing Bubble перемещает ограничивающую сферу, в то время как мой алгоритм увеличивает ее только в одном направлении за раз. Кроме того, я понимаю, что Bouncing Bubble использует приближение и может не содержать все точки; однако, похоже, есть параметр допуска, которым вы можете управлять.

Я думаю, что мой алгоритм все еще имеет некоторые хорошие свойства:

1. It's O(n).
2. Это пошаговое.
3. Он гарантированно содержит все точки.
4. Это чертовски просто. Единственное, что проще, это найти ограничивающую рамку, выровненную по осям, а затем поместить вокруг нее сферу, но это, безусловно, будет больший объем, чем дает этот алгоритм.

Я придумал инкрементный алгоритм O (n) самостоятельно, который, судя по описанию в Википедии, должен быть алгоритмом прыгающего пузыря. Я понятия не имею, насколько близко это подходит к минимально ограничивающему кругу. Сначала я попытаюсь объяснить свои рассуждения об алгоритме. Это может помочь нарисовать это на бумаге.

Представьте, что у вас есть ограничивающий круг с центром C и радиусом r. Вы хотите расширить круг, чтобы он содержал новую точку P. Вы вычисляете расстояние d от C до P и обнаруживаете, что оно больше, чем r, поэтому оно должно быть вне круга. Теперь, если провести луч из P через C, он выйдет из задней части круга в точке A.

Теперь представьте, что вы закрепляете круг в точке A, а затем расширяете его вдоль луча обратно к P. Он по-прежнему будет содержать все предыдущие точки, и теперь он содержит P. Диаметр этого нового круга равен r + d, поэтому его радиус равен r' = (r + d) / 2. Пройдите по лучу от P до A на расстояние нового радиуса, и вы окажетесь в центре нового круга.

Вот набросок алгоритма:

1. Инициализируйте ограничивающую сферу с центром в первой точке C с радиусом r равным 0.
2. For each remaining point P:
   1. If the distance d from P to C is <= r, skip it.
   2. If d > r:
      1. The new radius r' = (r + d) / 2.
      2. Новый центр C' = P + r' (C - P) / d.

А вот код на Objective-C, который я написал только сегодня:

- (void)updateBoundingVolume {
    self.boundingSphereCenter = GLKVector3Make(0, 0, 0);
    self.boundingSphereRadius = 0;

    GLfloat* vertex = self.vertices;
    GLfloat* endV = vertex + 3 * self.vertNum;

    if (vertex < endV) {
        // initialize to zero radius sphere centered on the first point
        self.boundingSphereCenter = GLKVector3Make(vertex[0], vertex[1], vertex[2]);
        vertex += 3;
    }

    while (vertex < endV) {
        GLKVector3 p = GLKVector3Make(vertex[0], vertex[1], vertex[2]);
        vertex += 3;

        float d = GLKVector3Distance(self.boundingSphereCenter, p);
        if (d <= self.boundingSphereRadius) {
            // already inside the sphere
            continue;
        }

        // length of line from other side of sphere through the center to the new point
        // divide by 2 for radius
        self.boundingSphereRadius = (self.boundingSphereRadius + d) / 2;

        // unit vector from new point to old center
        GLKVector3 u = GLKVector3DivideScalar(GLKVector3Subtract(self.boundingSphereCenter, p), d);
        // step back from the new point by the new radius to get the new center
        self.boundingSphereCenter = GLKVector3Add(p, GLKVector3MultiplyScalar(u, self.boundingSphereRadius));
    }
}