Vectors.cpp

//---------------------------------------------------------------------------

#pragma hdrstop

#include "Vectors.h"
#include "math.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)


IntersectionResult Intersect(Segment3D a, Segment3D b, Vector3D& result)
{
// вначале проверяем вектор на параллельность – 
// векторное умножение даст нуль-вектор
	Vector3D n = CrossProduct(a.GetVector(), b.GetVector());
	if (fabs(n.X) < null_vector_e && fabs(n.Y) < null_vector_e && fabs(n.Z) < null_vector_e) {
		return IntersectionResult::Parallel;
	}
// определяем координаты начальной точки отрезка В в новой системе координат
// Z координата даст показатель непланарности, что позволит сразу определить
// отсутствие пересечения
	Vector3D sb_start = a.TransformSystem(b,b.start);
	if (sb_start.Z > off_palne_tolerance)
	{
		return IntersectionResult::NotInPlane;
	}
// определяем координаты конца отрезка в новой СК
	Vector3D sb_end = a.TransformSystem(b,b.end);
// если координаты Y начала и конца имеют один знак, значит обе точки лежат по 
// одну сторону от отрезка A и пересечения нет
	if (std::signbit(sb_start.Y) == std::signbit(sb_end.Y) ) {

		return IntersectionResult::NoIntersection;
	}
// находим пересечение с осью X отрезка 
	double x_intersection_2d = Intersect2D(sb_start.X, sb_start.Y, sb_end.X, sb_end.Y);
	double length = a.GetVector().Length();
// если пересечение оказалось за пределами отрезка A, пересечение не засчитывается
	if (x_intersection_2d < 0 || x_intersection_2d > length) {
		return IntersectionResult::NoIntersection;
	}
// если пересечение в пределах отрезка А, находим точку на отрезке А соответветсвующую пересечению
// результат будет уже в исходной системе координат  
	result = a.start + (a.GetVector() * (x_intersection_2d / length));
	return IntersectionResult::Ok;
}

//эта функция определяет ортонормированную систему координат, где x параллельна текущему отрезку
//z перпендикулярна текущему отрезку и отрезку segment, y вычисляется на основе z и x
//функция возвращает координаты точки vector в полученной системе координат
//подразумевается, что проверка на параллельность выполена ранее

Vector3D Segment3D::TransformSystem(Segment3D segment, Vector3D vector)
{

	Vector3D n = CrossProduct(this->GetVector(), segment.GetVector());
    Vector3D new_z = n.Normalize();
	Vector3D new_x = this->GetVector().Normalize();
	Vector3D new_y = CrossProduct (new_z, new_x);

    //           А - матрица перехода       АТ - транспонированная матрица перехода 
	//      |  new_x.X new_y.X new_z.X |        |new_x.X new_x.Y new_x.Z|
	// A =  |  new_x.Y new_y.Y new_z.Y |   AT = |new_y.X new_y.Y new_y.Z|
	//      |  new_x.Z new_y.Z new_z.Z |        |new_z.X new_z.Y new_z.Z|

    // получаем вектор с учетом нуля новой системы координат (начало отрезка А)
	Vector3D p = vector - start;
	Vector3D result;

    // Умножаем вектор на матрицу АТ
	result.X = new_x.X*p.X + new_x.Y*p.Y + new_x.Z*p.Z;
	result.Y = new_y.X*p.X + new_y.Y*p.Y + new_y.Z*p.Z;
	result.Z = new_z.X*p.X + new_z.Y*p.Y + new_z.Z*p.Z;


	return result;
}

// эта функция ищет пересечение бесконечной кривой, проведенной через точки А и В с осью Х
//подразумевается, что проверка на параллельность выполена ранее
double Intersect2D (double ax, double ay, double bx, double by)
{
	double k = (bx-ax)/(by-ay);
	return ax-k*ay;
}


// Далее стандартные функции, необходимые для работы с векторами
double Vector3D::Length()
{
return sqrt(pow(X,2)+pow(Y,2)+pow(Z,2));
}

Vector3D Vector3D::Normalize()
{
	Vector3D vector;
	double l = this->Length();
	vector.X = X/l;
	vector.Y = Y/l;
	vector.Z = Z/l;
   return vector;
}

Vector3D CrossProduct(Vector3D a, Vector3D b)
{
	Vector3D v;
	v.X = a.Y*b.Z - a.Z*b.Y;
	v.Y = a.Z*b.X - a.X*b.Z;
	v.Z = a.X*b.Y - a.Y*b.X;
	return v;
}

Vector3D Vector3D::operator-(Vector3D other)
{
	Vector3D v;
	v.X = X - other.X;
	v.Y = Y - other.Y;
	v.Z = Z - other.Z;
	return v;
}

Vector3D Vector3D::operator+(Vector3D other)
{
	Vector3D v;
	v.X = X + other.X;
	v.Y = Y + other.Y;
	v.Z = Z + other.Z;
	return v;
}

Vector3D Vector3D::operator*(double value)
{
	Vector3D v;
	v.X = X*value;
	v.Y = Y*value;
	v.Z = Z*value;
	return v;
}

Vector3D Segment3D::GetVector()
{
	Vector3D v;
	v.X = end.X - start.X;
	v.Y = end.Y - start.Y;
	v.Z = end.Z - start.Z;
	return v;
}