[Из песочницы] Сплайны в 3d графике, максимально автоматизированный вариант30.10.2019 15:20

С месяц назад начал учить Python по книге Доусона и очнулся уже глубоко в процессе написания своей игры под pygame. ТЗ было таково, что наиболее перспективным показалось сделать игру с псевдо-трехмерной графикой, запихнув в спрайты сохраненные поверхности 3d-сплайнов. О последних и напишу.

Итак, имеются полигоны (проще всего работать с четырехугольниками), на которые мы хотим натянуть кубические поверхности так, чтобы они стыковались достаточно плавно — эти поверхности и есть сплайны.

Удобнее всего сплайн для одного полигона представить в виде функции

$$display$$g (u, v) = v\cdot g_1 (u) + (1 — v)\cdot g_2 (u) + u\cdot g_3 (v) + (1-u)\cdot g_4(v) (1) $$display$$

Здесь

$$display$$g_1, g_2, g_3, g_4$$display$$

— кубические многочлены, удовлетворяющие каким-то граничным условиям (на рисунке ниже — салатовые и рыжие кривые, а производные-начальные условия — сиреневые и голубые вектора); u и v меняются в пределах от 0 до 1.

В такой интерпретации теряются некоторые степени (произведения 2-х и 3-х степеней параметров), зато коэффициенты многочлена можно найти, задав всего 12 векторов начальных условий (4 точки полигона и вектора производных по u и по v для каждой точки). На стыках сплайны совпадают, если задаются аналогичные начальные условия для соседних полигонов (вектора производных должны быть коллинеарны, поверхность, проходит через те же точи полигона).

Одна беда — производная при такой постановке задачи на всей границе может не совпадать — будут небольшие артефакты на стыках. Можно выдумать еще 4 условия для исправления этого и аккуратно добавить в формулу еще один член

$$display$$u\cdot v\cdot (u-1)\cdot (v-1)\cdot g_5(u, v)$$display$$

который не испортит границы, но это тема для отдельной статьи.

Альтернативный вариант — Поверхность Безье, но в нем предлагается задавать 16 параметров непонятного (мне) математического смысла, т.е. предполагается, что будет работать своими руками художник. Мне это не очень подходит, поэтому далее следует велосипед с костылями.

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

NumPy отлично с этой задачей справляется.

Остался один вопрос — откуда взять векторы производных. Предполагается, что надо их как-то выбирать исходя из положения соседних (для полигона) точек и из соображения гладкости.
Лично для меня еще было совершенно контринтуитивно, как быть с длиной вектора производной. Направление — понятно, а вот длина?

В итоге родился такой алгоритм:

1. На первом этапе происходит некоторая классификация точек. В графе (который задают точки полигонов и их связи) ищутся и запоминаются циклы длиной 4, а кроме того записываются соседи, наиболее подходящие на роль продолжения ребра полигона (заранее определяется, какие ребра соответствуют изменению параметра u, а какие — v). Вот кусок кода, который ищет, сортирует и запоминает соседей для 0-й точки какого-то цикла:

            """   схематичное расположение точек:
                                        cykle[5]       cykle[7]
                              cykle[4]--cykle[0] --u-- cykle[1]-cykle[6]
                                            |v              |v
                              cykle[10]-cykle[3] --u-- cykle[2]-cykle[8]
                                       cykle[11]      cykle[9]
            """
                sosed = []
                for i in range(0, N):
                    if self.connects[i][ind] == 1 and i != cykle[0] and i != cykle[1] and i != cykle[3]:
                        sosed += [i] #в нашем доме поселился замечательный...
                if(len(sosed) < 2):
                    continue #пока такое не обрабатываем 
                else:
                    #вычисляем длины векторов между соседями точки cykle[0]
                    vs0c3 = MinusVecs(self.dots[sosed[0]], self.dots[cykle[3]])
                    vs1c1 = MinusVecs(self.dots[sosed[1]], self.dots[cykle[1]])
                    #есть два варианта и один из них длиннее
                    vs0c1 = MinusVecs(self.dots[sosed[0]], self.dots[cykle[1]])
                    vs1c3 = MinusVecs(self.dots[sosed[1]], self.dots[cykle[3]])
                    Rvsc = [ScalMult(vs0c3,vs0c3),ScalMult(vs1c1,vs1c1),ScalMult(vs0c1,vs0c1),ScalMult(vs1c3,vs1c3) ]
                    n1 = Rvsc.index(max(Rvsc))
                    if n1 < 2:
                        cykle += [sosed[1],sosed[0]] #сначала по u потом по v
                    else:
                        cykle += [sosed[0],sosed[1]]

Дальше, при построении сплайна, производная вдоль ребра (например по параметру u) в точке полигона выбирается исходя из расположения двух точек ребра и одной соседней точки (пусть это будут точки vec1, vec2 и vec3; точка в которой ищется производная — вторая).

