5

我有一个 3D 体积和一个 2D 图像以及两者之间的近似映射（没有倾斜的仿射变换、已知的缩放、旋转和平移近似已知并且需要拟合）。因为这个映射有错误，我想进一步将 2D 图像注册到 3D 体积。我以前没有为注册目的编写代码，但是因为我找不到任何程序或代码来解决这个问题，所以我想尝试这样做。我相信注册的标准是优化互信息。我认为这在这里也很合适，因为两个图像之间的强度不相等。所以我想我应该做一个转换函数，一个互信息函数和一个优化函数。

根据一篇文章，我确实在两年前的 mathworks线程上找到了一些 Matlab 代码。OP 报告说她设法让代码正常工作，但我不知道她是如何做到这一点的。在 matlab 的 IP 包中也有一个实现，但我没有那个包，而且似乎没有octave的等效项。SPM是一个使用 matlab 并实现了注册的程序，但不处理 2d 到 3d 的注册。在文件交换中，有一种使用互信息注册两个 2D 图像的蛮力方法。

她所做的是将多平面重建函数和相似性/误差函数传递到最小化算法中。但细节我不太明白。也许重新开始会更好：

load mri; volume = squeeze(D);

phi = 3; theta = 2; psi = 5; %some small angles
tx = 1; ty = 1; tz = 1; % some small translation
dx = 0.25, dy = 0.25, dz = 2; %different scales
t = [tx; ty; tz];
r = [phi, theta, psi]; r = r*(pi/180);
dims = size(volume);
p0 = [round(dims(1)/2);round(dims(2)/2);round(dims(3)/2)]; %image center
S = eye(4); S(1,1) = dx; S(2,2) = dy; S(3,3) = dz;

Rx=[1 0 0 0;
    0 cos(r(1)) sin(r(1)) 0;
    0 -sin(r(1)) cos(r(1)) 0;
    0 0 0 1];
Ry=[cos(r(2)) 0 -sin(r(2)) 0;
    0 1 0 0;
    sin(r(2)) 0 cos(r(2)) 0;
    0 0 0 1];
Rz=[cos(r(3)) sin(r(3)) 0 0;
    -sin(r(3)) cos(r(3)) 0 0;
    0 0 1 0;
    0 0 0 1];
R = S*Rz*Ry*Rx;
%make affine matrix to rotate about center of image
T1 = ( eye(3)-R(1:3,1:3) ) * p0(1:3);
T = T1 + t; %add translation
A = R;
A(1:3,4) = T;
Rold2new = A;
Rnew2old = inv(Rold2new);

%the transformation
[xx yy zz] = meshgrid(1:dims(1),1:dims(2),1:1);
coordinates_axes_new = [xx(:)';yy(:)';zz(:)'; ones(size(zz(:)))'];
coordinates_axes_old = Rnew2old*coordinates_axes_new;
Xcoordinates = reshape(coordinates_axes_old(1,:), dims(1), dims(2), dims(3));
Ycoordinates = reshape(coordinates_axes_old(2,:), dims(1), dims(2), dims(3));
Zcoordinates = reshape(coordinates_axes_old(3,:), dims(1), dims(2), dims(3));

%interpolation/reslicing
method = 'cubic'; 
slice= interp3(volume, Xcoordinates, Ycoordinates, Zcoordinates, method);
%so now I have my slice for which I would like to find the correct position

% first guess for A
A0 = eye(4); A0(1:3,4) = T1; A0(1,1) = dx; A0(2,2) = dy; A0(3,3) = dz; 
% this is pretty close to A
% now how would I fit the slice to the volume by changing A0 and examining some similarity measure?
% probably maximize mutual information?
% http://www.mathworks.com/matlabcentral/fileexchange/14888-mutual-information-computation/content//mi/mutualinfo.m
4

1 回答 1

2

好的，我希望其他人的方法，这可能会比我的更好，因为我以前从未做过任何优化或注册。所以我等待 Knedlsepps 的赏金几乎完成。但我确实有一些代码现在可以工作。它将找到一个局部最优值，因此初始猜测必须是好的。我自己编写了一些函数，按原样从文件交换中获取了一些函数，并从文件交换中广泛编辑了一些其他函数。现在我将所有代码放在一起作为示例工作，旋转已关闭，将尝试纠正它。我不确定示例和我的原始代码之间的代码差异在哪里，一定是在替换一些变量和数据加载方案时打错了。

