最近在看關于數字影像的知識點，目前在影像金字塔部分，實在是懶得用手作筆記了，就以其中比較出名的“高斯金字塔”和“拉普拉斯金字塔”為例，基于OpenCV的源代碼作決議存個檔；畢竟屬于基礎部分，以后有需要就當介面呼叫吧；有寫的不對需要改正的地方，還請大家指出，謝謝，

一、金字塔

1.何為影像金字塔？

影像金字塔是影像中多尺度表達形式的之一，最主要用于影像的分割，是一種以多解析度來解釋影像的有效但概念簡單的結構，影像金字塔最初用于機器視覺和影像壓縮，一幅影像的金字塔是一系列以金字塔形狀從低(下)到高(上)排列的，解析度逐步降低的影像集合，

具體以下示例圖和流程圖：

2.何為高斯金字塔？

高斯金字塔顧名思義，基于模板金字塔的構建的方式引入指定的低通高斯濾波（空間域）并結合下采樣的迭代方式進行層次之間的計算；（高斯金字塔底部是第0級）

而拉普拉斯金字塔恰好與高斯金字塔相反，拉普拉斯金字塔依賴于高斯金字塔構建的圖層；當高斯金字塔形成之時，我們需要對高斯金字塔的每一層進行上采樣，形成該圖層對應的近似預測影像，求去它們之間的差值，得到就是拉普拉斯金字塔該層次的影像；也就是說，高斯金字塔由下至上構建，而拉普拉斯金字塔又上至下構建（高斯金塔要比拉普拉斯金字塔多構建一次，用于頂層差值求取）；得到的拉普拉斯金字塔起到預測殘差的作用，

具體以下示例圖和流程圖：

二、金字塔構建流程

1. 高斯金字塔的構建

將原始影像與低通高斯濾波矩陣做卷積處理
利用基數為2的下采樣（洗掉偶數行和偶數列）得到維數減半的上一級影像

1.1 低通高斯濾波

"濾波"一詞借用于頻率域處理，"濾波"是指接受（通過）或拒絕一定的頻率的成分，低頻的濾波器成為低通濾波器，其最終的效果是模糊（平滑）一幅影像，

數值影像處理中，低通高斯濾波可以以不同的形式作用于空間域（線性）和頻率域，第一種是屬于空間濾波器（也被稱為高斯模板、高斯核、高斯掩模、高斯視窗），第二種方法是通過傅里葉變換后進行操作，本文涉及到的是第一種平滑空間濾波器，

常用的平滑線性空間濾波器有均值濾波以及高斯濾波等，均值濾波使用模板內所有像素的平均值代替模板中心像素灰度值，這種方法易收到噪聲的干擾，不能完全消除噪聲，只能相對減弱噪聲，且存在著不希望有的邊緣模糊負面效應；為了減少平滑處理中的模糊效應，得到更加自然的平滑效果，需要適當加大模板中心點的權重；隨著距離中心點的距離來增減控制權重占比大小，基于這樣的考慮形成的模板即為高斯模板，

常用的高斯模板核通常是3×3，5×5的奇數矩陣，根據OpenCV提供的代碼為例，采用5×5的模板如下（令右側矩陣為M，左側為歸一化系數 $1/\sum M(i,j)$ ）：

$kernel=\frac{1}{256}\begin{bmatrix} 1 & 4 & 6 & 4 & 1\\ 4& 16& 24 & 16& 4\\ 6& 24 &36 & 24 & 6\\ 4& 16& 24 & 16& 4 \\ 1 & 4 & 6 & 4 & 1 \end{bmatrix}_{5\times 5}$

該高斯模板核矩陣中的引數是通過二維高斯函式，即二維正態分布密度函式求得的，回憶一下（由一維延伸到二維），具體如下，

一維高斯函式，形式如下：

二維高斯函式，形式如下：

1.2 高斯模板的實作

高斯模板正是將連續的二維高斯函式的離散化表示，因此任意大小的高斯模板都可以通過建立一個 $(2k+1)*(2k+1)$ 的矩陣M(m×n)，得到其 $(i,j)$ 位置的引數可如下確定：

