[Insight-developers] [Insight-users] Convert Rigid3DTransform To VersorRigid3DTransform
Hans Johnson
hans-johnson at uiowa.edu
Wed Dec 17 23:41:09 EST 2008
All,
I wanted to follow up a bit on what I've determined with regards to
converting between itkVersor and Rotation Matrices (and back again). The
itk implmentation had the threshold for using alternate computations to
compute the conversions near the singularities that arise near roations of
pi as 1e-30, but at that very small value, the numerical stability of the
standard calculation was not sufficient to compute proper matricies. Upon
reading
http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternio
nToEuler/ I determined that the threshold needed to be changed.
A bug has been submitted to http://public.kitware.com/Bug/view.php?id=8312
I've changed this threshold to vcl_numeric_limits<double>::epsilon(), and
wrote a test to confirm that the proper desired behavior was occuring.
/cvsroot/Insight/Insight/Code/Common/itkVersor.txx,v <--
Code/Common/itkVersor.txx
new revision: 1.29; previous revision: 1.28
/cvsroot/Insight/Insight/Testing/Code/Common/itkVersorTest.cxx,v <--
Testing/Code/Common/itkVersorTest.cxx
new revision: 1.20; previous revision: 1.19
============================================
I believe that this change (and the test case that demonstrates it is stable
and reliable) allows itkVersorRigid3DTransform.h SetRotationMatrix function
to be declared as a safe way to initialize a the Versor.
Any comments would be greatly appreciated.
Thanks,
Hans
On 12/17/08 7:32 AM, "Hans Johnson" <hans-johnson at uiowa.edu> wrote:
> Bill,
>
> Thanks. I have that page printed on my desk and have been referencing it.
>
> The code I have supplied in my previous message is producing a
> VersorRigid3DTransformation that produces different resampling results than
> the original Rigid3DTransform.
>
> ========= From itkVersor.h ============
> 205 /** Set the versor using an orthogonal matrix.
> 206 * Based on code from:
> 207 * http://www.euclideanspace.com/maths/geometry/rotations/
> 208 * conversions/matrixToQuaternion/index.htm
> 209 */
> 210 void Set( const MatrixType & m );
>
> ======== From itkVersorRigid3DTransform.h ================
> 116
> 117 /** This method must be made protected here because it is not a safe
> way of
> 118 * initializing the Versor */
> 119 virtual void SetRotationMatrix(const MatrixType & matrix)
> 120 { this->Superclass::SetRotationMatrix( matrix ); }
>
>
> The comment in line 117 of itkVersorRigid3DTransform.h is concerning to me.
> I would think that this would be a safe way to initialize the rotation part
> of the itkVersor3DTransform, as long as the input matrix was tested to be
> orthogonal.
>
> ======== From itkVersorRigid3DTransform.h ================
> 70 // Check if input matrix is orthogonal to within tolerance
> 71 template<class TScalarType>
> 72 bool
> 73 Rigid3DTransform<TScalarType>
> 74 ::MatrixIsOrthogonal(
> 75 const MatrixType & matrix,
> 76 double tolerance )
> 77 {
> 78 typename MatrixType::InternalMatrixType test =
> 79 matrix.GetVnlMatrix() * matrix.GetTranspose();
> 80
> 81 if( !test.is_identity( tolerance ) )
> 82 {
> 83 return false;
> 84 }
> 85
> 86 return true;
> 87 }
>
>
> The itkVersorRigid3DTransform has 6 total parameters. The last three are
> the translation, and the first 3 are the rotation. The versor is a unit
> quaternion, so one of the 4 quaternion parameters can be computed from the
> other 3. I am not sure if I am placing the proper 3 parameters in the first
> three locations of the itk::VersorRigid3DTransform.
>
> I would like to see the following code work, or an explanation of why this
> will not work:
>
> itk::VersorRigid3DTransform<double>::Pointer
> ConvertToVersorRigid3D( itk::Rigid3DTransform<double>::Pointer RT)
> {
> if(! RT->MatrixIsOrthogonal(RT->GetRotationMatrix()))
> {
> throw exception...
> }
> itk::VersorRigid3DTransform<double>::Pointer returnVersorRigid3D=
> itk::VersorRigid3DTransform<double>::New();
> returnVersorRigid3D->SetFixedParameters(RT->GetFixedParameters());
> returnVersorRigid3D->SetRotationMatrix( RT->GetRotationMatrix() );
> returnVersorRigid3D->SetTranslation( RT->GetTranslation() );
> return returnVersorRigid3D;
> }
>
>
>
> Thanks,
> Hans
>
>
>
>
> On 12/16/08 11:02 PM, "Bill Lorensen" <bill.lorensen at gmail.com> wrote:
>
>> Hans,
>>
>> Here's a wikipedia entry:
>> http://en.wikipedia.org/wiki/Conversion_between_quaternions_and_Euler_angles
>>
>> Bill
>>
>> On Tue, Dec 16, 2008 at 10:20 PM, Hans Johnson <hans-johnson at uiowa.edu>
>> wrote:
>>> Hello All,
>>>
>>> I have a program that generates rho, phi, theta rotation euler angles, and I
>>> can easily generate the rotation matrix for creating the Rigid3DTransform.
>>>
>>> I now need to convert Rigid3DTransform to VersorRigid3DTransform so that I
>>> can initialize a mutual information registration with it.
>>>
>>> Could anyone provide some advice on how to perform this conversion?
>>>
>>> Thanks,
>>> Hans
>>>
>>> PS: The following code does not work.
>>> ===========================
>>> typedef itk::Rigid3DTransform<double> RigidTransformType;
>>> typedef itk::VersorRigid3DTransform<double> VersorTransformType;
>>>
>>> VersorTransformType::Pointer
>>> ConvertToVersorRigid3D(RigidTransformType::Pointer RT)
>>> {
>>> VersorTransformType::Pointer VT=VersorTransformType::New();
>>> VT->SetFixedParameters(RT->GetFixedParameters());
>>>
>>> itk::Matrix<double,3,3> R=RT->GetRotationMatrix();
>>> RigidTransformType::TranslationType T=RT->GetTranslation();
>>>
>>> VersorTransformType::ParametersType p;
>>> p.SetSize(6);
>>> itk::Versor<double> v;
>>> v.Set(R);
>>> //Get the first 3 elements of the versor;
>>> p[0]=v.GetRight()[0];
>>> p[1]=v.GetRight()[1];
>>> p[2]=v.GetRight()[2];
>>> p[3]=T[0];
>>> p[4]=T[1];
>>> p[5]=T[2];
>>> VT->SetParameters(p);
>>> return VT;
>>> }
>>>
>>>
>>> Notice: This UI Health Care e-mail (including attachments) is covered by the
>>> Electronic Communications Privacy Act, 18 U.S.C. 2510-2521, is confidential
>>> and may be legally privileged. If you are not the intended recipient, you
>>> are hereby notified that any retention, dissemination, distribution, or
>>> copying of this communication is strictly prohibited. Please reply to the
>>> sender that you have received the message in error, then delete it. Thank
>>> you.
>>>
>>> _______________________________________________
>>> Insight-users mailing list
>>> Insight-users at itk.org
>>> http://www.itk.org/mailman/listinfo/insight-users
>>>
>>>
>
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
More information about the Insight-developers
mailing list