使用Eigen函数库做球面四边形插值的方法

最新推荐文章于 2024-04-25 11:27:18 发布

牧童027

最新推荐文章于 2024-04-25 11:27:18 发布

阅读量2.9k

点赞数

文章标签： 3D

本文链接：https://blog.csdn.net/zgce/article/details/8856968

版权

正在将一段D3D的代码改造成OpenGL的代码，在原来D3D中使用了D3DXQuaternionSquadSetup与D3DXQuaternionSquad两个D3D的数学函数。

其中D3DXQuaternionSquadSetup这个函数的MSDN解释如下：

void D3DXQuaternionSquadSetup(
     __in  D3DXQUATERNION *pAOut, 
     __in  D3DXQUATERNION *pBOut,  
     __in  D3DXQUATERNION *pCOut, 
     __in  const D3DXQUATERNION *pQ0, 
     __in  const D3DXQUATERNION *pQ1, 
     __in  const D3DXQUATERNION *pQ2, 
     __in  const D3DXQUATERNION *pQ3);

q0 = |Q0 + Q1| < |Q0 - Q1| ? -Q0 : Q0

q2 = |Q1 + Q2| < |Q1 - Q2| ? -Q2 : Q2

q3 = |Q2 + Q3| < |Q2 - Q3| ? -Q3 : Q3

AOut = q1 * e^{[-0.25 *( Ln[Exp(q1)*q2] + Ln[Exp(q1)*q0] ) ]}

BOut = q2 * e^{[-0.25 *( Ln[Exp(q2)*q3] + Ln[Exp(q2)*q1] ) ]}

COut = q2

但通过代码测试，发现按上面的公式实现出来的三个控制点与D3DXQuaternionSquadSetup函数输出的三个点中qA与qB完全不一样。

后面又经过查询了OGRE的代码，发现另有一个计算公式：

let p = point[i], pInv = p.Inverse
tangent[i] = p * exp( -0.25 * ( log(pInv * point[i+1]) + log(pInv * point[i-1]) ) )

按照这个公式编写代码，输出的结果与D3DXQuaternionSquadSetup输出一致，代码如下：

#include "stdafx.h"
#include <iostream>
#include <Eigen/Eigen.hpp>
#include <D3DX10math.h>

static float msEpsilon = 1e-04f;
void QuaternionExp(Eigen::Quaternionf & kResult, Eigen::Quaternionf const & q) {
	// If q = A*(x*i+y*j+z*k) where (x,y,z) is unit length, then
	// exp(q) = cos(A)+sin(A)*(x*i+y*j+z*k).  If sin(A) is near zero,
	// use exp(q) = cos(A)+A*(x*i+y*j+z*k) since A/sin(A) has limit 1.
	float fAngle = Eigen::internal::sqrt(q.x() * q.x() + q.y() * q.y() + q.z() * q.z());
	float fSin = Eigen::internal::sin(fAngle);

	kResult.w() = Eigen::internal::cos(fAngle);

	if(Eigen::internal::abs(fSin) > msEpsilon) {
		float fCoeff = fSin / (fAngle);
		kResult.x() = fCoeff * q.x();
		kResult.y() = fCoeff * q.y();
		kResult.z() = fCoeff * q.z();
	}
	else {
		kResult.x() = q.x();
		kResult.y() = q.y();
		kResult.z() = q.z();
	}
}

void QuaternionLog(Eigen::Quaternionf & kResult, Eigen::Quaternionf const & q) {
	// If q = cos(A)+sin(A)*(x*i+y*j+z*k) where (x,y,z) is unit length, then
	// log(q) = A*(x*i+y*j+z*k).  If sin(A) is near zero, use log(q) =
	// sin(A)*(x*i+y*j+z*k) since sin(A)/A has limit 1.
	kResult.w() = 0.0f;

	if(Eigen::internal::abs(q.w()) < 1.0f) {
		float fAngle = Eigen::internal::acos(q.w());
		float fSin = Eigen::internal::sin(fAngle);
		if(Eigen::internal::abs(fSin) >= msEpsilon) {
			float fCoeff = fAngle / fSin;
			kResult.x() = fCoeff * q.x();
			kResult.y() = fCoeff * q.y();
			kResult.z() = fCoeff * q.z();
			return;
		}
	}

	kResult.x() = q.x();
	kResult.y() = q.y();
	kResult.z() = q.z();
}

