Допустим, вы используете небольшую структуру для описания координат, а другую для выровненных по оси ограничительных рамок:
typedef struct {
double x;
double y;
double z;
} vec3d;
typedef struct {
vec3d min;
vec3d max;
} box3d; /* Axis-aligned bounding box */
Допустим, вы моделируете N сферические объекты с радиусами, описанными в double radius[], текущими координатами в vec3d curr[] (или vec3d *curr) и предыдущими координатами в vec3d prev[] (или vec3d *prev).
Мы хотим проверить, не сталкиваются ли или не пересекаются ли какие-либо объекты, но делаем это эффективно.
Чтобы избежать ненужной работы, используйте массив для выровненных по оси ограничивающих рамок, чтобы в нем находился сферический объект с предыдущими и текущими координатами. Два объекта могут тогда сталкиваться или пересекаться только в том случае, если их ограничивающие оси ограничивающие рамки пересекаются:
static inline double dmin(const double d1, const double d2)
{
return (d1 <= d2) ? d1 : d2;
}
static inline double dmax(const double d1, const double d2)
{
return (d1 >= d2) ? d1 : d2;
}
void update_boxes(box3d *box, const size_t count,
const vec3d *curr, const vec3d *prev,
const double *radius)
{
size_t i;
for (i = 0; i < count; i++) {
box[i].min.x = dmin(curr[i].x, prev[i].x) - radius[i];
box[i].max.x = dmax(curr[i].x, prev[i].x) + radius[i];
box[i].min.y = dmin(curr[i].y, prev[i].y) - radius[i];
box[i].max.y = dmax(curr[i].y, prev[i].y) + radius[i];
box[i].min.z = dmin(curr[i].z, prev[i].z) - radius[i];
box[i].max.z = dmax(curr[i].z, prev[i].z) + radius[i];
}
}
Чтобы вычислить прямоугольники, мы берем координаты центра и расширяем прямоугольник на радиус. (Ограничивающие рамки с выравниванием по оси не обязательно должны быть точными; им нужно только покрыть сферы в обоих положениях. Если они больше, это просто означает, что мы сделаем некоторую ненужную работу. Мы используем ограничивающие рамки с выравниванием по оси, потому что проверка, пересекаются ли два из них, очень быстро.)
Если предположить, что два сферических объекта имеют постоянную скорость в течение каждого временного шага, то положение во время t (0 <= t && t <= 1) в течение временного шага равно
xi(t) = (1-t)*prev[i].x + t*curr[i].x;
yi(t) = (1-t)*prev[i].y + t*curr[i].y;
zi(t) = (1-t)*prev[i].z + t*curr[i].z;
Это простая линейная интерполяция в 3D между двумя позициями. t = 0 на предыдущем временном шаге, t = 1 на текущем временном шаге.
Если мы сформируем уравнение квадрата расстояния между двумя такими объектами, при индексах i и k мы получим
L(t) = SQUARE( ((1-t)*prev[i].x + t*curr[i].x) - ((1-t)*prev[k].x + t*prev[k].x) )
+ SQUARE( ((1-t)*prev[i].y + t*curr[i].y) - ((1-t)*prev[k].y + t*prev[k].y) )
+ SQUARE( ((1-t)*prev[i].z + t*curr[i].z) - ((1-t)*prev[k].z + t*prev[k].z) )
, где SQUARE(expr) = (expr)*(expr). Он достигает минимума, когда его производная равна нулю. Если мы решим это для t, мы обнаружим, что существует ровно один реальный корень, то есть два объекта вдоль этих линейных траекторий с постоянными скоростями ближе всего друг к другу в момент времени t:
t = ( (prev[i].x - prev[k].x) * ( (prev[i].x - prev[j].x) - (curr[i].x - curr[k].x) )
+ (prev[i].y - prev[k].y) * ( (prev[i].y - prev[j].y) - (curr[i].y - curr[k].y) )
+ (prev[i].z - prev[k].z) * ( (prev[i].z - prev[j].z) - (curr[i].z - curr[k].z) )
) / ( SQUARE( (prev[i].x - prev[k].x) - (curr[i].x - curr[k].x) )
+ SQUARE( (prev[i].y - prev[k].y) - (curr[i].y - curr[k].y) )
+ SQUARE( (prev[i].z - prev[k].z) - (curr[i].z - curr[k].z) ) )
Это действительно только в том случае, если делитель ненулевой (он никогда не может быть отрицательным, потому что это сумма квадратов).
Нас интересуют только случаи, когда t >= 0 и t <= 1; то есть случаи, когда два объекта находились ближе всего друг к другу между предыдущим и текущим временными шагами.
Если это произойдет, все, что нам нужно сделать, это подключить t обратно к уравнению и сравнить с SQUARE(radius[i] + radius[k]), чтобы выяснить, не столкнулись ли два объекта.
Давайте рассмотрим пример функции, которая также вычисляет ограничивающие прямоугольники по оси, использует их для быстрого выбора и должна правильно определять столкновения:
void handle_collisions(const size_t count,
box3d *const box,
vec3d *const curr,
const vec3d *const prev,
const double *const radius)
{
size_t i, k;
for (k = 0; k < count; k++) {
box[k].min.x = dmin(prev[k].x, curr[k].x) - radius[k];
box[k].max.x = dmax(prev[k].x, curr[k].x) + radius[k];
box[k].min.y = dmin(prev[k].y, curr[k].y) - radius[k];
box[k].max.y = dmax(prev[k].y, curr[k].y) + radius[k];
box[k].min.z = dmin(prev[k].z, curr[k].z) - radius[k];
box[k].max.z = dmax(prev[k].z, curr[k].z) + radius[k];
for (i = 0; i < k; i++) {
if (box[k].min.x <= box[i].max.x &&
box[k].min.y <= box[i].max.y &&
box[k].min.z <= box[i].max.z &&
box[k].max.x >= box[i].min.x &&
box[k].max.y >= box[i].max.y &&
box[k].max.z >= box[i].max.z) {
/* A collision is possible, since the axis-aligned
bounding boxes intersect. Check. */
const vec3d p = { prev[i].x - prev[k].x,
prev[i].y - prev[k].y,
prev[i].z - prev[k].z };
const vec3d d = { p.x - curr[i].x + curr[k].x,
p.y - curr[i].y + curr[k].y,
p.z - curr[i].z + curr[k].z };
const double tn = p.x * d.x + p.y * d.y + p.z * d.z;
const double td = d.x * d.x + d.y * d.y + d.z * d.z;
if (tn >= 0.0 && tn <= td) {
const double t1 = tn / td;
const double t0 = 1.0 - t1;
const vec3d loc_k = { t0*prev[k].x + t1*curr[k].x,
t0*prev[k].y + t1*curr[k].y,
t0*prev[k].z + t1*curr[k].z };
const vec3d loc_i = { t0*prev[i].x + t1*curr[i].x,
t0*prev[i].y + t1*curr[i].y,
t0*prev[i].z + t1*curr[i].z };
const vec3d delta = { loc_i.x - loc_k.x,
loc_i.y - loc_k.y,
loc_i.z - loc_k.z };
if (delta.x*delta.x + delta.y*delta.y + delta.z*delta.z <=
(radius[i] + radius[k])*(radius[i] + radius[k])) {
/* Collision occurs at time t (0 <= t && t <= 1),
between object k (at loc_k) and object i (at loc_i).
*/
}
}
}
}
}
}
Технически, более двух объектов могут столкнуться в течение одного временного шага, хотя это чрезвычайно редко. Дело в том, что если вы «удалите» любой из сталкивающихся объектов из рассмотрения во время вышеуказанного цикла, вы можете пропустить их. (Я подозреваю, что это достаточно редко, чтобы не волноваться, но я параноик и не люблю не беспокоиться о вещах.:)
Чтобы избежать этого, я просто сохранил бы коллизии в вышеупомянутом цикле в некоторый массив или в структуру данных с несвязанным множеством . Затем, после приведенного выше кода, объедините все сталкивающиеся объекты. (Обратите внимание, что структура данных несвязанного набора делает это проще, потому что, если у вас есть коллизии AB, BC и CD, непересекающийся набор фактически разрешает это в AB, AC и AD. В противном случае, если вы присоединяете B к A и уничтожаете B, соединение C с B. невозможно. Опять же, это угловой случай, но не стоит беспокоиться о таких угловых случаях, потому что существует разница между наличием надежного симулятора и симулятора, который дает сбой в каждой тысяче симуляций или, возможно, в миллиардных шагах. по невидимой непонятной причине.)