Comments and Ratings

Hi Steffen,
Thanks for bringing up some bugs in the cpd code.
1) In cpd_rigid, you right, it should be with abs: ntol=abs((L-L_old)/L). It was my mistake. Fix that please. Now your code should work fine when using scaling.
2) Without scaling estimation (rotation and translation only)
(opt.normalize=0; opt.scale=0;), It get stuck into a local minimum of rotation. This is partially because your data set is very sparse and co-planar. Try using a slightly different initialization of one of your point-sets, let's say rotate Y a bit, For instance:

angle=0.1;
ca=cos(angle); sa=sin(angle);
R=[0 ca -sa;
   0 sa ca ;
   1 0 0];
Y=Y*R;

and then run the registration with your settings. You should get an accurate registration result. Perhaps, more systematic solution would be to run the rigid registration with several different rotation initializations.
Please, email my if you find more bugs or have more questions.
I'll fix these bugs and update my code shortly.
-Andriy

To candelabro
1) As for the zero-padding, modifying fft(x,n) is incorrect.
Just pad the input 3D vector x before taking the dct.
For instance, use :

newsize=[M N K]; % Define padded size of the 3D array;
tmp=zeros([M N K]); % padded data buffer
cutsize=min([size(x); newsize]);
tmp(1:cutsize(1),1:cutsize(2),1:cutsize(3))=x(1:cutsize(1),1:cutsize(2),1:cutsize(3));
x=tmp;
y=mirt_dctn(x); take 3D DCT (padded or reduced to the M N K size)

2) The positive and negative frequencies is only appropriate for fft. If you really want to have the analogy with fft, then remember that DCT gives an output for positive frequencies. And negative frequencies can be obtained by mirror flip.

to Paul Matthews,
 thanks for the feedback,

I was surprised by your results at first, then I find out that your results in 1D are correct only for ROW vectors. Try the same with COLUMN vectors, please:

tic; for i=1:3000;x=rand(512,1);dct(x);end; toc % 0.81s
tic; for i=1:3000;x=rand(512,1);idct(x);end; toc % 0.91s

tic; for i=1:3000;x=rand(512,1);mirt_dctn(x);end; toc % 0.66s
tic; for i=1:3000;x=rand(512,1);mirt_idctn(x);end; toc % 0.70s

For the larger 1D vectors, mirt_dctn (and mirt_idctn) performance improvement is even bigger
 tic; for i=1:3000;x=rand(4096,1);mirt_dctn(x);end; toc % 2.1s
 tic; for i=1:3000;x=rand(4096,1);dct(x);end; toc % 4.55s.

PS: I'll take a look why my 'row' vector processing takes longer, and correct it shortly.