圖像紋理合成及紋理傳輸算法學習（附源碼）。

時間 2019-11-18

標籤圖像紋理合成傳輸算法學習源碼简体版

原文原文鏈接

有2到3年沒有逛CodeProject了，上班一時無聊，就翻翻這個比較有名的國外網站，在其Articles » Multimedia » General Graphics » Graphics一欄看到exture Transfer using Efros & Freeman's Image Quilting Algorithm頓感興趣，其效果頗有特點，又激起了我很久沒寫代碼的心，想一想在這個算法上大概前先後後思索了1個多星期，雖然最終仍是沒有獲得我想要的結果，但並不全無收穫，至此國慶佳節加班之際抽空總結一下，以便慰藉我孤獨的心靈。html

紋理圖像的合成算法在早期的Photoshop中我記得是有一個單獨的功能的，在後來的版本中不知道爲何被取消了，印象中他能將只有幾顆小草的圖片生長成不少小草，而且基本看不出什麼瑕疵和不天然的地方，那麼CodeProject上的這個Quilting算法也有相似的能力，首先貼幾張處理的結果圖欣賞一下。web

　　　　原始紋理圖像　　　　　　　　　　　　　　　　　　　　　　　　由小紋理合成的大紋理圖像算法

那麼簡單的描述下這篇文章所使用的算法的過程吧。編程

算法須要的輸入：原始的紋理圖像(W * H)，塊的大小TileSize，重疊部分的大小Overlap。函數

　第一步：咱們從原始的紋理圖像中一個隨機的抽取一個小塊，放到目標圖像的左上角。post

無需解釋，其中黑色的部分表示目標圖像還沒有處理的部分。大數據

第二步：按照從左到右，從上到下，抽取出塊重疊的部分的數據，並計算這部分數據和原始紋理圖像中各塊的類似度。優化

水平重疊垂直重疊水平和垂直部分具備重疊網站

　　要達到合成天然，在有了第一個塊以後，其周邊的塊應該從原圖中選擇和其最類似的部位，所以，咱們從已經合成的塊中抽取部分重疊的數據，如上圖所示，紅色方框內，沒有紅色網格線的部分即爲重疊的內容，在整個的處理過程當中，會有三種狀況出現，分別如上圖所示,即（1）只有垂直方向有部分重疊；（2）只有水平發方向有部分重疊；（3）水平和垂直方向都有重疊。ui

爲了選擇合適的下一個塊，咱們結算這些重疊的塊和原圖中的全部塊的類似度，塊和塊之間的類似度能夠有不少準則計算，最經常使用的莫過於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;
}