图像处理基础（十四）插值算法 - 网站源码_资源分享

在采用 OpenCV 的 resize 表达式时，存有对数类别的ZWG，如 INTER_NEAREST（前段时间邻对数）、INTER_LINEAR（双非线性对数）、INTER_CUBIC（双四次对数）等。在此做个历史记录。

双非线性对数

以水准路径 $y\to$ \vec{y} 路径为例，做一场非线性对数，非线性满足用户上式

$f(y2)-f(y1)y2-y1=f(y)-f(y1)y-y1$ \frac{f(y_2) – f(y_1)}{y_2 – y_1} = \frac{f(y) – f(y_1)}{y – y_1}

则能获得最终目标点 $y$ y 的对数结论，

$f(y)=(y2-y)\cdotf(y1)+(y-y1)\cdotf(y2)$ f(y) = \color{red}{(y_2 – y)} \cdot f(y_1) + \color{red}{(y – y_1)} \cdot f(y_2)

一般满足用户 $y1\leqy\leqy2,y2=y1+1$ y_1 \leq y \leq y_2,\,\,\,y_2 = y_1 + 1 ，

上面的对数直观理解， $y$ y 离 $y2$ y_2 越近， $(y2-y)$ (y_2 – y) 越小， $f(y1)$ f(y_1) 的权重越小， $y$ y 和 $y1$ y_1 的 $f$ f的表达式值越不相近；而 $(y-y1)=(y-(y2-1))=1-(y2-y)$ (y – y_1) = (y – (y_2 – 1))= 1 – (y_2 – y) 越大， $f(y2)$ f(y_2) 的权重越大， $y$ y 和 $y2$ y_2 的 $f$ f 的表达式值越相近，合理。

同理，能获得竖直路径 $x\to$ \vec{x} 的非线性对数

$f(x2)-f(x1)x2-x1=f(x)-f(x1)x-x1$ \frac{f(x_2) – f(x_1)}{x_2 – x_1} = \frac{f(x) – f(x_1)}{x – x_1} ,

双非线性对数能看作是先在水准路径做两次非线性对数，再在竖直路径上对前面的两个对数结论做一场非线性对数，具体如下图，

上图中， $Q12$ Q_{12} 和 $Q22$ Q_{22} 做一场非线性对数获得 $R2$ R_2 ，然后 $Q11$ Q_{11} 和 $Q21$ Q_{21} 做一场非线性对数获得 $R1$ R_1 ，最后 $R2$ R_2 与 $R1$ R_1 做一场非线性对数获得最终目标坐标 $P$ P的值（ $P$ P 是小数的坐标，不是离散的）。

（注：先在竖直路径上做，再水准路径做，二者是等价的）。

下采样

先以下采样为例

假如要将 $5\times5$ 5\times 5 的图像 $I$ I 对数成 $3\times3$ 3 \times 3 的图像 $O$ O，以水准方向 $y\to$ \vec{y} 为例，按照水准的宽比 $53=1.67$ \frac{5}{3} = 1.67 ，求最终目标图像 $O$ O 中 $(0,1)$ (0, 1) 的点，对应输入图像 $I$ I 中坐标 $(0,1.67)$ (0, 1.67)，表示为 $P$ P，则找到 $P$ P 周围四个点—— $(0,1),(0,2),(1,1),(1,2)$ (0,1),(0,2), (1,1), (1,2) ，根据距离的大小做加权（之前的非线性对数），获得坐标 $P$ P 处的值，即

$f(P)=0.33\cdotf(0,1)+0.67\cdotf(0,2)+0\cdotf(1,1)+0\cdotf(1,2)$ f(P) = 0.33 \cdot f_{(0,1)} + 0.67 \cdot f_{(0,2)} + 0 \cdot f_{(1,1)} + 0 \cdot f_{(1,2)} 。

同理，如果求最终目标图像中点 $Q$ Q $(1,2)$ (1,2) 的值，对应原图 $5\times5$ 5\times 5 中的坐标为 $(1.67,3.33)$ (1.67, 3.33) ，找到周围四个点 $(1,3),(1,4),(2,3),(2,4)$ (1, 3),(1, 4), (2, 3), (2, 4) ，做四次非线性对数——先水准路径做两次非线性对数，分别获得

$f(R2)=0.67\cdotf(1,3)+0.33\cdotf(1,4)$ f(R_2) = 0.67\cdot f_{(1,3)} + 0.33 \cdot f_{(1,4)}