Eigen::Quaternionf operator+ (Eigen::Quaternionf const & lq, Eigen::Quaternionf const & rq) {
	return Eigen::Quaternionf(lq.w()+rq.w(), lq.x()+rq.x(), lq.y()+rq.y(), lq.z()+rq.z());
}

Eigen::Quaternionf operator- (Eigen::Quaternionf const & lq, Eigen::Quaternionf const & rq) {
	return Eigen::Quaternionf(lq.w()-rq.w(), lq.x()-rq.x(), lq.y()-rq.y(), lq.z()-rq.z());
}

Eigen::Quaternionf operator* (Eigen::Quaternionf const & lq, float const & rf) {
	return Eigen::Quaternionf(lq.w()*rf, lq.x()*rf, lq.y()*rf, lq.z()*rf);
}

Eigen::Quaternionf operator* (float const & rf, Eigen::Quaternionf const & lq) {
	return Eigen::Quaternionf(lq.w()*rf, lq.x()*rf, lq.y()*rf, lq.z()*rf);
}

void QuaternionSquadSetup(
						  Eigen::Quaternionf & qa,
						  Eigen::Quaternionf & qb,
						  Eigen::Quaternionf & qc,
						  Eigen::Quaternionf const & Q0,
						  Eigen::Quaternionf const & Q1,
						  Eigen::Quaternionf const & Q2,
						  Eigen::Quaternionf const & Q3)
{
	// let p = point[i], pInv = p.Inverse
	// tangent[i] = p * exp( -0.25 * ( log(pInv * point[i+1]) + log(pInv * point[i-1]) ) )
	Eigen::Quaternionf q0 = 
		Eigen::internal::abs((Q0 + Q1).norm()) < Eigen::internal::abs((Q0 - Q1).norm()) 
		? (Eigen::Quaternionf(Q0.w()*(-1), Q0.x()*(-1), Q0.y()*(-1), Q0.z()*(-1))) : Q0;

	Eigen::Quaternionf q2 = 
		Eigen::internal::abs((Q1 + Q2).norm()) < Eigen::internal::abs((Q1 - Q2).norm()) 
		? (Eigen::Quaternionf(Q2.w()*(-1), Q2.x()*(-1), Q2.y()*(-1), Q2.z()*(-1))) : Q2;

	Eigen::Quaternionf q3 = 
		Eigen::internal::abs((Q2 + Q3).norm()) < Eigen::internal::abs((Q2 - Q3).norm()) 
		? (Eigen::Quaternionf(Q3.w()*(-1), Q3.x()*(-1), Q3.y()*(-1), Q3.z()*(-1))) : Q3;
	Eigen::Quaternionf q1 = Q1;

	Eigen::Quaternionf q1Inv = q1.inverse();
	Eigen::Quaternionf logt1, logt2, expt;
	QuaternionLog(logt1, (q1Inv * q2));
	QuaternionLog(logt2, (q1Inv * q0));
	QuaternionExp(expt, (-0.25) * (logt1 + logt2));
	qa = q1 * expt;

	Eigen::Quaternionf q2Inv = q2.inverse();
	QuaternionLog(logt1, (q2Inv * q3));
	QuaternionLog(logt2, (q2Inv * q1));
	QuaternionExp(expt, (-0.25) * (logt1 + logt2));
	qb = q2 * expt;
	qc = q2;
}

void QuaternionSquad(
		   Eigen::Quaternionf & kResult, 
		   Eigen::Quaternionf const & q1, 
		   Eigen::Quaternionf const & qa, 
		   Eigen::Quaternionf const & qb,
		   Eigen::Quaternionf const & qc,
		   float const t)
{
	Eigen::Quaternionf tq1 = q1.slerp(t, qc);
	Eigen::Quaternionf tq2 = qa.slerp(t, qb);
	kResult = tq1.slerp((2.0f * t * (1.0f - t)), tq2);
}

