ABSOR is a tool for least squares estimation of the rotation  and optionally also the
scaling and translation  that maps one collection of point coordinates to
another. It is based on Horn's quaternionbased algorithm. The function works for both 2D and 3D coordinates, and also gives the option of weighting the coordinates nonuniformly. The code avoids forloops so as to maximize speed.
DESCRIPTION:
As input data, one has
A: a 2xN or 3xN matrix whos columns are the coordinates of N source points.
B: a 2xN or 3xN matrix whos columns are the coordinates of N target points.
The basic syntax
[regParams,Bfit,ErrorStats]=absor(A,B)
solves the unweighted/unscaled registration problem
min. sum_i R*A(:,i) + t  B(:,i)^2
for unknown rotation matrix R and unknown translation vector t.
ABSOR can also solve the more general problem
min. sum_i w(i)*s*R*A(:,i) + t  B(:,i)^2
where s>=0 is an unknown global scale factor to be estimated, along with R and t,
and w is a usersupplied Nvector of weights. One can include/exclude any
combination of s, w, and translation t in the problem formulation. Which
parameters participate is controlled using the syntax,
[regParams,Bfit,ErrorStats]=absor(A,B,'param1',value1,'param2',value2,...)
with parameter/value pair options,
'doScale'  Boolean flag. If TRUE, the global scale factor, s, is included.
Otherwise, it is assumed that s=1. Default=FALSE.
'doTrans'  Boolean flag. If TRUE, the translation, t, is included. Otherwise,
zero translation is assumed. Default=TRUE.
'weights'  The length Nvector of weights, w. Default, no weighting.
OUTPUTS:
regParams: structure output with estimated registration parameters,
regParams.R: The estimated rotation matrix, R
regParams.t: The estimated translation vector, t
regParams.s: The estimated scale factor.
regParams.M: Homogenous coordinate transform matrix [s*R,t;[0 0 ... 1]].
For 3D problems, the structure includes
regParams.q: A unit quaternion [q0 qx qy qz] corresponding to R and
signed to satisfy max(q)=max(abs(q))>0
For 2D problems, it includes
regParams.theta: the counterclockwise rotation angle about the
2D origin
Bfit: The rotation, translation, and scaling (as applicable) of A that
best matches B.
ErrorStats: structure output with error statistics. In particular,
defining err(i)=sqrt(w(i))*norm( Bfit(:,i)B(:,i) ),
it contains
ErrorStats.errlsq = norm(err)
ErrorStats.errmax = max(err)
1.4.0.0  Small edits to the file description. 

1.4.0.0  Added the option of constraining translation to zero (see the 'doTrans' parameter). 

1.4.0.0  *The many previous files in the distribution have been consolidated into a single file absor.m. *New capability:weighted least squares registration. *Fixed minor bug in the 2D registration routines, occurring when rotation angle was 0. 

1.3.0.0  Added tools for purely 2D registration 

1.2.0.0  Added absorientParams. See README.txt for details 

1.1.0.0  1. Important typo fix in help doc. Model is s*R*A+t not s*R*At