$f(R1)=0.67\cdotf(2,3)+0.33\cdotf(2,4)$ f(R_1) = 0.67\cdot f_{(2,3)} + 0.33 \cdot f_{(2,4)}

获得的两个值，再在竖直路径做一场非线性对数

$fQ=0.33\cdotf(R2)+0.67\cdotf(R1)$ f_{Q} = 0.33\cdot f(R_2) + 0.67 \cdot f(R_1)

即可。

居中

下采样演算法按照上面的演算法，能获得一个下采样结论，但对数时会忽略右下角的内容，还是上面 $5\times5$ 5\times 5的图像 $I$ I 对数成 $3\times3$ 3 \times 3 的图像 $O$ O，在水准 $y\to$ \vec{y} 路径上，求解 $I$ I 每一个点对应原图 $O$ O 的 $y$ y 坐标如下

$0,1.67,3.33$ 0,\,1.67,\,3.33

竖直 $x\to$ \vec{x}方向同理，坐标是 $0,1.67,3.33$ 0,\,1.67,\,3.33

但原图 $I$ I 的坐标范围是 $y\to$ \vec{y} ~ $[0,4]$ [0, 4] 和 $x\to$ \vec{x} ~ $[0,4]$ [0,4]

能发现，对数时只会利用到原图 $O$ O 左上角的像素，右下角的内容被忽略掉了。

如果换成更大图像的下采样，比如 $120\times120$ 120 \times 120 的图像下采样为 $3\times3$ 3\times 3 ，两个路径都只会利用到 $0,40,80$ 0, 40, 80 这些坐标的值，而忽略掉了 $81-119$ 81 -119 坐标的内容，这样的下采样有点不合格。

因此，后面提出将两幅图像的中心对齐，如下图

求解 $h\timesw$ h\times w的小图，因为坐标都是从 0 开始的，所以在水准 $y\to$ \vec{y} 路径上最远能到达 $(w-1)\cdotW0w$ (w – 1) \cdot \frac{W_0 }{w} ，则右侧有

$(W0-1)-((w-1)\cdotW0w)=W0w-1$ \begin{aligned} &(W_0 – 1) – \big((w – 1)\cdot \frac{W_0}{w}\big) = \frac{W_0}{w} – 1 \end{aligned}

空出来了，将上图黄区域往右偏移 $12$ \frac{1}{2} 个空出来的距离，即可尽量利用到中间的信息。

原来从小图映射到大图的坐标计算是

$y\timesW0w$ y \times \frac{W_0}{w} ，

现在往右偏移

$y\timesW0w+12\times(W0w-1)$ y \times \frac{W_0}{w} + \frac{1}{2} \times \big( \frac{W_0}{w} – 1 \big)

化简一下，

$(y+0.5)\timesW0w-0.5$ (y + 0.5) \times \frac{W_0}{w} – 0.5 ，

同理 $x$ x 路径上小图映射到大图的坐标计算就是

$(x+0.5)\timesH0h-0.5$ (x + 0.5) \times \frac{H_0}{h} – 0.5

上采样

下采样没什么问题，但上采样有所不同（以下只是我个人的看法，没看源码也看权威实现）。

举个例子， $3\times3$ 3 \times 3 的图像放大为 $5\times5$ 5 \times 5的图像，如果按照下采样的思路（先不考虑中心对齐），获得的一个路径上的坐标有

$0,0.6,1.2,1.8,2.4$ 0, 0.6, 1.2, 1.8, 2.4

一共五个坐标，但最后的坐标 $2.4$ 2.4 超出了原图 $3\times3$ 3 \times 3在一个路径的坐标范围—— $(0,0),(0,1),(0,2)$ (0, 0), (0, 1), (0, 2) ；

再假如是， $3\times3$ 3 \times 3 的图像放大为 $7\times7$ 7 \times 7的图像，一个路径的坐标就有

$0,0.429,0.857,1.285,1.714,2.143,2.571$ 0 \color{red}{ , }0.429\color{red}{ , }0.857\color{red}{ , }1.285\color{red}{ , } 1.714\color{red}{ , }2.143\color{red}{ , }2.571

有两个坐标 $2.143,2.571$ 2.143\color{red}{ , }2.571 超出了原图的坐标范围，也不可取。