我所做的是获取起始仿射变换矩阵，将其分解为正交矩阵和上三角矩阵。然后我假设正交矩阵是我的旋转矩阵，所以我从中计算欧拉角。我直接从仿射矩阵中获取平移，并且如问题中所述，我假设我知道缩放矩阵并且没有剪切。所以我有仿射变换的所有自由度，我的优化函数会改变并从中构造一个新的仿射矩阵，将其应用于体积并计算互信息。matlab 优化函数 patternsearch 然后最小化 1-MI/MI_max。

在我的多模态大脑图像的真实数据上使用它时，我注意到它在大脑提取图像上的效果要好得多，因此在去除颅骨和颅骨外的组织后。

%data
load mri; volume = double(squeeze(D));

%transformation parameters
phi = 3; theta = 1; psi = 5; %some small angles
tx = 1; ty = 1; tz = 3; % some small translation
dx = 0.25; dy = 0.25; dz = 2; %different scales
t = [tx; ty; tz];
r = [phi, theta, psi]; r = r*(pi/180);
%image center and size
dims = size(volume);
p0 = [round(dims(1)/2);round(dims(2)/2);round(dims(3)/2)]; 
%slice coordinate ranges
range_x = 1:dims(1)/dx;
range_y = 1:dims(2)/dy;
range_z = 1;

%rotation 
R = dofaffine([0;0;0], r, [1,1,1]);
T1 = ( eye(3)-R(1:3,1:3) ) * p0(1:3); %rotate about p0
%scaling 
S = eye(4); S(1,1) = dx; S(2,2) = dy; S(3,3) = dz;
%translation
T = [[eye(3), T1 + t]; [0 0 0 1]];
%affine 
A = T*R*S;

