【影像縮放】雙立方（三次）卷積插值

dailc發表於2019-03-04

卷積

前言

影像處理中有三種常用的插值演算法：

最鄰近插值
雙線性插值
雙立方（三次卷積）插值

其中效果最好的是雙立方（三次卷積）插值，本文介紹它的原理以及使用

如果想先看效果和原始碼，可以拉到最底部

本文的契機是某次基於canvas做影像處理時，發現canvas自帶的縮放功能不盡人意，於是重溫了下幾種影像插值演算法，並整理出來。

為何要進行雙立方插值

對影像進行插值的目的是為了獲取縮小或放大後的圖片
常用的插值演算法中，雙立方插值效果最好
本文中介紹雙立方插值的一些數學理論以及實現

雙立方和三次卷積只是這個插值演算法的兩種不同叫法而已，可以自行推導，會發現最終可以將求值轉化為卷積公式

另外，像Photoshop等影像處理軟體中也有這三種演算法的實現

數學理論

雙立方插值計算涉及到16個畫素點，如下圖

簡單分析如下：

其中P00代表目標插值圖中的某畫素點(x, y)在原圖中最接近的對映點
- 譬如對映到原圖中的座標為(1.1, 1.1)，那麼P00就是(1, 1)
而最終插值後的影像中的(x, y)處的值即為以上16個畫素點的權重卷積之和

下圖進一步分析

如下是對圖的一些簡單分析

譬如計算插值圖中(distI, distJ)處畫素的值
首先計算它對映到原圖中的座標(i + v, j + u)
也就是說，卷積計算時，p00點對應(i, j)座標
最終，插值後的圖中(distI, distJ)座標點對應的值是原圖中(i, j)處鄰近16個畫素點的權重卷積之和
- i, j的範圍是[i - 1, i + 2]，[j - 1, j + 2]

卷積公式

設取樣公式為S(x)
原圖中每一個(i, j)座標點的值得表示式為f(i, j)
插值後對應座標的值為F(i + v, j + u)（這個值會作為(distI, distJ)座標點的值）

那麼公式為：

等價於(可自行推導)

提示

一定要區分本文中v, u和row, col的對應關係，v代表行數偏差，u代表列數偏差（如果混淆了，會造成最終的影像偏差很大）

如何理解卷積？

這是大學數學內容，推薦看看這個答案如何通俗易懂的解釋卷積-知乎

取樣公式

在卷積公式中有一個S(x)，它就是關鍵的卷積插值公式

不同的公式，插值效果會有所差異（會導致加權值不一樣）

本文中採用WIKI-Bicubic interpolation中給出的插值公式：

公式中的特點是：

S(0) = 1
S(n) = 0(當n為整數時)
當x超出範圍時，S(x)為0
當a取不同值時可以用來逼近不同的樣條函式（常用值-0.5, -0.75）

當a取值為-1

公式如下：

此時，逼近的函式是y = sin(x*PI)/(x*PI)，如圖

當a取值為-0.5

公式如下：

此時對應三次Hermite樣條

不同a的簡單對比

推導

可參考：

關於網上的一些推導公式奇怪實現

在網上查詢了不少相關資料，發現有不少文章中都用到了以下這個奇怪的公式（譬如百度搜尋雙立方插值）

一般這些文章中都聲稱這個公式是用來近似y = sin(x*PI)/(x)

但事實上，進過驗證，它與y = sin(x*PI)/(x)相差甚遠（如上圖中是將sin函式縮放到合理係數後比對）

由於類似的文章較多，年代都比較久遠，無從得知最初的來源

可能是某文中漏掉了分母的PI，亦或是這個公式只是某文自己實現的一個取樣公式，與sin無關，然後被誤傳了。

這裡都無從考據，僅此記錄，避免疑惑。

另一種基於係數的實現

可以參考：影像處理（一）bicubic解釋推導

像這類的實現就是直接計算最原始的係數，然後通過16個畫素點計算不同係數值，最終計算出目標畫素

本質是一樣的，只不過是沒有基於最終的卷積方程計算而已（也就是說在原始理論階段沒有推成插值公式，而是直接解出係數並計算）。

程式碼實現在github專案中可看到，參考最後的開源專案

程式碼實現

以下是JavaScript程式碼實現的插值核心方程

/**
 * 取樣公式的常數A取值,調整銳化與模糊
 * -0.5 三次Hermite樣條
 * -0.75 常用值之一
 * -1 逼近y = sin(x*PI)/(x*PI)
 * -2 常用值之一
 */
const A = -0.5;

