[Insight-developers] itk::Vector Normalize method
Bradley Lowekamp
blowekamp at mail.nih.gov
Thu Jan 14 20:22:52 EST 2010
Hello Arnaud,
The behavior of a divide by zero is defined by the IEEE floating point (754-1985). Here a related issue:
http://standards.ieee.org/reading/ieee/interp/754-1985.html
Specifically a 0/0 should produce a nan, but a (non-zero)/0 should produce an inf.
The issue at hand is the changing ( redefining or defining) of this special case for this methods. As far as I can tell I don't think one approach is more correct then the other. That is using IEEE floating point error handling or returning a zero vector (null vector). As there may be code relying on the current behavior, I don't see a reason to change it for what appears to be an arbitrary expected behavior
In many cases where I have thought about using the Normalize methods, I have instead used a GetNorm, comparison, and then division myself. This is because a small 2-norm can be very erroneous and produce rather inexact results. If the the 2-norm is below some small eps threshold, then I try to handle it based on the context of the calculation. It sounds like this approach should work for your calculation as well.
Thank you for running the performance comparison. My concern with introducing the conditional, is that is will prevent vectorization of loops. Also it would not surprise me if many implementations of sqrt already do a similar check. As for your performance test, you said that you built it in debug? That would not be a very good way to compare, as these performance comparisons must be done with optimizations turned on. Also the time of the random is much greater then the zero case, I presume that the drand32 is outside the loop. If that is the case, then it would be an indication that some optimization is already done inside the sqrt method for the zero case. It is an interesting results. Would you mind sharing the test code?
Thanks,
Brad
On Jan 14, 2010, at 5:25 PM, Arnaud Gelas wrote:
> I have compared the current implementation of GetNorm and the following implementation:
>
> double sq_norm = v.GetSquaredNorm();
> return ( sq_norm > 0. ) ? static_cast< RealValueType >( vcl_sqrt( sq_norm ) ) : 0.;
>
> I called 10^9 times both methods (in Debug mode, compiled with gcc 4.4.2, on linux 64 bits machine)
> * all vectors are NULL vector
> * current implementation -> 42.670 s
> * minor optimization -> 29.983 s
> * all vectors are randomly taken using vnl_random for each component, i.e. using drand32( -1., 1. )
> * current implementation -> 2 min 50.855 s
> * minor optimization -> 2 min 50.857 s
>
> In the general case (where vectors are non null), both implementation achieve similar performances.
> However whenever there are null vectors, this little modification really improves the performance.
>
> I have looked at the implementation of Normalize method in vnl, vtk.
> In both case they return a null vector whenever you provide a null vector as input.
>
> I have a question: is it normal that 0 divided by 0 returns "inf" ? Are there any examples of this implementation?
> or where is it implemented?
>
> Arnaud
>
> On Jan 14, 2010, at 1:36 PM, Bradley Lowekamp wrote:
>
>> I object to these changes.
>>
>> Why do you expect the normalize value of a zero vector to be 0 as apposed to "inf"? Can you not check for the "inf" value with vnl_math_isinf() function?
>>
>> Why do you think you GetNorm methods improves performance? Did you actually check?
>>
>> I also strongly object to exception being thrown in low level numerical methods.
>>
>>
>> Brad
>>
>> On Jan 14, 2010, at 1:17 PM, Arnaud Gelas wrote:
>>
>>> I would agree with you about throwing an exception in the case of NULL
>>> vector.
>>> However, when computing mesh normals on a mesh with boundaries it
>>> sometimes
>>> happens on some point of the boundary. In such a case, I am not sure
>>> an exception
>>> should be thrown.
>>>
>>> But how should we define normalization for NULL vectors?
>>> (I have just checked in VTK, the normalization of a NULL vector is a
>>> NULL vector and it returns the norm).
>>>
>>> Can I proceed and test these changes? any suggestions / thoughts?
>>>
>>> Arnaud
>>>
>>> On Jan 14, 2010, at 12:43 PM, Luis Ibanez wrote:
>>>
>>>> Hi Arnaud,
>>>>
>>>> All your suggestions sound very reasonable to me.
>>>>
>>>> I would just add that when checking for NULL in the
>>>> Normalize() method we probably should use the
>>>> Assert and Throw Exception macro (or a variant of it).
>>>>
>>>> In your current suggested modification, calling Normalize
>>>> in a Null vector will still return a null vector.... well... I guess
>>>> it comes down to picking a mathematical definition of what
>>>> should be the normalization of a null vector...
>>>>
>>>>
>>>> I would vote for throwing an Exception in that case....
>>>>
>>>>
>>>> Any thoughts ?
>>>>
>>>>
>>>> Luis
>>>>
>>>>
>>>> --------------------------------
>>>> On Thu, Jan 14, 2010 at 10:47 AM, Gelas, Arnaud Joel Florent
>>>> <Arnaud_Gelas at hms.harvard.edu> wrote:
>>>>> Hi guys,
>>>>>
>>>>> I was going through the implementation of the normalize method of
>>>>> itk::Vector, and I was wondering if there is any reason not to test
>>>>> if the norm is null (see the implementation below)?
>>>>> I guess it can create some problems when trying to normalize null
>>>>> vector...
>>>>>
>>>>> /**
>>>>> * Divide vector's components by vector's norm
>>>>> */
>>>>> template<class T, unsigned int TVectorDimension>
>>>>> void
>>>>> Vector<T, TVectorDimension>
>>>>> ::Normalize( void )
>>>>> {
>>>>> const RealValueType norm = this->GetNorm();
>>>>> for( unsigned int i=0; i<TVectorDimension; i++)
>>>>> {
>>>>> (*this)[i] = static_cast<T> (static_cast<RealValueType>((*this)
>>>>> [i]) / norm);
>>>>> }
>>>>> }
>>>>>
>>>>> What about the following modification?
>>>>>
>>>>> /**
>>>>> * Divide vector's components by vector's norm
>>>>> */
>>>>> template<class T, unsigned int TVectorDimension>
>>>>> void
>>>>> Vector<T, TVectorDimension>
>>>>> ::Normalize( void )
>>>>> {
>>>>> const RealValueType norm = this->GetNorm();
>>>>> if( norm > 0 )
>>>>> {
>>>>> for( unsigned int i=0; i<TVectorDimension; i++)
>>>>> {
>>>>> (*this)[i] = static_cast<T> (static_cast<RealValueType>((*this)
>>>>> [i]) / norm);
>>>>> }
>>>>> }
>>>>> }
>>>>>
>>>>> I also think that the method GetNorm could be slightly optimized in
>>>>> the case where the norm is null.
>>>>>
>>>>> /**
>>>>> * Returns vector's Euclidean Norm
>>>>> */
>>>>> template<class T, unsigned int TVectorDimension>
>>>>> typename Vector<T, TVectorDimension>::RealValueType
>>>>> Vector<T, TVectorDimension>
>>>>> ::GetNorm( void ) const
>>>>> {
>>>>> return RealValueType( vcl_sqrt(double(this->GetSquaredNorm()) ));
>>>>> }
>>>>>
>>>>> modification
>>>>> /**
>>>>> * Returns vector's Euclidean Norm
>>>>> */
>>>>> template<class T, unsigned int TVectorDimension>
>>>>> typename Vector<T, TVectorDimension>::RealValueType
>>>>> Vector<T, TVectorDimension>
>>>>> ::GetNorm( void ) const
>>>>> {
>>>>> double sq_norm = this->GetSquaredNorm();
>>>>> if( sq_norm > 0. )
>>>>> {
>>>>> return RealValueType( vcl_sqrt( sq_norm ) );
>>>>> }
>>>>> else
>>>>> {
>>>>> return static_cast<RealValueType >( sq_norm );
>>>>> }
>>>>> }
>>>>>
>>>>> or
>>>>>
>>>>> double sq_norm = this->GetSquaredNorm();
>>>>> return ( sq_norm > 0. ) : static_cast<RealValueType
>>>>>> ( vcl_sqrt( sq_norm ) ), sq_norm;
>>>>>
>>>>> Any objection?
>>>>>
>>>>> Arnaud
>>>>> _______________________________________________
>>>>> Powered by www.kitware.com
>>>>>
>>>>> Visit other Kitware open-source projects at
>>>>> http://www.kitware.com/opensource/opensource.html
>>>>>
>>>>> Kitware offers ITK Training Courses, for more information visit:
>>>>> http://kitware.com/products/protraining.html
>>>>>
>>>>> Please keep messages on-topic and check the ITK FAQ at:
>>>>> http://www.itk.org/Wiki/ITK_FAQ
>>>>>
>>>>> Follow this link to subscribe/unsubscribe:
>>>>> http://www.itk.org/mailman/listinfo/insight-developers
>>>>>
>>>
>>> _______________________________________________
>>> Powered by www.kitware.com
>>>
>>> Visit other Kitware open-source projects at
>>> http://www.kitware.com/opensource/opensource.html
>>>
>>> Kitware offers ITK Training Courses, for more information visit:
>>> http://kitware.com/products/protraining.html
>>>
>>> Please keep messages on-topic and check the ITK FAQ at:
>>> http://www.itk.org/Wiki/ITK_FAQ
>>>
>>> Follow this link to subscribe/unsubscribe:
>>> http://www.itk.org/mailman/listinfo/insight-developers
>>
>> ========================================================
>> 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/20100114/e7538bf5/attachment.htm>
More information about the Insight-developers
mailing list