The London Perl and Raku Workshop takes place on 26th Oct 2024. If your company depends on Perl, please consider sponsoring and/or attending.
 / 
Math-Vector-BestRotation-0.009
River stage zero No dependents
/ Math::Vector::BestRotation

NAME

Math::Vector::BestRotation - best rotation to match two vector sets

VERSION

Version 0.009

SYNOPSIS

    use Math::Vector::BestRotation;

    my $best = Math::Vector::BestRotation->new();

    $best->add_pair([1, 2, 3], [4, 5, 6]);
    $best->add_pair([0, -1, 5], [-3, 6, 0]);
    .
    .
    .
    $best->add_pair([3, -1, 5], [0.1, -0.7, 4]);

    my $ortho = $best->best_orthogonal;
    my $rot   = $best->best_rotation;
    my $flip  = $best->best_improper_rotation;

    my $axis  = $best->rotation_axis;
    my $angle = $best->rotation_angle;

    # start over
    $best->clear;

DESCRIPTION

Introduction

Assume that you have a list of vectors v_1, v_2, v_3, ..., v_n and an equally sized list of vectors w_1, w_2, ..., w_n. A way to quantify how similar these lists are to each other is to compute the sum of the squared distances between the vectors: sum((w_1 - v_1)**2 + ... + (w_n - v_n)**2). In the literature, this sum is sometimes divided by 2 or divided by n or divided by n and the square root is taken ("root mean square" or RMS deviation).

In some situations, one data set can be arbitrarily rotated with respect to the other one. In this case, one of them has to be rotated in order to calculate the RMS deviation in a meaningful way. Math::Vector::BestRotation solves this problem. It calculates the best orthogonal map U between the v_i and w_i. "Best" means here that the RMS deviation between Uv and w as calculated above is minimized.

An orthogonal map can be a (proper) rotation or a rotation combined with a reflection (improper rotation). This module enables you to find the best orthogonal map, the best proper rotation, or the best improper rotation between two given vector sets.

Analysis

Once you have obtained your optimal map you might be interested in what was actually needed to optimize the match. Currently, the method offers to calculate the rotation axis and angle for a proper rotation. Support for improper rotations is planned. It might also be interesting to know how much (in terms of RMS deviation) is gained by applying the map. Right now you have to do this yourself, but support for this is also planned.

Outlook and Limitations

The algorithm implemented here is based on two research papers listed in the ACKNOWLEDGEMENTS section. It works for higher dimensional vector spaces as well, but the current implementation supports only three-dimensional vectors. This limitation is going to be remedied in a future version of this module.

The two data sets could not only be rotated with respect to each other, but also translated. This translation can be removed prior to the determination of the rotation by aligning the centers of mass of the two vector sets. However, this procedure is not offered by Math::Vector::BestRotation and possibly will never be, because this would require to store the full data sets in memory which is not necessary now.

The underlying algorithm supports to assign different weights to the vector pairs to reflect that it might be more important to align some pairs then others (e.g. because there measurement had a smaller error). This is currently not implemented but planned for the future.

INTERFACE

Constructors

new

  Usage   : Math::Vector::BestRotation->new(%args)
  Function: creates a new Math::Vector::BestRotation object
  Returns : a Math::Vector::BestRotation object
  Args    : initial attribute values as named parameters

Creates a new Math::Vector::BestRotation object and calls init(%args). If you subclass Math::Vector::BestRotation overload init, not new.

init

  Usage   : only called by new
  Function: initializes attributes
  Returns : nothing
  Args    : initial attribute values as named parameters

If you overload init, your method should also call this one. It provides the following functions:

  • For each given argument it calls the accessor with the same name to initialize the attribute. If such an accessor does not exist a warning is printed and the argument is ignored.

Public Attributes

Each of the following attributes has an accessor method of the same name which can be used to set or retrieve the stored value (it is mentioned below if an attribute is readonly). The accessors always return the current value (the new one in case it has been updated).

matrix_r

This attributes holds a matrix built from the pairs of vectors and used to compute the desired orthogonal map. It is called R in the documentation and the underlying Research papers. The accessor is readonly. At startup, matrix_r is initialized with zeros.

Note that the matrix is stored internally as an array to speed up data acquisition. When you call the accessor a Math::MatrixReal object is created. This implies that such an object is not updated as you add more vector pairs. You have to call the accessor again to get a new object. Accordingly, changing of your retrieved matrix does not alter the underlying matrix stored in the Math::Vector::BestRotation object.

