[Insight-developers] itk::Vector Normalize method
Arnaud Gelas
Arnaud_Gelas at hms.harvard.edu
Thu Jan 14 17:25:12 EST 2010
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/f57981eb/attachment.htm>
More information about the Insight-developers
mailing list