gvdb.SetTransform() does not appear to calculate the inverse transform correctly

[digbeta]

I was debugging some issues getting a volume to spin around its center of mass and noted that when rendering in GVDB, the volume was orbiting the origin instead of spinning. As a sanity check, I passed the same computed mXform to a mesh object in the scene and it spun correctly and it did not orbit.

The orbiting behavior of the volume looked like what you might see if rotation was happening after the translation (assuming in the inverse transform because the mXform appeared to be correct). As I began to look into the transform more closely, I wanted to double check the mInvXform that was being calculated in SetTransform(), so I added this to the end of the function:

	// Check inverse
	Matrix4F should_be_ident = mXform;
	should_be_ident *= mInvXform;

And found that I was not getting an identity matrix from the result. So, I modified SetTransform as follows. Note I left the mInvXrot unchanged as there were issues I couldn't immediately diagnose, but I changed the mInvXform to calculate the inverse matrix using the RotateTZYXS function and a new tmp2 variable to make things readable below. As a result, I was now getting an identity matrix from the newly calculated mInvXform and the spinning behavior works correctly for volumes:

	mPretrans = pretrans;
	mScale = scal;
	mAngs = rotate;
	mTrans = trans;

	// p' = T R S PT p				pretrans -> scale -> rotate -> translate
	mXform.RotateTZYXS(mAngs, mTrans, mScale);
	mXform.PreTranslate(mPretrans);

	// p = PT^-1 S^-1 R^-1 T^-1 p'	inv trans -> inv rot -> inv scale -> inv pretrans
	mInvXform.Translate(-mPretrans.x, -mPretrans.y, -mPretrans.z);	// PT^-1
	Matrix4F tmp, tmp2;
	tmp.RotateTZYX(mAngs, Vector3DF(0, 0, 0));
	tmp.InvertTRS();
	tmp2.RotateTZYXS(mAngs, mTrans, mScale); // <----------------new 
	tmp2.InvertTRS();    // <----------------new 
	mInvXrot.Scale(1.0 / mScale.x, 1.0 / mScale.y, 1.0 / mScale.z);		// S ^-1
	mInvXrot *= tmp;												// R ^-1
	mInvXform *= tmp2;   // <----------------new 
	//mInvXform += Vector3DF(-mTrans.x, -mTrans.y, -mTrans.z); // <---------------- new

	// Check inverse
	Matrix4F should_be_ident = mXform;
	should_be_ident *= mInvXform;

This seems to have broken some things in the Optix rendering that I am still investigating.

Comments and thoughts are welcomed here. Thanks!

[ramakarl]

Ah ha, nice work!

I was investigating this as well and came to the same conclusion. The SetTransform was not computing those correctly. Here is my updated code for comparison:

void VolumeGVDB::SetTransform(Vector3DF pretrans, Vector3DF scal, Vector3DF angs, Vector3DF trans ) 
{
	mPretrans = pretrans;
	mScale = scal;
	mAngs = angs;
	mTrans = trans;

	// rotation only
	mXrot.RotateTZYX(mAngs, Vector3DF(0, 0, 0));		
	mInvXrot = mXrot;
	mInvXrot.InvertTRS();

	// p' = T S R PT p	APPLY_ORDER: pretrans -> rotate -> scale -> translate
	mXform.Identity();
	mXform.Translate ( mTrans );	
	mXform.Scale ( mScale );
	mXform.Rotate ( mXrot );
	mXform.Translate ( mPretrans );

	// p = PT^-1 R^-1 S^-1  T^-1 p'  APPLY_ORDER: inv trans -> inv scale -> inv rot -> inv pretrans
	mInvXform.Identity();
	mInvXform.InvTranslate ( mPretrans );	
	mInvXform.InvRotate ( mXrot );		
	mInvXform.InvScale ( mScale );	
	mInvXform.InvTranslate ( mTrans );
}

