As mentioned in the comment above, Jherico's answer really helped me solve this problem.
This answer builts uppon his post, just to make things a bit more complete, for future reference.
Please find the GLM library here, it's a header only library, so you don't need to build this, just link to it and include the necessary header files. It has some great tools, as can be seen below.
#define GLM_FORCE_RADIANS // for compatibility between GLM and OpenGL API legacy functions
#include <glm/glm.hpp>
#include <glm/gtc/type_ptr.hpp> // for glm::mat4_cast // casting quaternion->mat4
#include <glm/gtx/quaternion.hpp> // for glm::rotation
#include <glm/gtx/string_cast.hpp> // for glm::to_string
...
...
...
// _startBoneVec3 used below is Eigen::Vector3d
// _endddBoneVec3 used below is Eigen::Vector3d
...
...
...
//
// for JOINTS
//
glPushMatrix();
glTranslated( _startBoneVec3(0),_startBoneVec3(1),_startBoneVec3(2) );
glutWireSphere(1.0,30,30);
glPopMatrix();
...
...
//
// for BONES
//
glPushMatrix();
float boneLength = ( _endddBoneVec3 - _startBoneVec3 ).norm();
glTranslated( _startBoneVec3(0),_startBoneVec3(1),_startBoneVec3(2) );
glm::vec3 _glmStart( _startBoneVec3(0),_startBoneVec3(1),_startBoneVec3(2) );
glm::vec3 _glmEnddd( _endddBoneVec3(0),_endddBoneVec3(1),_endddBoneVec3(2) );
glm::vec3 _glmDirrr = glm::normalize( _glmEnddd - _glmStart ); // super important to normalize!!!
glm::quat _glmRotQuat = glm::rotation( glm::vec3(0,0,1),_glmDirrr); // calculates rotation quaternion between 2 normalized vectors
//glm::mat4 _glmRotMat = glm::mat4_cast(_glmRotQuat); // quaternion -> mat4
glm::mat4 _glmRotMat = glm::toMat4 (_glmRotQuat); // quaternion -> mat4
//std::cout << glm::to_string(_glmDirrr) << std::endl;
glMultMatrixf(glm::value_ptr(_glmRotMat));
glutWireCone(4.2,boneLength,4,20); // cone with 4 slices = pyramid-like
glPopMatrix();
Result: