基于像素的二维图像配准简要介绍

function [refimgout] = warpimgrigid2d(refimg, rigidmatrix, movimg, interptype)
% 将movimg图像根据刚性配准矩阵warp到refimg, 计算ref img上个像素点的值

if nargin < 3
    error('input params are too small');
end
if nargin == 3
    interptype = 'cubic';
end
    
size2dref = [refimg.dim(1), refimg.dim(2)];
% 插值的点 movimg变换到refimg
xi = zeros(size2dref);  %height方向
yi = zeros(size2dref);  %width方向
size2dmov = [movimg.dim(1), movimg.dim(2)];
% xmov = zeros(size2dmov);   % 原始的点 movimg
% ymov = zeros(size2dmov);

x = refimg.origin(1) :  refimg.pixelspacing(1) : ( refimg.origin(1) + (refimg.dim(1)-1) * refimg.pixelspacing(1) );
y = refimg.origin(2) :  refimg.pixelspacing(2) : ( refimg.origin(2) + (refimg.dim(2)-1) * refimg.pixelspacing(2) );

% xx = repmat(x', [1, size2dref(2)]);
% yy = repmat(y, [size2dref(1), 1]);
[yy, xx] = meshgrid(y, x);

for i = 1 : size2dref(1)
    for j = 1 : size2dref(2)
        xi(i, j) = rigidmatrix(1, :) * [xx(i, j), yy(i, j), 1]';
        yi(i, j) = rigidmatrix(2, :) * [xx(i, j), yy(i, j), 1]';
    end
end

xl = movimg.origin(1) :  movimg.pixelspacing(1) : movimg.origin(1) + (movimg.dim(1)-1) * movimg.pixelspacing(1);
yl = movimg.origin(2) :  movimg.pixelspacing(2) : movimg.origin(2) + (movimg.dim(2)-1) * movimg.pixelspacing(2);
% xmov = repmat(xl', [1,size2dmov(2)]);
% ymov = repmat(yl, [size2dmov(1), 1]);
[ymov, xmov] = meshgrid(yl, xl);

imgslice = interp2(ymov, xmov, movimg.img, yi, xi, interptype, 0);

refimgout.img = max(imgslice,0);
refimgout.dim = refimg.dim;
refimgout.origin = refimg.origin;
refimgout.pixelspacing = refimg.pixelspacing;