Версия для печати


Использование методов аналитической геометрии для просчёта столкновений

---

http://www.delphikingdom.com/asp/viewitem.asp?catalogID=922

Ушаков Юрий
дата публикации 16-02-2004 10:52

В своей статье (Вариант реализации простейших костных деформаций) я начал описывать один из своих проектов - библиотеку для трёхмерных игр. И в прошлой статье я описал первый шаг - я ввёл формат модели, описал способ, которым предлагается строить модели в пространстве.

Просматривая материалы Руннета, определяя тему своей очередной статьи, я обнаружил, что в Руннете очень мало пишут о том, о чем будет говориться в статье. Руннет переполнен сайтами, авторы которых помешаны на создании игр, и некоторые даже заявляют, что движки их игр по возможностям сравнимы в движком Quake, или чего покруче. Однако тут же делают поправку - нет просчёта столкновений.

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

Математические процедуры помогают реализовать передвижение персонажа в пространстве, помогают распознать препятствия на его пути и должным образом отреагировать на них. В этой статье я расскажу, как можно выполнять проверку столкновений. Но не спешите, создавать замысловатые лабиринты аля-Quake ещё рано.

Задача

Пусть материальная точка движется из точки (x₀;y₀;z₀) пространства на вектор (Δ_x;Δ_y;Δ_z). Имеется так же множество треугольников в пространстве, которые заданы координатами вершин:

Δ1(x_1a;y_1a;z_1a;x_1b;y_1b;z_1b;x_1c;y_1c;z_1c),
Δ2(x_2a;y_2a;z_2a;x_2b;y_2b;z_2b;x_2c;y_2c;z_2c),
Δ3(x_3a;y_3a;z_3a;x_3b;y_3b;z_3b;x_3c;y_3c;z_3c),
...,
Δi(x_ia;y_ia;z_ia;x_ib;y_ib;z_ib;x_ic;y_ic;z_ic).

