Problem mit Arcball in OpenGL

[Ein ehemaliger Benutzer]

Hallo Leute,
versuche seit geraumer Zeit nen Arcball zu implementieren.

Zunächst zur Theorie:
Wenn ich die rechte Maustaste drücke wird von der Clickposition aus ein Strahl in die Szene geschickt und der Schnittpunkt mit dem Arcball bestimmt. Bewege ich nun mit gedrückter Taste die Maus, so wird ein weiterer Strahl-Arcball-Schnittpunkt bestimmt.

Der Winkel zwischen den beiden Vektoren, die vom Mittelpunkt des Arcballs zu den beiden Schnittpunkten führen, ist mein Rotationswinkel und das Kreuzprodukt aus den beiden Vektoren liefert die Rotationsachse.
Mit Hilfe des Winkels und der Achse bestimme ich ein Quaternion mit dessen Hilfe ich wiederum eine Rotationsmatrix aufstelle.

Ich verwende gluLookat für die Kamera.

gluLookAt(eye_x,eye_y,eye_z,		pos_x, pos_y, pos_z,	0, 1, 0);

Deshalb dachte ich mir ich könnte einfach den Augenpunkt mit der Rotationsmatrix multiplizieren.
Es funktioniert auch wunderbar bei horizontalen Bewegungen.

Das Problem ist, das bei vertikalen Bewegungen mit der Maus, die erstellte Matrix zusätzlich zur Rotation eine Skalierung zur Folge hat.

Was mache ich falsch ? Kann mal bitte jemand im Quellcode schauen was da schief läuft ? habe so gut es geht kommentiert.
Ich finde den Fehler einfach nicht.

// Hier wird die Kamera-Rotation mit Quaternionen umgesetzt
	if(right_mousedown && hitSphere){

	GLdouble model_view[16];
	glGetDoublev(GL_MODELVIEW_MATRIX, model_view);

	GLdouble projection[16];
	glGetDoublev(GL_PROJECTION_MATRIX, projection);

	GLint viewport[4];
	glGetIntegerv(GL_VIEWPORT, viewport);

	//Liefert die 3d-Koordinaten der Maus-Bildschirm-Koordinaten
	gluUnProject(mousex, mousey, 0.01,
	model_view, projection, viewport,
	pos3D_x, pos3D_y, pos3D_z);

	//Aufstellen des Strahls
	Ray r;
	r.direction.x = *pos3D_x - eye_x;
	r.direction.y = *pos3D_y - eye_y;
	r.direction.z = *pos3D_z - eye_z;

	r.direction.normalizeVector();

	r.origin.x = eye_x;
	r.origin.y = eye_x;
	r.origin.z = eye_x;

	//Nachfolgend Bestimmung des Schnittpunktes
	Vec temp_a(*pos3D_x,*pos3D_y,*pos3D_z);
	Vec temp_b(trackball->center.x, trackball->center.y, trackball->center.z);
	Vec dst = Vec::offsetVector(temp_a,temp_b);
	float B = Vec::dot(dst, r.direction);
	float C = Vec::dot(dst, dst) - pow(trackball->radius,2);
	float D = B*B - C;

	real intersection;
	if(D > 0){
	intersection = -B - sqrt(D); 
	inSphere = true;
	}
	else{
	intersection = std::numeric_limits<float>::infinity();
	inSphere = false;
	}

	if(inSphere){
	current_intersection->x = *pos3D_x + intersection * r.direction.x;
	current_intersection->y = *pos3D_y + intersection * r.direction.y;
	current_intersection->z = *pos3D_z + intersection * r.direction.z;

	//Ermittle die Verbindungsvektoren zwischen Pivot = (pos_x, pos_y, pos_z) und den beiden Mausvektoren
	//current_intersection und init_intersection
	Vec m_a, m_b;

	m_a.x = current_intersection->x - pos_x;
	m_a.y = current_intersection->y - pos_y;
	m_a.z = current_intersection->z - pos_z;

	m_b.x = init_intersection->x - pos_x;
	m_b.y = init_intersection->y - pos_y;
	m_b.z = init_intersection->z - pos_z;

	m_a.normalizeVector();
	m_b.normalizeVector();

	//Berechne das Kreuzprodukt aus den beiden Verbindungsvektoren um die Rotationsachse zu erhalten
	rotationAxis = Vec::cross(m_a,m_b);

	//Berechne den Winkel zwischen den beiden Verbindungsvektoren.
	rotationAngle = acos(Vec::dot(m_a,m_b)/(Vec::magnitude(m_a)*Vec::magnitude(m_b)));

	//Stelle das Einheitsquaternion auf. Falls Winkel zu klein stelle das Null-Quaternion auf um eine Identitätsmatrix zu erhalten
	if(rotationAngle>1.0e-8){
	quat.W = cos(rotationAngle/2);
	quat.X = rotationAxis.x * sin(rotationAngle/2);
	quat.Y = rotationAxis.y * sin(rotationAngle/2);
	quat.Z = rotationAxis.z * sin(rotationAngle/2);
	}
	else{
	quat.W = 0;
	quat.X = 0;
	quat.Y = 0;
	quat.Z = 0;
	}

	quat.Normalize();

	//Bilde aus dem Einheitsquaternion die Rotationsmatrix
	rot_eye[0] = pow(quat.W,2) + pow(quat.X,2) - pow(quat.Y,2) - pow(quat.Z,2);
	rot_eye[1] = 2 * quat.X * quat.Y + 2 * quat.W * quat.Z;
	rot_eye[2] = 2 * quat.X * quat.Z + 2 * quat.W * quat.Y;

	rot_eye[3] = 2 * quat.W * quat.Z + 2 * quat.X * quat.Y;
	rot_eye[4] = pow(quat.W,2) - pow(quat.X,2) + pow(quat.Y,2) - pow(quat.Z,2);
	rot_eye[5] = 2 * quat.Y * quat.Z - 2 * quat.W * quat.X;

	rot_eye[6] = 2 * quat.X * quat.Z - 2 * quat.W * quat.Y;
	rot_eye[7] = 2 * quat.W * quat.X + 2 * quat.Y * quat.Z;
	rot_eye[8] = pow(quat.W,2) - pow(quat.X,2) - pow(quat.Y,2) + pow(quat.Z,2);

	//Wende die Rotationsmatrix auf den Augenpunkt an////////////////////
	real temp_x, temp_y, temp_z;

	temp_x = eye_x;
	temp_y = eye_y;
	temp_z = eye_z;

	eye_x = rot_eye[0] * temp_x + rot_eye[1] * temp_y + rot_eye[2] * temp_z;
	eye_y = rot_eye[3] * temp_x + rot_eye[4] * temp_y + rot_eye[5] * temp_z;
	eye_z = rot_eye[6] * temp_x + rot_eye[7] * temp_y + rot_eye[8] * temp_z;
}

Danke

Blue

[Ein ehemaliger Benutzer]

Kann mir eventuell jemand seine eigene ArcBall-Implementierung schicken ?
Die Sachen, die ich im Netz gefunden habe, helfen mir leider nicht weiter ...