[数字图像处理]简单的几何学图像变换与图像配准

       1.图像的几何学变换

       之前的博文里，我简单的介绍了图像的放大与缩小。放大与缩小也算是图像的几何学变换，本文介绍了其他的几何学变换，包括旋转、水平倾斜和垂直倾斜（当然，还有水平移动与垂直移动。这些变换很简单，不需要插值，所以这里就不着重介绍了）。

       假设输入图像为G(u,v)，其变换后的图像为F(x,y)。其变化的方法，如下所示。

       图像的几何学变换，主要有两种向前映射与向后映射。

       1.1向前映射

       所谓向前映射，就是从输入图像为g(u,v)的(0,0)点开始，将g(u,v)遍历一遍，依次计算g(u,v)变换后的坐标。当然，计算出来的坐标不会是整数（很大程度上不会是整数），就像下图所示那样。

      借由上图，我们重新理解一些向前映射。假设图像的一部分这5个点，经过变换，我们得到了右边的5个点。其实，这个5个点经过变换之后，计算出来的坐标并不是整数。 假设现在遍历到 这个点（图中以蓝色表示），通过计算得到的点是。我们将的坐标进行四舍五入，将的值赋值给（这种处理相当于选择最近的点，属于最初级的插值方法）。

       使用上述的向前映射的思想，将图像旋转，我们能得到如下效果。

       可以看出来，所得到的结果非常“斑驳”。其原因是，我们将g(u,v)进行变换，计算之后，所得到的坐标四舍五入之后，有的点没有被赋值，而有点被赋值了多次。这就产生了“空穴”，所以让变换后的图像非常斑驳。

       我看了一些文章，基本到这里都会说，“由于向前映射会产生空穴，所以向前映射一般不使用”。诸如此类的话，其实这样理解不太对，产生空穴的根本原因是由于插值法的不恰当所产生的。使用向前映射，也可以不产生“空穴”，所使用的方法如下所示。

       我们可以通过最接近的四个点，去确定的值。这样的话，就可以不产生空穴。但是，这个方法实现起来很困难，所以，向前映射才一般不被使用。

       1.2 向后映射

       向后映射，就是将输出图像f(x,y)遍历一遍，然后计算输出(x,y)的时候，所需要g(x,y)的坐标。当然，这个数不一定是整数。为了方便理解，还是看下图。

       假设，我们遍历到，通过计算，我们得到这样一个点。也就意味着，如果需要确定的灰度值，那么，我们需要点处的灰度值。这样，几何变换的问题又归结到了插值问题。最方便的办法就是选择最近的点，入上图的例子，的灰度值等于的灰度值。当然，还有更加“高级”点的方法，精度也还算可以。的灰度值的确定方法如下所示。

       如上图，我们使用四个点（），去确定的值。首先与之间使用线性插值，确定出的值。同样的方法，确定出的值。与之间，进行线性插值，就可以得到的值了。（图画的非常清楚了）

       使用向后映射，将图像旋转，我们可以得到如下结果。

       1.3 上述两部分的Matlab代码

close all;
clear all;
clc;

%% -----------Geometric_spatial_transform------------------
f = imread('./letter_T.tif');
f = mat2gray(f,[0 255]);

[M,N] = size(f);

%-----by forward mapping-----%
g_fm = zeros(M,N);
seta = -pi/8;

for v = (-M/2):1:(M/2)-1
    for w = (-N/2):1:(N/2)-1
        x = round(v*cos(seta) - w*sin(seta));
        y = round(v*sin(seta) + w*cos(seta));
        if (((y>=(-N/2))&&(y<=(N/2)))&&((x>=(-M/2))&&(x<=(M/2))))
            g_fm(x+(M/2)+1,y+(N/2)+1) = f(v+(M/2)+1,w+(N/2)+1);
        end
    end
end

figure();
subplot(1,2,1);
imshow(f,[0 1]);
xlabel('a).Original Image');

subplot(1,2,2);
imshow(g_fm,[0 1]);
xlabel('b).Ruselt of Geometric spatial transform');

%-----by inverse mapping----
g_im = zeros(M,N);
seta = -pi/8;