Требуется определить, столкнётся ли при движении на заданный вектор материальная точка с одним из треугольников, и, если столкнётся, то определить, на какой вектор (Δx`;Δy`;Δz`) она будет двигаться после столкновения.

Решение

Разобьём нашу задачу на несколько частей. Первая часть будет определять, произойдёт ли столкновение. Следующая часть будет просчитывать координаты нового вектора, и основная часть алгоритма будут осуществлять перебор всех поверхностей и проверять, не произойдёт ли с одной из них столкновение. Нужно сказать, что такая задача требует решения методами аналитической геометрии. Поэтому для тех, кто по какой-либо причине не владеет этими методами, я буду приводить некоторые теоретические сведения в материале.

Для этих целей организуем отдельный модуль Collision.pas, в котором будет находиться всего две функции в разделе интерфейса, а в разделе реализации будет присутствовать ещё несколько математических функций.

Проверка столкновения

Начнём с того, что определим, будет ли столкновение. Такую проверку удобно реализовать в виде отдельной функции. Эта функция поможет программе не только проверять столкновения при их просчёте, но и выполнять, например, просчёт попадания при стрельбе. В общем, запишем в разделе Interface ссылку на функцию DetectCollision. Она будет просчитывать столкновение материальной точки с отдельным треугольником, поэтому в качестве параметров ей необходимо передать:

• Положение материальной точки
• Вектор, на который движется материальная точка
• Координаты вершин треугольника
• Нормаль к поверхности треугольника.

Стоит объяснить последний пункт. На самом деле, нам не потребуется нормаль при просчёте столкновения. Просто нормаль потребуется нам в другой функции, а, раз она уже будет к этому моменту найдёна, так почему бы её не использовать для получения уравнения плоскости, в которую может врезаться материальная точка? Ведь коэффициенты при переменных в уравнении плоскости пропорциональны одноимённым координатам вектора нормали к поверхности. Видимо, нам придётся находить нормаль, и дальше я напомню, как это делается.

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

Найдём углы следующим образом:


Учтём так же тот случай, когда Δx =0, в этом случае получим π=±0.5? в зависимости от знака Δy.

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

Подставим формулы координат x, y, z в первое уравнение и получим расстояние до точки пересечения.

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

Итак, мы нашли расстояние до точки пересечения, и его можно сравнить с длиной вектора и с нулём. Очевидно, что если r меньше нуля, то материальная точка не пересечёт треугольник - она будет двигаться в противоположном направлении. Однако задача не решена. Нужно теперь проверить, будет ли точка пересечения лежать внутри треугольника.

Естественно, что такую проверку проще всего выполнить следующим образом. При попадании должны выполняться условия для углов:

CAM+MAB=CAB
MCA+MCB=ACB

Однако при реализации такой проверки перед нами становится новая проблема. И дело даже не в том, чтобы найти углы. Их можно найти по теореме косинусов, предварительно получив по формуле расстояния между двумя точками длины всех сторон. Дело в том, что при просчете чисел в типе Single, да ещё и с использованием тригонометрических функций, значительную роль начинает играть погрешность. Поэтому нужно либо учитывать её, либо искать другой способ.

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

Введём плоскую систему координат в плоскости треугольника следующим образом. Точка А будет лежать в начале отсчёта, точка B будет находиться где-то на оси абсцисс, а точка C - в первой или второй четверти. Для того, чтобы воплотить эту идею, необходимо найти длины сторон AB, AC, BC, AM, MC, MB - это можно сделать по формуле расстояния между точками в пространстве. Затем по теореме косинусов найдём углы BAM, CAM и BAC. После этого можно найти координаты точек в новой системе координат.

Когда треугольник и точка столкновения лежат в плоскости, можно проверить, лежит ли эта точка внутри треугольника. Этот процесс проще всего реализовать в виде отдельной функции Belong (см. текст модуля).

Текст функции DetectCollision

Перед тем, как ввести новую функцию, необходимо добавить новый тип - вектор.

Type
TVector=record
  x:single;
  y:single;
  z:single;
end;

Далее я привожу текст неописанной функции определения принадлежности точки на плоскости треугольнику. В параметрах, думаю, всё очень понятно кроме последней переменной. Дело в том, что не всегда программа работает в дробных числах так, как должна работать. Очень часто возникают большие погрешности, особенно при работе с тригонометрическими и обратными тригонометрическими функциями. Всвязи с этим программа должна учитывать погрешность измерения длины сторон. Эта погрешность будет определена в функции проверки столкновения, поэтому в данную подпрограмму она передаётся в готовом виде.

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

Function belong(x1,z1,x2,z2,x3,z3,x,z,e:extended):boolean; //Функция возвращает true, если точка лежит внутри треугольника var a, b:double; //Коэффициенты в уравнении сторон треугольника begin belong := True; //Если сторона параллельна оси ординат if x1=x2 then begin if not(((x3>=x1-e)and(x>=x1-e))or((x3<=x1+e)and(x<=x1+e))) then belong:=False; //Иначе, получить уравнение этой прямой end else begin a:=(z2-z1)/(x2-x1); b:=-(a*x1-z1); //И произвести проверку if not(((z3<=a*x3+b+e)and(z<=a*x+b+e))or((z3>=a*x3+b-e)and(z>=a*x+b-e))) then belong:=False; end; if x2=x3 then begin if not(((x1>=x3-e)and(x>=x3-e))or((x1<=x3+e)and(x<=x3+e))) then belong:=False; end else begin a:=(z3-z2)/(x3-x2); b:=-(a*x2-z2); if not(((z1<=a*x1+b+e)and(z<=a*x+b+e))or((z1>=a*x1+b-e)and(z>=a*x+b-e))) then belong:=False; end; if x1=x3 then begin if not(((x2>=x1-e)and(x>=x1-e))or((x2<=x1+e)and(x<=x1+e))) then belong:=False; end else begin a:=(z3-z1)/(x3-x1); b:=-(a*x1-z1); if not(((z2<=a*x2+b+e)and(z<=a*x+b+e))or((z2>=a*x2+b-e)and(z>=a*x+b-e))) then belong:=False; end; end;

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

Второе отличие - это проверка попадания в треугольник. Не смотря на все проблемы с погрешностью, мы вынуждены проверять углы хотя бы при одной вершине, чтобы удостовериться, что точка лежит в первой или второй четвертях. Очевидно, угол, определённый по теореме косинусов, не даёт такой информации. Однако проверять этот факт равенством я так же не собираюсь. Известно, что треугольник целиком находится в первой и второй четвертях. Значит, угол треугольника с вершиной в начале отсчёта, считая против часовой стрелки, не может привысить ста восьмидесяти градусов. Этот факт и проверяет программа.

Function DetectCollision(Normal,Vec0:TVector; vx,vy,vz,x1,y1,z1,x2,y2,z2,x3,y3,z3:single):TVertex; var alpha, beta, //Направление вектора движения l, //Длина вектора движения a,b,c,d, //Коэффициенты в уравнении плоскости x,y,z, //Координаты точки столкновения r, //Расстояние от данной точки до точки столкновения l1,l2,l3,l4,l20, //Длины сторон (см. комментарии ниже) an310,an210,an213, //Углы в треугольнике //Координаты точки в новой системе координат xpl1,xpl2,xpl3,ypl1,ypl2,ypl3,xpl,ypl, da:extended; //Погрешность TV:TVertex; //Точка столкновения Begin //Длина вектора движения l:=sqrt(sqr(vec0.x)+sqr(vec0.y)+sqr(vec0.z)); DetectCollision.x:=0; DetectCollision.y:=0; DetectCollision.z:=0; //Определим длины сторон треугольника l1 :=sqrt(sqr(x3-x2)+sqr(y3-y2)+sqr(z3-z2)); //BC l2 :=sqrt(sqr(x3-x1)+sqr(y3-y1)+sqr(z3-z1)); //AC l3 :=sqrt(sqr(x2-x1)+sqr(y2-y1)+sqr(z2-z1)); //AB //Погрешность da:=((1E-7)/l1+(1E-7)/l2+(1E-7)/l3); //Определим направление вектора движения beta:=arcsin(Vec0.z/(sqrt(sqr(Vec0.x)+sqr(Vec0.y)+sqr(Vec0.z)))); if Vec0.y=0 then begin if Vec0.x>0 then alpha:=pi/2 else alpha:=-pi/2; end else alpha:=arctan(Vec0.x/Vec0.y); if Vec0.x<0 then alpha:=-alpha; //Используя нормаль, найдём коэффициенты в уравнении плоскости a:=Normal.x;b:=Normal.y;c:=Normal.z;d:=-(a*x3)-(b*y3)-(c*z3); //Проверим, будут ли пересекаться продолжение вектора и плоскость if abs((a*cos(beta)*sin(alpha)+b*cos(beta)*cos(alpha)+c*sin(beta)))then exit; //Если да, то определить расстояние до точки пересечения r:=-(d+a*vx+b*vy+c*vz)/(a*cos(beta)*sin(alpha)+ b*cos(beta)*cos(alpha)+c*sin(beta)); //и координаты этой точки x:=vx+r*cos(beta)*sin(alpha); y:=vy+r*cos(beta)*cos(alpha); z:=vz+r*sin(beta); //Если точка пересечения достаточно далеко, // значит столкновения не произойдёт if (r>l) or (r<=da) then exit; //Определим длины оставшихся сторон l4 :=sqrt(sqr( x-x3)+sqr( y-y3)+sqr( z-z3)); //OC l20:=sqrt(sqr( x-x2)+sqr( y-y2)+sqr( z-z2)); //OB l :=sqrt(sqr( x-x1)+sqr( y-y1)+sqr( z-z1)); //OA //И углы при одной вершине по теореме косинусов an310:=arccos((sqr(l)+sqr(l2)-sqr(l4))/(2*l*l2)); an210:=arccos((sqr(l)+sqr(l3)-sqr(l20))/(2*l*l3)); an213:=arccos((sqr(l3)+sqr(l2)-sqr(l1))/(2*l3*l2)); //Если определённый косвенно угол в начале отсчёта больше 180, if an310+an210>pi then exit; //то столкновения не могло быть. xpl1:=0; ypl1:=0; xpl2:=l2; ypl2:=0; xpl3:=l3*cos(an213); ypl3:=l3*sin(an213); xpl:=l*cos(an310); ypl:=l*sin(an310); if not Belong(xpl1,ypl1,xpl2,ypl2,xpl3,ypl3,xpl,ypl,da) then exit; TV.x:=x; TV.y:=y; TV.z:=z; DetectCollision:=TV; end;

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

Отражение от поверхности

Пусть нам известно, что материальная точка, которая движется на данный вектор, сталкивается с поверхностью, нормаль к которой {x_n;y_n;z_n}. Необходимо определить, на какой вектор будет двигаться материальная точка после столкновения. В самом деле, для такого расчета необходимо знать только нормаль к поверхности и данный вектор.

Для тех, кто не знаком с теоретическим материалом, напомню алгоритм нахождения нормали к плоскости методом произведения векторов. Пусть заданы координаты трех точек этой плоскости x₁;y₁;z₁;x₂;y₂;z₂;x₃;y₃;z₃. Введём два вектора, которые будут соединять одну из этих точек с двумя другими. Например,

Векторы, лежащие в плоскости, перпендикулярны вектору нормали, то есть их скалярные произведения равны нулю. Имеем:

Третью координату найдём из оставшегося уравнения.

Примечание

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

Можно, однако, слегка доработать алгоритм нахождения нормали, чтобы с его помощью можно было бы настраивать освещение. Поскольку нам нет разницы, с какой стороны мы смотрим на треугольник, то требуется, чтобы нормаль, проводимая из плоскости, была обращена к нам. Поэтому нужно ввести вектор из точки наблюдения к какой-нибудь точке плоскости (например, к одной из вершин треугольника, или к центру этого треугольника). Затем нужно найти угол между этим вектором и имеющейся нормалью. Если этот угол окажется меньше прямого, то нормаль обращена от нас, и её координаты нужно умножить на -1, а уж затем использовать в процедуре glNormal.

Вообще, удобно находить нормаль единичной длины. Это часто приходит на помощь в самых разных ситуациях, в том числе и в нашей. Найдём длину найденной нормали, а затем поделим её координаты на эту длину, тем самым получив координаты нормали единичной длины.

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

Вектор, который на рисунке обозначен n, на самом деле, не обычная нормаль. Как известно, метод нахождения нормали с помощью векторного произведения требует задания одной координаты нормали произвольно, а получения двух других исходя из заданной. Поэтому совершенно не обязательно, что длина нормали будет так соотноситься с длиной и положением данного вектора, как показано на рисунке. Для этого нам опять же послужит угол α. Найдём угол 180-α, а затем, используя косинус угла, найдём координаты нужной нормали, умножив координаты найденной нормали единичной длины на необходимую величину.

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

Текст функции GetNewV

В качестве параметров функция принимает координаты нормали и координаты данного вектора.

Function GetNewV(Normal,Vec1:TVector):TVector; var nX, nY, nZ, //Координаты нового вектора нормали alpha, //Угол между вектором и нормалью l:single; //Длина данного вектора begin nX:=Normal.x; nY:=Normal.y; nZ:=Normal.z; //Найти длину данного вектора l:=sqrt(sqr(Vec1.x)+sqr(Vec1.y)+sqr(Vec1.z)); //И угол между нормалью и данным вектором alpha:=arccos((nX*Vec1.x+nY*Vec1.y+nZ*Vec1.z)/ (sqrt(sqr(nx)+sqr(ny)+sqr(nz))*l)); //Если угол меньше прямого, перевернуть нормаль if alpha2 then begin nx:=-nx; ny:=-ny; nz:=-nz; end else //В противном случае, перевернём угол alpha:=pi-alpha; //Новые координаты вектора нормали nX:=2*nX*l*cos(alpha); nY:=2*nY*l*cos(alpha); nZ:=2*nZ*l*Cos(alpha); //Сложим вектор нормали с данным вектором GetNewV.x:=Vec1.x+nx; GetNewV.y:=Vec1.y+ny; GetNewV.z:=Vec1.z+nz; end;

Текст функции GetNewVec.

Эта последняя процедура будет осуществлять перебор всех полигонов и проверять, произойдёт ли с ними столкновение, а, если столкновение произойдёт, то она будет пытаться найти новый вектор. Объясню, почему пытаться. Может случиться такая ситуация, когда, отразившись от одной стенки, материальная точка сразу же ударится в другую стенку. Однако нашей процедуры такая ситуация не касается. Я объясню, как можно рассмотреть этот случай при использовании составляемого нами модуля, позже, а сейчас приведу текст функции GetNewVec. Однако предварительно нужно ввести несколько новых типов.

Это ВЕРШИНА (Vertex), ПОЛИГОН (Face), а так же их массивы.

//Вершина, как и вектор имеет три координаты, поэтому определение //вектора очень похоже на определение вершины, однако эти типы нельзя совмещать, //поэтому я ставлю перед TVector слово type Type TVertex= type TVector; PVertexArray=^TVertexArray; TVertexArray=array[word] of TVertex; TFace=array[0..2] of word; PFaceArray=^TFaceArray; TFaceArray=array[word] of TFace; //Скопировано из калькулятора Const pi=3.1415926535897932384626433832795; Procedure GetNewVec(C0:TVertex; //Начальное положение точки Vector:TVector; //Вектор, на который движется точка VertexCount:word; //Количество вершин препятствия Vertices:PVertexArray; //Массив вершин препятствия FacesCount:word; //Количество полигонов препятствия Faces:PFaceArray; //Массив полигонов препятствия Var NewVec:TVector; //Возвращаемый новый вектор движения var NewCoord:TVertex //Возвращаемые координаты точки столкновения );

Текст функции GetNewVec Вы можете увидеть в модуле Collision. Поскольку эта функция очень проста, я не стану приводить её с комментариями.

Теперь пришло время использовать наш модуль для написания какой-нибудь программы.

Использование модуля Collision

Я подготовил один пример использования модуля Collision. В нём на каменную поверхность падают непонятные тела и, отталкиваясь, улетают вниз. Остановлюсь на этом примере более подробно.

Во-первых, в качестве каменной поверхности используется сетка, загружаемая из файла PrevMesh.eye. Этот файл создан с помощью макроса 3dStudioMax MeshExport, описание которого Вы можете найти в моей статье (Вариант реализации простейших костных деформаций).

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

repeat k:=k+1; oldvec:=Drops[i].DropVector; Drops[i].DropCoord.x:=(tvr.x+Drops[i].DropCoord.x)/2; Drops[i].DropCoord.y:=(tvr.y+Drops[i].DropCoord.y)/2; Drops[i].DropCoord.z:=(tvr.z+Drops[i].DropCoord.z)/2; GetNewVec(Drops[i].DropCoord, Drops[i].DropVector, PrevMesh.VertexCount, PrevMesh.Vertices, PrevMesh.FacesCount, PrevMesh.Faces, Drops[i].DropVector,tvr); until ((Drops[i].DropVector.x=0) and (Drops[i].DropVector.x=0) and (Drops[i].DropVector.y=0)) or (k=10);

Этот цикл каждый раз помещает данную материальную точку (Drop) между её текущим положением и точкой столкновения и проверяет, произойдёт ли столкновение, если материальная точка будет двигаться из текущего положения на полученный вектор. Если при движении после столкновения материальная точка вновь встречает препятствие, то программа снова находит новый вектор после столкновения с этим препятствием. Такой цикл повторяется до тех пор, пока точка не находит свободное направление, или пока не зашкаливает счетчик. После десяти повторений почти очевидно, что точка не сможет найти свободное направление и продолжать движение, поэтому она должна быть остановлена.

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

Теперь может возникнуть вопрос, как все эти падающие камешки могут быть связаны с движком? Я отвечаю - непосредственно. Одно из применений данному модулю - просчёт попадания в цель. Можно так же использовать модуль для просчёта трайектории движения, скажем, брошенноу гранаты.

Ещё одно применение, которое хотелось бы придумать модулю - движение камеры в ограниченном полигонами пространстве. В самом деле, раз материальные точки можно заставить прыгать по поверхности, то почему бы к одной из них не приделать камеру? Такое движение будет отвечать всем законам физики, а так же позволит реализовать ряд других возможностей.

В заключении повторю, что организация всех приведённых здесь алгоритмов, а тем более написание движка, немыслимы без методов аналитической геометрии, которая повсеместно применяется при создании компьютерных игр.

К материалу прилагаются файлы:

• Примеры к статье (235 K) обновление от 2/16/2004 10:52:00 AM