% first guess for A
r00 = [1,1,1]*pi/180;
R00 = dofaffine([0;0;0], r00, [1 1 1]);
t00 = T1 + t + ( eye(3) - R00(1:3,1:3) ) * p0;
A0 = dofaffine( t00, r00, [dx, dy, dz] );
[ t0, r0, s0 ] = dofaffine( A0 );
x0 = [ t0.', r0, s0 ];

%the transformation
slice = affine3d(volume, A, range_x, range_y, range_z, 'cubic');
guess = affine3d(volume, A0, range_x, range_y, range_z, 'cubic');

%initialisation
Dt = [1; 1; 1]; 
Dr = [2 2 2].*pi/180;
Ds = [0 0 0];
Dx = [Dt', Dr, Ds];
%limits
LB = x0-Dx;
UB = x0+Dx;
%other inputs
ref_levels = length(unique(slice));
Qref = imquantize(slice, ref_levels);
MI_max = MI_GG(Qref, Qref);
%patternsearch options
options = psoptimset('InitialMeshSize',0.03,'MaxIter',20,'TolCon',1e-5,'TolMesh',5e-5,'TolX',1e-6,'PlotFcns',{@psplotbestf,@psplotbestx});

%optimise
[x2, MI_norm_neg, exitflag_len] = patternsearch(@(x) AffRegOptFunc(x, slice, volume, MI_max, x0), x0,[],[],[],[],LB(:),UB(:),options);

%check
p0 = [round(size(volume)/2).'];
R0 = dofaffine([0;0;0], x2(4:6)-x0(4:6), [1 1 1]);
t1 = ( eye(3) - R0(1:3,1:3) ) * p0;
A2 = dofaffine( x2(1:3).'+t1, x2(4:6), x2(7:9) ) ;
fitted = affine3d(volume, A2, range_x, range_y, range_z, 'cubic');
overlay1 = imfuse(slice, guess);
overlay2 = imfuse(slice, fitted);
figure(101); 
ax(1) = subplot(1,2,1); imshow(overlay1, []); title('pre-reg')
ax(2) = subplot(1,2,2); imshow(overlay2, []); title('post-reg');
linkaxes(ax);

function normed_score = AffRegOptFunc( x, ref_im, reg_im, MI_max, x0 )
    t = x(1:3).';
    r = x(4:6);
    s = x(7:9);

    rangx = 1:size(ref_im,1);
    rangy = 1:size(ref_im,2);
    rangz = 1:size(ref_im,3);

    ref_levels =  length(unique(ref_im));
    reg_levels =  length(unique(reg_im));

    t0 = x0(1:3).';
    r0 = x0(4:6);
    s0 = x0(7:9);
    p0 = [round(size(reg_im)/2).'];
    R = dofaffine([0;0;0], r-r0, [1 1 1]);
    t1 = ( eye(3) - R(1:3,1:3) ) * p0;
    t = t + t1;
    Ap = dofaffine( t, r, s );

    reg_im_t = affine3d(reg_im, A, rangx, rangy, rangz, 'cubic');

    Qref = imquantize(ref_im, ref_levels);
    Qreg = imquantize(reg_im_t, reg_levels);

    MI = MI_GG(Qref, Qreg);

    normed_score = 1-MI/MI_max;
end

function [ varargout ] = dofaffine( varargin )
% [ t, r, s ] = dofaffine( A )
% [ A ] = dofaffine( t, r, s )
if nargin == 1
    %affine to degrees of freedom (no shear)
    A = varargin{1};

    [T, R, S] = decompaffine(A);
    r = GetEulerAngles(R(1:3,1:3));
    s = [S(1,1), S(2,2), S(3,3)];
    t = T(1:3,4);

    varargout{1} = t;
    varargout{2} = r;
    varargout{3} = s;
elseif nargin == 3
    %degrees of freedom to affine (no shear)
    t = varargin{1};
    r = varargin{2};
    s = varargin{3};

    R = GetEulerAngles(r); R(4,4) = 1;
    S(1,1) = s(1); S(2,2) = s(2); S(3,3) = s(3); S(4,4) = 1;
    T = eye(4); T(1,4) = t(1); T(2,4) = t(2); T(3,4) = t(3);
    A = T*R*S;

    varargout{1} = A;
else
    error('incorrect number of input arguments');
end
end

function [ T, R, S ] = decompaffine( A )
    %I assume A = T * R * S
    T = eye(4);
    R = eye(4);
    S = eye(4);
    %decompose in orthogonal matrix q and upper triangular matrix r
    %I assume q is a rotation matrix and r is a scale and shear matrix
    %matlab 2014 can force real solution
    [q r] = qr(A(1:3,1:3));
    R(1:3,1:3) = q;
    S(1:3,1:3) = r;

    % A*S^-1*R^-1 = T*R*S*S^-1*R^-1 = T*R*I*R^-1 = T*R*R^-1 = T*I = T
    T = A*inv(S)*inv(R);
    t = T(1:3,4);
    T = [eye(4) + [[0 0 0;0 0 0;0 0 0;0 0 0],[t;0]]];
end

function [varargout]= GetEulerAngles(R)

assert(length(R)==3)

dims = size(R);

    if min(dims)==1
        rx = R(1); ry = R(2); rz = R(3);
        R = [[                           cos(ry)*cos(rz),                          -cos(ry)*sin(rz),          sin(ry)];...
             [ cos(rx)*sin(rz) + cos(rz)*sin(rx)*sin(ry), cos(rx)*cos(rz) - sin(rx)*sin(ry)*sin(rz), -cos(ry)*sin(rx)];...
             [ sin(rx)*sin(rz) - cos(rx)*cos(rz)*sin(ry), cos(rz)*sin(rx) + cos(rx)*sin(ry)*sin(rz),  cos(rx)*cos(ry)]];
        varargout{1} = R;
    else
        ry=asin(R(1,3));
        rz=acos(R(1,1)/cos(ry));
        rx=acos(R(3,3)/cos(ry));

        if nargout > 1 && nargout < 4
            varargout{1} = rx;
            varargout{2} = ry;
            varargout{3} = rz;
        elseif nargout == 1
            varargout{1} = [rx ry rz];
        else
            error('wrong number of output arguments');
        end
    end
end
于 2015-02-05T11:11:57.090 回答