Качественный расчет TBN

struct Vertex
{
	Vector3  position;	// позиция вершины
	Vector2  texcoord;	// текстурные координаты
	Vector3  normal;	// нормаль к вершине
	Vector3  tangent;	// тангент
	Vector3  binormal;	// бинормаль
};

struct Triangle
{
	int a, b, c;	// индексы в массив вершин, образующие треугольник
};

void CalculateTangentsAndBinormals( int vertexCount, Vertex* vertexes, int numFaces, Triangle* faces )
{
	register int i, j;
	std::vector<Vector3> tangents, binormals;	// касательные вектора к треугольнику
 
	// проходимся по всем треугольникам, и для каждого считаем базис
	for ( i = 0; i < numFaces; i++ )
	{
		int ind0 = faces[ i ].a;
		int ind1 = faces[ i ].b;
		int ind2 = faces[ i ].c;
 
		Vector3 v1 = vertexes[ ind0 ].position;
		Vector3 v2 = vertexes[ ind1 ].position;
		Vector3 v3 = vertexes[ ind2 ].position;
		float s1      = vertexes[ ind0 ].texcoord.x;
		float t1      = vertexes[ ind0 ].texcoord.y;
		float s2      = vertexes[ ind1 ].texcoord.x;
		float t2      = vertexes[ ind1 ].texcoord.y;
		float s3      = vertexes[ ind2 ].texcoord.x;
		float t3      = vertexes[ ind2 ].texcoord.y;
 
		Vector3  t, b;
		CalcTriangleBasis( v1, v2, v3, s1, t1, s2, t2, s3, t3, t, b );
		tangents.push_back( t );
		binormals.push_back( b );
	}
 
	// теперь пройдемся по всем вершинам, для каждой из них найдем
	// грани ее содержащие, и запишем все это хозяйство на будущее =)
	for ( i = 0; i < vertexCount; i++ )
	{
		std::vector<Vector3> rt, rb;
		for ( j = 0; j < numFaces; j++ )
		{
			if ( faces[ j ].a == i || faces[ j ].b == i || faces[ j ].c == i )
			{
				// нашли грань которой эта вершина принадлежит.
				// добавим вектора этой грани в наш массив
				rt.push_back( tangents[ j ] );
				rb.push_back( binormals[ j ] );
			}
		}
 
		// все прилежащие вектора нашли, теперь просуммируем их
		// и разделим на их количество, т.е. усредним.
		Vector3 tangentRes( 0, 0, 0 );
		Vector3 binormalRes( 0, 0, 0 );
		for ( j = 0; j < rt.size(); j++ )
		{
			tangentRes += rt[ j ];
			binormalRes += rb[ j ];
		}
		tangentRes /= float( rt.size() );
		binormalRes /= float( rb.size() );
 
		// а теперь то, о чем многие забывают. Как мы помним,
		// TBN базис представляет собой всего-навсего систему координат.
		// поэтому все три направляющие вектора этого базиса
		// обязаны быть попарно перпендикулярны. Вот об этом мы
		// сейчас и побеспокоимся, выполнив ортогонализацию
		// векторов методом Грама-Шмидта
 
		tangentRes = Ortogonalize( vertexes[ i ].normal, tangentRes );
		binormalRes = Ortogonalize( vertexes[ i ].normal, binormalRes );
 
		// вот и все, теперь просто запишем результат
 
		vertexes[ i ].tangent =  tangentRes;
		vertexes[ i ].binormal =  binormalRes;
	}
}

void CalcTriangleBasis( const Vector3& E, const Vector3& F, const Vector3& G, float sE,
		float tE, float sF, float tF, float sG, float tG, Vector3& tangentX,
		Vector3& tangentY )
{
	Vector3 P = F - E;
	Vector3 Q = G - E;
	float s1 = sF - sE;
	float t1 = tF - tE;
	float s2 = sG - sE;
	float t2 = tG - tE;
	float pqMatrix[2][3];
	pqMatrix[0][0] = P[0];
	pqMatrix[0][1] = P[1];
	pqMatrix[0][2] = P[2];
	pqMatrix[1][0] = Q[0];
	pqMatrix[1][1] = Q[1];
	pqMatrix[1][2] = Q[2];
	float temp = 1.0f / ( s1 * t2 - s2 * t1);
	float stMatrix[2][2];
	stMatrix[0][0] =  t2 * temp;
	stMatrix[0][1] = -t1 * temp;
	stMatrix[1][0] = -s2 * temp;
	stMatrix[1][1] =  s1 * temp;
	float tbMatrix[2][3];
	// stMatrix * pqMatrix
	tbMatrix[0][0] = stMatrix[0][0] * pqMatrix[0][0] + stMatrix[0][1] * pqMatrix[1][0];
	tbMatrix[0][1] = stMatrix[0][0] * pqMatrix[0][1] + stMatrix[0][1] * pqMatrix[1][1];
	tbMatrix[0][2] = stMatrix[0][0] * pqMatrix[0][2] + stMatrix[0][1] * pqMatrix[1][2];
	tbMatrix[1][0] = stMatrix[1][0] * pqMatrix[0][0] + stMatrix[1][1] * pqMatrix[1][0];
	tbMatrix[1][1] = stMatrix[1][0] * pqMatrix[0][1] + stMatrix[1][1] * pqMatrix[1][1];
	tbMatrix[1][2] = stMatrix[1][0] * pqMatrix[0][2] + stMatrix[1][1] * pqMatrix[1][2];
	tangentX.Set( tbMatrix[0][0], tbMatrix[0][1], tbMatrix[0][2] );
	tangentY.Set( tbMatrix[1][0], tbMatrix[1][1], tbMatrix[1][2] );
	tangentX.Normalize();
	tangentY.Normalize();
}

Vector3 ClosestPointOnLine( const Vector3& a, const Vector3& b, const Vector3& p )
{
	Vector3 c = p - a;
	Vector3 V = b - a;
	float d = V.Length();
	V.Normalize();
	float t = DotProduct( V, c );		// скалярное произведение векторов
 
	// проверка на выход за границы отрезка
	if ( t < 0.0f )
		return a;
	if ( t > d )
		return b;
 
	// Вернем точку между a и b
	V *= t;
	return ( a + V );
}

Vector3 Ortogonalize( const Vector3& v1, const Vector3& v2 )
{
	Vector3 v2ProjV1 = ClosestPointOnLine( v1, -v1, v2 );
	Vector3 res = v2 - v2ProjV1;
	res.Normalize();
	return res;
}