For large N most of the time is spent in the two for loops.
To remove them bsxfun is useful. Replace the lines in Kabsch.m:
------------------------------------
Pdm = zeros(D,N) ;
for i=1:N
Pdm(:,i) = m(i)*P(:,i) ;
end
C = Pdm*Q' ;
------------------------------------
by
------------------------------------
C = bsxfun(@times,m,P)*Q';
------------------------------------
and the lines
------------------------------------
lrms = 0 ;
for i=1:N
lrms = lrms + m(i)*Diff(:,i)'*Diff(:,i) ;
end
lrms = sqrt(lrms) ;
------------------------------------
by
------------------------------------
lrms = sqrt(sum(sum(bsxfun(@times,m,Diff).*Diff))) ;
------------------------------------
and the execution speed will be increased by ~70 for N=3'000'000

22 Oct 2012

Kabsch algorithm
Find the rigid transformation & Least Root Mean Square distance between two paired sets of points
Author: Ehud Schreiber

For large N most of the time is spent in the two for loops.
To remove them bsxfun is useful. Replace the lines in Kabsch.m:
------------------------------------
Pdm = zeros(D,N) ;
for i=1:N
Pdm(:,i) = m(i)*P(:,i) ;
end
C = Pdm*Q' ;
------------------------------------
by
------------------------------------
C = bsxfun(@times,m,P)*Q';
------------------------------------
and the lines
------------------------------------
lrms = 0 ;
for i=1:N
lrms = lrms + m(i)*Diff(:,i)'*Diff(:,i) ;
end
lrms = sqrt(lrms) ;
------------------------------------
by
------------------------------------
lrms = sqrt(sum(sum(bsxfun(@times,m,Diff).*Diff))) ;
------------------------------------
and the execution speed will be increased by ~70 for N=3'000'000

5

22 Oct 2012

Kabsch algorithm
Find the rigid transformation & Least Root Mean Square distance between two paired sets of points

I found a little mistake that results in numerical jitter under some circumstances.
If you change the line:
"if (det(C) < 0)" into "if (det(W*V') < 0)"
all works fine.
Cheers, Andreas