運動模糊的製造&去除（c++&opencv）
本文是看了以下文章，有感而發，特地做一些補充記錄。
Wang Hawk：40。 如何消除攝影中的運動模糊？
如果消除攝影中的運動模糊？ 在影象質量評價中，物體是否清晰是最重要的指標之一。而受限拍攝條件和場景，成片常常會存在各種各樣原因引起的模糊。這篇文章討論的是，
相機固定
，由
運動物體
引起的模糊，我們應該怎麼去除。最常見的場景就是
監控攝像頭
啦，拍攝街道，路口，怎麼把運動車輛的車牌拍清楚呢？
1。 如何用程式碼模擬運動模糊
首先理解一下拍攝運動物體是如何引起模糊的呢？這張圖詮釋得很好：
在曝光時間內，不同位移的影象記錄在感測器上，相當於是一張靜止物體的影象卷積上矩形波，可以把這個矩形波看成是一維的卷積核函式：
![[1,1,1,1,1,1]](https://img.heatask.com/upload/website_attach/vtkQ3SODNTPQ2R1wnvKE0vbF5DhkPSbZ=+yQ13C-bEok3WgX_b=rsX9PjYm+7JOCje3mtMvWOY8eK+kX_boU+.jpeg)
，最終的成片相當於是這個靜止影象卷積上這個核函式。在空域上我們這麼操作，上程式碼：
//實現單方向卷積：空域卷積
Mat img = imread（“D：\\data\\lenna。jpg”， 0）；
Mat outimg = img。clone（）；
int h = img。rows；
int w = img。cols；
int wk = 24； //卷積核長度
Mat borderimg；
copyMakeBorder（img， borderimg， 0， 0， wk， 0， BORDER_CONSTANT， Scalar：：all（0））；
auto src = borderimg。data + wk；
auto dst = outimg。data；
for（int i = 0； i < h； i++）
{
for（int j = 0； j < w； j++）
{
int sum = 0；
for（int k = 0； k < wk； k++）
{
sum += src［j - k］；//卷積核加權，sum型別的範圍不要溢位
}
dst［j］ = sum / （wk）；
}
src += w + wk；
dst += w；
}
imwrite（“outimg。jpg”， outimg）；
在這裡我們要注意的是：
1。 理論上曝光時間是連續的，但我們在計算機的實現中已經把卷積核進行了離散化這麼一個過程。
2。 上面的kernel模擬的是
勻速直線運動
的物體，離散化就是直接對這個方波進行了均勻取樣；但是實際情況中物體的
速度和方向都可能發生變化
，也就是實際產生的模糊是一個
blind
的過程，那麼這個blur kernel其實是沒辦法去猜的。
繼續模擬勻速直線運動物體的blur，在
頻域
上我們這麼操作，用到了opencv的dft函式（據說內部也是fft實現的（fft：快速傅立葉變換））：
enum ConvolutionType // 函式 conv2 卷積時引數的型別
{
CONVOLUTION_FULL， // 卷積時的引數，和 matlab 的 full 一致
CONVOLUTION_SAME， // 卷積時的引數，和 matlab 的 same 一致
CONVOLUTION_VALID // 卷積時的引數，和 matlab 的 valid 一致
}；
void Conv2DFT（Mat& convRes， const Mat& img， const Mat& kernel， ConvolutionType type， int ddepth）
{//為了保持和空域卷積輸出尺寸的一致，以及更快的計算
int dft_M = getOptimalDFTSize（img。rows + kernel。rows - 1）； // 行數
int dft_N = getOptimalDFTSize（img。cols + kernel。cols - 1）； // 列數
Mat imagePad（dft_M， dft_N， CV_32FC1， Scalar（0））；
Mat imagePadROI = imagePad（Rect（0， 0， img。cols， img。rows））；
img。convertTo（imagePadROI， CV_32FC1， 1， 0）；
Mat kernelPad（dft_M， dft_N， CV_32FC1， Scalar（0））；
Mat kernelPadROI = kernelPad（Rect（0， 0， kernel。cols， kernel。rows））；
kernel。convertTo（kernelPadROI， CV_32FC1， 1， 0）；
dft（imagePad， imagePad， 0， imagePad。rows）；//最後一個引數， imagePad。rows，每一行都做卷積
dft（kernelPad， kernelPad， 0， kernelPad。rows）；
// set the last parameter to false to compute convolution instead of correlation
mulSpectrums（imagePad， kernelPad， imagePad， DFT_COMPLEX_OUTPUT， false）； // false： A。*B；true：xf。*conj（yf）
idft（imagePad， imagePad， cv：：DFT_SCALE | cv：：DFT_REAL_OUTPUT， imagePad。rows）；
Rect r；
switch（type） //因為卷積最終的長度大小是img。rows + kernel。rows - 1，所以要決定是否把輸出圖裁剪回原來大小
{
case CONVOLUTION_FULL： // full
r = Rect（0， 0， img。cols + kernel。cols - 1， img。rows + kernel。rows - 1）；
break；
case CONVOLUTION_SAME： // same
r = Rect（（kernel。cols + 0。5） / 2， （kernel。rows + 0。5） / 2， img。cols， img。rows）；
break；
case CONVOLUTION_VALID： // valid
r = Rect（（kernel。cols + 0。5） / 2， （kernel。rows + 0。5） / 2， img。cols - kernel。cols + 1， img。rows - kernel。rows + 1）；
break；
default： // same
r = Rect（（kernel。cols + 0。5） / 2， （kernel。rows + 0。5） / 2， img。cols， img。rows）；
break；
}
imagePad（r）。convertTo（convRes， ddepth， 1， 0）；
}
void test（）{
//實現單方向卷積：頻域乘積
Mat img = imread（“D：\\data\\lenna。jpg”， 0）；
Mat outimgdft = img。clone（）；
int wk = 24； //卷積核長度
Mat kernelx = Mat：：ones（1， wk， CV_32FC1） / wk；
Conv2DFT（outimgdft， img， kernelx， CONVOLUTION_SAME， 0）；
imwrite（“outimgdft。jpg”， outimgdft）；
timer。get_time_ms（“time2：”）；
}
上圖中，左邊是頻域乘積結果，右邊是時域卷積結果，兩者輸出是一致的。滑步到模糊的的lenna！
ps： 頻域操作會比時域快1/2-1/3左右哦。
如果在頻域處理deblur
在現實中，我們只有已經blur了的圖，卻不知道是什麼kernel讓它blur的，我們的目標是把blur的圖還原到清晰的圖。所以選擇對的kernel是很重要的一個步驟。這裡先假設我們已經知道blur kernel是什麼了，那麼我們怎麼在空域把影象還原成原圖呢？好像很少見到。這種時候，在頻域處理的便利性就大大地展現了出來。時域反捲積等於頻域除法。
F（cleanimg）*F（blurkernel）->F（blurimg）->blurimg
blurimg->F（blurimg）->/F（blurkernel）->F（cleanimg）->cleanimg
Mat deblur_DFTi（const Mat& blurimg， const Mat& kernel， ConvolutionType type， int ddepth）
{
Mat convRes；
int dft_M = getOptimalDFTSize（blurimg。rows + kernel。rows - 1）； // 行數
int dft_N = getOptimalDFTSize（blurimg。cols + kernel。cols - 1）； // 列數
Mat imagePad（dft_M， dft_N， CV_32FC1， Scalar（0））；
Mat imagePadROI = imagePad（Rect（0， 0， blurimg。cols， blurimg。rows））；
blurimg。convertTo（imagePadROI， CV_32FC1， 1， 0）；
//要分別用到實部和虛部計算，所以都要保留
Mat planes［］ = { Mat_
Mat complexI；
merge（planes， 2， complexI）；
Mat kernelPad（dft_M， dft_N， CV_32FC1， Scalar（0））；
Mat kernelPadROI = kernelPad（Rect（0， 0， kernel。cols， kernel。rows））；
kernel。convertTo（kernelPadROI， CV_32FC1， 1， 0）；
Mat planesk［］ = { Mat_
Mat complexk；
merge（planesk， 2， complexk）；
dft（complexI， complexI）；
dft（complexk， complexk）；
split（complexI， planes）；
split（complexk， planesk）；
Mat mulout［2］；
//和上面的公式對應的操作
Mat sub = planesk［0］。mul（planesk［0］） + planesk［1］。mul（planesk［1］）；
mulout［0］ = planes［0］。mul（planesk［0］） + planes［1］。mul（planesk［1］）；
divide（mulout［0］， sub， mulout［0］）；//divide（I1，I2，dst，scale，int dtype=-1）；
mulout［1］ = planesk［0］。mul（planes［1］） - planesk［1］。mul（planes［0］）；
divide（mulout［1］， sub， mulout［1］）；
Mat out；
merge（mulout， 2， out）；
idft（out， out， DFT_REAL_OUTPUT）；
normalize（out， out， 0， 1， CV_MINMAX）；
out。convertTo（out， CV_8UC1， 255）；
Rect r；
switch（type）
{
case CONVOLUTION_FULL： // full
r = Rect（0， 0， blurimg。cols + kernel。cols - 1， blurimg。rows + kernel。rows - 1）；
break；
case CONVOLUTION_SAME： // same
r = Rect（（kernel。cols + 0。5） / 2， （kernel。rows + 0。5） / 2， blurimg。cols， blurimg。rows）；
break；
case CONVOLUTION_VALID： // valid
r = Rect（（kernel。cols + 0。5） / 2， （kernel。rows + 0。5） / 2， blurimg。cols - kernel。cols + 1， blurimg。rows - kernel。rows + 1）；
break；
default： // same
r = Rect（（kernel。cols + 0。5） / 2， （kernel。rows + 0。5） / 2， blurimg。cols， blurimg。rows）；
break；
}
out（r）。convertTo（convRes， CV_8UC1， 1， 0）；
return convRes；
}
在這裡要注意一個問題就是blurimg要加合適的border才能得到正確的deblur結果，因為平時我們得到的blur圖片是
，
和
是感測器出圖的長寬，但它大小原本應該是
，因為它是卷積某個blur kernel得到的，所以跟kernel的尺寸有關，在計算的時候需要補全。
不同blurkernel的可逆性
最後，不是所有blur都是可逆的。
blur圖和原圖：
當blurkernel = （1， 0， 1， 0， 1， 1， 0， 1， 0， 0） / 5，deblur的效果，基本可以很好地恢復。
當blurkernel = （1， 1， 1， 1， 1， 1， 0， 1， 0， 1） / 8， deblur的效果，會有輕微水波紋：
當blurkernel = （0， 0， 1， 1， 1， 1， 1， 1， 0， 0） / 6，完全恢復不了。
當blurkernel = （1， 1， 1， 1， 1， 1， 1， 1， 1， 1） / 10，完全恢復不了。
當blurkernel的傅立葉變換有很多0點時，就意味著很難恢復原圖了。