$M(i,j)=\frac{1}{2\pi \sigma ^{2}}\exp (-\frac{(i-k-1)^{2}+(j-k-1)^{2}}{2\sigma ^{2}})$

此公式基于Matlab的實作，因此i和j的起始為1非0（陣列索引必須為正整數或邏輯值）；基于VS的實作，i和j的起始為0（不要帶入-1的引數）；k為每個方向距離該方向中心的距離，即 $k=\frac{m-1}{2}=\frac{n-1}{2}$ ，

獲取OpenCV中smooth.dispatch.cpp的getGaussianKernel函式的原始碼：

Mat getGaussianKernel(int n, double sigma, int ktype)
{
    CV_Assert(n > 0);
    const int SMALL_GAUSSIAN_SIZE = 7;
    static const float small_gaussian_tab[][SMALL_GAUSSIAN_SIZE] =
    {
        {1.f},                                                                     //1×1
        {0.25f, 0.5f, 0.25f},                                                      //3×3
        {0.0625f, 0.25f, 0.375f, 0.25f, 0.0625f},                                  //5×5
        {0.03125f, 0.109375f, 0.21875f, 0.28125f, 0.21875f, 0.109375f, 0.03125f}   //7×7
    };

    const float* fixed_kernel = n % 2 == 1 && n <= SMALL_GAUSSIAN_SIZE && sigma <= 0 ?
        small_gaussian_tab[n>>1] : 0;

    CV_Assert( ktype == CV_32F || ktype == CV_64F );
    Mat kernel(n, 1, ktype);
    float* cf = kernel.ptr<float>();
    double* cd = kernel.ptr<double>();

    double sigmaX = sigma > 0 ? sigma : ((n-1)*0.5 - 1)*0.3 + 0.8;
    double scale2X = -0.5/(sigmaX*sigmaX);
    double sum = 0;

    int i;
    for( i = 0; i < n; i++ )
    {
        double x = i - (n-1)*0.5;
        double t = fixed_kernel ? (double)fixed_kernel[i] : std::exp(scale2X*x*x);
        if( ktype == CV_32F )
        {
            cf[i] = (float)t;
            sum += cf[i];
        }
        else
        {
            cd[i] = t;
            sum += cd[i];
        }
    }

    CV_DbgAssert(fabs(sum) > 0);
    sum = 1./sum;
    for( i = 0; i < n; i++ )                            //歸一化
    {
        if( ktype == CV_32F )
            cf[i] = (float)(cf[i]*sum);
        else
            cd[i] *= sum;
    }

    return kernel;
}

代碼實作矩陣：

#include<opencv2/opencv.hpp>
#include<iostream>

using namespace cv;
using namespace std;

int main()
{
	Mat g, g1, g2;
	g= getGaussianKernel(5, 0, CV_32F);    //size,sigma,type
	g1 = g * g.t();                        //g * g的轉置得到二維高斯卷積核
	cout << g1 << endl;
	g2 = g * g.t() * 256;                  //歸一化右側整數矩陣
	cout << g2 << endl;
	cin.get();
	return 0;
}

計算出來的高斯模板引數為：

當第二個引數sigma的取值為0時，getGaussianKernel函式中已經指定了1，3，5，7這四個模板的引數，為的是呼叫常用模板取固定整型引數，去除小數點也方便運算；
當第二個引數sigma的取值不為0時，getGaussianKernel函式將按照指定sigma值代入下兩行代碼運算：

double x = i - (n-1)*0.5;          //相當于式中 i-k=i-(n-1)/2
double t = fixed_kernel ? (double)fixed_kernel[i] : std::exp(scale2X*x*x);

上兩行代碼對應了上述二維高斯函式（不帶引數-1）的實作，其中濾去了系數 $1/2 \pi \sigma ^{2}$ ，是因為在下面代碼作歸一化處理時，可以消除該權重，因此可以省去加快計算速度，

獲取Matlab中的fspecial函式：

 filter=fspecial('gaussian',5,0.5);    //mood,size,sigma

*關于標準差 $\sigma$ (sigma)的取值和歸一化處理做個兩個注釋：

1. 標準差 $\sigma$ (sigma)