This should be a lot cleaner/clearer. One thing that bothered me for a while was that the Matrix/Vector classes didn't have composing operators. I have added those and used them here. They apply the correct matrix multiplications to clean up the code a lot.

//------------------------------------------- Composing operations
// M' = M*T
Matrix4F& Matrix4F::Translate ( const Vector3DF& t )
{
	data[12] += (VTYPE) t.x;
	data[13] += (VTYPE) t.y;
	data[14] += (VTYPE) t.z;
	return *this;
}
// M' = M*R
Matrix4F& Matrix4F::Rotate ( const Matrix4F& rot )
{
	Multiply ( rot );
	return *this;
}
// M' = M*S    -- actual a matrix-vector multiplication
Matrix4F& Matrix4F::Scale ( const Vector3DF& s )
{
	data[0] *= (VTYPE) s.x;	data[1] *= (VTYPE) s.y;	data[2] *= (VTYPE) s.z;	
	data[4] *= (VTYPE) s.x;	data[5] *= (VTYPE) s.y;	data[6] *= (VTYPE) s.z;	
	data[8] *= (VTYPE) s.x;	data[9] *= (VTYPE) s.y;	data[10] *= (VTYPE) s.z;
	data[12] *= (VTYPE) s.x; data[13] *= (VTYPE) s.y; data[14] *= (VTYPE) s.z;
	return *this;
}
// M' = M * T^-1
Matrix4F& Matrix4F::InvTranslate ( const Vector3DF& t )
{
	data[12] -= (VTYPE) t.x;
	data[13] -= (VTYPE) t.y;
	data[14] -= (VTYPE) t.z;
	return *this;
}
Matrix4F& Matrix4F::InvRotate ( Matrix4F mtx )
{
	mtx.InvertTRS();
	Multiply ( mtx );
	return *this;
}
Matrix4F& Matrix4F::InvScale ( const Vector3DF& s )
{	
	data[0] /= (VTYPE) s.x;	data[1] /= (VTYPE) s.y;	data[2] /= (VTYPE) s.z;	
	data[4] /= (VTYPE) s.x;	data[5] /= (VTYPE) s.y;	data[6] /= (VTYPE) s.z;	
	data[8] /= (VTYPE) s.x;	data[9] /= (VTYPE) s.y;	data[10] /= (VTYPE) s.z;
	data[12] /= (VTYPE) s.x; data[13] /= (VTYPE) s.y; data[14] /= (VTYPE) s.z;
	return *this;
}

Regarding the OptiX rendering, what I have now is:

RT_PROGRAM void vol_intersect( int primIdx )
{
	float3 hit = make_float3(NOHIT,NOHIT,NOHIT);	
	float3 norm = make_float3(0,0,0);	
	float4 hclr = make_float4(1,1,1,1);
	float t;
	if (ray.ray_type == MESH_RAY ) return;

	//-- Ray march	
	float3 orig = mmult(ray.origin, SCN_INVXFORM);
	float3 dir = mmult(normalize(ray.direction), SCN_INVXROT);	
	rayCast ( &gvdbObj, gvdbChan, orig, dir, hit, norm, hclr, raySurfaceTrilinearBrickF4 );	
	if ( hit.z >= NOHIT ) return;
	t = length ( hit - ray.origin );

	// report intersection to optix
	if ( rtPotentialIntersection( t ) ) {	

		shading_normal = norm;
		geometric_normal = shading_normal;
		front_hit_point = mmult(hit, SCN_XFORM); 
		back_hit_point  = mmult(hit, SCN_XFORM); 
		deep_color = hclr;		

		rtReportIntersection( mat_id );
	}
}

This works correctly now, including shadows, with OptiX.

Nice find on SetTransform! Glad we found that independently as confirmation.

[digbeta]