因此，我修改一下坐标映射

$y\timesW0-1w$ y \times \frac{W_0 – 1}{w}

然后以上采样为 $7\times7$ 7 \times 7 为例， $y\times3-17$ y \times \frac{3 – 1}{7} ，一个路径上从 0, 1, 2, 3, 4, 5, 6的七个坐标有

$0,0.285,0.571,0.857,1.1428,1.428,1.714$ 0\color{red}{,}0.285\color{red}{,}0.571\color{red}{,}0.857\color{red}{,}1.1428\color{red}{,}1.428\color{red}{,}1.714

右下角也没利用到，再和之前一样做个中心对齐，

$y\timesW0-1w+12\times(W0-1w-1)=(y+0.5)\timesW0-1w-0.5$ y \times \frac{W_0 – 1}{w} + \frac{1}{2} \times \big( \frac{W_0 – 1}{w} – 1 \big) = (y + 0.5) \times \frac{W_0 – 1}{w} – 0.5

除了上面写的一个可能可行的方案，还有另一种办法，不做中心对齐，

直接坐标映射改成 $y\timesW0-1w-1$ y \times \frac{W_0 – 1}{w – 1} ，以 $3\times3$ 3 \times 3 的图像放大为 $7\times7$ 7 \times 7 的图像为例，拉伸比变成了

$37\to3-17-1=26=0.33$ \frac{3}{7} \rightarrow \frac{3 – 1}{7 – 1} = \frac{2}{6} = 0.33 ,

一个路径上从 0, 1, 2, 3, 4, 5, 6的七个坐标有

$0,0.33,0.67,1.0,1.33,1.67,2$ 0\color{red}{,}0.33\color{red}{,}0.67\color{red}{,}1.0\color{red}{,}1.33\color{red}{,}1.67\color{red}{,}2

恰好布满原图像 $3\times3$ 3 \times 3在一个路径上的区间范围！而且不需要做中心对齐。

上面都是两个路径都是同时放大，或者同时缩小，如果一个路径拉伸，另一个路径缩小，就分别采用上采样的映射和下采样的映射即可。

实验

341 x 512

下采样

140 x 200

上采样

1600 x 2400

代码

