Наверное, все знают этот алгоритм, но от меня «власти скрывали». Нашел его словесное описание на третьей странице поисковика в архиве автопереводов англоязычного форума. Мне кажется, его подробное описание (и с кодом) достойно хабростатьи.
Итак, например вам надо генерировать мобов для игрушки и где-то в процессе отсеивать тех, кто не стоит на ногах. Для этого нужно найти центр масс моба (а это почти то же самое, что найти его объем) и убедиться, что он находится где-то над ногами моба.
Моб — это многогранник, для простоты считаем, что многогранник состоит только из треугольников (в алгоритме внутри сидит формула площади Гаусса, так что можно расширить его для любого многогранника, но зачем...). Кроме того, многогранник должен не иметь самопересечений и ограничивать замкнутый объем, как и положено приличным многогранникам.
(ну типа такого)
Маленький UPD, поясняющий, почему на КДПВ правый моб Не Ок, а левый Ок:
Правая картинка не Ок потому что моб завалится вперед, т.к. его центр масс вынесен за площадь опоры. Площадь опоры стоящего на поверхности многоугольника определяется как минимальный многоугольник, внутри которого окажутся все точки, находящиеся на поверхности. В левом случае площадь опоры сдвинута под центр масс и больше (т.к. динозаврские лапы больше), а на правой картинке сама площадь меньше и ближе к хвосту.
Соотношение опорной площади и центра масс будет примерно такое:
Сразу начну с кода поиска объема (Python, входные данные — список точек и матрица переходов):
Суть алгоритма — считаем объемы фигур, которые образуют «падающие» на плоскость xy грани многогранника. Для этого надо знать площадь проекции треугольника и знак, с которым надо суммировать объем фигуры (усеченнной призмы). На самом деле, если заранее упорядочить треугольники, и объем и знак сводятся к одному вычислению.
Поэтому первым делом рекурсивная функция собирает треугольники из входных данных. Собирает таким образом, чтобы при взгляде «снаружи» на многогранник, направления обхода треугольников были одинаковыми (в идеале против часовой стрелки; если взять направления по часовой стрелке, то результат получится правильным, но отрицательным — поэтому в ретурн отдается модуль объема).
Добиться этого очень просто — берем какой-то треугольник (точки a1, a2, a3), ищем его соседей и перечисляем две совпавшие вершины в обратном порядке (например, так: a2, a1, b1).
Получается что-то вроде этого:
Теперь, если мы спроецируем такой треугольник на плоскость xy, то порядок обхода для проекции «верхнего» треугольника будет совпадать с изначально выбранным, а порядок обхода для проекци «нижнего» треугольника поменяет свое направление. Как следствие, поменяет знак и площадь этого треугольника, вычисленная по формуле Гаусса. Здесь «нижний» треугольник — понятие условное — имеется ввиду, что объем непосредственно под ним не входит в объем многогранника. «Нижний» треугольник у невыпуклого многогранника может быть выше «верхнего».
После этих предварительных действий, чтобы вычислить полный объем многогранника, надо просто сложить (с учетом знака, который получается «сам собой») все объемы усеченных призм, собранных из граней и проекций этих граней на плоскость xy. А объемы призм считаются как произведение площади (по Гауссу, со знаком) и среднего арифметического z-координат вершин треугольника.
Если многогранник пересекается плоскостью xy, то при вычислении объема, все знаки скомпенсируют друг друга и результат остается правильным (надо только брать высоты призмы без модуля).
(как-то так выглядит «верхняя» усеченная призма)
С поиском центра масс все приблизительно также. Аналогично надо найти центры масс для каждой усеченной призмы и просуммировать покоординатно, умножая на объем призмы (предполагается, что масса распределена равномерно по объему и можно одно заменить другим). Чтобы найти центр масс усеченной призмы, придется посчитать центры масс двух тетраэдеров (+1 функция) и одной обычной призмы. Алгоритм так же «не портится», если многогранник пересекает плоскость xy (а здесь могла бы быть репродукция Магритта).
(вот эти два тетраэдера, обозначены красным и рыжим, вместе с треугольной призмой (ниже красного тетраэдера) образуют искомую усеченную призму. Нам надо найти центры масс и объемы всех трех фигур. Обозначения примерно соответствуют обозначениям в коде)
Код, который считает то и то:
Кусок алгоритма, где считаются направления треугольников и используются для понимания внешнего и внутреннего объема — это очень сильный ход, его много как можно применить при работе с многогранниками. Например, если надо посчитать направление нормалей «наружу» — достаточно знать направление «против часовой стрелки» для одной грани — и вуаля!
(угадай фильм!)
Итак, например вам надо генерировать мобов для игрушки и где-то в процессе отсеивать тех, кто не стоит на ногах. Для этого нужно найти центр масс моба (а это почти то же самое, что найти его объем) и убедиться, что он находится где-то над ногами моба.
Моб — это многогранник, для простоты считаем, что многогранник состоит только из треугольников (в алгоритме внутри сидит формула площади Гаусса, так что можно расширить его для любого многогранника, но зачем...). Кроме того, многогранник должен не иметь самопересечений и ограничивать замкнутый объем, как и положено приличным многогранникам.
(ну типа такого)
Маленький 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)
Кусок алгоритма, где считаются направления треугольников и используются для понимания внешнего и внутреннего объема — это очень сильный ход, его много как можно применить при работе с многогранниками. Например, если надо посчитать направление нормалей «наружу» — достаточно знать направление «против часовой стрелки» для одной грани — и вуаля!
(угадай фильм!)
