[Insight-developers] itk::Vector Normalize method
Arnaud Gelas
Arnaud_Gelas at hms.harvard.edu
Fri Jan 15 11:46:46 EST 2010
Hi Brad,
On Jan 14, 2010, at 8:22 PM, Bradley Lowekamp wrote:
> 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.
This is what I have been doing until now, but I was wondering if it
was not better to move this piece of code (which is finally quite
redundant over the the toolkit) into itk::Vector.
>
> 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?
At first, I was executing the same piece of code within the main
(reducing the calls), but to be more realistic I created a template
function ComputeNormTest.
Here is the latest version of the code.
Arnaud
>
> 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/20100115/df830aa0/attachment-0003.htm>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: CMakeLists-3.txt
URL: <http://www.itk.org/mailman/private/insight-developers/attachments/20100115/df830aa0/attachment-0001.txt>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/mailman/private/insight-developers/attachments/20100115/df830aa0/attachment-0004.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: VectorNormalization.cxx
Type: application/octet-stream
Size: 1713 bytes
Desc: not available
URL: <http://www.itk.org/mailman/private/insight-developers/attachments/20100115/df830aa0/attachment-0001.obj>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/mailman/private/insight-developers/attachments/20100115/df830aa0/attachment-0005.htm>
More information about the Insight-developers
mailing list