inline float _min(const float x, const float y) { return x > y ? y : x; } inline float _max(const float x, const float y) { return x < y ? y : x; } cv::Mat bilinear_interpolate(const cv::Mat& origin, const std::pair<int, int>& _size) { // 获取信息 const int H0 = origin.rows; const int W0 = origin.cols; const int C = origin.channels(); const int h = _size.first; const int w = _size.second; // 计算 x 路径和 y 路径上的比率 const float x_rate = H0 > h ? H0 * 1.f / h : (H0 – 1) * 1.f / h; const float y_rate = W0 > w ? W0 * 1.f / w : (W0 – 1) * 1.f / w; // 准备一个结论 cv::Mat result(h, w, CV_8UC3); uchar* const res_ptr = result.data; // 用来计算在结论 .data 存放的位置 int cnt = 0; // 一共要对数 h X w 次 for(int x = 0; x < h; ++x) { // 找到结论中 x, 对应原图中的 x 坐标 float x_pos = x_rate * (x + 0.5f) – 0.5f; // 找到这个 x_pos 的上下界 const int x_down = _max(std::floor(x_pos), 0); const int x_up = _min(x_down + 1, H0 – 1); // 计算 x 路径上的对数参数 float x_left = x_up – x_pos; float x_right = x_pos – x_down; // 1 – x_left // 原图中(x_down 和 x_up) 两行的指针 const uchar* const ori_ptr = origin.data + x_down * W0 * C; const uchar* const ori_ptr_2 = ori_ptr + W0 * C; // 填充结论的第 x 行的 y 个像素 for(int y = 0;y < w; ++y) { // 计算 y 对应原图中 y 的坐标, 放大或者缩小 float y_pos = y_rate * (y + 0.5f) – 0.5f; // 计算 y 的上下界, 此时映射到原图中 (x_pos, y_pos) 的周围四个点都找到了 const int y_down = std::floor(y_pos); const int y_up = _min(y_down + 1, W0 – 1); // 计算 y 路径上的对数参数 float y_left = y_up – y_pos; float y_right = y_pos – y_down; // 1 – y_left // 多个通道分别计算 for(int c = 0;c < C; ++c) { // y 路径上第一场非线性对数, 获得两个值 float f_E = y_left * ori_ptr[y_down * C + c] + y_right * ori_ptr[y_up * C + c]; float f_F = y_left * ori_ptr_2[y_down * C + c] + y_right * ori_ptr_2[y_up * C + c]; // x 路径上第二次非线性对数 float target = x_left * f_E + x_right * f_F; res_ptr[cnt++] = cv::saturate_cast<uchar>(target); } } } return result; }

BiCubic 对数

双立方对数，原理看得我有点混乱，我也是在网上胡乱搜，好多博客都是直接给一个定义好的局部 $4\times4$ 4 \times 4 加权表达式，比如

至于为什么，有论文，但我没去看。我选择了一种更简单的理解方式，源于这个博客，直接从一维的 cubic 出发，给定 4 个连续点，求解四次表达式 $f(x)=ax3+bx2+cx+d$ f(x) = ax^3 + bx^2 + cx + d 的四个参数 $a,b,c,d$ a,b,c,d ，就可估计出映射坐标下界从 [-1, 2] 内的任意点的值，原理如下：

原理

Cubic 对数

$f(0)=df(1)=a+b+c+df'(0)=f(1)-f(-1)2f'(1)=f(2)-f(0)2$ \begin{aligned} f(0) &= d \\ f(1) &= a + b + c + d \\ f(0) &= \frac{f(1) – f(-1)}{2} \\ f(1) &= \frac{f(2) – f(0)}{2} \end{aligned}

这里直接用 $f(x+1)-f(x-1)2$ \frac{f(x + 1) – f(x – 1)}{2} 近似获得 $f'(x)$ f(x)，当然也能采用更加精确的近似方法。上面四个方程 $f(-1),f(0),f(1),f(2)$ f(-1), f(0), f(1), f(2) 已知，求 $a,b,c,d$ a,b,c,d ，如下：

$a=-0.5*f(-1)+1.5*f(0)-1.5*f(1)+0.5*f(2)b=f(-1)-2.5*f(0)+2*f(1)-0.5*f(2)c=-0.5*f(-1)+0.5*f(1)d=f(0)$ \begin{aligned} a &= -0.5 * f(-1) + 1.5 * f(0) – 1.5 * f(1) + 0.5 * f(2) \\ b &=f(-1) – 2.5 * f(0) + 2 * f(1) – 0.5 * f(2) \\ c &= -0.5 * f(-1) + 0.5 * f(1) \\ d &= f(0) \end{aligned}

举个例子，

假如映射坐标是 $3.2$ 3.2 ，获得下界 $floor(3.2)=3$ floor(3.2) = 3 ，从坐标 $3$ 3 开始，偏移量 $[-1,2]$ [-1,2] 之内的对数表达式能用上面的 $f(x)$ f(x) 来近似，只要知道 $f(-1),f(0),f(1),f(2)$ f(-1), f(0), f(1), f(2) 四个偏移量上的灰度值，就能估算出偏移量 $[-1,2]$ [-1,2] 内每一个点的值，例如 $f(3.2)$ f(3.2)

拓展到二维，就是 bicubic，能先水准路径 cubic 对数获得四个值，获得的四个值在竖直路径上再做一场 cubic 对数获得映射坐标的灰度值，如下图

实验

原图 314 x 512

先双立方对数缩小图像

双立方对数 50 x 75

再对数放大图像

我写的双非线性对数 600 * 900

我写的双立方对数

OpenCV 的双立方对数

双立方对数，我暂时没看出来哪里更好，可能是我写错了，但是我采用 OpenCV 内置的 INTER_CUBIC 也获得了不是很好的效果，看样子这种情况下还是需要超分辨适合一点。

代码

C++, 双立方对数即使开了 O2 也很慢，毕竟是 16 个点参与对数计算。

实现细节上，和之前的双非线性对数一样，映射之后的坐标偏左上方，需要求右边、下面空余的部分，然后每次映射的坐标都分别加上0.5 倍的空余部分，即可中心化。

因为每次对数，都需要周围 4 \times 4个点，所以我对参考图像做了长度 1 的 padding。

inline float cubic(const float x) { return x * x * x; } inline float square(const float x) { return x * x; } template<typename T> float make_cubic_interpolation(const std::vector<T>& F, const float input) { // 根据这 4 个点计算四次表达式的参数 float a = –0.5 * F[0] + 1.5 * F[1] – 1.5 * F[2] + 0.5 * F[3]; float b = F[0] – 2.5 * F[1] + 2 * F[2] – 0.5 * F[3]; float c = –0.5 * F[0] + 0.5 * F[2]; float d = F[1]; // 根据四次表达式, 对数算这个点 input 的值 return a * cubic(input) + b * square(input) + c * input + d; } cv::Mat bicubic_interpolate(const cv::Mat& origin, const std::pair<int, int>& _size) { int H = origin.rows; int W = origin.cols; const int H2 = _size.first; const int W2 = _size.second; const int C = origin.channels(); // 计算纵向跟横向的缩放比 const float h_ratio = H2 > H ? (H – 1) * 1.f / H2 : H * 1.f / H2; const float w_ratio = W2 > W ? (W – 1) * 1.f / W2 : W * 1.f / W2; // 计算纵向跟横向需要偏移的距离 const float h_add = 0.5 * ((H – 1) – (H2 – 1) * h_ratio); const float w_add = 0.5 * ((W – 1) – (W2 – 1) * w_ratio); // 做 padding, 因为是周围 16 个点做对数 const int pad = 1; const auto padded_image = make_pad(origin, pad, pad); const uchar* const pad_ptr = padded_image.ptr<uchar>(); // 准备一个结论 cv::Mat result(H2, W2, origin.type()); uchar* const res_ptr = result.ptr<uchar>(); int cnt = 0; // 准备几个临时变量 std::vector<uchar> temp_Y(4); // 存储横向一场对数的结论 std::vector<float> temp_X(4); // 存储纵向一场对数的四个点 // 对数每一个行 for(int x = 0;x < H2; ++x) { // 算这一行在 3 x 3 中的位置, 下界和偏移 float x_pos = x * h_ratio + h_add; int x_down = std::floor(x_pos); const float x_offset = x_pos – x_down; // 对数每一个点 for(int y = 0;y < W2; ++y) { float y_pos = y * w_ratio + w_add; int y_down = std::floor(y_pos); const float y_offset = y_pos – y_down; // 多通道 for(int ch = 0;ch < 3; ++ch) { // 首先, 计算从第 x_down 行开始, [-1, 0, 1, 2] 的对数 for(int i = –1; i <= 2; ++i) { // x_down + i 是在 3 x 3 图像中的坐标, pad 是做了 padding 的偏移量 const int X = x_down + i + pad; // 找到 x_down + i 行的数据起始的指针, + pad * C 是因为有横向 pad const uchar* const X_ptr = pad_ptr + X * padded_image.cols * C + pad * C; // 对数 (x_down + i, y_pos), 需要找到 (x_down + i, y_pos) 的四个点, 存储在 temp_Y 中 for(int j = –1; j <= 2; ++j) temp_Y[j + 1] = X_ptr[(y_down + j) * C + ch]; // 这 4 个点做 cubic 对数, 作为 (x_down + i, y_pos) 的结论 temp_X[i + 1] = make_cubic_interpolation<uchar>(temp_Y, y_offset); } // 获得了 (x_down + i, y_pos) 四个点的对数结论, 做一场 cubic 对数, 作为 (x_pos, y_pos) 的对数结论 const float one = make_cubic_interpolation<float>(temp_X, x_offset); res_ptr[cnt++] = cv::saturate_cast<uchar>(one); } } } return result; } int main() { // 读取图像 const std::string image_path(“./images/input/a1016-050716_115658__I2E4159.png”); cv::Mat origin_image = cv::imread(image_path); assert(not origin_image.empty()); // 先把图变小 const auto small = bicubic_interpolate(origin_image, {50, 75}); // 双立方对数, 把图变大 const auto big = bicubic_interpolate(small, {600, 900}); // 和双非线性对数对比 const auto bilinear_big = bilinear_interpolate(small, {600, 900}); // 和 OpenCV 内置实现对比 cv::Mat cv_big; cv::resize(small, cv_big, {900, 600}, cv::INTER_CUBIC); // 展示 cv_show(small); cv::Mat concat; cv::vconcat(std::vector<cv::Mat>({big, cv_big}), concat); cv_show(concat); // 保存结论 const std::string output_path(“./images/output/”); cv_write(small, output_path + “small.png”); cv_write(big, output_path + “big.png”); return 0; }