图像纹理合成及纹理传输算法学习（附源码）。

有2到3年没有逛CodeProject了，上班一时无聊，就翻翻这个比较有名的国外网站，在其Articles » Multimedia » General Graphics » Graphics一栏看到exture Transfer using Efros & Freeman's Image Quilting Algorithm顿感兴趣，其效果很有特色，又激起了我好久没写代码的心，想想在这个算法上大概前前后后思索了1个多星期，虽然最终还是没有得到我想要的结果，但并不全无收获，至此国庆佳节加班之际抽空总结一下，以便慰藉我孤独的心灵。

纹理图像的合成算法在早期的Photoshop中我记得是有一个单独的功能的，在后来的版本中不知道为什么被取消了，印象中他能将只有几颗小草的图片生长成很多小草，并且基本看不出什么瑕疵和不自然的地方，那么CodeProject上的这个Quilting算法也有类似的能力，首先贴几张处理的结果图欣赏一下。

　　　　原始纹理图像　　　　　　　　　　　　　　　　　　　　　　　　由小纹理合成的大纹理图像

那么简单的描述下这篇文章所使用的算法的过程吧。

算法需要的输入：原始的纹理图像(W * H)，块的大小TileSize，重叠部分的大小Overlap。

　第一步：我们从原始的纹理图像中一个随机的抽取一个小块，放到目标图像的左上角。

无需解释，其中黑色的部分表示目标图像尚未处理的部分。

第二步：按照从左到右，从上到下，抽取出块重叠的部分的数据，并计算这部分数据和原始纹理图像中各块的相似度。

水平重叠垂直重叠水平和垂直部分具有重叠

　　要达到合成自然，在有了第一个块之后，其周边的块应该从原图中选择和其最相似的部位，因此，我们从已经合成的块中抽取部分重叠的数据，如上图所示，红色方框内，没有红色网格线的部分即为重叠的内容，在整个的处理过程中，会有三种情况出现，分别如上图所示,即（1）只有垂直方向有部分重叠；（2）只有水平发方向有部分重叠；（3）水平和垂直方向都有重叠。

为了选择合适的下一个块，我们结算这些重叠的块和原图中的所有块的相似度，块和块之间的相似度可以有很多准则计算，最常用的莫过于SSD(Sum of Squared Distance)，这里我们也采用这种准则。

　　对于第一种只有垂直方向有部分重叠的情况，块的可能性有（W - Overlap）* （H - TileSize）种，而只有水平发方向有部分重叠时，块的可能性有（W - TileSize)* (H - Overlap)种，当水平和垂直方向都有重叠时，只有（W - TileSize）* （H - TileSize）种块的选择方案。

对于水平和垂直有交叉的情况，需要分别计算水平重叠的相似性A，垂直重叠的相似性B以及交叉部位的相似性AB，最终的相似性由A + B - AB决定。

计算完相似性后，一般情况下，我们可取相似度最小的块为下一个块，当然为了随机性更强，也可以取相似性序列中前N个最小值种的某个块为选中的块。

第三步：如果对选中的块直接进行拼贴到目标图中，则很明显两个块之间由过渡不会太自然，一种较好的方式就是在选中的块和重叠的部分找到一条路径，在该路径的两侧像素的距离和最小，如果使用暴力的方式去寻找这条路径，则是非常耗时和不必要的，文章介绍使用了贪心算法来寻找这条路径。

我们拿水平重叠的块来说明问题，首先找到第一行最小距离和的像素的位置，假定为(x,1),接着我们搜索第二行的（x - 1, 2), （x, 2), （x + 1, 2)三个位置的距离和的最小值，假定是（x - 1, 2)是第二行的最小值，则第三行需要搜索的位置为（x - 2, 3), （x - 1, 3), （x, 3)，依类类推直到最后一行。

其实路径必然是连续的，因此这种只搜索向下的三领域的做法是完全合理的。

一般情况下，做图像算法会按照从左到右，从上到下的顺序进行处理，这样能够充分利用cache line的优势，如果按照从上到下，然后在从左到右的方式来处理，虽然逻辑和结果是一样的，通常会需要更多的时间，因此，在计算垂直的块的路径时，我们可以借用水平块的算法，只要对垂直的块的数据先进行转置，处理转置的数据得到对应的数据，然后在把这个数据转置就得到了垂直块的结果。注意这里是用的转置，而不是旋转，因为旋转对于该问题是的到的结果是镜像的，所以必须注意。

在水平和垂直部分具有重叠的块的计算时，我们是分别计算垂直和水平的路径，然后取两个路径的交集，计算过程分别如下图。

第四步：按照路径的位置贴入新的数据。

编程技巧：

