Comments and Ratings

Great Work.
I want to register pair of images with your toolbox. But the size of my images are 320 by 240. How can I change them into N by 2 point sets. Thanks.

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

great work, Andriy.

Is there any way to register with translation and rotation only without scaling? I set the .scale option to 0, but the scaling value in the computed transformation is not equal to 1.

2Jichao Zhao:
Looks like you haven't added all CPD subfolders to the Matlab path.
In Matlab, go to File-Set_Path-Add_with_subfolders, select the CPD main folder, and click save.

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.

These functions can be changed to allows for zero padding?
It is an important concern in FFT...
I'm not sure if it can be done by simply changing in your codes the term "fft(x)" by "fft(x,n)", where n is the size of the resulting padded output...

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.