function interpolationCalculate(x) {
    const absX = x > 0 ? x : -x;
    const x2 = x * x;
    const x3 = absX * x2;

    if (absX <= 1) {
        return 1 - (A + 3) * x2 + (A + 2) * x3;
    } else if (absX <= 2) {
        return -4 * A + 8 * A * absX - 5 * A * x2 + A * x3;
    }
    return 0;
}複製程式碼

以上是卷積方程的核心實現。下面則是一套完整的實現

/**
 * 取樣公式的常數A取值,調整銳化與模糊
 * -0.5 三次Hermite樣條
 * -0.75 常用值之一
 * -1 逼近y = sin(x*PI)/(x*PI)
 * -2 常用值之一
 */
const A = -1;

function interpolationCalculate(x) {
    const absX = x >= 0 ? x : -x;
    const x2 = x * x;
    const x3 = absX * x2;

    if (absX <= 1) {
        return 1 - (A + 3) * x2 + (A + 2) * x3;
    } else if (absX <= 2) {
        return -4 * A + 8 * A * absX - 5 * A * x2 + A * x3;
    }

    return 0;
}

function getPixelValue(pixelValue) {
    let newPixelValue = pixelValue;

    newPixelValue = Math.min(255, newPixelValue);
    newPixelValue = Math.max(0, newPixelValue);

    return newPixelValue;
}

/**
 * 獲取某行某列的畫素對於的rgba值
 * @param {Object} data 影像資料
 * @param {Number} srcWidth 寬度
 * @param {Number} srcHeight 高度
 * @param {Number} row 目標畫素的行
 * @param {Number} col 目標畫素的列
 */
function getRGBAValue(data, srcWidth, srcHeight, row, col) {
    let newRow = row;
    let newCol = col;

    if (newRow >= srcHeight) {
        newRow = srcHeight - 1;
    } else if (newRow < 0) {
        newRow = 0;
    }

    if (newCol >= srcWidth) {
        newCol = srcWidth - 1;
    } else if (newCol < 0) {
        newCol = 0;
    }

    let newIndex = (newRow * srcWidth) + newCol;

    newIndex *= 4;

    return [
        data[newIndex + 0],
        data[newIndex + 1],
        data[newIndex + 2],
        data[newIndex + 3],
    ];
}

function scale(data, width, height, newData, newWidth, newHeight) {
    const dstData = newData;

    // 計算壓縮後的縮放比
    const scaleW = newWidth / width;
    const scaleH = newHeight / height;

    const filter = (dstCol, dstRow) => {
        // 源影像中的座標（可能是一個浮點）
        const srcCol = Math.min(width - 1, dstCol / scaleW);
        const srcRow = Math.min(height - 1, dstRow / scaleH);
        const intCol = Math.floor(srcCol);
        const intRow = Math.floor(srcRow);
        // 計算u和v
        const u = srcCol - intCol;
        const v = srcRow - intRow;

        // 真實的index，因為陣列是一維的
        let dstI = (dstRow * newWidth) + dstCol;

        dstI *= 4;

        // 儲存灰度值的權重卷積和
        const rgbaData = [0, 0, 0, 0];
        // 根據數學推導，16個點的f1*f2加起來是趨近於1的（可能會有浮點誤差）
        // 因此就不再單獨先加權值，再除了
        // 16個鄰近點
        for (let m = -1; m <= 2; m += 1) {
            for (let n = -1; n <= 2; n += 1) {
                const rgba = getRGBAValue(
                    data,
                    width,
                    height,
                    intRow + m,
                    intCol + n,
                );
                // 一定要正確區分 m,n和u,v對應的關係，否則會造成影像嚴重偏差（譬如出現噪點等）
                // F(row + m, col + n)S(m - v)S(n - u)
                const f1 = interpolationCalculate(m - v);
                const f2 = interpolationCalculate(n - u);
                const weight = f1 * f2;

                rgbaData[0] += rgba[0] * weight;
                rgbaData[1] += rgba[1] * weight;
                rgbaData[2] += rgba[2] * weight;
                rgbaData[3] += rgba[3] * weight;
            }
        }

        dstData[dstI + 0] = getPixelValue(rgbaData[0]);
        dstData[dstI + 1] = getPixelValue(rgbaData[1]);
        dstData[dstI + 2] = getPixelValue(rgbaData[2]);
        dstData[dstI + 3] = getPixelValue(rgbaData[3]);
    };

    // 區塊
    for (let col = 0; col < newWidth; col += 1) {
        for (let row = 0; row < newHeight; row += 1) {
            filter(col, row);
        }
    }
}

export default function bicubicInterpolation(imgData, newImgData) {
    scale(imgData.data,
        imgData.width,
        imgData.height,
        newImgData.data,
        newImgData.width,
        newImgData.height);

    return newImgData;
}複製程式碼