Haha - thank you, Rama! I'm glad you could confirm it -- I wanted to make really sure I was interpreting things correctly. It looks like you fixed a bunch of things today; I can test all of these out in the next day or so and then close them out/confirm the fixes if that works for you (including depth merge, which I am also using). Thanks again for digging in and providing the updates!

[digbeta]

Hi, Rama - I am still getting orbiting behavior (for basic meshes/models and volumes that use the computed matrix using this code) and when checking the inverse:

			// Check inverse
			should_be_ident = mXform;
			should_be_ident *= mInvXform;

I don't get an identity matrix using the new code...

I also get different mXform values when compared against this code:

		mXform.RotateTZYXS(mRotation, mTranslation, mScale);
		mXform.PreTranslate(mPreTranslation);

Which appears to give the desired behavior. I'm still poking around, but wanted to share the results so you knew there seems to still be some issues with the transform...

Is InvertTRS() still valid when called using a matrix where translate and scale have not been set?

[ramakarl]

I will be traveling the next few weeks, so I won't be able to take a look for a while.

It is concerning that mXform * mInvXform doesn't result in identity, so that would be interesting to understand better. The difference between the new mXform and the older RotateTZYXS /w PreTrans may be due to the order of operations. My latest code is doing p'=T S R PT p, where as the original may have been doing p'=T R S PT p.. The ordering of the pretrans, scale, rotate, translate is a matter of API design, but of course in either case mXform * mInvXform should be identity.

[digbeta]

Not a problem, I appreciate the earlier confirmation as it helps me focus in on the problem spots. The process is still useful as it helps me understand the library better as I navigate through the design. I'll let you know what I find out! Thanks again and safe travels.

[ramakarl]

I think I have this fixed now..

Here is my latest SetTransform:

void VolumeGVDB::SetTransform(Vector3DF pretrans, Vector3DF scal, Vector3DF angs, Vector3DF trans ) 
{
	mPretrans = pretrans;
	mScale = scal;
	mAngs = angs;
	mTrans = trans;

	// rotation only
	Matrix4F mXrot;
	mXrot.RotateTZYX(mAngs, Vector3DF(0, 0, 0));		
	
	mInvXrot.Identity();	
	mInvXrot.InvRotate ( mXrot );
	mInvXrot.InvScale ( mScale );           // scale here supports non-uniform voxels

	// p' = PT S R T           RIGHT-TO-LEFT-APPLY: trans -> rotate -> scale -> pre-translate
	mXform.Identity();
	mXform.Translate ( mPretrans );	
	mXform.Scale ( mScale );
	mXform.Rotate ( mXrot );	
	mXform.Translate ( mTrans );

	// p = T^-1 R^-1 S^-1 PT^-1 p'		RIGHT-TO-LEFT-APPLY: inv pre-trans -> inv scale -> inv rot -> inv trans
	mInvXform.Identity();
	mInvXform.InvTranslate ( mTrans );		
	mInvXform.InvRotate ( mXrot );		
	mInvXform.InvScale ( mScale );	
	mInvXform.InvTranslate ( mPretrans );

	Matrix4F test;
	test = mXform;
	test *= mInvXform;
	test.Print ();
}

The test at the end confirms mXform * mInvXform is identity.

It was also necessary to modify the camera ray generation in cuda_gvdb_geom.cuh:

inline __device__ float3 getViewPos ()
{
	return mmult ( scn.campos, SCN_INVXFORM );  
}

inline __device__ float3 getViewRay ( float x, float y )
{
  #ifdef CUDA_PATHWAY
	float3 v = x*scn.camu + y*scn.camv + scn.cams;
  #else
	float3 v = make_float3(0,0,0);
  #endif
  return normalize( mmult(v, SCN_INVXROT) );
}

[NBickford-NV]

So, I'm looking at integrating this into GVDB at the moment, and I think there's confusion inside the library about whether Matrix4F objects are row-major vs. column-major, and whether vectors are column vectors or row vectors. I'm leaning towards column-major and column vectors, but does anyone have a preference?