for x = (-M/2):1:(M/2)-1
    for y = (-N/2):1:(N/2)-1
        v = x*cos(-seta) - y*sin(-seta);
        w = x*sin(-seta) + y*cos(-seta);
        if (((w>=(-N/2)+1)&&(w<=(N/2)-1))&&((v>=(-M/2)+1)&&(v<=(M/2)-1)))
            %g_im(x+(M/2)+1,y+(N/2)+1) = 1;%f(round(v+(M/2)+1),round(w+(N/2)+1));
            Q_11 = f(floor(v)+(M/2)+1,floor(w)+(N/2)+1);
            Q_21 =  f(floor(v)+(M/2)+1,ceil(w)+(N/2)+1);
            Q_12 =  f(ceil(v)+(M/2)+1,floor(w)+(N/2)+1);
            Q_22 =   f(ceil(v)+(M/2)+1,ceil(w)+(N/2)+1);
            
            R1 = (Q_21 - Q_11)*(w-floor(w)) + Q_11;
            R2 = (Q_22 - Q_12)*(w-floor(w)) + Q_12;

            g_im(x+(M/2)+1,y+(N/2)+1) = (R2-R1)*(v-floor(v)) + R1;
        end
    end
end

figure();
subplot(1,2,1);
imshow(f,[0 1]);
xlabel('a).Original Image');

subplot(1,2,2);
imshow(g_im,[0 1]);
xlabel('b).Ruselt of Geometric spatial transform');

      2.图像的配准

      图像的配准，常常用于超解析领域。当然了，作为基础学习，博主没学那么深。这次的图像配准，试图去还原被倾斜变换的图像。首先，先将图像进行两个方向的倾斜。

       我们这次的目的是，将图像b).还原为a).。这次变换相对简单，我们使用如下模型去拟合变换关系。

我们在原图与变换后的图像上，选择出四个标准点，然后带入方程，并求出系数。使用这个方程，将图像b).还原为a)。这个过程比较简单，下面是结果。另外，所选择的标准点，我已经在上图中用红圈标定出来了。

       可以看到，还原的效果非常的好（额，其实这也是应为畸变比较简单的缘故）。为了看出于原图的区别，我们做出了差分图像。可以看出来，还原不是完美的。

       至此，上个代码作为本文的结尾。额，由于特征点的选择是手工的，这个代码可能没有多少意义。（不要吐槽最后这一句话。一般的，在使用特征点做图象配准的时候，一般选择在某个实际的物体上放入某个实际的标志，然后通过检测这个实际的标志，去进行变换。【这里可以参考《Digital Image Processing》 Rafael C. Gonzalez / Richard E. Woods的第二章的相关内容，这里有叙述的！！】）

close all;
clear all;
clc;
%% --------------image registration---------------------------
f_Original = imread('./characters_test_pattern.tif');
f_Original = mat2gray(f_Original,[0 255]);
[M,N] = size(f_Original);

f = zeros(M+4,N+4);
for x = 1:M
    f(x+2,:) = [0 0 f_Original(x,:) 0 0];
end

[M,N] = size(f);

s_v = 0.3;
s_w = 0.02;

P = round(M+s_v*N);
Q = round(N+s_w*M);
g_1 = zeros(P,N);

for x = (-P/2):1:(P/2)-1
    for y = (-N/2):1:(N/2)-1
        v = x - s_v * y;
        w = y;
        if (((w>=(-N/2))&&(w<=(N/2)-1))&&((v>=(-M/2))&&(v<=(M/2)-1)))
            Q_11 = f(floor(v)+(M/2)+1,floor(w)+(N/2)+1);
            Q_21 =  f(floor(v)+(M/2)+1,ceil(w)+(N/2)+1);
            Q_12 =  f(ceil(v)+(M/2)+1,floor(w)+(N/2)+1);
            Q_22 =   f(ceil(v)+(M/2)+1,ceil(w)+(N/2)+1);
            
            R1 = (Q_21 - Q_11)*(w-floor(w)) + Q_11;
            R2 = (Q_22 - Q_12)*(w-floor(w)) + Q_12;

            g_1(x+(P/2)+1,y+(N/2)+1) = (R2-R1)*(v-floor(v)) + R1;
        end
    end
