function    [TParams, time] = rigidRegisterEngine3D_EMST(vol1, vol2,gamma,numOfLevels,TParams_init, angleStep, numOfMaxIters)
% USAGE:
% [TParams, time] = rigidRegisterEngine3D_EMST(vol1,vol2,gamma = 1.9,numOfLevels = ?,TParams_init = [0 0 0 0 0 0], angleStep = 2.0, numOfMaxIters = 1000)
% Mert Rory Sabuncu (c) 2005-2008 CSAIL, MIT

% A multi-resolution, multi-modal rigid registration algorithm that uses an
% EMST (Eculidean Minimum Spanning Tree) based similarity measure

% INPUTS:
% vol1: 3D volume (double precision)
% vol2: 3D volume (double precision) of the same size as vol1
% gamma: The exponential weight in the EMST weight function -- gamma = 2*(1 - alpha), where alpha is the alpha of the Renyi entropy.
% numOfLevels: Number of Levels in the image pyramid
% TParams_init: Initial guess of the Transformation parameters
% [tx, ty, tz, theta, omega, phi]
% tx, ty, tz: the amount of translation along the first, second and third
% dimensions (in voxels)
% theta, omega, phi: Euler angles in degrees(!)
%
% angleStep: The step you take in gradient descent: it is in degrees and
% for instance 1 would correspond to a 1 degree accuracy in registration

% numOfMaxIters: number of maximum iterations in gradient descent


% **********************************
% needs:
% 1) volumeRigidRegisterEMSTwStochSampling.dll (or mex) -- which is
% generated from volumeRigidRegisterEMSTwStochSampling.c
% **********************************

[m1, n1, k1] = size(vol1);
[m2, n2, k2] = size(vol2);

m = min(m1,m2);
n = min(n1,n2);
k = min(k1,k2);

% Default values:
if (nargin < 7)
    numOfMaxIters = 1000;
end

if (nargin < 6)
    angleStep = 2.0;
end;

if (nargin < 5)
    TParams_init = zeros(1,6);
end
if (nargin < 4)
    numOfLevels = floor(log(min(m,n,k))/log(2)) - 3;
end
if (nargin < 4)
    gamma = 1.9;
end
if (nargin < 2)
    error('Not enough input parameters');
end

Lmin = 1;
weight1 = 10;

TParams = TParams_init;

% Need to readjust initial transformation parameters to the lo-res in the
% pyramid
TParams(1) = TParams(1)/2^(numOfLevels-1);
TParams(2) = TParams(2)/2^(numOfLevels-1);
TParams(3) = TParams(3)/2^(numOfLevels-1);

gamma_max = gamma;
tic

%IMAGE PYRAMID -- work on each level

cur_gamma = gamma_max;
angleStep_cur = angleStep;

pyramid1 = cell(numOfLevels, 1);
pyramid2 = cell(numOfLevels, 1);
pyramid1{1,1} = weight1 * renorm(vol1);
pyramid2{1,1} = renorm(vol2);

for level = 2:1:numOfLevels
    size1 = size(pyramid1{level-1,1});
    size1 = size1 - mod(size1,2);
    pyramid1{level,1} = pyramid1{level-1,1}(1:2:size1(1),1:2:size1(2), 1:2:size1(3));
    pyramid1{level,1} = pyramid1{level,1} + pyramid1{level-1,1}(2:2:size1(1),1:2:size1(2),1:2:size1(3));
    pyramid1{level,1} = pyramid1{level,1} + pyramid1{level-1,1}(1:2:size1(1),2:2:size1(2),1:2:size1(3));
    pyramid1{level,1} = pyramid1{level,1} + pyramid1{level-1,1}(1:2:size1(1),1:2:size1(2),2:2:size1(3));
    pyramid1{level,1} = pyramid1{level,1} + pyramid1{level-1,1}(1:2:size1(1),2:2:size1(2),2:2:size1(3));
    pyramid1{level,1} = pyramid1{level,1} + pyramid1{level-1,1}(2:2:size1(1),1:2:size1(2),2:2:size1(3));
    pyramid1{level,1} = pyramid1{level,1} + pyramid1{level-1,1}(2:2:size1(1),2:2:size1(2),1:2:size1(3));
    pyramid1{level,1} = pyramid1{level,1} + pyramid1{level-1,1}(2:2:size1(1),2:2:size1(2),2:2:size1(3));
    pyramid1{level,1} = pyramid1{level,1}/8;

    size2 = size(pyramid2{level-1,1});
    size2 = size2 - mod(size2,2);
    pyramid2{level,1} = pyramid2{level-1,1}(1:2:size2(1),1:2:size2(2), 1:2:size2(3));
    pyramid2{level,1} = pyramid2{level,1} + pyramid2{level-1,1}(2:2:size2(1),1:2:size2(2),1:2:size2(3));
    pyramid2{level,1} = pyramid2{level,1} + pyramid2{level-1,1}(1:2:size2(1),2:2:size2(2),1:2:size2(3));
    pyramid2{level,1} = pyramid2{level,1} + pyramid2{level-1,1}(1:2:size2(1),1:2:size2(2),2:2:size2(3));
    pyramid2{level,1} = pyramid2{level,1} + pyramid2{level-1,1}(1:2:size2(1),2:2:size2(2),2:2:size2(3));
    pyramid2{level,1} = pyramid2{level,1} + pyramid2{level-1,1}(2:2:size2(1),1:2:size2(2),2:2:size2(3));
    pyramid2{level,1} = pyramid2{level,1} + pyramid2{level-1,1}(2:2:size2(1),2:2:size2(2),1:2:size2(3));
    pyramid2{level,1} = pyramid2{level,1} + pyramid2{level-1,1}(2:2:size2(1),2:2:size2(2),2:2:size2(3));
    pyramid2{level,1} = pyramid2{level,1}/8;
end
tic
for level = numOfLevels:-1:Lmin
    qLevels  =  10 * (numOfLevels - level + 1);
    qLevels = max(qLevels, 10);
    
    if (level == Lmin)
        numOfMaxIters = 10;
    end
    numOfSamples = max(round(0.0001*numel(pyramid1{level,1})), 250);
    TParams = volumeRigidRegisterEMSTwStochSampling(pyramid1{level,1}, pyramid2{level,1}, TParams,numOfSamples, angleStep_cur, qLevels, cur_gamma, numOfMaxIters);
    cur_gamma = cur_gamma - 0.1;
    cur_gamma = max(min(cur_gamma, 1.9), 0.1);
    angleStep_cur = angleStep_cur / 2.0;
    
    numOfMaxIters = round(numOfMaxIters*0.5);
    angleStep     = 0.9 * angleStep;
    if (level ~= 1)
        TParams(1) = 2*TParams(1);
        TParams(2) = 2*TParams(2);
        TParams(3) = 2*TParams(3);
    end;
end
if (Lmin ~= 1)
    TParams(1) = 2^(Lmin-1)*TParams(1);
    TParams(2) = 2^(Lmin-1)*TParams(2);
    TParams(3) = 2^(Lmin-1)*TParams(3);
end
time = toc;
return