To sum it up (no pun intended!), most GVDB functions and comments seem to refer to column-major matrices and column vectors. For instance, RotateTZYX represents rotation around the X axis, followed by rotation around the Y axis, followed by rotation around the Z axis, followed by translation - i.e. it returns the matrix

T * Z * Y * X

where T is a translation matrix, and Z, Y, and X are rotation matrices. This is natural for multiplying by a vector on the right. Additionally, comments throughout the codebase such as SetTransform above seem to point to multiplying by a vector on the right. If we multiplied by vectors on the left instead, then the matrix TZYX would actually be a translation, followed by rotation around Z, then Y, then X - so in this case, a better function if we were to multiply on the left instead would be RotateXYZT.

However, there are also parts in the codebase that seem to imply row-major matrices with row vectors. For instance, mmult(v, m) represents either multiplying a row vector v by a row-major matrix m, or multiplying a column-major matrix m with a column matrix v. The current semantics for vector-matrix multiplication for Vector3DF are v *= m, instead of v = m * v. Finally - and this would probably have the most impact - Matrix4F's GetF and GetRowVec both seem to assume row-major. (I don't know if either of these is used in external code; GetRowVec is used to fill camivprow0...camivprow3, but these seem to be unused).

I feel like the best thing to do would be to define that Matrix4Fs are column-major and that vectors are column vectors; this would match comments, mean that there's no need to rename functions such as RotateTZYX, and follow how OpenGL documents its matrices. Internally, the order of arguments for mmult would be reversed, and I'd probably remove GetF, GetRowVec, and camivprov0...camivprow3. (I could leave v *= m as is, to mean v = m * v.) However, I could also do row-major and row vectors instead (like DirectX), in which case I'd probably go through, rename functions such as RotateTZYX, and fix comments to specify that vectors are row vectors. Let me know if anyone would prefer one or the other.

Also, I'm looking at removing vec.h, vec.cpp, camera.h, and camera.cpp from sample_utils, since a version of them is part of the GVDB core and accessible to samples (e.g. gvdb_vec.h), but the two versions have diverged slightly in their list of features. This might also allow me to remove some preprocessor directives in vec.h.

[ramakarl]

Hi Nick,
You said: "I feel like the best thing to do would be to define that Matrix4Fs are column-major and that vectors are column vectors"

That would be my inclination as well. I generally conceived of gvdb following the opengl paradigm with column-major notation, with the vector on right, and matrix multiply order right-to-left.

As you say many funcs are documented this way, and the naming of funcs follows this.
I think it should be fine to remove GetF and GetRowVec as these are rarely, if ever, used.
We should also check that operator() and operator[ ] for single elements follow the desired notation.

You are probably already aware, but this is keeping in mind that column-major is only about notation w opengl, not storage ordering. From OpenGL docs: " rght-side vector, with pre-multiplied matrices, in column-major order has an identical memory layout to left-side vector, with post-multiplied matrices, in row-major order."

So, in column-major notation used by opengl and gvdb, the translation elements are still in the 13th, 14th and 15th memory slots of the 16 contiguous elements.

The operator v * = m was for performance reasons, as it avoids the extra temporary variable generated by compiler when one does v = m * v. I could see adding a new func: vector3df& matrix4f::operator* (const vector3df op), which follows the proper column-major ordering notation. The func v *= m could be kept around for performance with a clear note as to why.

Lately I favor functions being closer to the desired notation ordering (easier to read). Hence the matrix compositing operations I've added and made use of in SetTransform.

PS. Another important change imo for user-friendliness. Matrix4F constructor:
Matrix4F() { Identity(); } // not all zeros

[NBickford-NV]

Hi all,

I think this should be fixed in the latest version now! Please see the summary of changes at #89 (comment).

Thanks,

[digbeta]

This has been fixed with all the related changes - thank you both!