博客

(9条消息) 局部边(9条消息) 局部边缘保持滤波(LEP)高动态范围图像HDR压缩 matlab程序（一）_

 本文程序根据论文”Gu B, Li W, Zhu M, et al. Local edge-preserving multiscale decomposition for high dynamic range image tone mapping[J]. Image Processing, IEEE Transactions on, 2013, 22(1): 70-79.”的matlab程序改编而成，关于局部边缘保持滤波(LEP)高动态范围图像HDR压缩可以参考如下博文

https://blog.csdn.net/majinlei121/article/details/50420927/

https://blog.csdn.net/majinlei121/article/details/50420980/

本人实现的python程序基于python3，代码为：

import numpy as np
import cv2
from math import pi
import os
from tqdm import *
import warnings
warnings.filterwarnings("ignore")
from PIL import Image


def boxfilter(imSrc, r):
    '''
       BOXFILTER   O(1) time box filtering using cumulative sum

    - Definition imDst(x, y)=sum(sum(imSrc(x-r:x+r,y-r:y+r)));
    - Running time independent of r; 
    - Equivalent to the function: colfilt(imSrc, [2*r+1, 2*r+1], 'sliding', @sum);
    - But much faster.
    '''
    hei, wid = imSrc.shape
    imDst = np.zeros(imSrc.shape)
    # cumulative sum over Y axis
    imCum = np.cumsum(imSrc, axis=0)
    # difference over Y axis
    imDst[0:r+1, :] = imCum[r:2*r+1, :]
    imDst[r+1:hei-r, :] = imCum[2*r+1:hei, :] - imCum[0:hei-2*r-1, :]
    # m = np.tile(imCum[hei-1, :], (r, 1))
    # n = imCum[hei-2*r-1:hei-r-1, :]
    imDst[hei-r:hei, :] = np.tile(imCum[hei-1, :], (r, 1)) - imCum[hei-2*r-1:hei-r-1, :]

    # cumulative sum over X axis
    imCum = np.cumsum(imDst, axis=1)
    # difference over X axis
    imDst[:, 0:r+1] = imCum[:, r:2*r+1]
    imDst[:, r+1:wid-r] = imCum[:, 2*r+1:wid] - imCum[:, 0:wid-2*r-1]
    # m = np.tile(np.expand_dims(imCum[:, wid-1], 1), (1, r))
    # n = imCum[:, wid-2*r-1:wid-r-1]
    imDst[:, wid-r:wid] = np.tile(np.expand_dims(imCum[:, wid-1], 1), (1, r)) - imCum[:, wid-2*r-1:wid-r-1]
    return imDst





def LocalWls_HDR(I, alpha=0.1, beta=1, r=4):
    hei , wid = I.shape
    N = boxfilter(np.ones((hei, wid)), r)# the size of each local patch; N=(2r+1)^2 except for boundary pixels.

    mean_I = boxfilter(I, r) / N
    mean_II = boxfilter(I*I, r) / N
    var_I = mean_II - mean_I * mean_I

    I_x = np.diff(I, n=1, axis=1)
    I_x = np.pad(I_x,((0,0),(0,1)),'constant',constant_values=(0,0))#进行列方向差分，求dx
    I_y = np.diff(I, n=1, axis=0)
    I_y = np.pad(I_y,((0,1),(0,0)),'constant',constant_values=(0,0))#进行行方向差分，求dy
    I_gradient = abs(I_x+I_y)
    I_gradient = np.power(I_gradient,2-beta)
    I_gradient = alpha*boxfilter(I_gradient,r)/N

    a = var_I / (var_I + I_gradient+0.00000001)
    b = mean_I - a * mean_I

    mean_a = boxfilter(a, r) / N
    mean_b = boxfilter(b, r) / N

    return mean_a * I + mean_b


def lepHDR(src,dst,alpha = 0.1,beta = 1,r = 2,nLevel = 3):
    HDR = cv2.imread(src)
    m,n,_ = HDR.shape
    HDR = cv2.resize(HDR,(int(n/3),int(m/3)),interpolation=cv2.INTER_CUBIC)
    # HDR = HDR/255.0
    HDR=cv2.normalize(HDR.astype('float'), None, 0.0, 1.0, cv2.NORM_MINMAX)
    L_in = (HDR[:,:,0]+HDR[:,:,1]+HDR[:,:,2])/3
    L = np.log(L_in*1e6+1)
    L = L/np.max(L)  
    B = [[],[],[]]
    D = [[],[],[]]
    D_compression=D_compression = [[],[],[]]
    B[nLevel-1] = L
    for j in (2,1):
        B[j-1] = LocalWls_HDR(B[j], alpha, beta, r)
        D[j] = B[j]-B[j-1]
        r = 20

    B0=np.mean(B[0])*np.ones(L.shape)
    D[0] = B[0]-B0

    for j in (2,1,0):
        D_compression[j] = (2/pi)*np.arctan(20*D[j])
        D_compression[j] = cv2.normalize(D_compression[j].astype('float'), None, 0.0, 1.0, cv2.NORM_MINMAX)

    L_out=D_compression[0]*0.5+D_compression[1]+D_compression[2]

    Rmax_clip = np.percentile(L_out,99.5)
    Rmin_clip = np.percentile(L_out,0.5)
    L_out[L_out>Rmax_clip] = Rmax_clip
    L_out[L_out<Rmin_clip] = Rmin_clip
    L_out=cv2.normalize(L_out.astype('float'), None, 0.0, 1.0, cv2.NORM_MINMAX)
    out = np.zeros(HDR.shape)
    out[:,:,0] = (np.power((HDR[:,:,0]/L_in),0.6))*L_out
    out[:,:,1] = (np.power((HDR[:,:,1]/L_in),0.6))*L_out
    out[:,:,2] = (np.power((HDR[:,:,2]/L_in),0.6))*L_out
    out = out*255
    out[out>255] = 255
    out[out<0] = 0
    out = np.uint8(out)
    # out = np.uint8(cv2.normalize(out.astype('float'), None, 0.0, 255.0, cv2.NORM_MINMAX))
    out = cv2.resize(out,(n,m),interpolation=cv2.INTER_CUBIC)
    cv2.imwrite(dst,out)




if __name__ =="__main__":
    path = "/run/media/charlie/DATA/1.JPG"
    save_path = "/run/media/charlie/DATA/temp"
    for file in tqdm(os.listdir(path)):
        lepHDR(os.path.join(path,file),os.path.join(save_path,file[:-4]+".png"))
    # lepHDR("/raid/逆光数据/test.JPG","/raid/逆光数据/tests.png")