（1）在整个的计算中，计算SSD是最为耗时的，因此对其优化是调高程序效率的关键。

　　我们看一段matlab的代码，如下所示：

%function Z = ssd(X, Y)
%Computes the sum of squared distances between X and Y for each possible
% overlap of Y on X.  Y is thus smaller than X
%
%Inputs:
%   X - larger image
%   Y - smaller image
%
%Outputs:
%   Each pixel of Z contains the ssd for Y overlaid on X at that pixel

function Z = ssd(X, Y)

K = ones(size(Y,1), size(Y,2));

for k=1:size(X,3),
    A = X(:,:,k);
    B = Y(:,:,k);
    
    a2 = filter2(K, A.^2, 'valid');
    b2 = sum(sum(B.^2));
    ab = filter2(B, A, 'valid').*2;

    if( k == 1 )
        Z = ((a2 - ab) + b2);
    else
        Z = Z + ((a2 - ab) + b2);
    end;
end;

　　这里的优化时通过卷积实现的，因为SSD的计算就是(a-b)^2的累积和，展开为a^2 + b^2 - 2ab，其中a为固定的图像，那么a^2则为一定值，b^2是不断变化的，但是也可以用积分图之类的算法实现快速效果，唯一的耗时时2ab项，可以用卷积来实现。

快速的卷积算法可以通过FFT来实现，也可以借鉴我在图像处理中任意核卷积(matlab中conv2函数)的快速实现一文中提到的用SSE的方式实现，鉴于这里的卷积更有其特殊性，他是2幅图像之间的卷积，因此参与计算的都是byte类型，则可以通过选择更为合适的SSE函数来进一步提高效率，这里要感谢有关高手提供的一段SSE代码。

/// <summary>
/// 基于SSE的字节数据的乘法。
/// </summary>
/// <param name="Kernel">需要卷积的核矩阵。</param>
/// <param name="Conv">卷积矩阵。</param>
/// <param name="Length">矩阵所有元素的长度。</param>
/// <remarks> 1: 使用了SSE优化。
/// <remarks> 2: 感谢采石工（544617183）提供的SSE代码<。</remarks>
///    https://msdn.microsoft.com/zh-cn/library/t5h7783k(v=vs.90).aspx

int MultiplySSE(unsigned char* Kernel, unsigned char * Conv, int Length)
{
    int Y, Sum;
    __m128i vsum = _mm_set1_epi32(0);
    __m128i vk0 = _mm_set1_epi8(0);
    for (Y = 0; Y <= Length - 16; Y += 16)
    {
        __m128i v0 = _mm_loadu_si128((__m128i*)(Kernel + Y));                //    对应movdqu指令,不需要16字节对齐
        __m128i v0l = _mm_unpacklo_epi8(v0, vk0);                            
        __m128i v0h = _mm_unpackhi_epi8(v0, vk0);                            //    此两句的作用是把他们分别加载到两个128位寄存器中，供下面的_mm_madd_epi16的16位SSE函数调用（vk0的作用主要是把高8位置0）                
        __m128i v1 = _mm_loadu_si128((__m128i*)(Conv + Y));
        __m128i v1l = _mm_unpacklo_epi8(v1, vk0);
        __m128i v1h = _mm_unpackhi_epi8(v1, vk0);
        vsum = _mm_add_epi32(vsum, _mm_madd_epi16(v0l, v1l));                //    _mm_madd_epi16 可以一次性进行8个16位数的乘法，然后把两个16的结果在加起来，放到一个32数中， r0 := (a0 * b0) + (a1 * b1)，详见 https://msdn.microsoft.com/zh-cn/library/yht36sa6(v=vs.90).aspx
        vsum = _mm_add_epi32(vsum, _mm_madd_epi16(v0h, v1h));
    }
    for (; Y <= Length - 8; Y += 8)
    {
        __m128i v0 = _mm_loadl_epi64((__m128i*)(Kernel + Y));
        __m128i v0l = _mm_unpacklo_epi8(v0, vk0);
        __m128i v1 = _mm_loadl_epi64((__m128i*)(Conv + Y));
        __m128i v1l = _mm_unpacklo_epi8(v1, vk0);
        vsum = _mm_add_epi32(vsum, _mm_madd_epi16(v0l, v1l));
    }
    vsum = _mm_add_epi32(vsum, _mm_srli_si128(vsum, 8));
    vsum = _mm_add_epi32(vsum, _mm_srli_si128(vsum, 4));
    Sum = _mm_cvtsi128_si32(vsum);                                            //    MOVD函数，Moves the least significant 32 bits of a to a 32-bit integer：   r := a0

    for ( ; Y < Length; Y++)
    {
        Sum += Kernel[Y] * Conv[Y];
    }
    return Sum;
}