[Insight-developers] [Insight-users] Convert Rigid3DTransform To VersorRigid3DTransform
Bradley Lowekamp
blowekamp at mail.nih.gov
Thu Dec 18 11:46:36 EST 2008
Hello,
It's been a while a while since I used a versor/quaternion much so I
am not completely sure of some of how this will be used. So this is
the bit of code:
const T epsilon = vcl_numeric_limits<T>::epsilon();
double trace = m(0,0) + m(1,1) + m(2,2) + 1.0;
if( trace > epsilon)
{
const double s = 0.5 / vcl_sqrt(trace);
m_W = 0.25 / s;
m_X = (m(2,1) - m(1,2)) * s;
m_Y = (m(0,2) - m(2,0)) * s;
m_Z = (m(1,0) - m(0,1)) * s;
}
else ..
Let me see what I can remember from numerical analysis. lf we assume
that the trace is 2*eps then we will only have about 2 bits of
precision, since 1.0 is in the computation of trace. This does not
look like enough precision. The other methods may be better before
this point. Perhaps eps should be multiplied by the number of bits
needed of the computation. Hrmm.. I am uncertain, but I think some of
that computation is only occurring as floats too.
Just some thought for consideration :)
Brad
On Dec 17, 2008, at 11:41 PM, Hans Johnson wrote:
> 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
>
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
========================================================
Bradley Lowekamp
Lockheed Martin Contractor for
Office of High Performance Computing and Communications
National Library of Medicine
blowekamp at mail.nih.gov
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/mailman/private/insight-developers/attachments/20081218/913ea29f/attachment.htm>
More information about the Insight-developers
mailing list