Create scripts with code, output, and formatted text in a single executable document.
JeanNoel Saugy (view profile)
Matt J (view profile)
Hi Saeid,
It's been a while since I did that derivation, but I remember that in one of the appendices in Horn's paper, he gives the quaternion form for the rotation about a known axis. By taking that parametric form for the case when the rotation axis is the zaxis, and also by assuming that the zcoordinates of all the input data are zero, the equations in the paper simplify and the optimal quaternion calculation reduces from a 3D eigenvalue problem to a 2D eigenvalue problem.
Saeid Nooshabadi (view profile)
Hi Matt,
The formulation for quaternionbased method in the original Horn's reference is for 3Dimensions. Can you please send me the reference to the formulation for the case of 3Dimensions? Alternatively please explain how the formulation works for 2Dimensions. The maths for 2D eludes me!
Matt J (view profile)
@Conrade,
No, it doesn't sound like this would be the right tool.
Conrade Muyambo (view profile)
Can I use your code to find the orientation/major/symmetry axis of a 2D shape from its xy coordinates, then rotate it through this axis?The centroid should also pass through this major axis, if possible.
Matt J (view profile)
Hi Neil,
Thanks for the feedback. When I wrote this code, it was to serve my own projects in which essentially all the errors were in B and none in A (a fairly common case when A are the points in some calibration rig). In section 2E, he mentions that an asymmetric formulas are more appropriate when the errors are asymmetric. I somehow convinced myself that the section 2D asymmetric formula was optimal for the case where all errors lay in B. But perhaps I should consider adding support for symmetric scaling in the code.
Neil Weisenfeld (view profile)
Hey this is great  I just wanted to point about the scalefactor  you reference section 2D in the original paper, but if you continue on to 2E, he demonstrates that 2D is suboptimal in that it's asymmetric and depends on the rotation, while 2E derives a symmetric scalefactor that can be computed prior to or independent of the rotation.
Starnav include (view profile)
Good Job, Thanks
Lizhao Du (view profile)
good job！
puni (view profile)
Matt J (view profile)
Hi Gir,
It will work with either double or single. Undoubtedly, the error is occurring because you downloaded absor.m to a folder that Matlab can't see.
Gir Val (view profile)
Hi Matt,
I am trying to use ABSOR, but am not sure what data type to convert my coordinates to. I've tried double, single, int8, int 64. Could you provide some input? I'm reading in two 2x1105 matrices.
My error is <Undefined function 'absor' for input arguments of type 'double'.>.
Emre Keski (view profile)
Matt J (view profile)
Hi Carlos,
Why couldn't you just project your 3D source points onto the camera plane and then do 2D2D registration in that plane?
Carlos Ruiz (view profile)
Hi Matt, thanks so much for this great code!
I have a source set of points in 3D but a target that is projected to 2D (camera plane), do you have any ideas of how I could modify the code to allow for a 3Dto2D affine projection matrix?
Thanks again!!
Matt J (view profile)
Hi Mathieu,
A linear transformation T*X=Y with anisotropic scaling has no special structure. You would just do T=Y/X to estimate it.
Mathieu Alloing (view profile)
Thanks a lot for this great code !
I was wondering if you had an idea of how to modify your code to handle anisotropic scaling ?
Viktor Erdelyi (view profile)
Matt J (view profile)
Hi Ha,
[V,D] =eig(regParams.R) MUST work. The form of the transformation is y=R*x and a point on the axis of rotation must be unchanged by the rotation. In other words, it must satisfy y=x or equivalently x=R*x. From the last, you can see that x must be an eigenvalue of R with eigenvalue 1.
Ha HWANG (view profile)
Hi, Matt
Thanks for your answer.
I'm still testing your code. It works for a 2D model. But for a 3D model, it doesn't work at all although I get a very small error and Bfit is also perfect but regParams.t vectors and three euler angles I calculate from regParams.R are completely off. The origin of the rotation [V,D] =eig(regParams.R) does not work, either.
Does an object have to be oriented about zero for the code to work?
My 3D model is quite simple: a cone that has 51 data points. Will you be able to give me some pointers if I post two data sets and all necessary info?
Thanks
Matt J (view profile)
Hi Ha,
The rotation axis for a 3x3 rotation matrix R, can be found by finding the eigenvector of R that has eigenvalue 1, see <https://en.wikipedia.org/wiki/Euler%27s_rotation_theorem#Matrix_proof>
In Euler angle decomposition, there are always 3 axes/angles. If the code you have found offers the option of decomposing into an arbitrary sequence of rotations, it is doing something weird. In any case, if you have come across code that does not appear to be working, you should probably consult the author of that code.
Ha HWANG (view profile)
Hi Matt
Thanks a lot for your quick answer.
Your code works great. It gives me translational information that is correct. Although, I have some problems trying to get the rotation information;
 How do I find the origin of the rotations (a base point that my model is rotated about)?
 Is an order of the rotational axes (XYZ or XZY or ZYX, etc...) or a number of the axes (X or XY or XYZ) important when calculating the euler angles???? I found a coupled of codes here that can convert your 'regParams.R' to euler angles. When my model is rotated twice (about X and Z), I can get correct angles about both axes. But with the 3rd rotation about Y, all three angles are completely off. Does it have anything to do with an order or a number of rotations?
