Andriy Myronenko

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.

I checked the code and found out that one cannot use normalization and fixed scaling together. So, I set opt.normalize=0. however, the computed transformation seems to be far away from the optimal transformation. here is my test code:

function test_rigidTranform_noScaling( )

X = [121.3102 58.0957 46.3863
    83.3762 82.3011 69.5807
    88.7430 46.7092 84.9432
    51.3302 67.5140 109.3134];

Y = [36.7239 68.1019 52.6539
    61.0168 46.4281 88.7562
    76.5892 82.4971 84.0007
    99.3448 58.1160 121.3272];

% Set the options
opt.method='rigid'; % use rigid registration
opt.viz=0; % show every iteration
opt.outliers=0; % do not assume any noise

opt.normalize=0; % normalize to unit variance and zero mean before registering (default)
opt.scale=0; % fixed global scaling = 1
opt.rot=1; % estimate strictly rotational matrix (default)
opt.corresp=1; % compute the correspondence vector at the end of registration (default)

opt.max_it=100; % max number of iterations
opt.tol=1e-8; % tolerance

% registering Y to X
[Transform] = cpd_register(X,Y,opt);

figure,cpd_plot_iter(X, Y); title('Before');
figure,cpd_plot_iter(X, Transform.Y); title('After registering Y to X');

% compute error as mean minimum distance between points
error_before = computeError(X, Y)
error_after = computeError(X, Transform.Y)

    function [E, D] = computeError(X, Y)
        %compute for each x the euclidean distance to the closest y
        for i = 1 : size(X,1)
            D(i) = min(sqrt(sum((repmat(X(i,:), size(Y,1), 1) - Y).^2,2)));
        end
        E = mean(D);
    end
end

Any idea, what's going on? I noticed that in cpd_rigid.m the value of ntol can be negative: ntol=(L-L_old)/L; I'm not sure if this is intended. I changed it to ntol=abs((L-L_old)/L). The registration result is only slightly better, though.

Another minor thing. In cpd_register.m the line
if ~isfield(opt,'viz') || isempty(opt.tol), opt.viz = 1; end;
should be
if ~isfield(opt,'viz') || isempty(opt.viz), opt.viz = 1; end;

Best,
Steffen.