int _tmain(int argc, _TCHAR* argv[]) {
	float q1[4] = { -0.302491f, -0.138833f, 0.053131f, 0.988886f };
	float q2[4] = { 0.109162f, -0.195262f, 0.070018f, 0.978206f };
	float q3[4] = { -0.137615f, 0.797013f, 0.145604f, 0.959571f };
	float q4[4] = { -0.170421f, 0.268557f, 0.178246f, 0.931162f };
	float t = 0.5f;


	Eigen::Quaternionf qf1(q1), qf2(q2), qf3(q3), qf4(q4);
	Eigen::Quaternionf qa, qb, qc, qr;

	D3DXQUATERNION d1(q1), d2(q2), d3(q3), d4(q4);
	D3DXQUATERNION da, db, dc, dr;

	D3DXQuaternionLn(&da, &d1);
	Eigen::Quaternionf lr; QuaternionLog(lr, qf1);

	D3DXQuaternionExp(&db, &d4);
	Eigen::Quaternionf er; QuaternionExp(er, qf4);


	D3DXQuaternionSquadSetup(&da, &db, &dc, &d1, &d2, &d3, &d4);
	QuaternionSquadSetup(qa, qb, qc, qf1, qf2, qf3, qf4);
	std::cout << "D3D Qa: \n" << Eigen::Vector4f((float*)da).transpose() << std::endl;
	std::cout << "Eigen Qa: \n" << Eigen::Vector4f(qa.x(), qa.y(), qa.z(), qa.w()).transpose() << std::endl;

	std::cout << std::endl;
	std::cout << "D3D Qb: \n" << Eigen::Vector4f((float*)db).transpose() << std::endl;
	std::cout << "Eigen Qb: \n" << Eigen::Vector4f(qb.x(), qb.y(), qb.z(), qb.w()).transpose() << std::endl;

	std::cout << std::endl;
	std::cout << "D3D Qc: \n" << Eigen::Vector4f((float*)dc).transpose() << std::endl;
	std::cout << "Eigen Qc: \n" << Eigen::Vector4f(qc.x(), qc.y(), qc.z(), qc.w()).transpose() << std::endl;

	D3DXQuaternionSquad(&dr, &d1, &da, &db, &dc, t);
	QuaternionSquad(qr, qf1, qa, qb, qc, t);

	std::cout << std::endl;
	std::cout << "D3D Qr: \n" << Eigen::Vector4f((float*)dr).transpose() << std::endl;
	std::cout << "Eigen Qr: \n" << Eigen::Vector4f(qr.x(), qr.y(), qr.z(), qr.w()).transpose() << std::endl;


	return 0;
}

输出

D3D Qa:
0.264123 -0.449604 0.0453855 0.858993
Eigen Qa:
0.264123 -0.449604 0.0453855 0.858993

D3D Qb:
-0.166177 1.02695 0.120697 0.706634
Eigen Qb:
-0.166177 1.02695 0.120697 0.706634

D3D Qc:
-0.137615 0.797013 0.145604 0.959571
Eigen Qc:
-0.137615 0.797013 0.145604 0.959571

D3D Qr:
-0.0803391 0.363393 0.106954 1.02754
Eigen Qr:
-0.0803391 0.363393 0.106954 1.02754

前后找来找去花了俺四五个小时，又一次被MS给坑了呀~

牧童027

关注

0
点赞
踩
2

收藏

觉得还不错? 一键收藏
0
评论
使用Eigen函数库做球面四边形插值的方法

正在将一段D3D的代码改造成OpenGL的代码，在原来D3D中使用了D3DXQuaternionSquadSetup与D3DXQuaternionSquad两个D3D的数学函数。其中D3DXQuaternionSquadSetup这个函数的MSDN解释如下：void D3DXQuaternionSquadSetup( __in D3DXQUATERNION *pAOut,
复制链接

扫一扫