Thank you
Matt J (view profile)
Hi Ha Hwang,
1. Reconstructing B from regParams.R*A+regParams.t is only possible if B satisfies the model B=R*A+t exactly for some R and t. But your example data does not. Incidentally, you can call absor with more output arguments to obtain the refitted values of B, and the errors in the fitting: [regParams,Bfit,ErrorStats]=absor(A,B) .
2. In general, it takes 3 angles (not one) represent a 3x3 rotation matrix R. However, this representation is not unique, so you need to decide which one you want. if you have the Robotics System Toolbox, you can convert R matrices to Euler angles using https://www.mathworks.com/help/robotics/ref/rotm2eul.html. If you don't have the toolbox , there are various offerings for Euler angle decomposition elsewhere on the FEX, see https://www.mathworks.com/matlabcentral/fileexchange/?utf8=%E2%9C%93&term=Euler+angles.
Ha HWANG (view profile)
Hi,
I'm pretty new to this matrix transformation. I'd like to use your code but I don't think I'm understanding it well. Can I ask you a couple of questions?
1 *****
I just tested your code but I don't get the result I expected. Can you please tell me what's wrong?
A=[1 25 62;2 58 74] and B=[2 55 10;4 65 87]. Then I get;
regParams.R
0.90880 0.41722
0.41722 0.90880
To recalculate B from A and R, I did; regParams.R*A+regParams.t which gives me
14.3854 12.8322 39.7824
1.4030 62.3094 92.2876
Aren't I supposed to get exactly the B values from this????
2*****
Eventually, I want to use your code to find out a rotation angle and translation values in XYZ coordinates for a mesh model (I have two mesh models that are identical but one of them is rotated and translated to unknown degrees). So is it possible to calculate an angle (theta) from your rotational matrix R and extract translational values for XYZ??
Your help would be greatly appreciated.
Thank you
Matt J (view profile)
Hi sjp228,
Thanks for your feedback and kind rating. For conversion of R matrices to Euler angles, there is https://www.mathworks.com/help/robotics/ref/rotm2eul.html if you have the Robotics System Toolbox. If not, there are various offerings for Euler angle decomposition on the FEX, see https://www.mathworks.com/matlabcentral/fileexchange/?utf8=%E2%9C%93&term=Euler+angles.
sjp228 (view profile)
Hi Matt,
I was wondering if there is any way to get the 3D rotation angles from the R matrix?
a j (view profile)
Domenico Scaramozzino (view profile)
Natalie Larson (view profile)
Matt J (view profile)
@Brandon
Thanks for your feedback. To find the center of rotation, you must take the R,t matrices that the code gives you and use MATLAB's linear algebra tools to find the solutions of
x=R*x+t
Note that the space of solutions, if nonempty, will be nonunique, since the 3x3 matrix IR always has a nonempty null space.
Brandon Marshall (view profile)
Great function. I have a quick question, if I have two sets of 3D coordinates, it will give me the rotation and translation matrices. How can I calculate the center of rotation (3D) based on this?
Rob Campbell (view profile)
Just what I need: Thanks!
Matt J (view profile)
Hi Jonathan,
You can compute the transformation based on any selection of points that obey that same transformation rule. If all 3820 points transform rigidly, for example, you can use a 3x3820 matrices for A and B. However, to compute the transformation, you must be able to match the points, i.e., A(:,i) must correspond to B(:,i).
Jonathan Cyganik (view profile)
I have two complex 3D objects each with 3820 vertices of 3 points that join to make faces. I am looking for the transformation matrix of 4x4x3820 points by which I can multiply each point of one of these objects to transform it into the other. Would I have to run this function for each vertex? Or the entire matrices at once? Would I then use the M matrix and corresponding B vector to make each of the 4x4 transformation matrices? Any help would be greatly appreciated.
Ahmed Ayman (view profile)
Matt J (view profile)
Hi Chauvin,
Yes, the routine is only for finding transformations of points that have been rotated, translated, or scaled. They must be the same number of points and you must know which pairs of points correspond and feed them to ABSOR in corresponding order.
Chauvin Axel (view profile)
Hi Matt,
First, congrats for your work and thank you for the contribution.
I am using absor() function and I have a few issues and maybe misunderstandings.
Is the method made only for finding the transformation of a set of points that has been rotated and translated to its original ?
Does it handle only sets with the same number of points ? (I think the answer is yes according to your code but I prefer to be sure)
If I have two set of points representing two consecutive slices of an object but one is rotated, does your method fit this problem ?
Thank you
PS : Sorry for my english is bad
Mendi Barel (view profile)
Viktor Erdelyi (view profile)
Here is a suggestion to handle reflections (at least in the 2D case):
[regParamsOriginal,~,ErrorStatsNormal]=absor(X_est, X_true);
[regParamsReflected,~,ErrorStatsReflected]=absor([1 0; 0 1] * X_est, X_true);
if ErrorStatsReflected.errlsq < ErrorStatsNormal.errlsq
X_est = [1 0; 0 1] * X_est;
regParams = regParamsReflected;
disp('Choosing reflected');
else
regParams = regParamsOriginal;
disp('Choosing original');
end
X_est = regParams.R * X_est + repmat(regParams.t, 1, M);
ChaoMin Huang (view profile)
Good to use
Jonas Eberhard (view profile)
ADD (view profile)
Vali ollah Maraghi (view profile)
ADD (view profile)
@ Matt J,
That's Right. Essentially whenever matrix N has repeating eigenvalues that are larger than all the other eigenvalues, the machine precision determines what eigenvector gets picked and that's where different answers might be produced. Thanks for your help and for the wonderful implementation of Professor Horn's work.
Matt J (view profile)
@ADD,
This could occur if the solution to the registration problem is ambiguous. For example, on my machine absor(A,B) with
A =
1 1
0 0
B =
1 1
0 0
gives R=eye(2), t=[1;0]. However, an equally valid solution is R=eye(2), t=[1;0]
ADD (view profile)
@ Matt J:
Does this algorithm depend on a particular form type of eigenvector calculation? To be more clear, if I replace the results of line#194 with what other eigenvector computation algorithms give as the result for instance using numpy.linalg.eig() or scipy.linalg.eig() and run the rest of the absor command what you would get is not the same answer as what Matlab's eig() function produces. Note that the results of numpy and scipy are perfectly correct as the A.X=lambda.X is satisfied but for some reason absor algorithm is sensitive to what the Eigenvector is. If you could clarify why this is true and how it can be overcome. Many thanks.
Matt J (view profile)
@Mirit,
VARARGIN is a Matlab 'command', which you can read about here,
http://www.mathworks.com/help/matlab/ref/varargin.html
It might be useful to you when writing your own MATLAB functions, but it's nothing you need to know about for using ABSOR. Typing 'help absor' should be all you need.
Matt J (view profile)
@John,
I'm not sure you can.
Mirit Sharabi (view profile)
can you please explain what is the input "varargin"?
John Kuzack (view profile)
Could you provide an example of how i could use this to calculate the center of rotation of one fixed point and one moving point.
Matt J (view profile)
harman,
Are you actually calling the function with output arguments, like below. MATLAB cannot just return _nothing_ when outputs are specifically requested. It must at least return an error message. If you are getting error message I would need to see copy/pastes of them. In any case, I get a not unreasonablelooking registration with your data:
>> [regParams,Bfit,ErrorStats]=absor(A,B)
regParams =
R: [3x3 double]
t: [3x1 double]
s: 1
M: [4x4 double]
q: [4x1 double]
Bfit =
0.9294 3.3154 3.6553
2.1461 2.4824 3.2715
1.0077 1.5681 2.3242
ErrorStats =
errlsq: 0.9306
errmax: 0.6896
harman Litt (view profile)
Everything besides regParams is working
harman Litt (view profile)
I am having trouble getting any output from this function.
I am using the following as my two points :
A= [1.2,3.4,3.5; 1.4,2.2,3.1; 1.1,1.9,2.6];
B= [1.5,3.3,3.1; 2.4,2.3,3.2; 1.3,1.4,2.2 ];
Then I type absor(A,B);
but for some reason when I run the program I get no output. Can someone tell me what i am doing wrong?
Matt J (view profile)
Thanks for the feedback, Meade. I'm not sure what issues you were talking about with the dealr() command, though.
Meade (view profile)
A fantastic script; very useful.
I appreciate that you've written this in a way that ensures maximum compatibility with multiple versions (e.g. bsxfun checks).
With that in mind, solve the issues with the 'deal' command and drop the function call at 214 by changing the following:
M = num2cell(M);
[Sxx,Syx,Szx,...
Sxy,Syy,Szy,...
Sxz,Syz,Szz] = M{:};
It helped me to incorporate your code in a larger wrapper I'm working on.
Thanks again!
RyanG (view profile)
Hi Matt,
I understand now.
Thank you!
Matt J (view profile)
Hi Ryan,
M is the optimum transformation operator that best maps A to B.
As for Bfit, my explanation in my previous comment was in error. Bfit(:,i) is what you get if you apply the optimum transformation M to the original data A(:,i). It is the best fit to B(:,i). If you are not interested in fitting the B data, there is no need to call Bfit.
RyanG (view profile)
Hi Matt,
Thank you again for the great script and helpful response. I assumed that the output of M was already the best fit solution. Is this not the case? Do I need to call 'Bfit' when I call absorb for this to be applied?
Ryan
Matt J (view profile)
Hi Ryan,
The general way to transform an arbitrary point P is s*R*P+t, though s will simply be 1 when you call absor without the doScale option. To perform the transformation in homogeneous coordinates, you can do M*[P;1] padding P with a 1, as you mentioned.
If you call absor with the Bfit output, this transformation will be applied for you internally on all B(:,i), sparing you the effort of doing so yourself manually.
RyanG (view profile)
Hi Matt,
A follow up to my earlier question, what is the use and purpose of the Bfit? Should I be using Bfit to map points A to B versus the M transformation matrix?
RyanG (view profile)
Hi Matt,
Thank you for supplying this useful script. I had a question regarding applying the resultant M transformation matrix. Below you state that to map your A to B points you can perform the following: Pnew=M*P. With P being A basically.
But that it requires homogeneous coordinates. I padded the 4th row of my data with 1's so that the above matrix multiplication could be performed. What is the fourth row? Are 1's required there so that when scaling is applied that the matrix dimensions match?
Thank you,
Ryan
Matt J (view profile)
Hi Yingjuan,
I'm not sure you can do it if you mean to consider rotation only (no translation). If you do,
absor([A1,...,A15],[B1,...,B15])
then the code will give you the best combination of rotation matrix and translation vector simultaneously matching all A to the corresponding B. It cannot constrain the transformation to be a rotation only.
Yingjuan Wu (view profile)
Hi Matt,
Thanks for your code. I have 15 A vectors (A1 to A15) and 15 correspondign B vectors (B1 to B15)in 3D. Now I want to figure out if these vectors have the same rotation angle (or rotation matrix). For example, A1 to B1 a has rotation angle theta, I want to know if A2 to B2 or A3 to B3 has the same rotation angle. Could I use your code to solve this problem? Thanks!
Matt J (view profile)
Hi Alexander,
No, the algorithm does not handle refelections. Only rotations, translations, and (optionally) scalings.
Alexander Kraskov (view profile)
Hi Matt,
I probably should read the original paper but though it would be quicker to ask you.
It seems that your code (and maybe method) does not treat correctly reflection cases. Is it correct?
Here is an example my example is
A=
17.6607 53.2300 69.2854
14.5424 31.3154 50.5354
28.9079 5.8403 55.3454
43.7231 3.0544 78.9054
48.0822 25.4468 98.7854
34.5058 50.4406 92.8754
19.4431 37.4269 83.5354
B=
24.0000 24.3000 72.3000
5.4000 2.5000 78.2000
10.3000 25.4000 66.5000
34.0000 31.3000 51.7000
55.4000 8.4000 42.2000
48.2000 19.3000 56.0000
38.5000 7.8000 72.4000
It produces wrong (big error) transformation matrix.
regParams.M
0.7865 0.3739 0.4916 40.6255
0.5476 0.7903 0.2751 29.5564
0.2856 0.4855 0.8262 102.4473
0 0 0 1.0000
>> ErrorStats
ErrorStats =
errlsq: 20.7210
errmax: 18.1723
If I reflect the first coordinate of A (multiply by
1 0 0
0 1 0
0 0 1)
than I get more reasonable transformation matrix
0.0127 0.0020 0.9999 44.3410
0.1158 0.9933 0.0005 27.5493
0.9932 0.1158 0.0129 96.5021
0 0 0 1.0000
ErrorStats =
errlsq: 3.9922
errmax: 2.9032
My question is whether it is possible to deal with it in a automatic way?
cheers
Alexander
Matt J (view profile)
Hi Fox,
I'm not sure if you've been through the help documentation or not, but the input syntax is described there and in the Description section above, in particular where it says
A: a 2xN or 3xN matrix whos columns are the coordinates of N source points.
B: a 2xN or 3xN matrix whos columns are the coordinates of N target points.
So, put your points to be matched in the columns of matrices A and B. They must be in correspondng order, i.e., B(:,j) must be the rototranslation of A(:,j).
Fox (view profile)
I'd just like to know how to write the coordinates out. How would I input the 3 sets of A values and 3 sets of B values?
Thanks
Fox (view profile)
Thanks for the quick response. I've just sent you a message.
I've just realised this might not work for what I require.
Matt J (view profile)
@Fox,
You haven't shown your call to absor(), so we don't know precisely what you're doing. However, your data example,
A = [1,2,3]'
B = [2,1,2]'
only contains 1 point in each set. You need at least 3 points in each set to do a sensible registration. They are to be placed in the columns of A,B in corresponding order.
Fox (view profile)
I'm pretty new to this type of coding, could someone tell me where to put my two sets of 3D coordinates?
If I put for example:
A = [1,2,3]'
B = [2,1,2]'
It says that A and B are used in the function, so they get changed to  in the first line of the code.
It does run, however doesn't give me any values. Just gives like [3x1 double matrix] in writing.
Thanks.
ZZZZZ (view profile)
Zohar BarYehuda (view profile)
Matt J (view profile)
@Steve,
The points must be placed into the columns of a and b, not the rows.
Steve (view profile)
One question about the translation vector (regParams.t):
I have the following 3 points before rotation:
A1 = (400, 400, 200,)
A2 = (600, 600, 200,)
A3 = (2000, 2000, 2000)
and after rotation:
B1 = (400, 400, 200,)
B2 = (600, 600, 200,)
B3 = (2000, 2000, 2001)
so all point are the same, just B3z is 2001 instead of 2000. Then the computed translation vector using your code is: (162, 76,4001)
What doesn't make sense to me? Because just one single point changed it's zvalue from 2000 to 2001 what is just a difference of 1...?
Nevertheless, the result seems to be right, because:
b(:,1) = regParams.R * a(:,1) + regParams.t
is true...
Steve (view profile)
Thanks for this very nice tool!
Cong (view profile)
Matt J (view profile)
Hi Cong,
If you're asking how to apply the rotation matrix R to rotate a point P=[x;y;z], you would just use matrix multiplication Pnew=R*P.
Note, however, that ABSOR calculates a translation vector, t, as well (and sometimes also a scaling if you select the doScale option). The full transformation would be Pnew=R*P+t.
Finally, if you are working in homogeneous coordinates P=[x;y;z;1], ABSOR also returns a 4x4 total transformation matrix M. In that case, the transformation can then be done Pnew=M*P.
Cong (view profile)
Thank you to provide this very useful toolbox. I have a question: if I obtain the rotation matrix from your method, then how can I get the coordinates of any point in the new coordinate system (target).
monica (view profile)
Thank you very much.
Matt J (view profile)
Hi Monica. You can find the derivation of the scaling factor in the original paper by Horn, Section 2D
http://people.csail.mit.edu/bkph/papers/Absolute_Orientation.pdf
monica (view profile)
Thank you Matt for this submission. I want to know how you got scaling factor sss=summ(right.*(R*left))/summ(left.*left). I know t+RSpq=0 but i am not able to figure out the formula you have used for getting scaling factor. Could you please elaborate.
Georg Wiora (view profile)
Very nice work!
Matt J (view profile)
Thanks for the feedback, Georg. I'm not entirely sure why the case of 3 points has been giving you trouble, though. It is true that alternative SVDbased methods and orthogonal matrix methods had to be modified to handle coplanar point data, but I always understood that to be a fairly stable solution.
Georg Stillfried (view profile)
great work! (seems to work fine and also allows matching of three points)