很長時間不敲代碼，感覺一閑下來就忘了很多。想着把一些圖像算法自己實作一遍也好，一方面加深算法的學習和了解，另一方面又可以練練編碼能力。對于我這個非科班出身的，也挺有好處的。
不管那麼多，先把這個坑挖了。閑着想敲代碼了就慢慢填。
下面先把圖像進行中常見的三種插值算法實作了。先解釋一下什麼叫插值。老樣子先看看維基百科怎麼解釋。
數學的數值分析領域中，内插或稱插值(英語：interpolation)是一種通過已知的、離散的資料點，在範圍内推求新資料點的過程或方法。
450px-Splined_epitrochoid.svg.png
一組離散資料點在一個外延的插值。曲線中實際已知資料點是紅色的；連接配接它們的藍色曲線即為插值。
再舉個例子：
x1 = 1, y1 = 3
x2 = 3, y2 = 7
x3 = 5, y3 = 24
求x = 4時， y = ？
是以可以看到，可以這麼了解插值，在一個函數裡面，自變量是離散有間隔的，插值就是往自變量的間隔之間插入新的自變量，然後求解新的自變量函數值。這有什麼作用呢，從上圖可以看到，散點圖是可以利用插值來拟合曲線的，藍色的就是插入的密密麻麻的點。
然而在圖像處理上的應用可以展現在圖像的縮放上面。如放大一張圖檔，在像素點的層面上其實就是往像素點之間插入新的像素點進而增大圖像的分辨率。
插值算法有很多種，具體可以參考維基百科插值
這裡就實作圖像處理比較常用的三種：最近鄰域內插補點、雙線性插值、雙三次插值
最近鄰域插值
最近鄰域内插法(Nearest Neighbor interpolation)
最近鄰域是三種插值之中最簡單的一種，原理就是選取距離插入的像素點(x+u, y+v)【注：x,y為整數， u，v為小數】最近的一個像素點，用它的像素點的灰階值代替插入的像素點。說到距離，順便複習一下像素點之間的三種空間距離。
對于像素p、q和z，分别具有坐标(x,y)，(s,t)和(u,v)，如果
(1) D(p,q) ≥ 0 (當且僅當p=q時，D(p,q)=0)
(2) D(p,q) = D(q,p)
(3) D(p,z) ≤ D(p,q) + D(q,z)
則稱D是距離函數或度量
歐氏距離.gif
D4距離(城市距離)
D4距離.gif
D8距離(棋盤距離)
D8距離.gif
這裡代碼實作采用的是歐式距離
//最近鄰域插值
//根據目标圖像的像素點(浮點坐标)找到原始圖像中的4個像素點，取距離該像素點最近的一個原始像素值作為該點的值。
#include
#include
namespace mycv {
void nearestIntertoplation(cv::Mat& src, cv::Mat& dst, const int rows, const int cols);
}//mycv
double distance(const double x1, const double y1, const double x2, const double y2);//兩點之間距離，這裡用歐式距離
int main()
{
cv::Mat img = cv::imread("lena.jpg", 0);
if (img.empty()) return -1;
cv::Mat dst;
mycv::nearestIntertoplation(img, dst, 600, 600);
cv::imshow("img", img);
cv::imshow("dst", dst);
cv::waitKey(0);
return 0;
return 0;
}//main
void mycv::nearestIntertoplation(cv::Mat& src, cv::Mat& dst, const int rows, const int cols)
{
//比例尺
const double scale_row = static_cast(src.rows) / rows;
const double scale_col = static_cast(src.rows) / cols;
//擴充src到dst
dst = cv::Mat(rows, cols, src.type());
assert(src.channels() == 1 && dst.channels() == 1);
for (int i = 0; i < rows; ++i)//dst的行
for (int j = 0; j < cols; ++j)//dst的列
{
//求插值的四個點
double y = (i + 0.5) * scale_row + 0.5;
double x = (j + 0.5) * scale_col + 0.5;
int x1 = static_cast(x);//col對應x
if (x1 >= (src.cols - 2)) x1 = src.cols - 2;//防止越界
int x2 = x1 + 1;
int y1 = static_cast(y);//row對應y
if (y1 >= (src.rows - 2)) y1 = src.rows - 2;
int y2 = y1 + 1;
//根據目标圖像的像素點(浮點坐标)找到原始圖像中的4個像素點，取距離該像素點最近的一個原始像素值作為該點的值。
assert(0 < x2 && x2 < src.cols && 0 < y2 && y2 < src.rows);
std::vector dist(4);
dist[0] = distance(x, y, x1, y1);
dist[1] = distance(x, y, x2, y1);
dist[2] = distance(x, y, x1, y2);
dist[3] = distance(x, y, x2, y2);
int min_val = dist[0];
int min_index = 0;
for (int i = 1; i < dist.size(); ++i)
if (min_val > dist[i])
{
min_val = dist[i];
min_index = i;
}
switch (min_index)
{
case 0:
dst.at(i, j) = src.at(y1, x1);
break;
case 1:
dst.at(i, j) = src.at(y1, x2);
break;
case 2:
dst.at(i, j) = src.at(y2, x1);
break;
case 3:
dst.at(i, j) = src.at(y2, x2);
break;
default:
assert(false);
}
}
}
double distance(const double x1, const double y1, const double x2, const double y2)//兩點之間距離，這裡用歐式距離
{
return (x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2);//隻需比較大小，傳回距離平方即可
}
雙線性插值
雙線性插值，顧名思義，在像素點矩陣上面，x和y兩個方向的線性插值所得的結果。那麼先看看線性插值
Linear_interpolation.png
image.png
image.png
以上是一維的，接下來看看二維中的雙線性插值
Bilinear_interpolation.png
首先在x方向上面線性插值，得到R2、R1
image.png
然後以R2,R1在y方向上面再次線性插值
image.png
如果選擇一個坐标系統使得 f 的四個已知點坐标分别為 (0, 0)、(0, 1)、(1, 0) 和 (1, 1)，那麼插值公式就可以化簡為
image.png
用矩陣表示
image.png
下面貼上代碼實作
//線性内插(雙線性插值)
#include
#include
#include
namespace mycv {
void bilinearIntertpolatioin(cv::Mat& src, cv::Mat& dst, const int rows, const int cols);
}//mycv
int main()
{
cv::Mat img = cv::imread("lena.jpg", 0);
if (img.empty()) return -1;
cv::Mat dst;
mycv::bilinearIntertpolatioin(img, dst, 600, 600);
cv::imshow("img", img);
cv::imshow("dst", dst);
cv::waitKey(0);
return 0;
}//main
void mycv::bilinearIntertpolatioin(cv::Mat& src, cv::Mat& dst, const int rows, const int cols)
{
//比例尺
const double scale_row = static_cast(src.rows) / rows;
const double scale_col = static_cast(src.rows) / cols;
//擴充src到dst
dst = cv::Mat(rows, cols, src.type());
assert(src.channels() == 1 && dst.channels() == 1);
for(int i = 0; i < rows; ++i)//dst的行
for (int j = 0; j < cols; ++j)//dst的列
{
//求插值的四個點
double y = (i + 0.5) * scale_row + 0.5;
double x = (j + 0.5) * scale_col + 0.5;
int x1 = static_cast(x);//col對應x
if (x1 >= (src.cols - 2)) x1 = src.cols - 2;//防止越界
int x2 = x1 + 1;
int y1 = static_cast(y);//row對應y
if (y1 >= (src.rows - 2)) y1 = src.rows - 2;
int y2 = y1 + 1;
assert(0 < x2 && x2 < src.cols && 0 < y2 && y2 < src.rows);
//插值公式,參考維基百科矩陣相乘的公式https://zh.wikipedia.org/wiki/%E5%8F%8C%E7%BA%BF%E6%80%A7%E6%8F%92%E5%80%BC
cv::Matx12d matx = { x2 - x, x - x1 };
cv::Matx22d matf = { static_cast(src.at(y1, x1)), static_cast(src.at(y2, x1)),
static_cast(src.at(y1, x2)), static_cast(src.at(y2, x2)) };
cv::Matx21d maty = {
y2 - y,
y - y1
};
auto val = (matx * matf * maty);
dst.at(i, j) = val(0,0);
}
}
雙三次插值
雙三次插值就稍有點複雜了，先看維基百科
在數值分析這個數學分支中，雙三次插值(英語：Bicubic interpolation)是二維空間中最常用的插值方法。在這種方法中，函數 f 在點 (x, y) 的值可以通過矩形網格中最近的十六個采樣點的權重平均得到，在這裡需要使用兩個多項式插值三次函數，每個方向使用一個
雙三次插值計算公式
image.png
即
image.png
那麼這個a(i, j)便是介紹裡面所說的權重系數了，是以關鍵是要把它求解出來。
求解權重系數的公式如下
image.png
求解加強系數代碼如下
std::vector mycv::getW(double coor, double a)
{
std::vector w(4);
int base = static_cast(coor);//取整作為基準
double e = coor - static_cast(base);//多出基準的小數部分
std::vector tmp(4);//存放公式中 |x| <= 1, 1 < |x| < 2四個值
// 4 x 4的16個點，是以tmp[0]和tmp[4]距離較遠，值在[1, 2]區間
tmp[0] = 1.0 + e;// 1 < x < 2
tmp[1] = e;//x <= 1
tmp[2] = 1.0 - e;// x <= 1
tmp[3] = 2.0 - e;// 1 < x < 2
//按照bicubic公式計算系數w
// x <= 1
w[1] = (a + 2.0) * std::abs(std::pow(tmp[1], 3)) - (a + 3.0) * std::abs(std::pow(tmp[1], 2)) + 1;
w[2] = (a + 2.0) * std::abs(std::pow(tmp[2], 3)) - (a + 3.0) * std::abs(std::pow(tmp[2], 2)) + 1;
// 1 < x < 2
w[0] = a * std::abs(std::pow(tmp[0], 3)) - 5.0 * a * std::abs(std::pow(tmp[0], 2)) + 8.0*a*std::abs(tmp[0]) - 4.0*a;
w[3] = a * std::abs(std::pow(tmp[3], 3)) - 5.0 * a * std::abs(std::pow(tmp[3], 2)) + 8.0*a*std::abs(tmp[3]) - 4.0*a;
return w;
}
求解出來之後，wx和wy就是4x4組成的16個系數，與插值所需要用到4x4的16個點相比對。插值步驟的關鍵代碼
//4x4數量的點(rr, cc) -> (y, x)
std::vector<:vector> > src_arr = {
{ src.at(rr - 1, cc - 1), src.at(rr, cc - 1), src.at(rr + 1, cc - 1), src.at(rr + 2, cc - 1)},
{ src.at(rr - 1, cc), src.at(rr, cc), src.at(rr + 1, cc), src.at(rr + 2, cc)},
{ src.at(rr - 1, cc + 1), src.at(rr, cc + 1), src.at(rr + 1, cc + 1), src.at(rr + 2, cc + 1)},
{ src.at(rr - 1, cc + 2), src.at(rr, cc + 2), src.at(rr + 1, cc + 2), src.at(rr + 2, cc + 2)}
};
for(int p = 0; p < 3; ++p)
for (int q = 0; q < 3; ++q)
{
//val(p, q) = w(p,q) * src(p, q)
val += wr[p] * wc[q] * static_cast(src_arr[p][q]);
}
assert(i < dst.rows && j < dst.cols);
dst.at(i, j) = static_cast(val);
下面放出完整代碼
//雙三次插值
#include
#include
#include
#include
namespace mycv {
void bicubicInsterpolation(cv::Mat& src, cv::Mat& dst, const int rows, const int cols);
std::vector getW(double coor, double a = -0.5);//a預設-0.5
}//mycv
int main(void)
{
cv::Mat img = cv::imread("lena.jpg", 0);
if (img.empty()) return -1;
cv::Mat dst;
mycv::bicubicInsterpolation(img, dst, 600, 600);
cv::imshow("img", img);
cv::imshow("dst", dst);
cv::waitKey(0);
return 0;
}//main
void mycv::bicubicInsterpolation(cv::Mat& src, cv::Mat& dst, const int rows, const int cols)
{
dst = cv::Mat(rows, cols, src.type(), cv::Scalar::all(0));//初始化dst
//比例尺
double row_scale = static_cast(src.rows) / rows;
double col_scale = static_cast(src.cols) / cols;
switch (src.channels())
{
case 1://灰階
for(int i = 2; i < dst.rows - 2; ++i)
for (int j = 2; j < dst.cols - 2; ++j)
{
//計算系數w
double r = static_cast(i * row_scale);
double c = static_cast(j * col_scale);
//防止越界
if (r < 1.0) r += 1.0;
if (c < 1.0) c += 1.0;
std::vector wr = mycv::getW( r);
std::vector wc = mycv::getW( c);
//最後計算插值得到的灰階值
double val = 0;
int cc = static_cast(c);
int rr = static_cast(r);
//防止越界
if (cc > src.cols - 3)
{
cc = src.cols - 3;
}
if (rr > src.rows - 3) rr = src.rows - 3;
assert(0 <= (rr - 1) && 0 <= (cc - 1) && (rr + 2) < src.rows && (cc + 2) < src.cols);
//4x4數量的點(rr, cc) -> (y, x)
std::vector<:vector> > src_arr = {
{ src.at(rr - 1, cc - 1), src.at(rr, cc - 1), src.at(rr + 1, cc - 1), src.at(rr + 2, cc - 1)},
{ src.at(rr - 1, cc), src.at(rr, cc), src.at(rr + 1, cc), src.at(rr + 2, cc)},
{ src.at(rr - 1, cc + 1), src.at(rr, cc + 1), src.at(rr + 1, cc + 1), src.at(rr + 2, cc + 1)},
{ src.at(rr - 1, cc + 2), src.at(rr, cc + 2), src.at(rr + 1, cc + 2), src.at(rr + 2, cc + 2)}
};
for(int p = 0; p < 3; ++p)
for (int q = 0; q < 3; ++q)
{
//val(p, q) = w(p,q) * src(p, q)
val += wr[p] * wc[q] * static_cast(src_arr[p][q]);
}
assert(i < dst.rows && j < dst.cols);
dst.at(i, j) = static_cast(val);
}
break;
case 3://彩色(原理一樣多了兩個通道而已)
break;
default:
break;
}
}
std::vector mycv::getW(double coor, double a)
{
std::vector w(4);
int base = static_cast(coor);//取整作為基準
double e = coor - static_cast(base);//多出基準的小數部分
std::vector tmp(4);//存放公式中 |x| <= 1, 1 < |x| < 2四個值
// 4 x 4的16個點，是以tmp[0]和tmp[4]距離較遠，值在[1, 2]區間
tmp[0] = 1.0 + e;// 1 < x < 2
tmp[1] = e;//x <= 1
tmp[2] = 1.0 - e;// x <= 1
tmp[3] = 2.0 - e;// 1 < x < 2
//按照bicubic公式計算系數w
// x <= 1
w[1] = (a + 2.0) * std::abs(std::pow(tmp[1], 3)) - (a + 3.0) * std::abs(std::pow(tmp[1], 2)) + 1;
w[2] = (a + 2.0) * std::abs(std::pow(tmp[2], 3)) - (a + 3.0) * std::abs(std::pow(tmp[2], 2)) + 1;
// 1 < x < 2
w[0] = a * std::abs(std::pow(tmp[0], 3)) - 5.0 * a * std::abs(std::pow(tmp[0], 2)) + 8.0*a*std::abs(tmp[0]) - 4.0*a;
w[3] = a * std::abs(std::pow(tmp[3], 3)) - 5.0 * a * std::abs(std::pow(tmp[3], 2)) + 8.0*a*std::abs(tmp[3]) - 4.0*a;
return w;
}
三次插值的代碼計算量特别大，而且代碼沒什麼優化，速度是慢的可以。但感覺可讀性還是ok得，便于了解吧。