Алгоритмы поиска объема и центра масс многогранника

    Наверное, все знают этот алгоритм, но от меня «власти скрывали». Нашел его словесное описание на третьей странице поисковика в архиве автопереводов англоязычного форума. Мне кажется, его подробное описание (и с кодом) достойно хабростатьи.

    Итак, например вам надо генерировать мобов для игрушки и где-то в процессе отсеивать тех, кто не стоит на ногах. Для этого нужно найти центр масс моба (а это почти то же самое, что найти его объем) и убедиться, что он находится где-то над ногами моба.



    Моб — это многогранник, для простоты считаем, что многогранник состоит только из треугольников (в алгоритме внутри сидит формула площади Гаусса, так что можно расширить его для любого многогранника, но зачем...). Кроме того, многогранник должен не иметь самопересечений и ограничивать замкнутый объем, как и положено приличным многогранникам.


    (ну типа такого)

    Маленький UPD, поясняющий, почему на КДПВ правый моб Не Ок, а левый Ок:
    Правая картинка не Ок потому что моб завалится вперед, т.к. его центр масс вынесен за площадь опоры. Площадь опоры стоящего на поверхности многоугольника определяется как минимальный многоугольник, внутри которого окажутся все точки, находящиеся на поверхности. В левом случае площадь опоры сдвинута под центр масс и больше (т.к. динозаврские лапы больше), а на правой картинке сама площадь меньше и ближе к хвосту.
    Соотношение опорной площади и центра масс будет примерно такое:


    Сразу начну с кода поиска объема (Python, входные данные — список точек и матрица переходов):

    немного кода
    def RecSetDirsTriangles(para, Connects, TR):
        """рекурсивная функция, которая принимает на вход ребро и по нему выбирает направление треугольника"""
        #1.найти треугольник включающий пару, убедиться, что такого еще нет
        for i in range(0,len(Connects)):
            if i != para[0] and i != para[1] and Connects[i][para[0]] and Connects[i][para[1]]: #вот этот треугольник!
                fl = 1
                for T in TR:
                    if i in T and para[1] in T and para[0] in T:
                        fl = 0 #этот треугольник уже обработан
                        break
                if fl: #найден треугольник!
                    TR += [(para[1],para[0],i)]
                    Recc((para[0], i) , Connects, TR)
                    Recc((i, para[1]) , Connects, TR)
    
    def FindV(dots, Connects):
        """ищем объем. Входные данные - dots список вершин многогранника вида [x, y, z], Connects - квадратная матрица, Connects[i][j]=1 если есть связь между вершинами i, j, иначе =0 """
        #1. сделать треугольники с упорядоченными вершинами
        TR = []
        for i in range(1,len(Connects)):#выбираем первый треугольник с нулевой точкой и еще каким-то
            if Connects[i][0]:
                for j in range(i+1, len(Connects)):
                    if Connects[0][j] and Connects[i][j]:
                        TR += [(0,i,j)]
                        break
                RecSetDirsTriangles((0,i),Connects, TR)
                break
        print("найдено треугольников: ", len(TR))
        
        #2. посчитать площадь базы и объем усеченной призмы
        V = 0
        for T in TR:
            '''Гаус рулит:           x1y2                         x2y3                         x3y1                           x2y1                        x3y2                        x1y3'''
            S = 0.5 * (dots[T[0]][0]*dots[T[1]][1] + dots[T[1]][0]*dots[T[2]][1] + dots[T[2]][0]*dots[T[0]][1] - dots[T[1]][0]*dots[T[0]][1] - dots[T[2]][0]*dots[T[1]][1] - dots[T[0]][0]*dots[T[2]][1])
            #S может быть + или - в зависимости от того, как направлен треугольник
            V += S*(dots[T[0]][2] + dots[T[1]][2] + dots[T[2]][2])/3 #объем усеченной призмы считается просто...
    
        return math.fabs(V)
    


    Суть алгоритма — считаем объемы фигур, которые образуют «падающие» на плоскость xy грани многогранника. Для этого надо знать площадь проекции треугольника и знак, с которым надо суммировать объем фигуры (усеченнной призмы). На самом деле, если заранее упорядочить треугольники, и объем и знак сводятся к одному вычислению.

    Поэтому первым делом рекурсивная функция собирает треугольники из входных данных. Собирает таким образом, чтобы при взгляде «снаружи» на многогранник, направления обхода треугольников были одинаковыми (в идеале против часовой стрелки; если взять направления по часовой стрелке, то результат получится правильным, но отрицательным — поэтому в ретурн отдается модуль объема).

    Добиться этого очень просто — берем какой-то треугольник (точки a1, a2, a3), ищем его соседей и перечисляем две совпавшие вершины в обратном порядке (например, так: a2, a1, b1).
    Получается что-то вроде этого:



    Теперь, если мы спроецируем такой треугольник на плоскость xy, то порядок обхода для проекции «верхнего» треугольника будет совпадать с изначально выбранным, а порядок обхода для проекци «нижнего» треугольника поменяет свое направление. Как следствие, поменяет знак и площадь этого треугольника, вычисленная по формуле Гаусса. Здесь «нижний» треугольник — понятие условное — имеется ввиду, что объем непосредственно под ним не входит в объем многогранника. «Нижний» треугольник у невыпуклого многогранника может быть выше «верхнего».

    После этих предварительных действий, чтобы вычислить полный объем многогранника, надо просто сложить (с учетом знака, который получается «сам собой») все объемы усеченных призм, собранных из граней и проекций этих граней на плоскость xy. А объемы призм считаются как произведение площади (по Гауссу, со знаком) и среднего арифметического z-координат вершин треугольника.

    Если многогранник пересекается плоскостью xy, то при вычислении объема, все знаки скомпенсируют друг друга и результат остается правильным (надо только брать высоты призмы без модуля).


    (как-то так выглядит «верхняя» усеченная призма)

    С поиском центра масс все приблизительно также. Аналогично надо найти центры масс для каждой усеченной призмы и просуммировать покоординатно, умножая на объем призмы (предполагается, что масса распределена равномерно по объему и можно одно заменить другим). Чтобы найти центр масс усеченной призмы, придется посчитать центры масс двух тетраэдеров (+1 функция) и одной обычной призмы. Алгоритм так же «не портится», если многогранник пересекает плоскость xy (а здесь могла бы быть репродукция Магритта).


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

    Код, который считает то и то:

    чуть больше кода
    def RecSetDirsTriangles(para, Connects, TR):
        #1.найти треугольник включающий пару, убедиться, что такого еще нет
        for i in range(0,len(Connects)):
            if i != para[0] and i != para[1] and Connects[i][para[0]] and Connects[i][para[1]]: #вот этот треугольник!
                fl = 1
                for T in TR:
                    if i in T and para[1] in T and para[0] in T:
                        fl = 0
                        break
                if fl: #найден треугольник!
                    TR += [(para[1],para[0],i)]
                    Recc((para[0], i) , Connects, TR)
                    Recc((i, para[1]) , Connects, TR)
        
    def TetrV(mas):#dot1, dot2, dot3, dot4):
        """объем тетраэдера по вершинам"""
        M = np.zeros((3,3),float)
        for i in range(1,4):
            for j in range(0,3):
                M[i-1][j] = mas[i][j] - mas[0][j]
        #print(M)
        return math.fabs(np.linalg.det(M)/6)
    
    def FindVandCM(dots, Connects):
        """ищем объем и центр масс многогранника"""
        #1. сделать треугольники с упорядоченными вершинами
        TR = []
        for i in range(1,len(Connects)): #выбираем первый треугольник с нулевой точкой и еще каким-то
            if Connects[i][0]:
                for j in range(i+1, len(Connects)):
                    if Connects[0][j] and Connects[i][j]:
                        TR += [(0,i,j)]
                        break
                RecSetDirsTriangles((0,i),Connects, TR)
                break
        print("найдено треугольников: ", len(TR))
        
        #2. посчитать площадь базы, объем усеченной призмы и вклад в центр масс каждого
        V = 0
        CM = [0, 0, 0]
        for T in TR:
            '''Гаус рулит:           x1y2                         x2y3                         x3y1                           x2y1                        x3y2                        x1y3'''
            S = 0.5 * (dots[T[0]][0]*dots[T[1]][1] + dots[T[1]][0]*dots[T[2]][1] + dots[T[2]][0]*dots[T[0]][1] - dots[T[1]][0]*dots[T[0]][1] - dots[T[2]][0]*dots[T[1]][1] - dots[T[0]][0]*dots[T[2]][1])
            #S может быть + или - в зависимости от того, как направлен треугольник
            V += S*(dots[T[0]][2] + dots[T[1]][2] + dots[T[2]][2])/3 #объем усеченной призмы считается просто...
            
            #c центром масс я так просто не отделаюсь
            c1 = ((dots[T[0]][0] + dots[T[1]][0] + dots[T[2]][0])/3, (dots[T[0]][1]+ dots[T[1]][1]+ dots[T[2]][1])/3) #центральная точка проекции треугольника
            hm = min([dots[T[0]][2] , dots[T[1]][2] , dots[T[2]][2]])
            hM = max([dots[T[0]][2] , dots[T[1]][2] , dots[T[2]][2]])
            indM = [dots[T[0]][2] , dots[T[1]][2] , dots[T[2]][2]].index(hM)
            indm = [dots[T[0]][2] , dots[T[1]][2] , dots[T[2]][2]].index(hm)
            V3 = S * hm
            if indM == indm: #горизонтальный прямоугольник!
                CM[0] += V3*c1[0]
                CM[1] += V3*c1[1]
                CM[2] += V3*hm/2
                continue
            L = [0,1,2]
            
            L.remove(indM)
            L.remove(indm)
            indmidle = L[0]
            dots1 = [dots[T[0]], dots[T[1]], dots[T[2]], (dots[T[indM]][0], dots[T[indM]][1] , hm)] #верхний тетраэдер
            V1 = TetrV(dots1)
            if S < 0:
                V1 = -V1
            V2 = S * ( dots[T[indmidle]][2] - hm)/3
            #V3 = S * hm
            CM[0] += V1*((dots[T[0]][0] + dots[T[1]][0] + dots[T[2]][0] + dots[T[indM]][0])/4) + V2*((dots[T[0]][0] + dots[T[1]][0] + dots[T[2]][0] + dots[T[indmidle]][0])/4) + V3*c1[0]
            CM[1] += V1*((dots[T[0]][1] + dots[T[1]][1] + dots[T[2]][1] + dots[T[indM]][1])/4) + V2*((dots[T[0]][1] + dots[T[1]][1] + dots[T[2]][1] + dots[T[indmidle]][1])/4) + V3*c1[1]
            CM[2] += V1*((dots[T[0]][2] + dots[T[1]][2] + dots[T[2]][2] + hm)/4) + V2*((dots[T[0]][2] + dots[T[1]][2] + dots[T[2]][2] + hm)/4) + V3*hm/2
    
        CM[0] = CM[0]/V
        CM[1] = CM[1]/V
        CM[2] = CM[2]/V
    
        return (math.fabs(V), CM)
    


    Кусок алгоритма, где считаются направления треугольников и используются для понимания внешнего и внутреннего объема — это очень сильный ход, его много как можно применить при работе с многогранниками. Например, если надо посчитать направление нормалей «наружу» — достаточно знать направление «против часовой стрелки» для одной грани — и вуаля!


    (угадай фильм!)
    AdBlock похитил этот баннер, но баннеры не зубы — отрастут

    Подробнее
    Реклама

    Комментарии 22

      0

      Очень интересная тема, но есть несколько замечаний.


      Итак, например вам надо генерировать мобов для игрушки и где-то в процессе отсеивать тех, кто не стоит на ногах. Для этого нужно найти центр масс моба (а это почти то же самое, что найти его объем) и убедиться, что он находится где-то над ногами моба.

      Что именно значит "где-то над ногами"? Каким образом знание координат точки центра масс поможет понять, что модель "стоит на ногах"?


      Почему левая картинка ок, а правая — неок?


      Кроме того, многогранник должен не иметь самопересечений и ограничивать замкнутый объем, как и положено приличным многогранникам.

      Алгоритм поиска центра масс должен проверять исходные данные и не должен возвращать результат в случае многогранников с самопересечениями и в "незамкнутости".


      Код поиска объема настолько лаконичен

      … что с трудом не влезает на экран монитора


      входные данные — список точек и матрица переходов

      Не хватает точного описания схемы данных. Как именно хранятся x,y,z в точках? Как именно хранится матрица переходов?


      Замечания по коду первого фрагмента:


      • Магические числа в индексах массивов нужно заменить на человекопонятные идентификаторы или подготовить данные так, чтобы вместо массивов использовались АТД с понятными полями.
      • Названия функций невыразительны. Что вернет функция Recc, в случае когда не удалось найти треугольник, включающий пару?
      • Для проверки "содержится ли треугольник T во множестве TR" лучше воспользоваться соответствующей структурой данных — Set. Тогда можно будет избавиться от итерирования по TR.
      • Модифицировать аргумент TR в теле функции Recc на мой взгляд — плохая идея. Лучше реализовать эту функцию с возвращаемым значением.
      • TR нужно переименовать, этот идентификатор совершенно ничего не говорит о своем содержимом.
      • Выражение "if fl:" на самом деле значит "треугольник не найден", что противоречит комментарию возле него. Если бы он был найден, то fl было бы равно нулю.
      • Функцию FindV тоже нужно переименовать. Что такое V?
      • Почему в функции FinvV итерирование выполняется с 1, а не с 0?

      Из статьи не понятно, зачем нам вообще искать объем. Допустим мы его нашли. Это какое-то число 42. Как из него получить координаты центра масс?
      Напишите сначала алгоритм в общих чертах, а потом уже переходите к ньюансам каждого отдельного этапа.


      Чтобы найти центр масс усеченной призмы, придется посчитать центры масс двух тетраэдеров (+1 функция) и одной обычной призмы.

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


      Замечания по второму фрагменту кода:


      • Зачем в нем код из первого фрагмента? Весь листинг в тексте статьи лучше разбить на отдельные функции, а полный текст привести в конце статьи.
      • Число 6 в функции TetrV нужно заменить на константу и добавить объяснение, почему именно 6. Остальные магические числа тоже нужно заменить константами.
      • Портянка со сложением-умножением элементов массива в конце нечитаема. Разбейте ее на отдельные функции с понятными именами.

      Кусок алгоритма, где считаются направления треугольников и используются для понимания внешнего и внутреннего объема — это очень сильный ход, его много как можно применить при работе с многогранниками

      Его нигде нельзя применить, потому что вообще не понятно, какая часть этого кода считает "направления" треугольников. Как это вообще относится к алгоритму центра масс?


      Пожалуйста, уберите статью в черновики как только получите зачет.

        –6
        Есть замечания по делу, но хамить то зачем? Теперь исправлять желания нет — и для чего вы так длинно писали?
          0
          Где в комментарии вы увидели «хамство»? Все расписано детально и по делу, именно эти вопросы и возникают у читателя. Вы как-то больно близко к сердцу принимаете замечания. Передохните и перечитайте свой пост свежим взглядом.

          И да, меня тоже интересует, как из поста следует, что левая картинка ок, а правая — не ок?
            0
            Замечания нормальные. Хамство в последней строчке (она по мотивам недавнего «расследования» про студентов, короче меня школьником обозвали).

            Правая картинка не Ок потому что моб завалится вперед, т.к. его центр масс вынесен за площадь опоры. Площадь опоры стоящего на поверхности многоугольника определяется как минимальный многоугольник, внутри которого окажутся все точки, находящиеся на поверхности. В левом случае площадь опоры сдвинута под центр масс и больше (т.к. динозаврские лапы больше), а на правой картинке сама площадь меньше и ближе к хвосту.
            Ну как-то так:
              0
              Добавьте это в статью, как описание проблемы, которую алгоритм помогает решить.
                0
                А разве в статье этого нет?
                0
                Ну так-то я понял же. Но специально добавил уточнение:

                > как из поста следует, что левая картинка ок, а правая — не ок?

                В посте это никак не обыгрывается, ни про ноги, ни про еще что-то. МОжет, мне как раз это и было интересно, а из статьи непонятно…
                  0
                  а из статьи непонятно

                  Эмм… ИМХО, в контексте поиска центра масс, это по самой картинке очевидно.
                    +1
                    Как нам говорили в университете, «очевидно — значит, легко доказать». Ну и вот мне ни разу не очевидно ни по картинке ни из статьи.
                    Из статьи непонятно (потому что статья вообще не о том была). Ну и топ-коммент же как раз об этом, ну:

                    Что именно значит «где-то над ногами»? Каким образом знание координат точки центра масс поможет понять, что модель «стоит на ногах»?

                    Если мы из носа добавляем тонкую невесомую «ногу», то получим треногу при условно неизменном центре масс (поскольку «нога» может быть бесконечно тонкая и бесконечно легкая). Как знание центра масс помогает определить устойчивость?
                      0
                      Нужно знать центр масс и набор точек, в которых тело касается к поверхности. Причем будет куда проще, если эта поверхность задана как плоскость z = 0. И если «третья нога» достигает этой поверхности — нужно знать и в ней координаты «всех точек».
                      Но мы задаем вершины, так что нужно знать конкретные точки, у которых z = 0.
                        0
                        Та я все понимаю ж. Но почему это объяснение я читаю в комментах, а не в статье? Вопрос-то к автору был именно об этом. И еще раз напоминаю, что автор в статье искал центре масс\объем, а слово «опорная площадь» появилось только в «upd» части.

                        Ну и в контексте многогранников я так думаю, что многогранник задан точками\полигонами, и у него нет преопределенной ориентации «низ там». Соответственно, как вы определяете, где же та плоскость, где z=0? Почему ориентация нашего динозавра именно такая, а не другая? Если там снизу плоское копыто, это одно, а если внизу «лапа», то она же неоднородная и там много вариантов, как провести плоскость «пола».
                          0
                          Хотел бы я знать ответ на последний вопрос. Особенно с учетом того, что у меня там вместо лапы сплайн, натянутый на низкополигональный меш (в прошлой моей статье есть подробности).
                          Пока рабочее решение — уровень «пола» определяется нижней точкой, а границы «многоугольника опоры» определяются x,y координатами нижних 5% точек (т.е. точки, z координата которых находится в пределах min(z)..min(z)+(max(z)-min(z))*0.05). Такого определения должно хватить для псевдо-3D игры на спрайтах. Когда-нибудь я доползу до прототипа и буду знать точно.
                            0
                            z=0 я взял как наиболее простой вариант, в сочетании с предположением, что гравитация направлена по вектору (0,0,-1).
                            Тогда нам нужно считать только 2 координаты центра масс для ответа на вопрос, находится ли он над площадкой опоры.
                            И «площадку опоры» строго выбирать из точек, у которых z=0. А если вердикт «динозавр не стоит устойчиво», тогда сила тяжести должна сместить все его вершины, пока не будет выполняться это условие.
                          0
                          очевидно — значит, легко доказать

                          Почему бы нет?
                          Контекст статьи — поиск центра масс фигуры, составленной из полигонов.
                          На картинке показаны две фигуры, используемые, вероятно, в какой-то игре (или чём-то подобном).
                          Для чего в физике используется понятие центра масс? Ну, вообще для ряда вещей, но наиболее известное применение (известное ещё по школьной физике) — для определения устойчивости тела.
                          Из этого естественным образом вытекает гипотеза о том, что «не ОК» — это неустойчивая фигура.
                          Для проверки этой гипотезы мысленно дорисовываем картинку как в комменте выше.
                          Да, гипотеза подтвердилось, значит весьма вероятно, что автор именно это и имел в виду
                    +1
                    Пожалуйста, уберите статью в черновики как только получите зачет.

                    Я думаю, что это вполне себе оскорбление, хоть и завуалированное.
                  0
                  Число 6 в функции TetrV нужно заменить на константу и добавить объяснение, почему именно 6. Остальные магические числа тоже нужно заменить константами.

                  Если это объем тетраэдра через смешанное произведение, тогда все правильно.
                    0
                    Алгоритм поиска центра масс должен проверять исходные данные и не должен возвращать результат в случае многогранников с самопересечениями и в «незамкнутости».


                    Совершенно необязательно. Алгоритм поиска центра масс должен искать центр масс. Вполне нормально не проверять входные данные, если отсутствие такой проверки оговорено заранее, например, в комментарии перед функцией/документации к апи и т.д. Причин несколько. Первое — проверка может занимать больше времени, чем сам поиск центра масс. В данном случае скорее всего так и будет, корректная проверка многогранника на самопересечения штука не такая уж простая, особенно если координаты точек представлены числами с плавающей запятой. Второе — входной многогранник может быть проверен где-то сильно раньше, до вызова функции, может быть вообще по построению несамопересекающийся и «замкнутый». В таких случаях лишние проверки сделают алгоритм в пару раз медленнее, а безопасности особо не добавят. Разумеется, при этом желательно иметь также «безопасную» версию со всеми проверками для дебажных целей с возможностью быстро сменить версию, если нужно. Третье — результат может быть полезен и для некоторых самопересекающихся многогранников, если большая точность не нужна. Допустим, у вас есть меш того же динозавра, в нем подвинули ногу и хотят узнать, куда переместился центр масс. При этом одна нога может, условно, на пару сантиметров зайти в другую ногу. Да, формально меш самопересекающийся, но фактически вас это не волнует, потому что влияние на результат околонулевое, а у вас игра и идеальной точности не требуется.
                      0
                      Допустим, у вас есть меш того же динозавра, в нем подвинули ногу и хотят узнать, куда переместился центр масс.


                      Для конкретной задачи «понять, устойчиво ли стоит» можно вообще считать только 2 координаты ЦМ и потом определять, попадают ли они в многоугольник «ноги + площадь между ними». Но это нужно ввести условие, что у нас есть точки «ног», которые все в одной плоскости и ничего ниже нет.
                    0
                    формула площади Гаусса

                    Для выпуклого многоугольника эта формула может быть доказана через векторные произведения?
                      0
                      Честно говоря, понятия не имею :)
                        0
                        Сама по себе формула площади по идее может быть записана как
                        A = 1/2*|sum(i=1..n-1){[ri-r1,ri+1-r1]}|,
                        где под [x,y] имеется в виду векторное произведение векторов.
                      0
                      del

                      Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                      Самое читаемое