# ref[1] He, Kaiming, Jian Sun, and Xiaoou Tang. "Single image haze removal using dark channel prior."
# Pattern Analysis and Machine Intelligence, IEEE Transactions on 33.12 (2011): 2341-2353.
# ref[2] He, Kaiming, Jian Sun, and Xiaoou Tang. "Guided image filtering."
# Pattern Analysis and Machine Intelligence, IEEE Transactions on 35.6 (2013): 1397-1409.
# HAZY IMAGE MODEL:
# I(x) = J(x)t(x) + A(1 - t(x))
# I(x): observed intensity(hazy image)
# J(x): scene radiance
# t(x): atomsphere transmission
# A: global atomsphere light
# AIM: given observed intensity I(x), recover J(x) (scene radiance without haze)
import os
import numpy as np
from numpy.linalg import inv
from PIL import Image
%pylab inline
rcParams['figure.figsize'] = (10, 4) #wide graphs by default
from __future__ import print_function
from __future__ import division
from scipy.io import wavfile
%matplotlib inline
SRC_DIR = 'img/'
DEST_DIR = 'result'
IMG_NAMES = '7.JPG' # ,"2.jpg","3.jpg")
Patch_Size = 10
Top_atom_Rad = 0.05
Rad_Threshold = 170
Omega = 0.95
Min_Transmission = 0.1
Populating the interactive namespace from numpy and matplotlib
WARNING: pylab import has clobbered these variables: ['test'] `%matplotlib` prevents importing * from pylab and numpy
#w is patch size
# the dark channel of an non hazy outdoor image is typically very "dark" See ref1 section 3
# I oberserved intensity(image), W: window size
def get_dark_channel(I, w):
M,N,C = I.shape
darkchnl = zeros((M,N))
padded = np.pad(I, ((int(w / 2), int(w / 2)), (int(w / 2), int(w / 2)), (0, 0)), 'edge')
for i, j in np.ndindex(darkchnl.shape):
darkchnl[i, j] = np.min(padded[i:i + w, j:j + w, :]) # ref[1] tpami, eq.5
return darkchnl
# get atomsphere raidance
# I: oberved intensity. darkchnl: obtained darkchannel. p cut off percent
# ref[1] section 4.3 Get the p percent brightest pixel in darkchannel, find the maximum of the intensity
# for those pixels as atomsphere radiance.
def get_atom_Radiance (I, darkchnl, p):
M, N = darkchnl.shape
flatdark = darkchnl.ravel() #flatten array
SearchIndex = (-flatdark).argsort() # [indexof Maxiumum .... index of Min ]
SearchIndex = SearchIndex[:M*N*p] #cut off percent index
Iflat = I.reshape(M*N, 3)
IflatCut = Iflat.take(SearchIndex, axis = 0)
return np.max(IflatCut, axis = 0)
# get scene transimssion, I observed intensity, A global atomsphere radiance, w: window size
# this function will get raw transmission, which is constant on patch window size w, so the
# raw transmission will suffer from blocky artifacts.
def get_transmission(I, A , w):
return 1 - Omega * get_dark_channel(I/A, w) # ref[1] tpami, eq 12
#get scene radiance. I: observed intensity. t :atom transmission. A estimiated atomphere light
#
def get_Scene_radiance(I, A, t):
trans = np.zeros(I.shape)
trans[:,:,0] = trans[:,:,1] = trans[:,:,2] = t[:,:]
return (I - A) / np.maximum(trans, Min_Transmission) + A #ref1 tpami eq 22
#fast box sum routine
# X, 2d matrix need to be sumed, w: window width
# input x, output is a 2d matrix same shape as the input matrix, each cell is the sum of value of window size w around
# that cell.
def FastBoxSum(X, w):
M,N = X.shape
result = np.zeros((M,N))
prefixSumY = np.cumsum(X, axis=0) #prefix sum along y axis
result[:w+1] = prefixSumY[w:2*w+1]
result[w+1:M-w] = prefixSumY[2*w+1:] - prefixSumY[:M-2*w-1]
result[-w:] = np.tile(prefixSumY[-1], (w,1)) - prefixSumY[M-2*w-1:M-w-1]
prefixSumX = np.cumsum(result , axis=1) #prefix sum along x axis
result[:,:w+1] = prefixSumX[:, w:2*w+1]
result[:,w+1:N-w] = prefixSumX[:, 2*w+1:] - prefixSumX[:, :N-2*w-1]
result[:, -w:] = np.tile(prefixSumX[:, -1][:, None], (1, w)) - prefixSumX[:, N-2*w-1:N-w-1] #[:, None] transpose to column matrix
return result
#p input matrix to be filtered. I, guide image used to filter p, contains, 3 channel. W: filter width
# reference: ref[2] guided image filtering
# the raw extracted tramission is assumed to be constant within a window size w, However, the tranmission
# is not constant on an image. This function will filter the raw transmission to the refined transmission
# under the guidence of the original image
# involves a guidance image I, an filtering input image p, and an output image q
# Core idea: q is a linear transformation of I in the window size w
# Basically it filter the raw transmission so that the refined transmission is 'similiar' in shape of the
# guided image I. By this way, it would remove the blocky artifacts of raw transmission
def Guilded_Filter(p, I, eps, w):
M, N = p.shape
#get average p and norm factor
AveP = np.zeros((M,N))
Norm = np.ones((M,N))
Norm = FastBoxSum(Norm, w)
AveP = FastBoxSum(p, w) / Norm
#get A
AveA = np.zeros(I.shape)
AveU = np.zeros(I.shape)
afin = np.zeros(I.shape)
VarM = np.zeros((3,3,M,N))
AveU[:,:,0] = FastBoxSum(I[:,:,0], w)/Norm
AveU[:,:,1] = FastBoxSum(I[:,:,1], w)/Norm
AveU[:,:,2] = FastBoxSum(I[:,:,2], w)/Norm
AveA[:,:,0] = FastBoxSum(p*I[:,:,0],w)/Norm - AveU[:,:,0]*AveP
AveA[:,:,1] = FastBoxSum(p*I[:,:,1],w)/Norm - AveU[:,:,1]*AveP
AveA[:,:,2] = FastBoxSum(p*I[:,:,2],w)/Norm - AveU[:,:,2]*AveP
#ShowImage(AveA*256)
for i in xrange(3):
for j in xrange(3):
VarM[i,j,:,:] = FastBoxSum(I[:,:,i]*I[:,:,j], w)/Norm - AveU[:,:,i]*AveU[:,:,j]
for i in xrange(M):
for j in xrange(N):
Var = np.array([[VarM[0][0][i, j], VarM[0][1][i, j], VarM[0][2][i, j]],
[VarM[1][0][i, j], VarM[1][1][i, j], VarM[1][2][i, j]],
[VarM[2][0][i, j], VarM[2][1][i, j], VarM[2][2][i, j]]])
Var = Var + eps * np.eye(3)
IP3 = AveA[i,j,:]
afin[i,j,:] = np.dot(inv(Var), IP3)
bfin = AveP - afin[:,:,0]*AveU[:,:,0] - afin[:,:,1]*AveU[:,:,1] - afin[:,:,2]*AveU[:,:,2]
#q_i = ave(a[r])*I[r] + ave(a[g])*I[g] + ave(a[b])*I[b] + ave[b]
qfin = FastBoxSum(afin[:,:,0], w) * I[:,:,0] + FastBoxSum(afin[:,:,1], w) * I[:,:,1] + FastBoxSum(afin[:,:,2], w) * I[:,:,2]\
+ FastBoxSum(bfin, w)
qfin = qfin/Norm
return qfin
#ShowImage(AveU)
def ShowImage(I):
cut = np.maximum(np.minimum(I, 256 - 1), 0).astype(np.uint8)
TmpImg = Image.fromarray(cut)
TmpImg.show()
def SaveImage(I, Path_name):
cut = np.maximum(np.minimum(I, 256 - 1), 0).astype(np.uint8)
TmpImg = Image.fromarray(cut)
TmpImg.save(Path_name)
file_dir = os.path.dirname(os.path.realpath('__file__'))
parent_dir, _ = os.path.split(file_dir)
src_path = os.path.join(parent_dir, SRC_DIR)
dest_path = os.path.join(parent_dir, DEST_DIR)
#filenames = []
base, ext = os.path.splitext(IMG_NAMES)
tempname = base + ext
src_path = src_path + IMG_NAMES
#dest_path = dest_path + tempname
#print src_path
#print dest_path
im = imread(src_path)
figure(figsize = (20,30))
imshow(im)
title("Input Image")
<matplotlib.text.Text at 0x110a78690>
I = np.asarray(im, dtype=np.float64)
Darkchnl = get_dark_channel(I, Patch_Size) #get dark channel
SaveImage(Darkchnl, dest_path + base + 'Dark' +ext)
imshow(Darkchnl)
title("dark channel")
<matplotlib.text.Text at 0x10c2d3150>
atom_Radiance = get_atom_Radiance(I, Darkchnl, Top_atom_Rad)
atom_Radiance = np.minimum(atom_Radiance, Rad_Threshold)
transmission = get_transmission(I, atom_Radiance, Patch_Size)
Scene_radance = get_Scene_radiance(I, atom_Radiance, transmission)
SaveImage(Scene_radance, dest_path + base + 'Ran' + ext)
SaveImage(transmission*256, dest_path + base + 'Rawdepth' + ext)
/usr/local/lib/python2.7/site-packages/ipykernel/__main__.py:10: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
test = imread(dest_path + base + 'Ran' + ext);
figure(figsize = (20,30))
imshow(test)
title("raw scence randiance")
<matplotlib.text.Text at 0x110316310>
#print transmission.shape
#print atom_Radiance
Inorm = I / (I.max() - I.min())
tfine = Guilded_Filter(transmission, Inorm, 1e-3, 40)
#tfine = guided_filter(I, transmission, 40, 1e-3)
Scene_radance_new = get_Scene_radiance(I, atom_Radiance, tfine)
ShowImage(tfine*256) #refined depth map
imshow(tfine)
SaveImage(Scene_radance_new, dest_path + base + 'RadN' + ext)
SaveImage(tfine*256, dest_path + base + 'depth' + ext)
test = imread(dest_path + base + 'RadN' + ext);
figure(figsize = (20,30))
imshow(test)
<matplotlib.image.AxesImage at 0x110243cd0>
tt = maximum(tfine, Min_Transmission);
SaveImage(-log(tt)*100, dest_path + base + 'map' + ext)