matrix_u

Holds the result matrix after calling one of the best_... methods (undef before the first such call and after calling clear). Note that it will not be reset by calling add_pair or add_many_pairs. It will still hold the result of the last best_... call. The accessor is readonly.

Methods for Data Input

add_pair

  Usage   : $obj->add_pair([1, 2, 3], [0, 7, 5])
  Function: updates matrix_r
  Returns : nothing
  Args    : a pair of vectors, each as an ARRAY reference

The orthogonal map computed by this module will try to map the first vector of each pair onto the corresponding second vector. This method uses the new vector pair to update the matrix R which is later used to compute the best map. The vectors are discarded afterwards and can therefore not be removed once they have been added.

In some applications, very many vector pairs will be added making this the rate limiting step of the calculation. Therefore, no convenience functionality is offered. For example, the method strictly requires ARRAY references. If your vectors are stored e.g. as Math::VectorReal objects you have to turn them into ARRAY references yourself. Furthermore, no error checking whatsoever is performed. The method expects you to provide valid data. See also add_many_pairs.

add_many_pairs

  Usage   : $obj->add_many_pairs([[1, 2, 3], [3, 0, 6]],
                                 [[0, 7, 5], [-2, 1, -1]])
  Function: updates matrix_r
  Returns : nothing
  Args    : a pair of vectors, each as an ARRAY reference

An alternative to add_pair. It expects two lists of vectors. The first one contains the first vector of each pair, the second one contains the second vector of each pair (see add_pair). If you have many vector pairs to add it is probably faster to build these lists and then use this method since it saves you a lot of method calls.

For perfomance reasons, no checks are performed not even if the two lists have equal sizes. You are expected to provide valid data.

clear

  Usage   : $obj->clear
  Function: resets the object
  Returns : nothing
  Args    : none

This method resets matrix_r to the null matrix (all entries equal zero). This enables you to start from scratch with two new vector sets without destroying the object. Note that matrix_u is not reset.

Methods for Finding Maps

best_orthogonal

  Usage   : $matrix = $obj->best_orthogonal
  Function: computes the best orthogonal map between the vector sets
  Returns : a Math::MatrixReal object
  Args    : none

Computes the best orthogonal map between the two vector sets, i.e. the orthogonal map that minimizes the sum of the squared distances between the image of the first vector of each pair and the corresponding second vector. This map can be either a rotation or a rotation followed by a reflection (improper rotation).

The representing matrix of the found map is returned in form of a Math::MatrixReal object.

best_rotation

  Usage   : $matrix = $obj->best_rotation
  Function: computes the best rotation between the vector sets
  Returns : a Math::MatrixReal object
  Args    : none

This is identical to best_orthogonal except that it finds the best special orthogonal map (this means that the determinant is +1, i.e. the map is a true rotation).

The method computes the best rotation between the two vector sets, i.e. the rotation that minimizes the sum of the squared distances between the image of the first vector of each pair and the corresponding second vector.

The representing matrix of the found map is returned in form of a Math::MatrixReal object.

best_proper_rotation

Alias for best_rotation.

best_improper_rotation

  Usage   : $matrix = $obj->best_improper_rotation
  Function: computes the best rotation combined with a reflection
  Returns : a Math::MatrixReal object
  Args    : none

This is identical to best_orthogonal except that it finds the best orthogonal map with determinant -1. I do not know why one would want that, but the method is included for completeness.

The representing matrix of the found map is returned in form of a Math::MatrixReal object.

best_flipped_rotation

Alias for best_improper_rotation for backwards compatibility.

Methods for Result Analysis

rotation_axis

  Usage   : $axis = $obj->rotation_axis
  Function: computes the rotation axis of the last map found
  Returns : a unit vector in form of a Math::MatrixReal column
  Args    : optional a Math::MatrixReal object

In the three-dimensional case, a proper rotation (with an angle which is not a multiple of pi) leaves exactly one line fixed. This method takes the matrix stored in matrix_u and computes a unit vector along this axis and returns it in form of a Math::MatrixReal object (column vector). There are two special cases

Special cases

1. The rotation angle is a multiple of 2pi. This is the same as no rotation at all and an axis cannot be determined uniquely (in fact, each vector is as good as any other). In this case, a warning is printed and undef is returned. A warning is also printed if the rotation is very close to the identity map such that numerical instability a threat. See also Warnings.
2. The rotation angle is a odd multiple of pi. In this case, also lines in the plane perpendicular to the axis are mapped onto themselves. Still, a vector in the direction of the rotation axis is returned.

