What’s in this article?
I’d like to explain and implement the Single Image Haze Removal Using Dark Channel Prior by Kaiming He.
The paper is the CVPR 2009 Best Paper.
Everyone can do Computer Vision (CV). It’s easy to come up with an over-engineered and complicated solution for a CV problem. But it’s hard to have an elegant solution.
In the paper, Kaiming He designed a simple and elegant solution for the haze removal problem.
1. What is the Haze Removal?
A picture explains the haze removal problem.
image comes from Kaiming’s website.
image comes from Kaiming’s website.
To remove haze computationally, we need to know the mathematical model of haze.
2. Haze Model
A natural image can be modeled as the weighted sum of the scene radiance \(J\) and the global atmospheric light \(A\).
The weight is given by the transmission map \(t(x)\). It is a function of the distance to an object. i.e. \(t(x)\) is different for each pixel.
The scene radiance \(J\) is the light directly reflected or emitted by objects.
the global atmospheric light (or the ambient light) \(A\) is basically light reflected by particles in the air, or, in this case, the haze. The \(A\) is modeled as a constant image.
The goal of the dehaze is:
Given a image \(I(x)\), remove the haze \(A\) to recover the scene radiance \(J(x)\).
Concretely, the high-level ideas of the paper are,
- Estimate the haze \(A\) by an ad-hoc method.
- Using the dark-channel prior to estimation the weight of the haze model \(t(x)\).
- Remove the haze by \(J(x) = \frac{I(x) – A(1 – t(x))}{t(x)}\).
3. The Dark Prior
The author defined the dark channel of the image and statistically show that the dark channel of the scene radiance is close to zero.
The author defined the dark channel of a image as,
The \(J^{\text{dark}}\) is the minimum of all RBG intensities in a local window \(\Omega\). By this definition, we can compute the dark channel of an image by the sliding window method.
A key observation from the author was,
“on haze-free outdoor images: in most of the non-sky patches, at least one color channel has very low intensity at some pixels”
Kaiming He
In other word, the dark channel of an image with only scene radiance \(J(x)\) is close to 0.
For example, consider this haze-free image,
The dark channel of the Pittsburgh haze-free image is very dark!
Whereas the dark channel of a hazed image is bright.
The author showed the result qualitatively,
and quantitatively.
The conclusion is,
The dark channel of haze-free image is 0.
Since the haze-free is an image with only the scene radiance, we have this equation,
Why the intensity of the dark channel of the scene radiance is close to zero? Kaiming gave an argument,
Point (b) is interesting. The draw channel of a colorful image is close to 0. e.g. The RGB of green is (0, 0, 255). The dark channel of green is 0. The RGB of gray is (150, 150, 150). The dark channel of gray is 150.
In later section, the author will use this prior to derive an equation to compute the transmission map \(t(x)\).
Implementation
To compute the dark channel, we simply need to follow the definition.
cv::Mat find_dark_channel(const cv::Mat& m, const int window_half_size)
{
cv::Mat m_pad;
cv::copyMakeBorder(m, m_pad, window_half_size, window_half_size,
window_half_size, window_half_size, cv::BORDER_REPLICATE);
cv::Mat dark_chaneel = cv::Mat(m.rows, m.cols, CV_32FC1, 0.);
for (int row_idx = window_half_size; row_idx < m_pad.rows - window_half_size; ++row_idx) {
for (int col_idx = window_half_size; col_idx < m_pad.cols - window_half_size; ++col_idx) {
dark_chaneel.at<float>(row_idx - window_half_size, col_idx - window_half_size)
= find_min_in_rbg_patch(m_pad, row_idx, col_idx, window_half_size);
}
}
return dark_chaneel;
}
inline float find_min_in_rbg_patch(const cv::Mat& m, const int row, const int col, const int window_half_size)
{
float min_val = std::numeric_limits<float>::max();
for (int dy = -window_half_size; dy <= window_half_size; ++dy) {
for (int dx = -window_half_size; dx <= window_half_size; ++dx) {
const cv::Vec3f gbr = m.at<cv::Vec3f>(row + dy, col + dx);
min_val = std::min(gbr[0], min_val);
min_val = std::min(gbr[1], min_val);
min_val = std::min(gbr[2], min_val);
}
}
return min_val;
}Line 10~15: the sliding window.
Line 20~30: find the dark channel in each window.
4. The magic to estimate \(t(x)\)
The next step is estimating the transmission map \(t(x)\) by the dark channel prior.
Given,
The author did a magic to estimate the \(t(x)\).
Start with equation (1) and assuming \(t(x)\) is constant in a window, the author took 2 minimum operation to find the minimum values of all RBG intensities in a window \(\Omega\). Note, it is the definition of the dark channel.
By definition, the mid term of (8),
\(\min_c \min_{y \in \Omega(x)} \frac{J^c (y)}{A^c}\)
is exactly the dark channel of a image given by \(\frac{J(x)}{A}\).
In section 3, we know the dark channel of the scene radiance is 0. It leads to equation (9). Divide (9) by \(A\), which is positive by physics, we have equation (10).
Using (10) the author got equation (11). It is the equation to compute the transmission \(t(x)\).
Interestingly, the second term in (11) is the dark channel of \(\frac{I(y)}{A}\).
The implementation is straightforward, I will skip it.
Now we have a method to compute the transmission map \(t(x)\).
Recall the big picture that we need to remove haze by \(J(x) = \frac{I(x) – A(1 – t(x))}{t(x)}\). So the next step is to compute the atmosphere light \(A\).
5 Estimate the Atmosphere Light(Haze)
How to estimate the atmosphere light in a hazed image?
To solution is surprising simple. Find pixels corresponding to haze!
For example,
The red pixels in the above image corresponds to haze. We can save the value or take the mean of these pixels to compute the intensity of the atmosphere light (haze) \(A\).
Please read the section 4.4 of the paper for details.
6. The Soft Matting
In section 5, the authors computed the transmission map \(t(x)\) by compute the dark channel of \(\frac{I(y)}{A}\).
Now we should able to remove haze by \(J(x) = \frac{I(x) – A(1 – t(x))}{t(x)}\). However, there is a problem.
By definition of the dark channel, the output of the dark channel is the minimum of a local window. As a result, the dark channel has an artifact. The dark channel looks like a dilated version of the original image.
For example,
Right: original image.
The dark channel is thicker than the original image.
Because we will use \(J(x) = \frac{I(x) – A(1 – t(x))}{t(x)}\) per pixel to remove haze, we better make the shape of \(t(x)\) closer to \(I(x)\). A classical method to do it is the Soft Matting.
For example,
Right: The dark channel after the soft matting.
The goal of the soft matting is to find an image such that,
(1) Its shape is similar to the original image \(p\).
(2) Its intensity is similar to the dark channel \(t\).
It can be formulated as an optimization problem.
\(\text{cost}(x, a, b) = \sum_{j \in I} \sum_{i \in \Omega} ||a_i x_j + b_i – p_j||^2 + ||x – t||^2\)
Where \(j \in I\) means (2d) index \(j\) in image coordinate \(I\). \(i \in \Omega\) means index i in a local window \(\Omega\).
The first term \(\sum_{j \in I} \sum_{i \in \Omega} ||a_i x_j + b_i – p_j||^2\) wants the prior image \(p_j\) be to a linear combination of the optimization variable \(x_j\) (a pixel in image x). Intuitively, If we have \(p = a x + b\), the edge (gradient) satisfies \(\nabla p = a \nabla x\). Thus the shape of \(x\) is similar to the original image \(p\). In this case, we want the shape of \(x\) to be close to the image p (the original image).
shape is similar
The second \(||x – t||^2\) term wants the optimization variable \(x\) to be close to a prior. In this case, we want the value of \(x\) to be close to the transmission map \(t(x)\).
Intensity is similar
The problem is formulated. We can solve it.
The problem is a least squares. We can do the Gaussian Newton to solve it. But the number of \(x\) is the 5 times number of pixels (3 for a, 1 for b, 1 for t per pixel). Solving the normal equation can be slow. (Although I feel like using a sparse solver is feasible).
6.1 Gradient Descent
We can do a gradient descent to optimize,
\(\text{cost}(x, a, b) = \sum_{j \in I} \sum_{i in \Omega} ||a_i x_j + b_i – p_j||^2 + ||x – t||^2\).
The cost is a quadratic \(\text{cost}(x) = x^T W x\). the gradient of it is \(2 W x\).
Implementation
// goal: 1. t = a * target + b
// 2. t close to t_prior
// cost: sum_{j:patch idx} sum_{i:pixel idx in patch} (t_i - a_j * im_i - b_j)^2
// + weight * sum{i: pixel_index } (t_i - t_prior_i)^2
void compute_gradient()
{
grad_A_.setTo(cv::Vec3f(0., 0., 0.));
grad_B_.setTo(cv::Scalar(0.));
grad_T_.setTo(cv::Scalar(0.));
const int rows = t_prior_.rows;
const int cols = t_prior_.cols;
const int half_window_size = half_window_size_;
const int window_stride = 1;
// sum_{j:patch idx} sum_{i:pixel idx in patch} (t_i - a_j * im_i - b_j)^2
// patches
for (int row_idx = half_window_size; row_idx < rows - half_window_size; row_idx += window_stride) {
for (int col_idx = half_window_size; col_idx < cols - half_window_size; col_idx += window_stride) {
const cv::Point2i j_idx(col_idx, row_idx);
const cv::Vec3f a_j = A_.at<cv::Vec3f>(j_idx);
const float b_j = B_.at<float>(j_idx);
// pixels in a patch
for (int dy = -half_window_size; dy <= half_window_size; ++dy) {
for (int dx = -half_window_size; dx <= half_window_size; ++dx) {
// i, j are defined in the cost function
const cv::Point2i i_idx(col_idx + dx, row_idx + dy);
const float t_i = T_.at<float>(i_idx);
const cv::Vec3f im_i = im_.at<cv::Vec3f>(i_idx);
// For clarity, I keep all the intermidate steps.
// let (t_i - a_j * im_i - b_j) = k
const float dcost_dk = 2 * (t_i - vec3f_dot(a_j, im_i) - b_j);
// grad of k w.r.t variables
const float dk_dti = 1.;
const cv::Vec3f dk_daj = -im_i;
const float dk_dbj = -1.;
// const float weight = 1./((half_window_size * 2 + 1) * (half_window_size * 2 + 1));
const float weight = 0.1;
grad_T_.at<float>(i_idx) += weight * dcost_dk * dk_dti;
grad_A_.at<cv::Vec3f>(j_idx) += weight * dcost_dk * dk_daj;
grad_B_.at<float>(j_idx) += weight * dcost_dk * dk_dbj;
}
}
}
}
// prior cost: weight * sum{i: pixel_index } (t_i - t_prior_i)^2
for (int row_idx = 0; row_idx < rows; ++row_idx) {
for (int col_idx = 0; col_idx < cols; ++col_idx) {
const cv::Point2i pidx(col_idx, row_idx);
const float t_i = T_.at<float>(pidx);
const float t_prior_i = t_prior_.at<float>(pidx);
grad_T_.at<float>(pidx) += prior_weight_ * 2 * (t_i - t_prior_i) * 1.;
}
}
}
void gradient_descent()
{
A_ -= gd_step * grad_A_;
B_ -= gd_step * grad_B_;
T_ -= gd_step * grad_T_;
}
cv::Mat run() override
{
for (int iter = 0; iter < max_iter_; ++iter) {
compute_gradient();
gradient_descent();
...
}
return T_;
}I basically computed the gradient for each pixel and did the gradient descent.
However, it was VERY slow! It took me more than 30 minutes to get a good result.
6.2 Soft Matting Formulation
A Closed Form Solution to Natural Image Matting
Anat Levin, Dani Lischinski, Yair Weiss
This paper proposed another method to do the optimization for the soft matting problem,
\(\text{cost}(x, a, b) = \sum_{j \in I} \sum_{i \in \Omega} ||a_i x_j + b_i – p_j||^2 + ||x – t||^2\).
The idea is to eliminate \((a,b)\) in the quadratic cost function.
Here is a simple refresh of the variable elimination from Convex Optimization.
Or you can read this simple tutorial for it.
To explain the Soft Matting, start with the edge similarity term,
We can rewrite it in matrix form.
Then, we can eliminate (a,b). i.e. expressing the optimal \((a^*,b^*)\) as a function of \(\alpha\). The cost become only a function of \(\alpha\).
The gradient of it is \(2 \bar G_k^T \bar G_k \bar \alpha_k\).
The key why this method is faster is the first part of the gradient \(2 \bar G_k^T \bar G_k\) is constant. So, we can pre-compute it. (The author did the optimization by the conjugate gradient. I’d like to save it to another article.)
Implementation
The analytical expression of \(\bar G_k^T \bar G_k\) is confusing. To make my life easier, we can simply compute the \(\bar G_k^T \bar G_k\) by equation (8).
Then I used the gradient descent to optimize the cost.
The running was okay. It take a few minutes to get a good result.
void compute_G()
{
const int rows = t_prior_.rows;
const int cols = t_prior_.cols;
G_map_ = std::vector<std::vector<cv::Mat>>(rows, std::vector<cv::Mat>(cols, cv::Mat()));
const int half_window_size = half_window_size_;
const int window_size = half_window_size * 2 + 1;
const int num_elements_in_window = window_size * window_size;
const int window_stride = 1;
auto I = cv::Mat::eye(num_elements_in_window, num_elements_in_window, CV_32FC1);
// sum_{j:patch idx} sum_{i:pixel idx in patch} (t_i - a_j * im_i - b_j)^2
// patches
for (int row_idx = half_window_size; row_idx < rows - half_window_size; row_idx += window_stride) {
for (int col_idx = half_window_size; col_idx < cols - half_window_size; col_idx += window_stride) {
cv::Mat G(num_elements_in_window, 4, CV_32FC1, 0.);
// pixels in a patch
int count = 0;
for (int dy = -half_window_size; dy <= half_window_size; ++dy) {
for (int dx = -half_window_size; dx <= half_window_size; ++dx) {
// i, j are defined in the cost function
const cv::Point2i i_idx(col_idx + dx, row_idx + dy);
const cv::Vec3f im_i = im_.at<cv::Vec3f>(i_idx);
G.at<float>(count, 0) = im_i[0];
G.at<float>(count, 1) = im_i[1];
G.at<float>(count, 2) = im_i[2];
G.at<float>(count, 3) = 1.;
count++;
}
}
assert(count == num_elements_in_window);
cv::Mat G_bar = I - G * (G.t() * G + 1e-2 * cv::Mat::eye(4, 4, CV_32FC1)).inv() * G.t();
cv::Mat matting_grad = 2 * G_bar.t() * G_bar;
G_map_[row_idx][col_idx] = matting_grad;
}
}
}
// goal: 1. t = a * target + b
// 2. t close to t_prior
// cost: sum_{j:patch idx} sum_{i:pixel idx in patch} (t_i - a_j * im_i - b_j)^2
// + weight * sum{i: pixel_index } (t_i - t_prior_i)^2
void compute_gradient()
{
grad_T_.setTo(cv::Scalar(0.));
const int rows = t_prior_.rows;
const int cols = t_prior_.cols;
const int half_window_size = half_window_size_;
const int window_size = half_window_size * 2 + 1;
const int num_elements_in_window = window_size * window_size;
const int window_stride = 1;
auto I = cv::Mat::eye(num_elements_in_window, num_elements_in_window, CV_32FC1);
// sum_{j:patch idx} sum_{i:pixel idx in patch} (t_i - a_j * im_i - b_j)^2
// patches
for (int row_idx = half_window_size; row_idx < rows - half_window_size; row_idx += window_stride) {
for (int col_idx = half_window_size; col_idx < cols - half_window_size; col_idx += window_stride) {
cv::Mat t_vec(num_elements_in_window, 1, CV_32FC1, 0.);
// pixels in a patch
int count = 0;
for (int dy = -half_window_size; dy <= half_window_size; ++dy) {
for (int dx = -half_window_size; dx <= half_window_size; ++dx) {
// i, j are defined in the cost function
const cv::Point2i i_idx(col_idx + dx, row_idx + dy);
t_vec.at<float>(count, 0) = T_.at<float>(i_idx);
count++;
}
}
assert(count == num_elements_in_window);
cv::Mat matting_grad = G_map_[row_idx][col_idx] * t_vec;
// pixels in a patch
count = 0;
for (int dy = -half_window_size; dy <= half_window_size; ++dy) {
for (int dx = -half_window_size; dx <= half_window_size; ++dx) {
// i, j are defined in the cost function
const cv::Point2i i_idx(col_idx + dx, row_idx + dy);
grad_T_.at<float>(i_idx) += matting_grad.at<float>(count, 0);
count++;
}
}
assert(count == num_elements_in_window);
}
}
// prior cost: weight * sum{i: pixel_index } (t_i - t_prior_i)^2
for (int row_idx = 0; row_idx < rows; ++row_idx) {
for (int col_idx = 0; col_idx < cols; ++col_idx) {
const cv::Point2i pidx(col_idx, row_idx);
const float t_i = T_.at<float>(pidx);
const float t_prior_i = t_prior_.at<float>(pidx);
grad_T_.at<float>(pidx) += prior_weight_ * 2 * (t_i - t_prior_i);
}
}
}line 39: compute the constant term for the gradient. Note I replace the \(\epsilon\) in the paper by a simple L2 regularization.
line 67-95: compute the gradient. To get \(\alpha\) (or t in this case), we need to iterate the local window.
7. Result
Let me recap what we have done.
- In section 5, we estimated the haze \(A\) by an ad-hoc method.
- In section 3 & 4, we used the dark-channel prior to estimation the weight of the haze model \(t(x)\).
- In section 6, we applied the soft matting to make \(t(x)\) smooth.
Now we can do the last step: remove the haze by \(J(x) = \frac{I(x) – A(1 – t(x))}{t(x)}\).
(Note that the author did a robust version of this equation. Please read the paper for details.)
One small trick is to increase the exposure of the result image.
8. Follow up
Mr. Kaiming hated the Soft Matting! He had something faster.
Guided Image Filtering
Kaiming He, Jian Sun, and Xiaoou Tang
The soft matting optimization is,
It can be a little bit confusing. It’s easier to say the following.
In a local window, the cost is a function of \((p_i, a_k, b_k)\).
\(E(p_i, a_k, b_k)\)
Given \(p_i\), the optimal \((a_k^*, b_k^*)\) is,
\((a_k^*, b_k^*) = \text{argmin}_{a_k, b_k} E(p_i, a_k, b_k)\)
It can be computed by least squares. The results are equation (5) and (6). \((a_k^*, b_k^*)\) is a function of the given \(p_i\).
Since the goal of the soft matting is finding a linear combination of I,
\(p \approx a I + b\)
such that it’s close to a transmission map \(p\) (\(t\) in the previous section). Please go back to section 6 if you are confused.
Given a transmission map \(p\), the optimal (in terms of the soft matting cost) linear combination coefficients \((a^*, b^*)\) are given by (5) and (6). Using \(a^* I + b^*\), we have the optimal soft matting which (1) preserves the edge of I (2) and its intensity is close to \(p\).
The above operation only applies to a local window. A faster but a little bit incorrect way (compare to the global soft matting optimization) to do apply \(a^* I + b^*\) for all windows is to do a weighted average for overlapping results.
It was already fast. But Kaiming made it faster !!!! (I hope you get the joke 🙂 )
The idea is, for (5), (6), and (8), the summations are over a box region. We can pre-compute Integral Images (accumulated sum) and use the Integral Images to compute the summations in (5), (6), and (8).
A leetcode question! Matrix block sum. I guess I should stop writing blogs and work on Leetcode.
Implementation
The code link for my implementation: https://github.com/yimuw/yimu-blog/tree/master/on_the_shoulders/he_kaiming/dehaze
Here is my implementation. I didn’t implement the box filter since my wife needs some attention.
cv::Mat run() override
{
const int rows = t_prior_.rows;
const int cols = t_prior_.cols;
const int half_window_size = half_window_size_;
const int window_size = half_window_size * 2 + 1;
const int num_elements_in_window = window_size * window_size;
const int window_stride = 1;
auto I = cv::Mat::eye(num_elements_in_window, num_elements_in_window, CV_32FC1);
// To handle pixels near the broader.
cv::Mat ab_counter(rows, cols, CV_32FC1, 0.);
cv::Mat a(rows, cols, CV_32FC3, 0.);
cv::Mat b(rows, cols, CV_32FC1, 0.);
// Compute a&b
for (int row_idx = half_window_size; row_idx < rows - half_window_size; row_idx += window_stride) {
for (int col_idx = half_window_size; col_idx < cols - half_window_size; col_idx += window_stride) {
cv::Mat G(num_elements_in_window, 4, CV_32FC1, 0.);
cv::Mat t_vec(num_elements_in_window, 1, CV_32FC1, 0.);
// pixels in a patch
int count = 0;
for (int dy = -half_window_size; dy <= half_window_size; ++dy) {
for (int dx = -half_window_size; dx <= half_window_size; ++dx) {
// i, j are defined in the cost function
const cv::Point2i i_idx(col_idx + dx, row_idx + dy);
const float t_i = t_prior_.at<float>(i_idx);
const cv::Vec3f im_i = im_.at<cv::Vec3f>(i_idx);
// std::cout << "count: " << count << std::endl;
G.at<float>(count, 0) = im_i[0];
G.at<float>(count, 1) = im_i[1];
G.at<float>(count, 2) = im_i[2];
G.at<float>(count, 3) = 1.;
t_vec.at<float>(count, 0) = t_i;
count++;
}
}
assert(count == num_elements_in_window);
cv::Mat ab = (G.t() * G + 1e-2 * cv::Mat::eye(4, 4, CV_32FC1)).inv() * G.t() * t_vec;
for (int dy = -half_window_size; dy <= half_window_size; ++dy) {
for (int dx = -half_window_size; dx <= half_window_size; ++dx) {
const cv::Point2i i_idx(col_idx + dx, row_idx + dy);
a.at<cv::Vec3f>(i_idx)[0] += ab.at<float>(0, 0);
a.at<cv::Vec3f>(i_idx)[1] += ab.at<float>(1, 0);
a.at<cv::Vec3f>(i_idx)[2] += ab.at<float>(2, 0);
b.at<float>(i_idx) += ab.at<float>(3, 0);
ab_counter.at<float>(i_idx) += 1.;
}
}
}
}
cv::Mat filtered(rows, cols, CV_32FC1, 0.);
// apply the filter
for (int row_idx = 0; row_idx < rows; ++row_idx) {
for (int col_idx = 0; col_idx < cols; ++col_idx) {
const cv::Point2i i_idx(col_idx, row_idx);
// Compute the mean
a.at<cv::Vec3f>(i_idx) /= ab_counter.at<float>(i_idx);
b.at<float>(i_idx) /= ab_counter.at<float>(i_idx);
filtered.at<float>(i_idx) = im_.at<cv::Vec3f>(i_idx).dot(a.at<cv::Vec3f>(i_idx)) + b.at<float>(i_idx);
}
}
return filtered;
}I found out this guy had a great implementation. The link. I copied the core part here.
def Guidedfilter(im,p,r,eps):
mean_I = cv2.boxFilter(im,cv2.CV_64F,(r,r));
mean_p = cv2.boxFilter(p, cv2.CV_64F,(r,r));
mean_Ip = cv2.boxFilter(im*p,cv2.CV_64F,(r,r));
cov_Ip = mean_Ip - mean_I*mean_p;
mean_II = cv2.boxFilter(im*im,cv2.CV_64F,(r,r));
var_I = mean_II - mean_I*mean_I;
a = cov_Ip/(var_I + eps);
b = mean_p - a*mean_I;
mean_a = cv2.boxFilter(a,cv2.CV_64F,(r,r));
mean_b = cv2.boxFilter(b,cv2.CV_64F,(r,r));
q = mean_a*im + mean_b;
return q;Line 7: Compute the covariance by \(\text{Var}[x] = E[x^2] + E[x]^2\)
Assuming the cv2.boxFilter is implemented by the Integral Image technique, the guided filter should be very fast.
9 The End
The paper is very interesting to read and to play with.
I am very curious how did Kaiming discover the dark channel prior. There must be struggles. But I guess he simply loves what he is doing and devoted lots of time and effort.
Respect.
10 thoughts on “Haze Removal”