當標準差 $\sigma$ 取不同的值時，二維高斯函式的形狀會有很大的變化：如果標準差 $\sigma$ 選擇過小，偏離中心的所有像素的權重將會非常小，相當于加權和影響基本不考慮鄰域像素的作用，這樣濾波操作退化為影像的點運算，無法起到平滑噪聲的作用；相反如果標準差 $\sigma$ 選擇過大，而鄰域相對較小，這樣在鄰域內高斯模板將會退化為平均模板；因此在實際應用中選擇合適的標準差 $\sigma$ 非常重要，

2.歸一化

不難看出，在矩陣核的左邊存在一個系數，它是歸一化的象征 $1/\sum M(i,j)$ ;歸一化的目的：對灰度級為常數的影像區域，高斯模板的回應和必須為1，若小于1，像素值發生偏移，產生了誤差，鄰域像素之間的差值將減小；若大于1，存在超過像素上限（255）的可能，形成區域亮度，因此要對初始形成的模板進項歸一化處理，且也存在提高整體像素精度的作用，

2.1 影像卷積

我們依然采用延伸的思維從一維過渡到二維影像卷積，先來看下連續信號的卷積：

1.連續信號的卷積

對于任意波形的信號 $x(t)$ 都可以分割成許多相鄰的矩陣脈沖， $\Delta$ 代表了脈沖的寬度，對于 $t=n\Delta$ 時刻的矩形脈沖，其高度即 $x(t)$ 的值為 $x(n\Delta )$ ，

無窮多個矩形脈沖的疊加可用來近似原信號 $x(t)$ ，即

$x(t)\approx \sum_{n=-\infty }^{\infty }x(n\Delta )g_{\Delta }(t-n\Delta )\Delta$

顯然，當脈沖寬度越窄，近似程度就越高，就越逼近原信號 $x(t)$ （類似于高數中的經典積分思想）；當 $\Delta \rightarrow 0$ 極限的情況下，高度在上升，但面積始終保持為1，因此們函式可表示為由強度形式表達的單位沖激函式 $\delta (t)$ ,上式變換為

$x(t)=\lim_{\Delta \rightarrow 0}\sum_{n=-\infty }^{\infty }x(n\Delta )\delta (t-n\Delta )\Delta$

我們用 $d\tau$ 表示 $\Delta$ ， $\tau$ 表示 $n\Delta$ ，求和變成連續新變數 $\tau$ 的卷積積分

$x(t)=\int_{-\infty }^{\infty }x(\tau )\delta (t-\tau )d\tau$

表明任意波形的信號 $x(t)$ 可以表示為無限多個強度為 $x(\tau )d\tau$ 的單位沖激信號 $[x(\tau )d\tau]\delta (t-\tau)$ 的積分
表明任意波形的信號 $x(t)$ 都可以分解為連續的加權（延遲）單位沖激信號之和
對于連續信號而言，卷積是一種特殊的積分運算，它的整個程序就是一個函式固定不動，另一個函式先以y軸為對稱軸翻轉，然后不斷執行相乘，積分

2.離散信號的卷積

離散時間信號是連續時間信號經過離散化（即取樣）的結果，即連續卷積積分離散化為

$x(k)=\sum_{n=-\infty }^{\infty }x(n)\delta (k-n)=x(k)\ast h(k)$

表明任一離散信號 $x(k)$ 均可表示為單位函式 $\delta (k)$ 的延時加權和的形式

根據線性時不變系統的零狀態回應疊加性和時不變性，則離散系統對 $x(n)\delta (k-n)$ 零狀態回應為 $x(n)h(k-n)$ ,把得到的零狀態回應稱為卷積和或離散卷積，記為

$x(k)=\sum_{n=-\infty }^{\infty }x(n)h(k-n)=x(k)\ast h(k)$

對于上式離散信號求卷積和而言，它的整個程序就是一個一維序列點集合和另一個一維序列點集合反向（翻轉）對應乘積的和
在此基礎之上，拓展到二維（有限范圍k內），便有了二維離散卷積，即影像卷積 $g(x,y)=\sum_{m=-k}^{k}\sum_{n=-k}^{k}w(m,n)f(x-m,y-n)=w(x,y)\ast f(x,y)$

3.影像卷積