Сначала я пытался использовать для этой роли нормированный вектор vec3 — vec1 (типа применил теорему Лагранжа), но проблемы возникли именно с длиной вектора производных — делать его константой было плохой идеей.

Лирическое отступление:

Для небольшой иллюстрации, что такое вектор производной в параметрическом варианте, обратимся к простой двумерной аналогии — вот кусок дуги окружности:

$$display$$g (u) = (Rsin (u\pi/2) , Rcos (u\pi/2))$$display$$

$$display$$g'(u) =(R\pi/2\cdot cos (u\pi/2), -R\pi/2\cdot sin (u\pi/2))$$display$$

Т.е. модуль производной = R*pi/2 и, вообще говоря, зависит от размера куска дуги, который мы задали через параметр.

Что же с этим теперь делать? Лев Николаевич Толстой завещал нам, что все круглое = хорошее, поэтому, если мы зададим производную в точке так, как если бы хотели построить там дугу, — у нас получится гладкая красивая кривая.

Конец лирического отступления.

Второй этап поиска производной:

2. Через наши три точки vec1, vec2, vec3 строим плоскость, в этой плоскости ищем окружность, которая проходит через все три точки. Искомая производная будет направлена по касательной к окружности в точке vec2, а модуль вектора производной должен равняться произведению радиуса окружности и угла сектора, который образуют две точки грани полигона (аналогично нашей простой плоской аналогии из лирического отступления).

$$display$$g' = R\cdot \phi$$display$$

Вроде все просто, вот код функции, например (опять используется NumPy):

def create_dd(vec1, vec2, tuda, vec3):
    """делает вектор производной для отрезка 1-2 во второй точке в направлении точки 3 если туда == 1"""
    #0. убедиться что точки не лежат на одной прямой - в этом случае окружность через них не провести :(    
    vecmult = VecMult(MinusVecs(vec2,vec1),MinusVecs(vec3,vec1)) #векторное произведение vec2-vec1 x vec3 - vec1
    if vecmult[0]*vecmult[0] + vecmult[1]*vecmult[1] + vecmult[2]*vecmult[2] < 0.000001: #ну почти 0
        if tuda == 1:
            return NormVec(MinusVecs(vec3,vec2))
        else:
            return NormVec(MinusVecs(vec1,vec2))
        
    #1. найти центр окружности, проходящей через эти точки
    #сначала два условия про то что расстояния до центра одинаковые порождают 2 линейных уравнения
    M = np.zeros((3,3),float)
    C = np.zeros((3,1),float)
    M[0][0] = 2*vec2[0] - 2*vec1[0]
    M[0][1] = 2*vec2[1] - 2*vec1[1]
    M[0][2] = 2*vec2[2] - 2*vec1[2]
    C[0] = -vec1[0]*vec1[0] - vec1[1]*vec1[1] - vec1[2]*vec1[2] + vec2[0]*vec2[0] + vec2[1]*vec2[1] + vec2[2]*vec2[2]
    M[1][0] = 2*vec3[0] - 2*vec2[0]
    M[1][1] = 2*vec3[1] - 2*vec2[1]
    M[1][2] = 2*vec3[2] - 2*vec2[2]
    C[1] = -vec2[0]*vec2[0] - vec2[1]*vec2[1] - vec2[2]*vec2[2] + vec3[0]*vec3[0] + vec3[1]*vec3[1] + vec3[2]*vec3[2]
    #еще одно уравнение указывает, что все 4 вектора в одной плоскости
    M[2][0] = vecmult[0]
    M[2][1] = vecmult[1]
    M[2][2] = vecmult[2]
    C[2] = vecmult[0]*vec1[0] + vecmult[1]*vec1[1] + vecmult[2]*vec1[2]
    # M * вектор центра окружности = C
    Mobr = np.linalg.inv(M)
    Rvec = np.dot(Mobr,C)  
    res = NormVec(VecMult(MinusVecs(Rvec,vec2), vecmult)) #направление с точностью до знаka
    #R = math.sqrt((Rvec[0]-vec1[0])**2 + (Rvec[1]-vec1[1])**2 + (Rvec[2]-vec1[2])**2)
    #длина вектора производной = 3.14*2*asin(l/(2*R))/180 * R, но можно попытаться заменить на l с коэффициентом, взятым с потолка
    l = 1.2*math.sqrt((vec1[0]-vec2[0])**2 + (vec1[1]-vec2[1])**2 + (vec1[2]-vec2[2])**2)
    if ScalMult(res, MinusVecs(vec2, vec1)) > 0:  #вектор смотрит в ту сторону
        return MultVxC(res, tuda*l)
    else:
        return MultVxC(res, -tuda*l)

Ну и в общем, это все как-то работает. Для демонстрации я взял куб 5×5х5 и менял положение точек рандомайзером:

гифка, 8Мб

© Habrahabr.ru