最近在做圖像增強方面的算法，在參考了一些博客，論文和源代碼后 ，自己整理了Retinex相關算法的opencv實現，在這里貼出來供大家參考，有不當的地方歡迎大家指出！
一.Retinex算法原理
基礎理論：物體的顏色是由物體對長波（紅色），中波（綠色），短波（藍色）光線的反射能力來決定的，而不是由反射光強度的絕對值來決定的；物體的顏色不受光照非均勻性的影響，具有一致性，即Retinex算法是基于物體的顏色恒常性來實現的。
Retinex算法可以在圖像的動態顏色范圍壓縮，邊緣增強和顏色恒常三個方面達到平衡，在圖像除霧方面有著較好的效果（對正常圖像的增強效果不明顯）。
二.SSR（Single Scale Retinex）
根據Retinex算法的基礎理論，我們可以得到以下數學表達式：
S(x,y)=R(x,y)*L(x,y)　 （2-1）
圖示如下：

其中R（x,y）表示物體的反射性質，即圖像的內在屬性，應該最大程度的保留；L（x，y）表示入射光圖像，決定了圖像像素能夠達到的動態范圍，應該盡量去除；S（x,y）為人眼觀察到或者相機接收到的圖像。
對公式2-1的等式兩邊取對數，得到：
r(x,y)=Log[R(x,y)] = Log[S(x,y)]-Log[L(x,y)] (2-2)

算法的關鍵在于如何得到圖像的入射光圖像L（x,y），其中比較經典且效果較好的方法是通過對相機接收的圖像S（x,y）做高斯模糊來估計L（x,y）。
算法的具體步驟如下：
（1）輸入高斯模糊的高斯環繞尺度C（即二維高斯函數的標準差）；
（2）根據C對原始圖像數據S（x，y）做高斯模糊得到L（x,y）;
（3）根據公式2-2對S（x,y）和L（x,y）分別取對數并作差得到r(x,y);
（4）將r（x,y）的像素值量化到0到255的范圍內得到R（x,y），R（x,y）即我們想得到的增強圖像。量化的公式如下：
R(x,y) = ( Value - Min ) / (Max - Min) * (255-0) (2-3)
算法的源代碼如下：