在執行線性空間濾波時，存在兩個近似的概念：一個是"相關"，另一個是"卷積"，

"相關"是濾波器模板移過影像并計算每個位置乘積之和的處理
"卷積"是模板先旋轉180°，再將濾波器模板移過影像并計算每個位置乘積之和的處理

"卷積"的基本特性是 $x(k)$ 與單位函式 $\delta (k)$ 的卷積和仍然是 $x(k)$ 本身，即

$x(k)=x(k)\ast \delta (k)$

基于以上這一點，我們延伸到二維影像中作卷積(mode:same)，我們令模板與同模板尺寸大小，除中心點像素點為1，其余點為0的矩陣進行"相關"運算，得到的結果是在中心點位產生該模板的一個旋轉180°的版本，因此，如果我們預先旋轉模板，并執行相同的滑動乘積求和的操作，就能得到希望的結果（中心點為模版矩陣），也契合公式的求取，但如果濾波器模板是對稱的，那么作"相關"和"卷積"運算將得到相同的結果，高斯模板正是如此，

我們也可以從另一個角度舉例說明，影像卷積的模板需要先旋轉180°；將二維卷積公式展開代入求取特定的值 $g(0,0)$ ，有

$g(0,0)=w(1,1)f(-1,-1)+w(1,0)f(-1,0)+w(1,-1)f(-1,1)+...+w(-1,-1)f(1,1)$

通過展開式和上圖明顯能看出要想實作滑動乘積求和，要先將模板w旋轉180°即可；對于超出原影像邊界的像素值默認賦0，

3.1 下采樣

下采樣用于減半計算得到的近似及上一層空間維數的影像，下采樣操作可視為洗掉偶數行和偶數列的像素點，賦給新的矩陣序列

根據OpenCV官方提供的代碼，pyrDown()函式專門用于影像的下采樣計算（包含了高斯模糊的卷積運算，模板引數大小默認為5×5的）：

pyrDown()函式原型
1.  void cv::pyrDown(InputArray src,                       //待下采樣的影像 
2.                   OutputArray dst,                      //輸出下采樣后的影像
3.                   const Size & dstsize = Size(),        //輸出影像尺寸（限制），默認是N/2
4.                   int borderType = BORDER_DEFAULT)      //像素邊界外推方式，默認即可

至此，我們反復的迭代計算（一般金字塔4~5層），便形成了高斯金字塔（具體圖片與拉普拉斯金字塔一并給出），

2. 拉普拉斯金字塔的構建

利用基數為2的上采樣（在偶數行和偶數列補0）作用在高斯N+1級影像上（尺寸與N級高斯影像一致）
對上采樣后影像進行高斯模糊（高斯模板核*4）
將模糊后的影像與原N級高斯影像作差值運算，得到第N級的拉普拉斯影像

2.1 上采樣

上采樣用于翻倍計算得到的近似同下一層空間維數的影像，下采樣操作可視為在偶數行和偶數列的像素值賦0（與下采樣形成互補操作），賦給新的矩陣序列

根據OpenCV官方提供的代碼，pyrUp()函式專門用于影像的上采樣計算（包含了高斯模糊的卷積運算，5×5的模板引數*4）：

pyrUp()函式原型
1.  void cv::pyrUp(InputArray src,                       //待上采樣的影像 
2.                 OutputArray dst,                      //輸出上采樣后的影像
3.                 const Size & dstsize = Size(),        //輸出影像尺寸（限制），默認是N*2
4.                 int borderType = BORDER_DEFAULT)      //像素邊界外推方式，默認即可

至此，我們反復的迭代計算，記得作減法運算，便形成了拉普拉斯金字塔，

*關于模板核*4和模板插值濾波器做個兩個注釋：

1.模板核*4

對于上采樣后需要模糊的高斯模板核*4，很多博主都沒詳細說明，我的理解是：符合歸一化，采用5×5的模板引數落在對應的像素點上，其中存在大量賦值為0的像素點（這些點的權重相當于不作用），無論是對應原矩陣中奇0偶數、偶0奇數、奇數偶0、偶數奇0，非0像素點對應的權重和一定滿足 $\sum M(i,j)=64$ ，所以將原有歸一化系數*4=1/64，滿足歸一化的作用，