低照图/低亮度视频的修复与优化(ALTM_Retinex实现高HDR)——Java及C实现
项目中遇到问题：低亮度视频与低亮度图片，无法进行后续的处理。
分别测试了传统的变换方法，如：Gamma变换、Log变换、拉普拉斯提升等，效果均不好，
参照大佬的Github：https://github.com/IsaacChanghau/OptimizedImageEnhance
其中选取了ALTMRetinex算法，发现Github主的Java方法省略了源Matlab代码的局部提升及部分公式。
算法原理
源PDF：http://koasas.kaist.ac.kr/bitstream/10203/172985/1/73275.pdf
参照github主：https://github.com/IsaacChanghau/OptimizedImageEnhance/tree/master/matlab/ALTMRetinex
全局提升：取亮度，经过一系列变换，与原亮度相除，有个亮度提升比矩阵。
局部提升：没看，应用为主，局部提升会使一些图片过曝，实际使用中去掉了。
Matlab源码
原版带有局部提升算法
Matlab源码如下：
function outval = ALTM_Retinex(I)
II = im2double(I);
Ir=double(II(:,:,1)); Ig=double(II(:,:,2)); Ib=double(II(:,:,3));
% Global Adaptation
Lw = 0.299 * Ir + 0.587 * Ig + 0.114 * Ib;% input world luminance values
Lwmax = max(max(Lw));% the maximum luminance value
[m, n] = size(Lw);
Lwaver = exp(sum(sum(log(0.001 + Lw))) / (m * n));% log-average luminance
Lg = log(Lw / Lwaver + 1) / log(Lwmax / Lwaver + 1);
% Local Adaptation
kenlRatio = 0.01;
krnlsz = floor(max([3, m * kenlRatio, n * kenlRatio]));
Lg1 = maxfilt2(Lg, [krnlsz, krnlsz]);
Lg1 = imresize(Lg1, [m, n]);
Hg = guidedfilter(Lg, Lg1, 10, 0.01);
eta = 36;
alpha = 1 + eta * Lg / max(max(Lg));
alpha = alpha .* (alpha .^ (1 ./ alpha));
b = max(max(alpha));
a = 1.35;
alpha = 2 * atan(a * alpha / b) / pi * b;
Lgaver = exp(sum(sum(log(0.001 + Lg))) / (m * n));
lambda = 10;
beta = lambda * Lgaver;
Lout = alpha .* log(Lg ./ Hg + beta);
%Lout = normfun(Lout, 1);
Lout = SimplestColorBalance(Lout, 0.005, 0.001, 1);
gain = Lout ./ Lw;
gain(find(Lw == 0)) = 0;
outval = cat(3, gain .* Ir, gain .* Ig, gain .* Ib);
不带有局部提升算法
function outval = ALTM_Retinex2(I)
II = im2double(I);
Ir=double(II(:,:,1)); Ig=double(II(:,:,2)); Ib=double(II(:,:,3));
% Global Adaptation
Lw = 0.299 * Ir + 0.587 * Ig + 0.114 * Ib;% input world luminance values
Lwmax = max(max(Lw));% the maximum luminance value
[m, n] = size(Lw);
Lwaver = exp(sum(sum(log(0.001 + Lw))) / (m * n));% log-average luminance
Lg = log(Lw / Lwaver + 1) / log(Lwmax / Lwaver + 1);
% Global Adaptation
Lw = double(I(:,:,1));
Lw = Lw ./ 255.0;
%Lw = 0.299 * Ir + 0.587 * Ig + 0.114 * Ib;% input world luminance values
Lwmax = max(max(Lw));% the maximum luminance value
[m, n] = size(Lw);
Lwaver = exp(sum(sum(log(0.001 + Lw))) / (m * n));% log-average luminance
Lg = log(Lw / Lwaver + 1) / log(Lwmax / Lwaver + 1);
Lout = SimplestColorBalance(Lg, 0.005, 0.001, 1);
gain = Lout ./ Lw;
gain(find(Lw == 0)) = 0;
outval = cat(3, gain .* Ir, gain .* Ig, gain .* Ib);
RGB颜色提升优化算法
原版代码直接对rgb进行亮度的比例提升，会有色彩溢出的异常，采用公式：
R = (1/2) * (Lg / Lw * (R + Lw) + R - Lw)
都是元素对应位置相乘，GB两通道同理，方法也是参考某博主
function outval = ALTM_Retinex(I)
II = im2double(I);
Ir=double(II(:,:,1)); Ig=double(II(:,:,2)); Ib=double(II(:,:,3));
% Global Adaptation
Lw = 0.299 * Ir + 0.587 * Ig + 0.114 * Ib;% input world luminance values
Lwmax = max(max(Lw));% the maximum luminance value
[m, n] = size(Lw);
Lwaver = exp(sum(sum(log(0.001 + Lw))) / (m * n));% log-average luminance
Lg = log(Lw / Lwaver + 1) / log(Lwmax / Lwaver + 1);
% Local Adaptation
kenlRatio = 0.01;
krnlsz = floor(max([3, m * kenlRatio, n * kenlRatio]));
Lg1 = maxfilt2(Lg, [krnlsz, krnlsz]);
Lg1 = imresize(Lg1, [m, n]);
Hg = guidedfilter(Lg, Lg1, 10, 0.01);
eta = 36;
alpha = 1 + eta * Lg / max(max(Lg));
alpha = alpha .* (alpha .^ (1 ./ alpha));
b = max(max(alpha));
a = 1.35;
alpha = 2 * atan(a * alpha / b) / pi * b;
Lgaver = exp(sum(sum(log(0.001 + Lg))) / (m * n));
lambda = 10;
beta = lambda * Lgaver;
Lout = alpha .* log(Lg ./ Hg + beta);
%Lout = normfun(Lout, 1);
Lout = SimplestColorBalance(Lout, 0.005, 0.001, 1);
gain = Lout ./ Lw;
gain(find(Lw == 0)) = 0;
Ir = (1/2.0) .* (gain .* (Ir+Lw) + Ir - Lw);
Ig = (1/2.0) .* (gain .* (Ig+Lw) + Ig - Lw);
Ib = (1/2.0) .* (gain .* (Ib+Lw) + Ib - Lw);
outval = cat(3, Ir,Ig,Ib);
PS：以上代码.* ./都是对应位置乘除法。
依赖的代码请去Github自行下载，在matlab中的压缩包中。
Java代码
用Opencv4Java实现了Java版本的Matlab版本算法，代码如下：
ALTMRetinex 主类
package cn.uestc.image.function;
import org.opencv.core.Core;
import org.opencv.core.CvType;
import org.opencv.core.Mat;
import org.opencv.core.Point;
import org.opencv.core.Scalar;
import org.opencv.core.Size;
import org.opencv.imgproc.Imgproc;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import java.util.stream.Stream;
import cn.uestc.image.util.Filters;
/**
* 基于YUV的Y通道进行的图像提升
*/
public class ALTMRetinex {
private static final boolean isDebug = false;
private static long time;
// Global Adaptation Parameters -- fixed
private static final double rParam = 0.299;
private static final double gParam = 0.587;
private static final double bParam = 0.114;
public static Mat enhance(Mat src, int r, double eps, double eta, double lambda, double krnlRatio) {
Mat srcClone = src.clone();
srcClone.convertTo(srcClone, CvType.CV_32F);
// 分离通道为bgr
List<Mat> bgr = new ArrayList<>();
// 先归一化到(0-1)
Core.divide(srcClone, new Scalar(255, 255, 255), srcClone);
Core.split(srcClone, bgr);
Mat bChannel = bgr.get(0);
Mat gChannel = bgr.get(1);
Mat rChannel = bgr.get(2);
int m = rChannel.rows();
int n = rChannel.cols();
// Global Adaptation
if (isDebug) {
time = System.currentTimeMillis();
}
List<Mat> list = globalAdaptation(rChannel.clone(), gChannel.clone(), bChannel.clone(), m, n);
if (isDebug) {
System.out.println(String.format("globalAdaptation-耗时:%d", (System.currentTimeMillis() - time)));
time = System.currentTimeMillis();
}
Mat Lw = list.get(0);
Mat Lg = list.get(1);
// Local Adaptation
Mat Hg = localAdaptation(Lg, m, n, r, eps, krnlRatio);
if (isDebug) {
System.out.println(String.format("localAdaptation-耗时:%d", (System.currentTimeMillis() - time)));
time = System.currentTimeMillis();
}
Lg.convertTo(Lg, CvType.CV_32F);
// process
// alpha = 1 + eta * Lg / max(max(Lg));
Mat alpha = new Mat(m, n, rChannel.type());
Core.divide(Lg, new Scalar(Core.minMaxLoc(Lg).maxVal / eta), alpha);
Core.add(alpha, new Scalar(1.0), alpha);
if (isDebug) {
System.out.println(String.format("alpha = 1 + eta * Lg / max(max(Lg));-耗时:%d", (System.currentTimeMillis() - time)));
time = System.currentTimeMillis();
}
// alpha = alpha .* (alpha .^ (1 ./ alpha));
float[] data = new float[alpha.width() * alpha.height()];
alpha.get(0, 0, data);
for (int i = 0; i < data.length; i++) {
// data[i] = (float) (data[i] * (Math.pow(data[i], 1 / data[i])));
data[i] = (float) ((Math.pow(data[i], 1 + 1 / data[i])));
}
if (isDebug) {
System.out.println(String.format("alpha = alpha .* (alpha .^ (1 ./ alpha));-耗时:%d", (System.currentTimeMillis() - time)));
time = System.currentTimeMillis();
}
alpha.put(0, 0, data);
double b = Core.minMaxLoc(alpha).maxVal;
// alpha = 2 * atan(alpha * a / b) / pi * b;
Core.multiply(alpha, new Scalar(1.35 / b), alpha);
alpha.get(0, 0, data);
for (int i = 0; i < data.length; i++) {
// data[i] = (float) (2 * Math.atan(data[i] * 1.35 / b) / Math.PI * b);
data[i] = (float) (Math.atan(data[i]));
}
alpha.put(0, 0, data);
Core.multiply(alpha, new Scalar(2 / Math.PI * b), alpha);
if (isDebug) {
System.out.println(String.format("alpha = 2 * atan(alpha * a / b) / pi * b;-耗时:%d", (System.currentTimeMillis() - time)));
time = System.currentTimeMillis();
}
// beta = exp(sum(sum(log(0.001 + Lg))) / (m * n)) * lambda;
Mat Lg_ = new Mat(m, n, rChannel.type());
Core.add(Lg, new Scalar(0.001), Lg_);
Core.log(Lg_, Lg_);
double beta = Math.exp(Core.sumElems(Lg_).val[0] / (m * n)) * lambda;
if (isDebug) {
System.out.println(String.format("beta = exp(sum(sum(log(0.001 + Lg))) / (m * n)) * lambda;-耗时:%d", (System.currentTimeMillis() - time)));
time = System.currentTimeMillis();
}
// Lout = alpha .* log(Lg ./ Hg + beta);
Mat Lout = new Mat(m, n, rChannel.type());
Core.divide(Lg, Hg, Lout);
Core.add(Lout, new Scalar(beta), Lout);
Core.log(Lout, Lout);
Lout = alpha.mul(Lout);
if (isDebug) {
System.out.println(String.format("Lout = alpha .* log(Lg ./ Hg + beta);-耗时:%d", (System.currentTimeMillis() - time)));
time = System.currentTimeMillis();
}
// Lout = Filters.SimplestColorBalance(Lout, 1);
Lout = Filters.SimplestColorBalance(Lout, 0.005, 0.001);
if (isDebug) {
System.out.println(String.format("SimplestColorBalance-耗时:%d", (System.currentTimeMillis() - time)));
time = System.currentTimeMillis();
}
// gain = Lout ./ Lw;
Mat gain = new Mat();
Core.divide(Lout, Lw, gain);
Core.patchNaNs(gain, 0);
// Mat gain = obtainGain(Lout, Lw, m, n);
if (isDebug) {
System.out.println(String.format("gain = Lout ./ Lw;-耗时:%d", (System.currentTimeMillis() - time)));
time = System.currentTimeMillis();
}
// output
// R'=1/2[(Y'/Y)(R+Y)+R-Y]
Mat temp2 = new Mat(m, n, CvType.CV_32FC1);
Core.add(rChannel, Lw, temp2);
Core.multiply(gain, temp2, temp2);
Core.add(temp2, rChannel, temp2);
Core.subtract(temp2, Lw, temp2);
Core.multiply(temp2, new Scalar(0.5), rChannel);
Core.multiply(rChannel, new Scalar(255), rChannel);
if (isDebug) {
System.out.println(String.format("R'=1/2[(Y'/Y)(R+Y)+R-Y]-耗时:%d", (System.currentTimeMillis() - time)));
time = System.currentTimeMillis();
}
// G'=1/2[(Y'/Y)(G+Y)+G-Y]
Mat temp3 = new Mat(m, n, CvType.CV_32FC1);
Core.add(gChannel, Lw, temp3);
Core.multiply(gain, temp3, temp3);
Core.add(temp3, gChannel, temp3);
Core.subtract(temp3, Lw, temp3);
Core.multiply(temp3, new Scalar(0.5), gChannel);
Core.multiply(gChannel, new Scalar(255), gChannel);
if (isDebug) {
System.out.println(String.format("G'=1/2[(Y'/Y)(G+Y)+G-Y]-耗时:%d", (System.currentTimeMillis() - time)));
time = System.currentTimeMillis();
}
// B'=1/2[(Y'/Y)(B+Y)+B-Y]
Mat temp4 = new Mat(m, n, CvType.CV_32FC1);
Core.add(bChannel, Lw, temp4);
Core.multiply(gain, temp4, temp4);
Core.add(temp4, bChannel, temp4);
Core.subtract(temp4, Lw, temp4);
Core.multiply(temp4, new Scalar(0.5), temp4);
Core.multiply(temp4, new Scalar(255), bChannel);
if (isDebug) {
System.out.println(String.format("B'=1/2[(Y'/Y)(B+Y)+B-Y]-耗时:%d", (System.currentTimeMillis() - time)));
time = System.currentTimeMillis();
}
// merge three color channels to a image
Mat outval = new Mat();
Core.merge(new ArrayList<>(Arrays.asList(bChannel, gChannel, rChannel)), outval);
outval.convertTo(outval, CvType.CV_8U);
if (isDebug) {
System.out.println(String.format("merge-耗时:%d", (System.currentTimeMillis() - time)));
time = System.currentTimeMillis();
}
// 及时释放
temp2.release();
temp3.release();
temp4.release();
rChannel.release();
gChannel.release();
bChannel.release();
Lw.release();
Lg.release();
Hg.release();
alpha.release();
Lg_.release();
Lout.release();
gain.release();
return outval;
}
public static Mat enhanceWithoutLocal(Mat src) {
if (isDebug) {
time = System.currentTimeMillis();
}
Mat srcClone = src.clone();
srcClone.convertTo(srcClone, CvType.CV_32F);
// 分离通道为bgr
List<Mat> bgr = new ArrayList<>();
// 先归一化到(0-1)
Core.divide(srcClone, new Scalar(255, 255, 255), srcClone);
if (isDebug) {
System.out.println(String.format("归一化-耗时:%d", (System.currentTimeMillis() - time)));
time = System.currentTimeMillis();
}
Core.split(srcClone, bgr);
Mat bChannel = bgr.get(0);
Mat gChannel = bgr.get(1);
Mat rChannel = bgr.get(2);
int m = rChannel.rows();
int n = rChannel.cols();
if (isDebug) {
System.out.println(String.format("分解通道-耗时:%d", (System.currentTimeMillis() - time)));
time = System.currentTimeMillis();
}
// Global Adaptation
if (isDebug) {
time = System.currentTimeMillis();
}
List<Mat> list = globalAdaptation(rChannel.clone(), gChannel.clone(), bChannel.clone(), m, n);
if (isDebug) {
System.out.println(String.format("globalAdaptation-耗时:%d", (System.currentTimeMillis() - time)));
time = System.currentTimeMillis();
}
Mat Lw = list.get(0);
Mat Lg = list.get(1);
// gain = Lout ./ Lw;
Mat gain = new Mat();
Core.divide(Lg, Lw, gain);
Core.patchNaNs(gain, 0);
if (isDebug) {
System.out.println(String.format("gain = Lg ./ Lw;-耗时:%d", (System.currentTimeMillis() - time)));
time = System.currentTimeMillis();
}
// output
// R'=1/2[(Y'/Y)(R+Y)+R-Y]
Mat temp2 = new Mat(m, n, CvType.CV_32FC1);
Core.add(rChannel, Lw, temp2);
Core.multiply(gain, temp2, temp2);
Core.add(temp2, rChannel, temp2);
Core.subtract(temp2, Lw, temp2);
Core.multiply(temp2, new Scalar(0.5), rChannel);
Core.multiply(rChannel, new Scalar(255), rChannel);
if (isDebug) {
System.out.println(String.format("R'=1/2[(Y'/Y)(R+Y)+R-Y]-耗时:%d", (System.currentTimeMillis() - time)));
time = System.currentTimeMillis();
}
// G'=1/2[(Y'/Y)(G+Y)+G-Y]
Mat temp3 = new Mat(m, n, CvType.CV_32FC1);
Core.add(gChannel, Lw, temp3);
Core.multiply(gain, temp3, temp3);
Core.add(temp3, gChannel, temp3);
Core.subtract(temp3, Lw, temp3);
Core.multiply(temp3, new Scalar(0.5), gChannel);
Core.multiply(gChannel, new Scalar(255), gChannel);
if (isDebug) {
System.out.println(String.format("G'=1/2[(Y'/Y)(G+Y)+G-Y]-耗时:%d", (System.currentTimeMillis() - time)));
time = System.currentTimeMillis();
}
// B'=1/2[(Y'/Y)(B+Y)+B-Y]
Mat temp4 = new Mat(m, n, CvType.CV_32FC1);
Core.add(bChannel, Lw, temp4);
Core.multiply(gain, temp4, temp4);
Core.add(temp4, bChannel, temp4);
Core.subtract(temp4, Lw, temp4);
Core.multiply(temp4, new Scalar(0.5), temp4);
Core.multiply(temp4, new Scalar(255), bChannel);
if (isDebug) {
System.out.println(String.format("B'=1/2[(Y'/Y)(B+Y)+B-Y]-耗时:%d", (System.currentTimeMillis() - time)));
time = System.currentTimeMillis();
}
// merge three color channels to a image
Mat outval = new Mat();
Core.merge(new ArrayList<>(Arrays.asList(bChannel, gChannel, rChannel)), outval);
outval.convertTo(outval, CvType.CV_8U);
if (isDebug) {
System.out.println(String.format("merge-耗时:%d", (System.currentTimeMillis() - time)));
time = System.currentTimeMillis();
}
// 及时释放
temp2.release();
temp3.release();
temp4.release();
rChannel.release();
gChannel.release();
bChannel.release();
Lw.release();
Lg.release();
gain.release();
srcClone.release();
if (isDebug) {
System.out.println(String.format("释放资源-耗时:%d", (System.currentTimeMillis() - time)));
}
return outval;
}
private static Mat localAdaptation(Mat Lg, int rows, int cols, int r, double eps, double krnlRatio) {
int krnlSz = Stream.of(3.0, rows * krnlRatio, cols * krnlRatio).max(Double::compare).orElse(3.0).intValue();
// maximum filter: using dilate to extract the local maximum of a image block, which acts as the maximum filter
// Meanwhile, minimum filter can be achieved by using erode function
Mat Lg_ = new Mat();
Mat kernel = Imgproc.getStructuringElement(Imgproc.MORPH_RECT, new Size(krnlSz, krnlSz), new Point(-1, -1));
Imgproc.dilate(Lg, Lg_, kernel);
// guided image filter
Mat guidedImage = Filters.GuidedImageFilter(Lg, Lg_, r, eps);
// 及时释放
Lg_.release();
kernel.release();
return guidedImage;
}
private static List<Mat> globalAdaptation(Mat r, Mat g, Mat b, int rows, int cols) {
// Calculate Lw(初始亮度矩阵) & maximum of Lw(亮度的最大值)
// 下面的LW是通过参数计算的灰度图,本意是YUV的Y通道
// Lw = 0.299 * Ir + 0.587 * Ig + 0.114 * Ib
Mat Lw = new Mat(rows, cols, r.type());
Core.multiply(r, new Scalar(rParam), r);
Core.multiply(g, new Scalar(gParam), g);
Core.multiply(b, new Scalar(bParam), b);
Core.add(r, g, Lw);
Core.add(Lw, b, Lw);
// the maximum luminance value
// 亮度的最大值
double LwMax = Core.minMaxLoc(Lw).maxVal;
// Calculate log-average luminance and get global adaptation result
// 计算log平均亮度
// Lw_average(log平均亮度) = exp(sum(sum(log(0.001 + Lw))) / (m * n));
Mat Lw_ = Lw.clone();
Core.add(Lw_, new Scalar(0.001), Lw_);
Core.log(Lw_, Lw_);
double LwAver = Math.exp(Core.sumElems(Lw_).val[0] / (rows * cols));
// 并获取新的亮度矩阵(Lg)
// Lg = log(Lw / Lw_aver + 1) / log(Lw_max / Lw_aver + 1);
Mat Lg = Lw.clone();
Core.divide(Lg, new Scalar(LwAver), Lg);
Core.add(Lg, new Scalar(1.0), Lg);
Core.log(Lg, Lg);
Core.divide(Lg, new Scalar(Math.log(LwMax / LwAver + 1.0)), Lg); // Lg is the global adaptation
List<Mat> list = new ArrayList<>();
list.add(Lw);
list.add(Lg);
// 及时释放
b.release();
g.release();
r.release();
Lw_.release();
return list;
}
private static Mat obtainGain(Mat Lout, Mat Lw, int rows, int cols) {
Mat gain = new Mat(rows, cols, Lout.type());
for (int i = 0; i < rows; i++) {
for (int j = 0; j < cols; j++) {
if (Lw.get(i, j)[0] == 0) gain.put(i, j, Lout.get(i, j)[0]);
else gain.put(i, j, Lout.get(i, j)[0] / Lw.get(i, j)[0]);
}
}
return gain;
}
@SuppressWarnings("unused")
private static Mat adjustment(Mat alpha, double a) {
double b = Core.minMaxLoc(alpha).maxVal;
int rows = alpha.rows();
int cols = alpha.cols();
for (int i = 0; i < rows; i++) {
for (int j = 0; j < cols; j++) {
//double val = alpha.get(i, j)[0];
alpha.put(i, j, (2 * Math.atan(a * alpha.get(i, j)[0] / b) / Math.PI * b));
}
}
return alpha;
}
}
滤波类Filters
package cn.uestc.image.util;
import org.opencv.core.Core;
import org.opencv.core.CvType;
import org.opencv.core.Mat;
import org.opencv.core.Point;
import org.opencv.core.Scalar;
import org.opencv.core.Size;
import org.opencv.imgproc.Imgproc;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
public class Filters {
/**
* Simplest Color Balance. Performs color balancing via histogram
* normalization.
*
* @param img input color or gray scale image
* @param percent controls the percentage of pixels to clip to white and black. (normally, choose 1~10)
* @return Balanced image in CvType.CV_32F
*/
public static Mat SimplestColorBalance(Mat img, int percent) {
if (percent <= 0)
percent = 5;
img.convertTo(img, CvType.CV_32F);
List<Mat> channels = new ArrayList<>();
int rows = img.rows(); // number of rows of image
int cols = img.cols(); // number of columns of image
int chnls = img.channels(); // number of channels of image
double halfPercent = percent / 200.0;
if (chnls == 3) Core.split(img, channels);
else channels.add(img);
List<Mat> results = new ArrayList<>();
for (int i = 0; i < chnls; i++) {
// find the low and high precentile values (based on the input percentile)
Mat flat = new Mat();
channels.get(i).reshape(1, 1).copyTo(flat);
Core.sort(flat, flat, Core.SORT_ASCENDING);
double lowVal = flat.get(0, (int) Math.floor(flat.cols() * halfPercent))[0];
double topVal = flat.get(0, (int) Math.ceil(flat.cols() * (1.0 - halfPercent)))[0];
// saturate below the low percentile and above the high percentile
// 优化成数组操作
Mat channel = channels.get(i);
float[] data = new float[rows * cols];
channel.get(0, 0, data);
for (int j = 0; j < data.length; j++) {
if (data[j] < lowVal) data[j] = (float) lowVal;
if (data[j] > topVal) data[j] = (float) topVal;
}
channel.put(0, 0, data);
// for (int m = 0; m < rows; m++) {
// for (int n = 0; n < cols; n++) {
// if (channel.get(m, n)[0] < lowVal) channel.put(m, n, lowVal);
// if (channel.get(m, n)[0] > topVal) channel.put(m, n, topVal);
// }
// }
Core.normalize(channel, channel, 0.0, 255.0, Core.NORM_MINMAX);
channel.convertTo(channel, CvType.CV_32F);
results.add(channel);
}
Mat outval = new Mat();
Core.merge(results, outval);
return outval;
}
public static Mat SimplestColorBalance(Mat img, double bottom, double top) {
img.convertTo(img, CvType.CV_32F);
List<Mat> channels = new ArrayList<>();
int rows = img.rows(); // number of rows of image
int cols = img.cols(); // number of columns of image
int chnls = img.channels(); // number of channels of image
if (chnls == 3) Core.split(img, channels);
else channels.add(img);
List<Mat> results = new ArrayList<>();
for (int i = 0; i < chnls; i++) {
// find the low and high precentile values (based on the input percentile)
Mat flat = new Mat();
channels.get(i).reshape(1, 1).copyTo(flat);
Core.sort(flat, flat, Core.SORT_ASCENDING);
int lowIndex = (int) Math.floor(flat.cols() * bottom);
double lowVal = flat.get(0, lowIndex)[0];
int topIndex = (int) Math.ceil(flat.cols() * (1.0 - top));
double topVal = flat.get(0, topIndex)[0];
// saturate below the low percentile and above the high percentile
// 优化成数组操作
Mat channel = channels.get(i);
float[] data = new float[rows * cols];
channel.get(0, 0, data);
for (int j = 0; j < data.length; j++) {
if (data[j] < lowVal) {
data[j] = (float) lowVal;
}
if (data[j] > topVal) {
data[j] = (float) topVal;
}
}
channel.put(0, 0, data);
data = null;
Core.normalize(channel, channel, 1, 0, Core.NORM_MINMAX);
channel.convertTo(channel, CvType.CV_32F);
results.add(channel);
flat.release();
}
Mat outval = new Mat();
Core.merge(results, outval);
for (Mat m : channels) {
m.release();
}
for(Mat m : results){
m.release();
}
return outval;
}
/**
* Guided Image Filter for grayscale image, O(1) time implementation of guided filter
*
* @param I guidance image (should be a gray-scale/single channel image)
* @param p filtering input image (should be a gray-scale/single channel image)
* @param r local window radius
* @param eps regularization parameter
* @return filtered image
*/
public static Mat GuidedImageFilter(Mat I, Mat p, int r, double eps) {
I.convertTo(I, CvType.CV_64FC1);
p.convertTo(p, CvType.CV_64FC1);
//[hei, wid] = size(I);
int rows = I.rows();
int cols = I.cols();
// N = boxfilter(ones(hei, wid), r); % the size of each local patch; N=(2r+1)^2 except for boundary pixels.
Mat N = new Mat();
Mat temp1 = Mat.ones(rows, cols, I.type());
Imgproc.boxFilter(temp1, N, -1, new Size(r * 2 + 1, r * 2 + 1), new Point(-1, -1), false, Core.BORDER_CONSTANT);
temp1.release();
// mean_I = boxfilter(I, r) ./ N;
Mat mean_I = new Mat();
Imgproc.boxFilter(I, mean_I, -1, new Size(r * 2 + 1, r * 2 + 1), new Point(-1, -1), false, Core.BORDER_CONSTANT);
Core.divide(mean_I, N, mean_I);
// mean_p = boxfilter(p, r) ./ N
Mat mean_p = new Mat();
Imgproc.boxFilter(p, mean_p, -1, new Size(r * 2 + 1, r * 2 + 1), new Point(-1, -1), false, Core.BORDER_CONSTANT);
Core.divide(mean_p, N, mean_p);
// mean_Ip = boxfilter(I.*p, r) ./ N;
Mat mean_Ip = new Mat();
Mat temp2 = I.mul(p);
Imgproc.boxFilter(temp2, mean_Ip, -1, new Size(r * 2 + 1, r * 2 + 1), new Point(-1, -1), false, Core.BORDER_CONSTANT);
temp2.release();
Core.divide(mean_Ip, N, mean_Ip);
// cov_Ip = mean_Ip - mean_I .* mean_p; % this is the covariance of (I, p) in each local patch.
Mat cov_Ip = new Mat();
Mat temp3 = mean_I.mul(mean_p);
Core.subtract(mean_Ip, temp3, cov_Ip);
temp3.release();
// mean_II = boxfilter(I.*I, r) ./ N;
Mat mean_II = new Mat();
Mat temp4 = I.mul(I);
Imgproc.boxFilter(temp4, mean_II, -1, new Size(r * 2 + 1, r * 2 + 1), new Point(-1, -1), false, Core.BORDER_CONSTANT);
temp4.release();
Core.divide(mean_II, N, mean_II);
// var_I = mean_II - mean_I .* mean_I;
Mat var_I = new Mat();
Mat temp5 = mean_I.mul(mean_I);
Core.subtract(mean_II, temp5, var_I);
temp5.release();
// a = cov_Ip ./ (var_I + eps); % Eqn. (5) in the paper;
Mat a = new Mat();
Core.add(var_I, new Scalar(eps), a);
Core.divide(cov_Ip, a, a);
//b = mean_p - a .* mean_I; % Eqn. (6) in the paper;
Mat b = new Mat();
Mat temp6 = a.mul(mean_I);
Core.subtract(mean_p, temp6, b);
temp6.release();
// mean_a = boxfilter(a, r) ./ N;
Mat mean_a = new Mat();
Imgproc.boxFilter(a, mean_a, -1, new Size(r * 2 + 1, r * 2 + 1), new Point(-1, -1), false, Core.BORDER_CONSTANT);
Core.divide(mean_a, N, mean_a);
// mean_b = boxfilter(b, r) ./ N;
Mat mean_b = new Mat();
Imgproc.boxFilter(b, mean_b, -1, new Size(r * 2 + 1, r * 2 + 1), new Point(-1, -1), false, Core.BORDER_CONSTANT);
Core.divide(mean_b, N, mean_b);
// q = mean_a .* I + mean_b; % Eqn. (8) in the paper;
Mat q = new Mat();
Mat temp7 = mean_a.mul(I);
Core.add(temp7, mean_b, q);
temp7.release();
q.convertTo(q, CvType.CV_32F);
// 及时释放
N.release();
mean_I.release();
mean_p.release();
mean_Ip.release();
cov_Ip.release();
mean_II.release();
var_I.release();
a.release();
b.release();
mean_a.release();
mean_b.release();
return q;
}
/**
* Guided Image Filter for Color Image
*
* @param origI guidance image (should be a gray-scale/single channel image)
* @param p filtering input image (should be a gray-scale/single channel image)
* @param r local window radius
* @param eps regularization parameter
* @param s blocks number, s * s
* @param depth image depth
* @return filtered image
*/
@SuppressWarnings("unused")
public static Mat GuidedImageFilter_Color(Mat origI, Mat p, int r, double eps, double s, int depth) {
//Pre-defined parameters
ArrayList<Mat> Ichannels = new ArrayList<>();
ArrayList<Mat> Isubchannels = new ArrayList<>();
int Idepth;
double r_sub;
Mat mean_I_r;
Mat mean_I_g;
Mat mean_I_b;
Mat invrr = new Mat();
Mat invrg = new Mat();
Mat invrb = new Mat();
Mat invgg = new Mat();
Mat invgb = new Mat();
Mat invbb = new Mat();
// Process
Mat I;
if (origI.depth() == CvType.CV_32F || origI.depth() == CvType.CV_64F) {
I = origI.clone();
} else {
I = convertTo(origI, CvType.CV_32F);
}
Idepth = I.depth();
Core.split(I, Ichannels);
Mat I_sub = new Mat();
Imgproc.resize(I, I_sub, new Size(I.cols() / s, I.rows() / s), 0.0, 0.0, Imgproc.INTER_NEAREST);
Core.split(I_sub, Isubchannels);
r_sub = r / s;
mean_I_r = boxfilter(Isubchannels.get(0), (int) r_sub);
mean_I_g = boxfilter(Isubchannels.get(1), (int) r_sub);
mean_I_b = boxfilter(Isubchannels.get(2), (int) r_sub);
// variance of I in each local patch: the matrix Sigma in Eqn (14).
// Note the variance in each local patch is a 3x3 symmetric matrix:
// rr, rg, rb
// Sigma = rg, gg, gb
// rb, gb, bb
Mat var_I_rr = new Mat();
Mat var_I_rg = new Mat();
Mat var_I_rb = new Mat();
Mat var_I_gg = new Mat();
Mat var_I_gb = new Mat();
Mat var_I_bb = new Mat();
Mat temp1 = new Mat();
Core.subtract(boxfilter(Isubchannels.get(0).mul(Isubchannels.get(0)), (int) r_sub),
mean_I_r.mul(mean_I_r), temp1);
Core.add(temp1, new Scalar(eps), var_I_rr);
Core.subtract(boxfilter(Isubchannels.get(0).mul(Isubchannels.get(1)), (int) r_sub),
mean_I_r.mul(mean_I_g), var_I_rg);
Core.subtract(boxfilter(Isubchannels.get(0).mul(Isubchannels.get(2)), (int) r_sub),
mean_I_r.mul(mean_I_b), var_I_rb);
Core.subtract(boxfilter(Isubchannels.get(1).mul(Isubchannels.get(1)), (int) r_sub),
mean_I_g.mul(mean_I_g), temp1);
Core.add(temp1, new Scalar(eps), var_I_gg);
Core.subtract(boxfilter(Isubchannels.get(1).mul(Isubchannels.get(2)), (int) r_sub),
mean_I_g.mul(mean_I_b), var_I_gb);
Core.subtract(boxfilter(Isubchannels.get(2).mul(Isubchannels.get(2)), (int) r_sub),
mean_I_b.mul(mean_I_b), temp1);
Core.add(temp1, new Scalar(eps), var_I_bb);
// Inverse of Sigma + eps * I
Core.subtract(var_I_gg.mul(var_I_bb), var_I_gb.mul(var_I_gb), invrr);
Core.subtract(var_I_gb.mul(var_I_rb), var_I_rg.mul(var_I_bb), invrg);
Core.subtract(var_I_rg.mul(var_I_gb), var_I_gg.mul(var_I_rb), invrb);
Core.subtract(var_I_rr.mul(var_I_bb), var_I_rb.mul(var_I_rb), invgg);
Core.subtract(var_I_rb.mul(var_I_rg), var_I_rr.mul(var_I_gb), invgb);
Core.subtract(var_I_rr.mul(var_I_gg), var_I_rg.mul(var_I_rg), invbb);
Mat covDet = new Mat();
Core.add(invrr.mul(var_I_rr), invrg.mul(var_I_rg), temp1);
Core.add(temp1, invrb.mul(var_I_rb), covDet);
Core.divide(invrr, covDet, invrr);
Core.divide(invrg, covDet, invrg);
Core.divide(invrb, covDet, invrb);
Core.divide(invgg, covDet, invgg);
Core.divide(invgb, covDet, invgb);
Core.divide(invbb, covDet, invbb);
Mat p2 = convertTo(p, Idepth);
Mat result = new Mat();
if (p.channels() == 1) {
result = filterSingleChannel(p2, s, Isubchannels, Ichannels, mean_I_r, mean_I_g, mean_I_b, invrr, invrg,
invrb, invgg, invgb, invbb, r_sub);
} else {
ArrayList<Mat> pc = new ArrayList<>();
Core.split(p2, pc);
for (int i = 0; i < pc.size(); i++) {
pc.set(i, filterSingleChannel(pc.get(i), s, Isubchannels, Ichannels, mean_I_r, mean_I_g, mean_I_b, invrr,
invrg, invrb, invgg, invgb, invbb, r_sub));
}
Core.merge(pc, result);
}
return convertTo(result, depth == -1 ? p.depth() : depth);
}
private static Mat boxfilter(Mat I, int r) {
Mat result = new Mat();
Imgproc.blur(I, result, new Size(r, r));
return result;
}
private static Mat convertTo(Mat mat, int depth) {
if (mat.depth() == depth) {
return mat;
}
Mat result = new Mat();
mat.convertTo(result, depth);
return result;
}
private static Mat filterSingleChannel(Mat p, double s, ArrayList<Mat> Isubchannels, ArrayList<Mat> Ichannels,
Mat mean_I_r, Mat mean_I_g, Mat mean_I_b, Mat invrr, Mat invrg, Mat invrb, Mat invgg, Mat invgb,
Mat invbb, double r_sub) {
Mat p_sub = new Mat();
Imgproc.resize(p, p_sub, new Size(p.cols() / s, p.rows() / s), 0.0, 0.0, Imgproc.INTER_NEAREST);
Mat mean_p = boxfilter(p_sub, (int) r_sub);
Mat mean_Ip_r = boxfilter(Isubchannels.get(0).mul(p_sub), (int) r_sub);
Mat mean_Ip_g = boxfilter(Isubchannels.get(1).mul(p_sub), (int) r_sub);
Mat mean_Ip_b = boxfilter(Isubchannels.get(2).mul(p_sub), (int) r_sub);
// convariance of (I, p) in each local patch
Mat cov_Ip_r = new Mat();
Mat cov_Ip_g = new Mat();
Mat cov_Ip_b = new Mat();
Core.subtract(mean_Ip_r, mean_I_r.mul(mean_p), cov_Ip_r);
Core.subtract(mean_Ip_g, mean_I_g.mul(mean_p), cov_Ip_g);
Core.subtract(mean_Ip_b, mean_I_b.mul(mean_p), cov_Ip_b);
Mat temp1 = new Mat();
Mat a_r = new Mat();
Mat a_g = new Mat();
Mat a_b = new Mat();
Core.add(invrr.mul(cov_Ip_r), invrg.mul(cov_Ip_g), temp1);
Core.add(temp1, invrb.mul(cov_Ip_b), a_r);
Core.add(invrg.mul(cov_Ip_r), invgg.mul(cov_Ip_g), temp1);
Core.add(temp1, invgb.mul(cov_Ip_b), a_g);
Core.add(invrb.mul(cov_Ip_r), invgb.mul(cov_Ip_g), temp1);
Core.add(temp1, invbb.mul(cov_Ip_b), a_b);
Mat b = new Mat();
Core.subtract(mean_p, a_r.mul(mean_I_r), b);
Core.subtract(b, a_g.mul(mean_I_g), b);
Core.subtract(b, a_b.mul(mean_I_b), b);
Mat mean_a_r = boxfilter(a_r, (int) r_sub);
Mat mean_a_g = boxfilter(a_g, (int) r_sub);
Mat mean_a_b = boxfilter(a_b, (int) r_sub);
Mat mean_b = boxfilter(b, (int) r_sub);
Imgproc.resize(mean_a_r, mean_a_r,
new Size(Ichannels.get(0).cols(), Ichannels.get(0).rows()), 0.0, 0.0, Imgproc.INTER_LINEAR);
Imgproc.resize(mean_a_g, mean_a_g,
new Size(Ichannels.get(0).cols(), Ichannels.get(0).rows()), 0.0, 0.0, Imgproc.INTER_LINEAR);
Imgproc.resize(mean_a_b, mean_a_b,
new Size(Ichannels.get(0).cols(), Ichannels.get(0).rows()), 0.0, 0.0, Imgproc.INTER_LINEAR);
Imgproc.resize(mean_b, mean_b,
new Size(Ichannels.get(0).cols(), Ichannels.get(0).rows()), 0.0, 0.0, Imgproc.INTER_LINEAR);
Mat result = new Mat();
Core.add(mean_a_r.mul(Ichannels.get(0)), mean_a_g.mul(Ichannels.get(1)), temp1);
Core.add(temp1, mean_a_b.mul(Ichannels.get(2)), temp1);
Core.add(temp1, mean_b, result);
return result;
}
}
以上代码均可以在Opencv4Java下运行，速度是一方面，大概一帧需要0.5秒，处理视频可以采用多线程并发，相互帧不影响。
相对原算法，去掉了局部提升、优化了RGB的提升方式，实现了完整的Java版本的代码。
C版本
可移植DLL或SO，供Windows或Linux或Android使用，只实现了全局提升。
//
// Created by ShenHongBin on 2019/3/14.
//
#ifndef SMILETEETH_IMAGEENHANCE_H
#define SMILETEETH_IMAGEENHANCE_H
#include <stdio.h>
#include <opencv2/opencv.hpp>
using namespace cv;
using namespace std;
class ImageEnhanceHelper {
public:
ImageEnhanceHelper(int rows, int cols) {
this->cols = cols;
this->rows = rows;
}
Mat enhance(Mat src, int r, double eps, double eta, double lambda, double krnlRatio);
Mat enhanceWithoutLocal(Mat src);
vector<Mat> globalAdaptation(Mat r, Mat g, Mat b);
Mat localAdaptation(Mat Lg, int rows, int cols, int r, double eps, double krnlRatio);
void close() {
}
private:
int rows;
int cols;
const int r = 10;
const double eps = 0.01;
const double eta = 36.0;
const double lambda = 10.0;
const double krnlRatio = 0.01;
};
Mat ImageEnhanceHelper::enhance(Mat src, int r, double eps, double eta, double lambda,
double krnlRatio) {
return Mat();
}
Mat ImageEnhanceHelper::enhanceWithoutLocal(Mat src) {
Mat srcClone, gain, outMat, lw, lg;
vector<Mat> bgr;
srcClone = src.clone();
srcClone.convertTo(srcClone, CV_32F);
srcClone = srcClone / 255.0f;
split(srcClone, bgr);
float a = bgr[0].at<float>(0,0);
a = bgr[1].at<float>(0,0);
a = bgr[2].at<float>(0,0);
vector<Mat> vectorMat = globalAdaptation(bgr[0], bgr[1], bgr[2]);
lw = vectorMat[0];
lg = vectorMat[1];
gain = vectorMat[1] / vectorMat[0];
patchNaNs(gain, 0);
a = gain.at<float>(0,0);
a = lw.at<float>(0,0);
a = lg.at<float>(0,0);
// R'=1/2[(Y'/Y)(R+Y)+R-Y]
Mat temp;
// cvtColor(lw, lw, COLOR_GRAY2BGR);
// cvtColor(gain, gain, COLOR_GRAY2BGR);
//
//
// add(outMat, lw, temp);
// multiply(temp, gain, temp);
// add(temp, outMat, temp);
// subtract(temp, lw, temp);
// outMat = temp.mul(Scalar(0.5));
temp = bgr[0] + lw;
multiply(temp, gain, temp);
bgr[0] = (0.5f) * (temp + bgr[0] - lw);
temp = bgr[1] + lw;
multiply(temp, gain, temp);
bgr[1] = (0.5f) * (temp + bgr[1] - lw);
temp = bgr[2] + lw;
multiply(temp, gain, temp);
bgr[2] = (0.5f) * (temp + bgr[2] - lw);
merge(bgr, outMat);
multiply(outMat, Scalar(255, 255, 255), outMat);
outMat.convertTo(outMat, CV_8U);
// 释放资源
lw.release();
lg.release();
gain.release();
temp.release();
bgr[0].release();
bgr[1].release();
bgr[2].release();
srcClone.release();
return outMat;
}
vector<Mat> ImageEnhanceHelper::globalAdaptation(Mat r, Mat g, Mat b) {
vector<Mat> vectorMat;
Mat lw, lwClone, lg;
lw = 0.299 * r + 0.587 * g + 0.114 * b;
double min, max;
Point2i minLoc, maxLoc;
minMaxLoc(lw, &min, &max, &minLoc, &maxLoc);
// Lw_average(log平均亮度) = exp(sum(sum(log(0.001 + Lw))) / (m * n));
double lwAver;
lwClone = lw + 0.001;
log(lwClone, lwClone);
lwAver = exp(sum(sum(lwClone))[0] / (rows * cols));
// Lg = log(Lw / Lw_aver + 1) / log(Lw_max / Lw_aver + 1);
lg = lw / lwAver + 1;
log(lg, lg);
lg = lg / log(max / lwAver + 1);
vectorMat.push_back(lw);
vectorMat.push_back(lg);
// 释放资源
lwClone.release();
return vectorMat;
}
Mat ImageEnhanceHelper::localAdaptation(Mat Lg, int rows, int cols, int r, double eps,
double krnlRatio) {
return Mat();
}
#endif //SMILETEETH_IMAGEENHANCE_H