end

g = zeros(P,Q);

for x = (-P/2):1:(P/2)-1
    for y = (-Q/2):1:(Q/2)-1
        v = x;
        w = y - s_w * x;
        if (((w>=(-N/2))&&(w<=(N/2)-1))&&((v>=(-P/2))&&(v<=(P/2)-1)))
            Q_11 = g_1(floor(v)+(P/2)+1,floor(w)+(N/2)+1);
            Q_21 =  g_1(floor(v)+(P/2)+1,ceil(w)+(N/2)+1);
            Q_12 =  g_1(ceil(v)+(P/2)+1,floor(w)+(N/2)+1);
            Q_22 =   g_1(ceil(v)+(P/2)+1,ceil(w)+(N/2)+1);
            
            R1 = (Q_21 - Q_11)*(w-floor(w)) + Q_11;
            R2 = (Q_22 - Q_12)*(w-floor(w)) + Q_12;

            g(x+(P/2)+1,y+(Q/2)+1) = (R2-R1)*(v-floor(v)) + R1;
        end
    end
end

figure();
subplot(1,2,1);
imshow(f,[0 1]);
xlabel('a).Original Image');
hold on;
plot(621, 78,'ro','MarkerSize',7);
plot(116,113,'ro','MarkerSize',7);
plot( 85,649,'ro','MarkerSize',7);
plot(624,641,'ro','MarkerSize',7);

subplot(1,2,2);
imshow(g,[0 1]);
xlabel('b).Ruselt of Geometric spatial transform');
hold on;
plot(625,264,'ro','MarkerSize',7);
plot(116,147,'ro','MarkerSize',7);
plot( 96,674,'ro','MarkerSize',7);
plot(638,828,'ro','MarkerSize',7);
%% 
f_reg = zeros(M,N);

Original = [621-(N/2)  78-(M/2)  (621-(N/2))*(78-(M/2)) 1;
            116-(N/2) 113-(M/2) (116-(N/2))*(113-(M/2)) 1;
             85-(N/2) 649-(M/2)  (85-(N/2))*(649-(M/2)) 1;
            624-(N/2) 641-(M/2) (624-(N/2))*(641-(M/2)) 1];
        
Ruselt = [625-(Q/2);116-(Q/2);96-(Q/2);638-(Q/2)];

c = (Original)\Ruselt;  %C1~C4

Ruselt = [264-(P/2);147-(P/2);674-(P/2);828-(P/2)];

c = [c (Original)\Ruselt];  %C5~C8 wwwwww

for y = (-M/2):1:(M/2)-1
    for x = (-N/2):1:(N/2)-1
        w = c(1,1)*y + c(2,1)*x + c(3,1)*x*y + c(4,1);
        v = c(1,2)*y + c(2,2)*x + c(3,2)*x*y + c(4,2);
        if (((w>=(-Q/2))&&(w<=(Q/2)-1))&&((v>=(-P/2))&&(v<=(P/2)-1)))
            
            Q_11 = g(floor(v)+(P/2)+1,floor(w)+(Q/2)+1);
            Q_21 =  g(floor(v)+(P/2)+1,ceil(w)+(Q/2)+1);
            Q_12 =  g(ceil(v)+(P/2)+1,floor(w)+(Q/2)+1);
            Q_22 =   g(ceil(v)+(P/2)+1,ceil(w)+(Q/2)+1);
            
            R1 = (Q_21 - Q_11)*(w-floor(w)) + Q_11;
            R2 = (Q_22 - Q_12)*(w-floor(w)) + Q_12;
            
            f_reg(x+(N/2)+1,y+(M/2)+1) = (R2-R1)*(v-floor(v)) + R1;
            %g(round(v+(P/2)+1),round(w+(Q/2)+1)) = 0.5;
        end
    end
end

g_diff = abs(f - f_reg);

figure();
subplot(1,2,1);
imshow(f_reg,[0 1]);
xlabel('c).Ruselt of image registration');

subplot(1,2,2);
imshow(g_diff,[0 1]);
xlabel('d).Difference image');

      原文发于博客：http://blog.csdn.net/thnh169/