很多人都認為retinex和暗通道去霧是八桿子都打不着的增強算法。的確，二者的理論、計算方法都完全迥異，本人直接從二者的公式入手來簡單說明一下，有些部分全憑臆想，不對之處大家一起討論。
首先，為描述方便，后面所有的圖像都是歸一化到[0,1]的浮點數圖像。
Retinex的公式就是：
J=I/L (1)
其中，J是所求的圖像，I是觀測圖像，L是估計的光照圖像。注意，由於有J、I、L的值都在區間[0,1]內，則有L>=I成立。（這里使用符號J和I，而不是常規的R和S，主要是為了和暗通道公式保持一致，便於比較）
暗通道去霧的公式為：
J=(I-A)/t+A (2)
其中，A是光照值，t是透射率。如果我們定義遮罩圖像V1， 並令t=1-V1/A, 將其帶入上面公式，則有：
J = (I-V1)/(1-V1/A) (3)
一般有I>=V1, 由於光照A的值一般都偏大，接近於1，那么上面公式再次簡化為：
J = (I-V1)/(1-V1) (4)
仔細觀察公式（1）和（4），你發現相似之處了嗎？
在公式（1）中，I值介於0和L之間，其作用就是將I線性拉伸到[0,1]之間，公式（4）中，I值介於[V1,1]之間，其作用也是將其值線性拉伸到[0,1]之間。
所以，二者是類似的。
如果現在圖像I值介於[V1,L]之間，那么自然地恢復公式是：
J = (I-V1)/(L-V1) (5)
如果直接套用上面公式到普通圖像，效果很容易增強太過，畢竟難以找到又有較強霧霾有光照不足的場景。什么圖像合適呢？標題已經點明了，HDR圖像。V1也就是我前幾天說的暗邊界，L就是亮邊界，三個RGB通道可以共用暗邊界，但L要各自計算。此外，hdr的預處理也非常關鍵，比如進行對數操作，最后要采用非線性歸一化操作。
軟件EXE下載地址：http://pan.baidu.com/s/1cxKU0u
程序采用python實現，未經性能優化，exe中打包了python及numpy wxpython opencv等重量級模塊，故體積較大，如殺毒軟件誤報為病毒，請信任運行。
下面是python處理代碼（沒有讀取hdr文件部分，參見下一篇博客），只有98行哦
import cv2,wx
import numpy as np
from readHdr import readHdr #readHdr程序代碼參見下一遍博客
def zmMinFilterGray(src, r=7): #計算最小值濾波，r是濾波器半徑
'''if r <= 0:
return src
h, w = src.shape[:2]
I = src
res = np.minimum(I , I[[0]+range(h-1) , :])
res = np.minimum(res, I[range(1,h)+[h-1], :])
I = res
res = np.minimum(I , I[:, [0]+range(w-1)])
res = np.minimum(res, I[:, range(1,w)+[w-1]])
return zmMinFilterGray(res, r-1)'''
return cv2.erode(src, np.ones((2*r+1, 2*r+1)))
def guidedfilter(I, p, r, eps): #引導濾波
height, width = I.shape
m_I = cv2.boxFilter(I, -1, (r,r))
m_p = cv2.boxFilter(p, -1, (r,r))
m_Ip = cv2.boxFilter(I*p, -1, (r,r))
cov_Ip = m_Ip-m_I*m_p
m_II = cv2.boxFilter(I*I, -1, (r,r))
var_I = m_II-m_I*m_I
a = cov_Ip/(var_I+eps)
b = m_p-a*m_I
m_a = cv2.boxFilter(a, -1, (r,r))
m_b = cv2.boxFilter(b, -1, (r,r))
return m_a*I+m_b
def stretchImage2(data, vv = 10.0): #非線性拉伸
m = data-data.mean()
S = np.sign(m)
A = np.abs(m)
A = 1.0 - vv**(-A)
m = S*A
vmin, vmax = m.min(), m.max()
return (m-vmin)/(vmax-vmin)
def getV1(m, r, eps, ratio): #對所有通道求同樣暗邊界
tmp = np.min(m,2)
V1 = cv2.blur(zmMinFilterGray(tmp, 7), (7,7))
V1 = guidedfilter(tmp, V1, r, eps)
V1 = np.minimum(V1, tmp)
V1 = np.minimum(V1*ratio, 1.0)
return V1
def getV2(m, r, eps, ratio):
Y = np.zeros(m.shape)
for k in range(3): #對每個通道單獨求亮邊界
v2 = 1 - cv2.blur(zmMinFilterGray(1-m[:,:,k],7), (7,7))
v2 = guidedfilter(m[:,:,k], v2, r, eps)
v2 = np.maximum(v2, m[:,:,k])
Y[:,:,k] = np.maximum(1-(1-v2)*ratio, 0.0)
return Y
def ProcessHdr(m, r, eps, ratio, para1): #單尺度處理
V1 = getV1(m, r, eps, ratio) #計算暗邊界
V2 = getV2(m, r, eps, ratio) #計算亮邊界
Y = np.zeros(m.shape)
for k in range(3):
Y[:,:,k] = ((m[:,:,k]-V1)/(V2[:,:,k]-V1))
Y = stretchImage2(Y,para1) #非線性拉伸
return Y
def ProcessHdrMs(m, r=[161], eps=[0.005,0.001, 0.01], ratio=[0.98, 0.98, 0.92], para1=10.0): #多尺度處理
Y = []
for k in range(len(r)):
Y.append(ProcessHdr(m, r[k], eps[k], ratio[k], para1))
return sum(Y)/len(r)
if __name__ == '__main__':
import glob,os.path
for d in ['auto.hdr',]:
m = readHdr(d) #讀取dhr文件， readHdr程序代碼參見下一遍博客
m1,m2 = m.max(), m.min()
m = (m-m2)/(m1-m2) *255 #數據拉伸到[0,255]
m1 = m[:,:,0].copy(); m[:,:,0] = m[:,:,2]; m[:,:,2]=m1 #顏色通道調整，opencv里R和B反了
m = np.log(m+1)/np.log(256) #log處理
for i in range(2): #如果圖像還是很暗，則需要多次log處理
tmp = np.max(m,2)
tmp = guidedfilter(tmp, tmp, 301, 0.01)
th = np.mean(tmp<0.05)
if th < 0.3:
break
m1 = np.log(m*255+1)/np.log(256)
tmp = np.clip(tmp, 0.0, 1.0) ** (0.05) #tmp是權重參數
for k in range(3): #取加權平均
m[:,:,k] = tmp*m[:,:,k] + (1-tmp)*m1[:,:,k]
m2 = ProcessHdrMs(m)*255
cv2.imwrite('%s.jpg' % d.split('.')[0], m2)
下面上圖。左邊是Luminance-HDR軟件的結果，右邊是我的增強結果。