void ssr(Mat src, Mat& dst, double sigma) {Mat src_log,gauss,gauss_log,dst_log;src_log= Mat(src.size(), CV_32FC3);gauss_log= Mat(src.size(), CV_32FC3);dst_log = Mat(src.size(), CV_32FC3);dst = Mat(src.size(), CV_32FC3);int height = dst_log.rows;int width = dst_log.cols;int ksize = (int)(sigma * 3/2);ksize = ksize * 2 + 1;//求Log(S(x,y)for (int i = 0;i < height;i++) {for (int j = 0;j < width;j++) {for (int k = 0;k < 3;k++) {float value = src.at<Vec3b>(i, j)[k];if (value <= 0.01) value = 0.01;src_log.at<Vec3f>(i, j)[k] = log10(value);}}}GaussianBlur(src, gauss, Size(ksize,ksize),sigma,sigma,4);//求Log（L(x,y)）for (int i = 0;i < height;i++) {for (int j = 0;j < width;j++) {for (int k = 0;k < 3;k++) {float value = gauss.at<Vec3b>(i, j)[k];if (value <= 0.01) value = 0.01;gauss_log.at<Vec3f>(i, j)[k] = log10(value);}}}//求Log（S（x,y）)-Log(L(x,y))for (int i = 0;i < height;i++) {for (int j = 0;j < width;j++) {for (int k = 0;k < 3;k++) {float value1 = src_log.at<Vec3f>(i, j)[k];float value2 = gauss_log.at<Vec3f>(i, j)[k];dst_log.at<Vec3f>(i, j)[k] = value1-value2;}}}float min[3] = { dst_log.at<Vec3f>(0, 0)[0], dst_log.at<Vec3f>(0, 0)[1],dst_log.at<Vec3f>(0, 0)[2] };float max[3] = { dst_log.at<Vec3f>(0, 0)[0], dst_log.at<Vec3f>(0, 0)[1],dst_log.at<Vec3f>(0, 0)[2] };//求R/G/B三通道的min,maxfor (int i = 0;i < height;i++) {for (int j = 0;j < width;j++) {for (int k = 0;k < 3;k++) {float value = dst_log.at<Vec3f>(i, j)[k];if (value > max[k]) max[k] = value;if (value < min[k]) min[k] = value;}}}//量化處理cout << min[0] << " " << min[1] << " " << min[2] << endl;cout << max[0] << " " << max[1] << " " << max[2] << endl;for (int i = 0;i < height;i++) {for (int j = 0;j < width;j++) {for (int k = 0;k < 3;k++) {float value = dst_log.at<Vec3f>(i, j)[k];dst.at<Vec3f>(i, j)[k] = (saturate_cast<float>(255 * (value - min[k]) / (max[k] - min[k])));}}}dst.convertTo(dst, CV_8UC3);return; }

效果圖如下：
原圖：

SSR（C=300）：

原圖：

SSR（C=300）：

三.MSR（Multi Scale Retinex）
MSR是在SSR算法的基礎上提出的，用不同的尺度C來估計L（x，y）。數學表達如下：
r(x,y)=∑k Wk*{logS(x,y)?log[Fk(x,y)?S(x,y)]} （3-1）
為了兼有SSR高，中，低三個尺度的優點的考慮，K通常取3，且有W1=W2=W3=1/3.
算法步驟：
（1）輸入權值矩陣Wk和K個高斯環繞尺度；
（2）根據公式3-1得到r（x,y）;
（3）對r（x，y）做量化處理得到增強圖像R（x,y）.
算法的源代碼如下：

void msr(Mat src, Mat& dst, vector<float>weight, vector<float>sigmas) {Mat src_log, gauss, gauss_log, dst_log;src_log = Mat(src.size(), CV_32FC3);gauss_log = Mat(src.size(), CV_32FC3);dst_log = Mat::zeros(src.size(), CV_32FC3);dst = Mat::zeros(src.size(), CV_32FC3);int height = dst_log.rows;int width = dst_log.cols;//求Log(S(x,y)for (int i = 0;i < height;i++) {for (int j = 0;j < width;j++) {for (int k = 0;k < 3;k++) {float value = src.at<Vec3b>(i, j)[k];if (value < 0.01) value = 0.01;src_log.at<Vec3f>(i, j)[k] = log10(value);}}}int scale = weight.size();for (int t = 0;t < scale;t++) {int ksize = (int)(sigmas[t] * 3 / 2);ksize = ksize * 2 + 1;GaussianBlur(src, gauss, Size(ksize, ksize), sigmas[t], sigmas[t], 4);//求Log（L(x,y)）for (int i = 0;i < height;i++) {for (int j = 0;j < width;j++) {for (int k = 0;k < 3;k++) {float value = gauss.at<Vec3b>(i, j)[k];if (value < 0.01) value = 0.01;gauss_log.at<Vec3f>(i, j)[k] = log10(value);}}}//求Log（S（x,y）)-Log(L(x,y))for (int i = 0;i < height;i++) {for (int j = 0;j < width;j++) {for (int k = 0;k < 3;k++) {float value1 = src_log.at<Vec3f>(i, j)[k];float value2 = gauss_log.at<Vec3f>(i, j)[k];dst_log.at<Vec3f>(i, j)[k] +=weight[t] * (value1 - value2);}}}}float min[3] = { dst_log.at<Vec3f>(0, 0)[0], dst_log.at<Vec3f>(0, 0)[1],dst_log.at<Vec3f>(0, 0)[2] };float max[3] = { dst_log.at<Vec3f>(0, 0)[0], dst_log.at<Vec3f>(0, 0)[1],dst_log.at<Vec3f>(0, 0)[2] };//求R/G/B三通道的min,maxfor (int i = 0;i < height;i++) {for (int j = 0;j < width;j++) {for (int k = 0;k < 3;k++) {float value = dst_log.at<Vec3f>(i, j)[k];if (value > max[k]) max[k] = value;if (value < min[k]) min[k] = value;}}}//量化處理for (int i = 0;i < height;i++) {for (int j = 0;j < width;j++) {for (int k = 0;k < 3;k++) {float value = dst_log.at<Vec3f>(i, j)[k];dst.at<Vec3f>(i, j)[k] =saturate_cast<float> (255 * (value - min[k]) / (max[k] - min[k]));}}}dst.convertTo(dst, CV_8UC3); }

算法的效果圖如下：
原圖：

MSR(C1=40,C2=100,C3=300):

原圖：

MSR（C1=40,C2=100,C3=300）：

四.MSRCR（帶顏色恢復的MSR）
由SSR和MSR算法的效果圖我們可以看到一般的Retinex算法在圖像去霧時可能會導致圖像失真。一般的Retinex算法處理圖像時，都會假設初始圖像灰度值是緩慢變化的，即圖像是平滑的，在實際情況下，Retinex算法在亮度差異大的區域會產生光暈。
為了解決上述的圖像顏色失真的情況，有人提出了一種帶有顏色恢復的MSR算法。MSRCR算法在MSR算法的基礎上，加入了色彩恢復因子來調節由于圖像局部區域對比度增強導致的顏色失真。
算法核心如下：

算法源代碼如下：

void scr(Mat &src, float low_clip, float high_clip) {int height = src.rows;int width = src.cols;int total = height*width;int u[3][256];//統計每個通道每個像素值出現的次數memset(u, 0, sizeof(u));for (int i = 0;i < height;i++) {for (int j = 0;j < width;j++) {for (int k = 0;k < 3;k++) {int v = src.at<Vec3b>(i, j)[k];u[k][v]++;}}}for (int i = 0; i < 3; i++) {for (int j = 1;j < 256;j++) {u[i][j] = u[i][j - 1] + u[i][j];}}int low_val[3], high_val[3];float low_rate = 0, high_rate = 0;for (int i = 0; i < 3; i++) {for (int j = 0;j < 256;j++) {float rate = u[i][j] / total;if (rate< low_clip&&rate>low_rate) {low_val[i] = j;low_rate = rate;}if (rate < high_clip&&rate>high_rate) { high_val[i] = j;high_rate = rate;}}}for (int i = 0;i < height;i++) {for (int j = 0;j < width;j++) {for (int k = 0;k < 3;k++) {int value = src.at<Vec3b>(i, j)[k];src.at<Vec3b>(i, j)[k] = max(min(value, high_val[k]), low_val[k]);}}} } void colorRestoration(Mat src, Mat& dst, float alpha, float beta) {dst = Mat(src.size(), CV_32FC3);for (int i = 0; i < src.rows; i++){for (int j = 0;j < src.cols;j++){ float sum = 0;for (int t = 0;t < src.channels();t++) sum += src.at<Vec3b>(i, j)[t];for (int k=0;k<3;k++){dst.at<Vec3f>(i, j)[k] = beta*(log10(alpha*src.at<Vec3b>(i, j)[k]) - log10(sum));}}} } void msrcr(Mat src, Mat& dst, vector<float>weight, vector<float>sigmas, float alpha, float beta, float low_clip, float high_clip) {Mat src_log, gauss, gauss_log, dst_log,dst_ci;src_log = Mat(src.size(), CV_32FC3);gauss_log = Mat(src.size(), CV_32FC3);dst_log = Mat::zeros(src.size(), CV_32FC3);dst_ci = Mat(src.size(), CV_32FC3);dst = Mat::zeros(src.size(), CV_32FC3);int height = dst_log.rows;int width = dst_log.cols;//求Log(S(x,y)for (int i = 0;i < height;i++) {for (int j = 0;j < width;j++) {for (int k = 0;k < 3;k++) {float value = src.at<Vec3b>(i, j)[k];if (value < 0.01) value = 0.01;src_log.at<Vec3f>(i, j)[k] = log10(value);}}}int scale = weight.size();for (int t = 0;t < scale;t++) {int ksize = (int)(sigmas[t] * 3 / 2);ksize = ksize * 2 + 1;GaussianBlur(src, gauss, Size(ksize, ksize), sigmas[t], sigmas[t], 4);//求Log（L(x,y)）for (int i = 0;i < height;i++) {for (int j = 0;j < width;j++) {for (int k = 0;k < 3;k++) {float value = gauss.at<Vec3b>(i, j)[k];if (value < 0.01) value = 0.01;gauss_log.at<Vec3f>(i, j)[k] = log10(value);}}}//求Log（S（x,y）)-Log(L(x,y))for (int i = 0;i < height;i++) {for (int j = 0;j < width;j++) {for (int k = 0;k < 3;k++) {float value1 = src_log.at<Vec3f>(i, j)[k];float value2 = gauss_log.at<Vec3f>(i, j)[k];dst_log.at<Vec3f>(i, j)[k] += weight[t] * (value1 - value2);}}}}//求色彩恢復因子CicolorRestoration(src, dst_ci, alpha, beta);for (int i = 0;i < height;i++) {for (int j = 0;j < width;j++) {for (int k = 0;k < 3;k++) {float v1 = dst_log.at<Vec3f>(i, j)[k];float v2 = dst_ci.at<Vec3f>(i, j)[k];dst_log.at<Vec3f>(i, j)[k] = v1*v2 ;}}}float min[3] = { dst_log.at<Vec3f>(0, 0)[0], dst_log.at<Vec3f>(0, 0)[1],dst_log.at<Vec3f>(0, 0)[2] };float max[3] = { dst_log.at<Vec3f>(0, 0)[0], dst_log.at<Vec3f>(0, 0)[1],dst_log.at<Vec3f>(0, 0)[2] };//求R/G/B三通道的min,maxfor (int i = 0;i < height;i++) {for (int j = 0;j < width;j++) {for (int k = 0;k < 3;k++) {float value = dst_log.at<Vec3f>(i, j)[k];if (value > max[k]) max[k] = value;if (value < min[k]) min[k] = value;}}}//量化處理for (int i = 0;i < height;i++) {for (int j = 0;j < width;j++) {for (int k = 0;k < 3;k++) {float value = dst_log.at<Vec3f>(i, j)[k];dst.at<Vec3f>(i, j)[k] = saturate_cast<float> (255 * (value - min[k]) / (max[k] - min[k]));}}}dst.convertTo(dst, CV_8UC3);scr(dst, low_clip, high_clip); }

上面的算法有著很多的經驗參數，不利于手動實現，有人提出了一種基于GIMP的MSR算法，用參數Dynamic來控制圖像的動態范圍實現去除圖像的色偏問題（一般取值2或3效果比較好）。這種算法的效果很好，能夠一定程度解決顏色失真的問題。
算法的簡要描述如下：
1.計算出 log[R(x,y)]中R/G/B各通道數據的均值Mean和均方差Var（注意是均方差）。
2.類似下述公式計算各通道的Min和Max值。
Min = Mean - Dynamic * Var;
Max = Mean + Dynamic * Var;
3.對Log[R(x,y)]的每一個值Value，進行線性映射：
R(x,y) = ( Value - Min ) / (Max - Min) * (255 - 0)
同時要注意增加一個溢出判斷，即：
if (R(x, y) > 255) R(x,y) = 255;
else if (R(x,y) < 0) R(x,y) = 0;

算法的源碼如下：

void msrcr_GIMP(Mat src, Mat& dst, vector<float>weight, vector<float>sigmas, int Dynamic) {Mat src_log, gauss, gauss_log, dst_log;src_log = Mat(src.size(), CV_32FC3);gauss_log = Mat(src.size(), CV_32FC3);dst_log = Mat(src.size(), CV_32FC3);dst = Mat(src.size(), CV_32FC3);float min[3] = { 0 };float max[3] = { 0 };float mean[3] = { 0 };float var[3] = { 0 };int height = dst_log.rows;int width = dst_log.cols;//求Log(S(x,y)for (int i = 0;i < height;i++) {for (int j = 0;j < width;j++) {for (int k = 0;k < 3;k++) {float value = src.at<Vec3b>(i, j)[k];if (value < 0.01) value = 0.01;src_log.at<Vec3f>(i, j)[k] = log10(value);}}}int scale = weight.size();for (int t = 0;t < scale;t++) {int ksize = (int)(sigmas[t] * 3 / 2);ksize = ksize * 2 + 1;GaussianBlur(src, gauss, Size(ksize, ksize), sigmas[t], sigmas[t], 4);//求Log（L(x,y)）for (int i = 0;i < height;i++) {for (int j = 0;j < width;j++) {for (int k = 0;k < 3;k++) {float value = gauss.at<Vec3b>(i, j)[k];if (value < 0.01) value = 0.01;gauss_log.at<Vec3f>(i, j)[k] = log10(value);}}}//求Log（S（x,y）)-Log(L(x,y))for (int i = 0;i < height;i++) {for (int j = 0;j < width;j++) {for (int k = 0;k < 3;k++) {float value1 = src_log.at<Vec3f>(i, j)[k];float value2 = gauss_log.at<Vec3f>(i, j)[k];dst_log.at<Vec3f>(i, j)[k] += weight[t] * (value1 - value2);}}}}//求R/G/B各通道的均值mean和均方差varfloat sum[3] = {0};for (int i = 0;i < height;i++) {for (int j = 0;j < width;j++) {for (int k = 0;k < 3;k++) {float value = dst_log.at<Vec3f>(i, j)[k];sum[k] += value;}}}for (int i = 0; i < 3; i++){mean[i] = sum[i] / (height*width);}for (int i = 0;i < height;i++) {for (int j = 0;j < width;j++) {for (int k = 0;k < 3;k++) {float value = dst_log.at<Vec3f>(i, j)[k];var[k] += (value - mean[k])*(value - mean[k]);}}}for (int i = 0; i < 3; i++){var[i] = sqrt(var[i] / (height*width));}//計算R/G/B各通道的Min，Maxfor (int i = 0; i < 3; i++){min[i] = mean[i] - Dynamic*var[i];max[i] = mean[i] + Dynamic*var[i];}//量化處理for (int i = 0;i < height;i++) {for (int j = 0;j < width;j++) {for (int k = 0;k < 3;k++) {float value = dst_log.at<Vec3f>(i, j)[k];dst.at<Vec3f>(i, j)[k] = saturate_cast<float>(255 * (value - min[k]) / (max[k] - min[k]));}}}dst.convertTo(dst, CV_8UC3); }

算法的效果圖如下：
MSR(C1=40,C2=100,C3=300)：

MSRCR_GIMP(C1=40,C2=100,C3=300,Dynamic=2):

原圖：

MSRCR_GIMP(C1=40,C2=100,C3=300,Dynamic=2):

以上就是我的Retinex算法的總結過程，算法本人已經實現過，效果圖也貼出來供大家參考，有不當的地方歡迎大家指出來，樓主會及時修正，最后希望我的總結能夠幫到大家！

總結

以上是生活随笔為你收集整理的Retinex算法的C++/opencv实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。