#include "precomp.hpp"
#include "opencv2/photo.hpp"
#include "opencv2/imgproc.hpp"
-//#include "opencv2/highgui.hpp"
#include "hdr_common.hpp"
namespace cv
samples(_samples),
lambda(_lambda),
random(_random),
- w(tringleWeights())
+ w(triangleWeights())
{
}
{
CV_INSTRUMENT_REGION()
+ // check inputs
std::vector<Mat> images;
src.getMatVector(images);
Mat times = _times.getMat();
CV_Assert(images.size() == times.total());
checkImageDimensions(images);
CV_Assert(images[0].depth() == CV_8U);
+ CV_Assert(times.type() == CV_32FC1);
+ // create output
int channels = images[0].channels();
int CV_32FCC = CV_MAKETYPE(CV_32F, channels);
+ int rows = images[0].rows;
+ int cols = images[0].cols;
dst.create(LDR_SIZE, 1, CV_32FCC);
Mat result = dst.getMat();
- std::vector<Point> sample_points;
+ // pick pixel locations (either random or in a rectangular grid)
+ std::vector<Point> points;
+ points.reserve(samples);
if(random) {
for(int i = 0; i < samples; i++) {
- sample_points.push_back(Point(rand() % images[0].cols, rand() % images[0].rows));
+ points.push_back(Point(rand() % cols, rand() % rows));
}
} else {
- int x_points = static_cast<int>(sqrt(static_cast<double>(samples) * images[0].cols / images[0].rows));
+ int x_points = static_cast<int>(sqrt(static_cast<double>(samples) * cols / rows));
+ CV_Assert(0 < x_points && x_points <= cols);
int y_points = samples / x_points;
- int step_x = images[0].cols / x_points;
- int step_y = images[0].rows / y_points;
-
+ CV_Assert(0 < y_points && y_points <= rows);
+ int step_x = cols / x_points;
+ int step_y = rows / y_points;
for(int i = 0, x = step_x / 2; i < x_points; i++, x += step_x) {
for(int j = 0, y = step_y / 2; j < y_points; j++, y += step_y) {
- if( 0 <= x && x < images[0].cols &&
- 0 <= y && y < images[0].rows )
- sample_points.push_back(Point(x, y));
+ if( 0 <= x && x < cols && 0 <= y && y < rows ) {
+ points.push_back(Point(x, y));
+ }
}
}
+ // we can have slightly less grid points than specified
+ //samples = static_cast<int>(points.size());
}
+ // we need enough equations to ensure a sufficiently overdetermined system
+ // (maybe only as a warning)
+ //CV_Assert(points.size() * (images.size() - 1) >= LDR_SIZE);
+
+ // solve for imaging system response function, over each channel separately
std::vector<Mat> result_split(channels);
- for(int channel = 0; channel < channels; channel++) {
- Mat A = Mat::zeros((int)sample_points.size() * (int)images.size() + LDR_SIZE + 1, LDR_SIZE + (int)sample_points.size(), CV_32F);
+ for(int ch = 0; ch < channels; ch++) {
+ // initialize system of linear equations
+ Mat A = Mat::zeros((int)points.size() * (int)images.size() + LDR_SIZE + 1,
+ LDR_SIZE + (int)points.size(), CV_32F);
Mat B = Mat::zeros(A.rows, 1, CV_32F);
- int eq = 0;
- for(size_t i = 0; i < sample_points.size(); i++) {
+ // include the data−fitting equations
+ int k = 0;
+ for(size_t i = 0; i < points.size(); i++) {
for(size_t j = 0; j < images.size(); j++) {
-
- int val = images[j].ptr()[3*(sample_points[i].y * images[j].cols + sample_points[i].x) + channel];
- A.at<float>(eq, val) = w.at<float>(val);
- A.at<float>(eq, LDR_SIZE + (int)i) = -w.at<float>(val);
- B.at<float>(eq, 0) = w.at<float>(val) * log(times.at<float>((int)j));
- eq++;
+ // val = images[j].at<Vec3b>(points[i].y, points[i].x)[ch]
+ int val = images[j].ptr()[channels*(points[i].y * cols + points[i].x) + ch];
+ float wij = w.at<float>(val);
+ A.at<float>(k, val) = wij;
+ A.at<float>(k, LDR_SIZE + (int)i) = -wij;
+ B.at<float>(k, 0) = wij * log(times.at<float>((int)j));
+ k++;
}
}
- A.at<float>(eq, LDR_SIZE / 2) = 1;
- eq++;
-
- for(int i = 0; i < 254; i++) {
- A.at<float>(eq, i) = lambda * w.at<float>(i + 1);
- A.at<float>(eq, i + 1) = -2 * lambda * w.at<float>(i + 1);
- A.at<float>(eq, i + 2) = lambda * w.at<float>(i + 1);
- eq++;
+
+ // fix the curve by setting its middle value to 0
+ A.at<float>(k, LDR_SIZE / 2) = 1;
+ k++;
+
+ // include the smoothness equations
+ for(int i = 0; i < (LDR_SIZE - 2); i++) {
+ float wi = w.at<float>(i + 1);
+ A.at<float>(k, i) = lambda * wi;
+ A.at<float>(k, i + 1) = -2 * lambda * wi;
+ A.at<float>(k, i + 2) = lambda * wi;
+ k++;
}
+
+ // solve the overdetermined system using SVD (least-squares problem)
Mat solution;
solve(A, B, solution, DECOMP_SVD);
- solution.rowRange(0, LDR_SIZE).copyTo(result_split[channel]);
+ solution.rowRange(0, LDR_SIZE).copyTo(result_split[ch]);
}
+
+ // combine log-exposures and take its exponent
merge(result_split, result);
exp(result, result);
}
}
protected:
- String name;
- int samples;
- float lambda;
- bool random;
- Mat w;
+ String name; // calibration algorithm identifier
+ int samples; // number of pixel locations to sample
+ float lambda; // constant that determines the amount of smoothness
+ bool random; // whether to sample locations randomly or in a grid shape
+ Mat w; // weighting function for corresponding pixel values
};
Ptr<CalibrateDebevec> createCalibrateDebevec(int samples, float lambda, bool random)