Orientation of the vector

Even if the rotation axis is unique, the orientation of the vector is not (a unit vector along the axis multiplied with -1 is as good). One could determine it in a way that the rotation angle in mathematical positive direction is less or equal than pi. This might come in the future. Currently, the orientation of the vector that is returned has to be considered as arbitrary (and might change between module versions).

Improper rotations

Improper rotations are currently not supported and lead to an exception. However, this is planned for a future version of this module.

Higher-dimensional vector spaces

Spaces with more than three dimensions are not supported. I do not know if they ever will be. In higher dimensions the eigenspace to the eigenvalue 1 can have more dimensions than one. I do not know what one would want to do with this. Moreover, the algorithm described below cannot be trivially extended to higher dimensions. There are other solutions, though. Please contact me if you would like to have some kind of this functionality.

Algorithm

In the case of a proper rotation, we look for an eigenvector of the rotation matrix with the eigenvalue 1. The canonical way is to solve the system of linear equations given by (U - I)v = 0 where I denotes the unity matrix. However, rounding errors can lead to the case where U is not strictly orthogonal and the system has only the trivial solution v = 0.

Instead, the method calculates the matrix (U + ~U) - (trU - 1)I where ~U is the transpose of U, trU is the trace of U, and I is the unity matrix. This matrix is a multiple of v * ~v where v is an axis vector. See reference [3] in the ACKNOWLEDGEMENTS for details. v is then extracted as a non-zero column of this matrix.

rotation_angle

  Usage   : $angle = $obj->rotation_angle
  Function: computes the rotation angle of the last map found
  Returns : the angle between 0 and pi
  Args    : none

Approach

The angle is calculated as acos((trU - 1) / 2). See reference [3] in the ACKNOWLEDGEMENTS for details. Currently, the sign of the angle is uncorrelated to the orientation of the axis vector returned by rotation_axis. This will be remedied in a future version, but for now, the returned value is the absolute value of the rotation angle.

Improper rotations

Improper rotations are currently not supported and raise an execption. However, this is planned for a future version of this module.

Higher-dimensional vector spaces

Spaces with more than three dimensions are not supported. I do not know if they ever will be. In higher dimensions the eigenspace to the eigenvalue 1 can have more dimensions than one. Currently, I do not know if a rotation angle can be defined in a meaningful way. Anyway, the algorithm mentioned above cannot be trivially extended to higher dimensions. Please contact me if you would like to have some kind of this functionality.

DIAGNOSTICS

Exceptions

Sorry, not documented, yet. Exceptions are thrown using croak (see Carp) in the case of user "misconduct".

Warnings

Sorry, not documented in detail, yet. Warnings are printed using carp (see Carp) when numerical instabilities have to be expected. This cannot be switched off at the moment, but in the future, there might be a verbose attribute.

BUGS AND LIMITATIONS

Bugs

No bugs have been reported, but the code should be considered as beta quality.

Please report any bugs or feature requests to bug-math-vector-bestrotation at rt.cpan.org, or through the web interface at http://rt.cpan.org/NoAuth/ReportBug.html?Queue=Math-Vector-BestRotation. I will be notified, and then you will automatically be notified of progress on your bug as I make changes.

Limitations

See DESCRIPTION.

AUTHOR

Lutz Gehlen, <perl at lutzgehlen.de>

ACKNOWLEDGEMENTS

The algorithm implemented here is based on two research papers by Wolfgang Kabsch, Max-Planck-Institut fuer Medizinische Forschung, Heidelberg, Germany:

[1] Kabsch, W. (1976). A solution for the best rotation to relate two sets of vectors. Acta Cryst., A32, 922
[2] Kabsch, W. (1978). A discussion of the solution for the best rotation to relate two sets of vectors. Acta Cryst., A34, 827-828

The determination of rotation axis and angle follows derivations laid out in

[3] Fillmore, J. P. (1984). A Note on Rotation Matrices. IEEE Comput. Graph. Appl., vol. 4, no. 2, pp. 30-33

LICENSE AND COPYRIGHT

Copyright 2010 Lutz Gehlen.

This program is free software; you can redistribute it and/or modify it under the terms of either: the GNU General Public License as published by the Free Software Foundation; or the Artistic License.

See http://dev.perl.org/licenses/ for more information.