View License

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video

Highlights from
Absolute Orientation - Horn's method

Join the 15-year community celebration.

Play games and win prizes!

» Learn more

4.85714
4.9 | 14 ratings Rate this file 49 Downloads (last 30 days) File Size: 2.95 KB File ID: #26186 Version: 1.4

Absolute Orientation - Horn's method

by

Matt J (view profile)

 

23 Dec 2009 (Updated )

Solves weighted absolute orientation problem using Horn's quaternion-based method.

| Watch this File

File Information
Description

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 quaternion-based algorithm. The function works for both 2D and 3D coordinates, and also gives the option of weighting the coordinates non-uniformly. The code avoids for-loops 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 user-supplied N-vector 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 N-vector 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 counter-clockwise 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)

Acknowledgements

Absolute Orientation inspired this file.

This file inspired Microwae Engineering.

MATLAB release MATLAB 7.9 (R2009b)
MATLAB Search Path
/
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (46)
02 Dec 2016 Viktor Erdelyi

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);

13 Oct 2016 Chao-Min Huang

Good to use

17 Aug 2016 Jonas Eberhard  
25 Mar 2016 ADD

ADD (view profile)

 
05 Mar 2016 Vali ollah Maraghi  
04 Mar 2016 ADD

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.

Comment only
04 Mar 2016 Matt J

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]

Comment only
03 Mar 2016 ADD

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.

Comment only
11 Feb 2016 Matt J

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.

Comment only
11 Feb 2016 Matt J

Matt J (view profile)

@John,
I'm not sure you can.

Comment only
11 Feb 2016 Mirit Sharabi

can you please explain what is the input "varargin"?

Comment only
09 Feb 2016 John Kuzack

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.

Comment only
31 Jan 2016 Matt J

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 unreasonable-looking 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

Comment only
31 Jan 2016 harman Litt

Everything besides regParams is working

Comment only
31 Jan 2016 harman Litt

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?

Comment only
15 Sep 2015 Matt J

Matt J (view profile)

Thanks for the feedback, Meade. I'm not sure what issues you were talking about with the dealr() command, though.

Comment only
14 Sep 2015 Meade

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!

31 Jul 2015 RyanG

RyanG (view profile)

Hi Matt,
I understand now.
Thank you!

Comment only
30 Jul 2015 Matt J

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.

Comment only
30 Jul 2015 RyanG

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

Comment only
25 Jul 2015 Matt J

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.

Comment only
24 Jul 2015 RyanG

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?

Comment only
24 Jul 2015 RyanG

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

10 Jul 2015 Matt J

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.

Comment only
10 Jul 2015 Yingjuan Wu

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!

Comment only
09 Jun 2015 Matt J

Matt J (view profile)

Hi Alexander,
No, the algorithm does not handle refelections. Only rotations, translations, and (optionally) scalings.

Comment only
09 Jun 2015 Alexander Kraskov

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

Comment only
15 Apr 2015 Matt J

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

Comment only
13 Apr 2015 Fox

Fox (view profile)

I'd just like to know how to write the co-ordinates out. How would I input the 3 sets of A values and 3 sets of B values?

Thanks

Comment only
12 Apr 2015 Fox

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.

Comment only
12 Apr 2015 Matt J

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.

Comment only
11 Apr 2015 Fox

Fox (view profile)

I'm pretty new to this type of coding, could someone tell me where to put my two sets of 3D co-ordinates?

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.

Comment only
06 Nov 2014 ZZZZZ

ZZZZZ (view profile)

 
23 Jun 2014 Zohar Bar-Yehuda  
14 Apr 2014 Matt J

Matt J (view profile)

@Steve,
The points must be placed into the columns of a and b, not the rows.

Comment only
14 Apr 2014 Steve

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 z-value 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...

Comment only
14 Apr 2014 Steve

Steve (view profile)

Thanks for this very nice tool!

14 Nov 2013 Cong

Cong (view profile)

 
14 Nov 2013 Matt J

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.

Comment only
14 Nov 2013 Cong

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

Comment only
18 Apr 2013 monica

monica (view profile)

Thank you very much.

Comment only
16 Apr 2013 Matt J

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

Comment only
16 Apr 2013 monica

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+RSp-q=0 but i am not able to figure out the formula you have used for getting scaling factor. Could you please elaborate.

06 Jun 2012 Georg Wiora

Very nice work!

10 May 2012 Matt J

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 SVD-based 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.

Comment only
10 May 2012 Georg Stillfried

great work! (seems to work fine and also allows matching of three points)

Updates
25 Dec 2009 1.1

1. Important typo fix in help doc. Model is s*R*A+t not s*R*A-t
2. Added version for earlier MATLAB versions lacking bsxfun

22 Jan 2010 1.2

Added absorientParams. See README.txt for details

09 Apr 2010 1.3

Added tools for purely 2D registration

29 Sep 2010 1.4

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

02 Sep 2015 1.4

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

02 Sep 2015 1.4

Small edits to the file description.

Contact us