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;