and 1 contributors

# NAME

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

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

``````  Usage   : \$obj->add_pair([1, 2, 3], [0, 7, 5])
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.

``````  Usage   : \$obj->add_many_pairs([[1, 2, 3], [3, 0, 6]],
[[0, 7, 5], [-2, 1, -1]])
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.

``````  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  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  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.

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:

 Kabsch, W. (1976). A solution for the best rotation to relate two sets of vectors. Acta Cryst., A32, 922
 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

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