А.Панов. Вычисление ОКРУЖНОСТИ в 3D системе координат

avs.chat.ru/sphere3d.htm
avs.chat.ru/sphere3d.avs [AVS WinAMP]

Цель: Построение модели 3D сцены для кругового обзора.
Задача: Вычисление окружности, которая образуется пересечением произвольной плоскости со сферой.

Сначала по заданному вектору точки вращения P(xp,yp,zp) вычисляется радиус опорной сферы rs и параметр секущей плоскости pc для заданного вектора нормали N(xn,yn,zn). "Неопределённость" может возникнуть, если вектор нормали имеет нулевое значение (nr=0), поэтому хотя бы одна из составляющих вектора нормали всегда должна оставаться не равной нулю. Затем вычисляется радиус секущей окружности rc и координаты её центра R(xr,yr,zr). Ось абсцисс аффинной системы координат (fi=0) задаётся вектором I(xi-xr,yi-yr,zi-zr). Точка отсчёта начала движения (fi=fp) задаётся по окружности (траектории) вектором P(xp-xr,yp-yr,zp-zr). Скалярное произведение двух радиус-векторов P и I позволяет получить модуль угла начала движения |fp| через арккосинус. Смешанное произведение трех векторов I, N и P даёт возможность полностью определить fp. Векторное произведение векторов I и N даёт направление оси ординат - вектор J(xj,yj,zj). Скалярное произведение векторов P и J определяет проекцию рассчётного вектора P на ось ординат, и, соответственно, знак fp. В процессе вычислений может происходить деление на ноль, когда нулю равен радиус секущей окружности (rc=0). Практически это означает, что точка объекта находится непосредственно на оси вращения, и поэтому не перемещается в пространстве. Такие точки необходимо исключать из пересчёта вращения. При вычислении вращения конкретного 3D объекта, траектории вычисляются отдельно для его каждой узловой точки. Все плоскости вращения удерживаются единым вектором нормали. Алгоритм обеспечивает корректный пересчёт углов вращения (собственного вращения, прецессии и нутации) и позволяет избежать "шарнирного замка", возникающего при использовании эйлеровых углов в классической 3D матрице вращения. Тем самым отменяется необходимость в использовании кватернионов при построении 3D сцены кругового движения и круговом обзоре трёхмерного объекта.

Сфера (0,0,0)
  x^2 + y^2 + z^2 = rs^2
Определение секущей плоскости по точке на плоскости и вектору нормали
  xn*(x-xp) + yn*(y-zp) + zn*(z-zp) = 0
Общее (полное) уравнение плоскости
  xn*x + yn*y + zn*z + ( - xn*xp - yn*zp - zn*zp) = 0
Нормированное уравнение секущей плоскости
  x*cos(oxn) + y*cos(oyn) + z*cos(ozn) - рc = 0
Направляющие косинусы
  nx = cos(oxn); ny = cos(oyn); nz = cos(ozn);
Единичный нормальный вектор секущей плоскости
  nx^2 + ny^2 + nz^2 = 1
Точка перемещается по секущей окружности
  rs^2 = rc^2 + pc^2

Скалярное произведение векторов P(xp-xr,yp-yr,zp-zr) и I(xi-xr,yi-yr,zi-zr)
  rc^2 * cos(fc) = (xp-xr)*(xi-xr) + (yp-yr)*(yi-yr) + (zp-zr)*(zi-zr)
   ip = rc * cos(fc) ; x=>ip ;
Смешанное произведение трех векторов
  Векторное произведение I(xi-xr,yi-yr,zi-zr) и N(xn,yn,zn)
   xj=(yi-yr)*zn-(zi-zr)*yn; yj=(zi-zr)*xn-(xi-xr)*zn; zj=(xi-xr)*yn-(yi-yr)*xn;
    rg=sqrt(xj*xj+yj*yj+zj*zj);
  Скалярное произведение P(xp-xr,yp-yr,zp-zr) и J(xj,yj,zj)
   rc*rg*cos(pi/2-fp) = (xp-xr)*xj + (yp-yr)*yj + (zp-zr)*zj
    jp = rc * cos(pi/2-fp) ; y=>jp ;

Вращение точки вокруг оси вектора нормали ON
  fi=i*2*pi; x1=r*cos(fi); y1=r*sin(fi); z1=0;
Матрица вращения вокруг оси OY
  x2=x1*cos(tetta)-z1*sin(tetta);
  z2=x1*sin(tetta)+z1*cos(tetta);
  y2=y1;
Матрица вращения вокруг оси OZ
  x3=xr+x2*cos(omega)-y2*sin(omega);
  y3=yr+x2*sin(omega)+y2*cos(omega);
  z3=zr+z2;

*

xp=-0.5; yp=-0.5; zp=-0.5;     // rs = 0.87;
zn = -0.05; xn = getkbmouse(1)*km; yn = getkbmouse(2)*km;
xn=getkbmouse(2)*km; yn=-getkbmouse(1)*km; // xn=0; yn=0; zn=1;
zn=zn+(getkbmouse(3)-getkbmouse(4))*kd;

// Сфера проходит через точку (xp,yp,zp) //
rs = sqrt(xp*xp+yp*yp+zp*zp);
// Нормирующий множитель //
mu = 1/sqrt(xn*xn+yn*yn+zn*zn);
// Единичный нормальный вектор секущей плоскости // nr=1 (!)
nx = xn*mu; ny = yn*mu; nz = zn*mu;    // nr = nx*nx + ny*ny + nz*nz;
// Расстояние до секущей плоскости //
pc = xp*nx + yp*ny + zp*nz;
// Центр секущей окружности //
xr = pc*nx; yr = pc*ny; zr = pc*nz;
// Радиус секущей окружности //
rc = sqrt(rs*rs - pc*pc);
// Угол вектора нормали от оси OZ (ozn) //
tetta = - acos(nz);
// Угол проекции нормали от оси OX //
omega = atan2(ny,nx);

// Вычисление оси абсцисс fi=0 //
xi=xr+rc*cos(omega)*cos(tetta); yi=yr+rc*sin(omega)*cos(tetta); zi=zr+rc*sin(tetta);
// Скалярное произведение // rc|=0 (!)
ip = ((xp-xr)*(xi-xr)+(yp-yr)*(yi-yr)+(zp-zr)*(zi-zr))/rc;
// Смешанное произведение // rg|=0 (!)
xj =(yi-yr)*zn-(zi-zr)*yn; yj=(zi-zr)*xn-(xi-xr)*zn; zj=(xi-xr)*yn-(yi-yr)*xn;
  rg = sqrt(xj*xj+yj*yj+zj*zj);
jp = ((xp-xr)*xj + (yp-yr)*yj + (zp-zr)*zj)/rg;
fp = - atan2(jp,ip);

// Ограничители деления на ноль //
fp = if( below(rc,0.001),0,fp ); // fp = if( below(rg,0.001),0,fp );

fi=fp+i*2*pi;

// Вычисление поворота //
x=xr+rc*(cos(fi)*cos(omega)*cos(tetta)-sin(fi)*sin(omega));
y=yr+rc*(cos(fi)*sin(omega)*cos(tetta)+sin(fi)*cos(omega));
z=zr+rc*(cos(fi)*sin(tetta));