Commit 1e1364ba authored by Baptiste Delos's avatar Baptiste Delos
Browse files

A CPU version of the Uniform condition based on the first 9 SH coefficients

parent 240d7e40
......@@ -3,158 +3,170 @@
#include <map>
#include <functional>
namespace Preprocessing
{
cv::Mat matFromOptiXBuffer(optix::Buffer &buffer)
namespace CvOptixInterop
{
RTsize w, h;
buffer->getSize(w, h);
buffer->getElementSize();
RTformat format = buffer->getFormat();
void *data = buffer->map(0, RT_BUFFER_MAP_READ);
cv::Mat mat(h, w, CV_32FC4);
for (unsigned int i = 0; i < w * h; ++i)
cv::Mat matFromOptiXBuffer(optix::Buffer &buffer)
{
float *pixel = &(static_cast<float *>(data)[4 * i]);
mat.at<cv::Vec4f>(i) = cv::Vec4f(pixel[0], pixel[1], pixel[2], pixel[3]);
}
RTsize w, h;
buffer->getSize(w, h);
buffer->getElementSize();
RTformat format = buffer->getFormat();
buffer->unmap();
void *data = buffer->map(0, RT_BUFFER_MAP_READ);
cv::Mat mat(h, w, CV_32FC4);
return mat;
}
for (unsigned int i = 0; i < w * h; ++i)
{
float *pixel = &(static_cast<float *>(data)[4 * i]);
mat.at<cv::Vec4f>(i) = cv::Vec4f(pixel[0], pixel[1], pixel[2], pixel[3]);
}
void setOptiXBufferFromMat(optix::Buffer &buffer, cv::Mat &mat)
{
RTsize w, h;
buffer->getSize(w, h);
buffer->unmap();
void *data = buffer->map(0, RT_BUFFER_MAP_WRITE);
return mat;
}
for (unsigned int i = 0; i < w * h; ++i)
void setOptiXBufferFromMat(optix::Buffer &buffer, cv::Mat &mat)
{
float *pixel = &(static_cast<float *>(data)[4 * i]);
RTsize w, h;
buffer->getSize(w, h);
cv::Vec4f mat_pixel = mat.at<cv::Vec4f>(i);
pixel[0] = mat_pixel[0];
pixel[1] = mat_pixel[1];
pixel[2] = mat_pixel[2];
pixel[3] = mat_pixel[3];
}
void *data = buffer->map(0, RT_BUFFER_MAP_WRITE);
buffer->unmap();
}
for (unsigned int i = 0; i < w * h; ++i)
{
float *pixel = &(static_cast<float *>(data)[4 * i]);
template<>
std::vector<float> matToVector(cv::Mat &mat)
{
std::vector<float> elements;
cv::Vec4f mat_pixel = mat.at<cv::Vec4f>(i);
pixel[0] = mat_pixel[0];
pixel[1] = mat_pixel[1];
pixel[2] = mat_pixel[2];
pixel[3] = mat_pixel[3];
}
for (unsigned int i = 0; i < mat.rows * mat.cols; ++i)
{
elements.emplace_back(mat.at<float>(i));
buffer->unmap();
}
return elements;
}
}// namespace cvOptixInterop
void quantify(cv::Mat &mat, int nbits)
namespace Preprocessing
{
float nvalues = float(1 << nbits) - 1.0f;
for (unsigned int i = 0; i < mat.rows * mat.cols; ++i)
template<>
std::vector<float> matToVector(cv::Mat &mat)
{
cv::Vec4f &pixel = mat.at<cv::Vec4f>(i);
pixel[0] = std::round(pixel[0] * nvalues);
pixel[1] = std::round(pixel[1] * nvalues);
pixel[2] = std::round(pixel[2] * nvalues);
std::vector<float> elements;
pixel /= nvalues;
}
}
cv::Mat cdfHist(cv::Mat &hist/*, const float **ranges*/)
{
cv::Mat sum;
hist.copyTo(sum);
for (unsigned int i = 0; i < mat.rows * mat.cols; ++i)
{
elements.emplace_back(mat.at<float>(i));
}
for (unsigned int i = 1; i < hist.rows; ++i)
{
sum.row(i) += sum.row(i - 1);
return elements;
}
for (unsigned i = 0; i < hist.cols; ++i)
void quantify(cv::Mat &mat, int nbits)
{
double scale_factor = /*double(ranges[i][1] - ranges[i][0] - 1.f)*/1.0 / double(sum.at<float>(hist.rows - 1, i));
sum.col(i) *= scale_factor;
}
float nvalues = float(1 << nbits) - 1.0f;
return sum;
}
std::vector<cv::Mat> cdf(cv::Mat const &mat, std::vector<int> channels, std::vector<const float *> ranges, float alpha, cv::ColorConversionCodes code)
{
cv::Mat cvt_mat;
mat.copyTo(cvt_mat);
for (unsigned int i = 0; i < mat.rows * mat.cols; ++i)
{
cv::Vec4f &pixel = mat.at<cv::Vec4f>(i);
pixel[0] = std::round(pixel[0] * nvalues);
pixel[1] = std::round(pixel[1] * nvalues);
pixel[2] = std::round(pixel[2] * nvalues);
if (mat.type() == CV_32FC1 || mat.type() == CV_32FC2 || mat.type() == CV_32FC3 || mat.type() == CV_32FC4)
{
cvt_mat.convertTo(cvt_mat, CV_32FC4, double(alpha));
pixel /= nvalues;
}
}
if (code != cv::COLOR_COLORCVT_MAX)
cv::Mat cdfHist(cv::Mat &hist)
{
cvt_mat.convertTo(cvt_mat, code);
}
cv::Mat sum;
hist.copyTo(sum);
std::vector<int> hist_sizes;
std::vector<const float*> hist_ranges;
for (auto bounds : ranges)
{
hist_sizes.emplace_back(static_cast<int>(bounds[1] - bounds[0]));
for (unsigned int i = 1; i < hist.rows; ++i)
{
sum.row(i) += sum.row(i - 1);
}
for (unsigned i = 0; i < hist.cols; ++i)
{
double scale_factor = 1.0 / double(sum.at<float>(hist.rows - 1, i));
sum.col(i) *= scale_factor;
}
return sum;
}
std::vector<cv::Mat> hists;
for (int c : channels)
cv::Mat cdfs(cv::Mat const &mat, std::vector<int> channels, std::vector<const float *> ranges, float alpha, cv::ColorConversionCodes code)
{
cv::Mat hist;
cv::calcHist(&cvt_mat, 1, &c, cv::Mat(), hist, 1, hist_sizes.data(), ranges.data());
hists.emplace_back(hist);
cv::Mat cvt_mat;
mat.copyTo(cvt_mat);
if (mat.depth() == CV_32F)
{
cvt_mat.convertTo(cvt_mat, cvt_mat.type(), double(alpha));
}
if (code != cv::COLOR_COLORCVT_MAX)
{
cvt_mat.convertTo(cvt_mat, code);
}
std::vector<int> hist_sizes;
std::vector<const float*> hist_ranges;
for (auto bounds : ranges)
{
hist_sizes.emplace_back(static_cast<int>(bounds[1] - bounds[0]));
}
std::vector<cv::Mat> channel_cdfs;
for (int c : channels)
{
cv::Mat hist;
cv::calcHist(&cvt_mat, 1, &c, cv::Mat(), hist, 1, hist_sizes.data(), ranges.data());
channel_cdfs.emplace_back(cdfHist(hist));
}
cv::Mat cdfs;
cv::merge(channel_cdfs.data(), channel_cdfs.size(), cdfs);
return cdfs;
}
std::vector<cv::Mat> cdf;
for (cv::Mat h : hists)
void equalize(const cv::Mat &mat_in, cv::Mat &mat_out, std::vector<int> channel_indices, int cvt_code_in, int cvt_code_out)
{
cdf.emplace_back(cdfHist(h));
}
std::vector<cv::Mat> channels;
return cdf;
}
mat_in.convertTo(mat_out, CV_8U, 255.);
void equalize(cv::Mat &mat)
{
std::vector<cv::Mat> channels;
cv::Mat mat_uchar;
mat.convertTo(mat_uchar, CV_8U, 255.);
cv::cvtColor(mat_uchar, mat_uchar, cv::COLOR_RGBA2RGB);
cv::cvtColor(mat_uchar, mat_uchar, cv::COLOR_RGB2HLS);
if (cvt_code_in != cv::COLOR_COLORCVT_MAX)
{
cv::cvtColor(mat_out, mat_out, cvt_code_in);
}
else
{
cv::copyTo(mat_out, mat_out, cv::Mat());
}
cv::split(mat_out, channels);
cv::split(mat_uchar, channels);
for (int index : channel_indices)
{
cv::equalizeHist(channels[index], channels[index]);
}
// cv::Mat channels_equalized[4];
cv::equalizeHist(channels[1], channels[1]);
// cv::equalizeHist(channels[1], channels_equalized[1]);
// cv::equalizeHist(channels[2], channels_equalized[2]);
// channels_equalized[3] = channels[3];
cv::merge(channels.data(), channels.size(), mat_out);
cv::merge(channels.data(), channels.size(), mat_uchar);
if (cvt_code_out != cv::COLOR_COLORCVT_MAX)
{
cv::cvtColor(mat_out, mat_out, cvt_code_out);
}
cv::cvtColor(mat_uchar, mat_uchar, cv::COLOR_HLS2RGB);
cv::cvtColor(mat_uchar, mat_uchar, cv::COLOR_RGB2RGBA);
mat_uchar.convertTo(mat, mat.type(), 1./255.);
}
mat_out.convertTo(mat_out, mat_in.type(), 1./255.);
}
}
}// namespace Preprocessing
#ifndef IMAGE_PROCESSING_H
#define IMAGE_PROCESSING_H
#ifndef PREPROCESSING_H
#define PREPROCESSING_H
#include "optix_backend/optix_util.hpp"
......@@ -12,102 +12,172 @@
#include <functional>
namespace Preprocessing
namespace CvOptixInterop
{
cv::Mat matFromOptiXBuffer(optix::Buffer &buffer);
void setOptiXBufferFromMat(optix::Buffer &buffer, cv::Mat &mat);
}
namespace Preprocessing
{
template<typename T>
std::vector<T> matToVector(cv::Mat &mat);
void quantify(cv::Mat &mat, int nbits);
void equalize(cv::Mat &mat);
void equalize(const cv::Mat &mat_in, cv::Mat &mat_out, std::vector<int> channel_indices, int cvt_code_in=cv::COLOR_COLORCVT_MAX, int cvt_code_out=cv::COLOR_COLORCVT_MAX);
cv::Mat cdfHist(cv::Mat &hist/*, const float **ranges*/);
std::vector<cv::Mat> cdf(cv::Mat const &mat, std::vector<int> channels, std::vector<const float *> ranges, float alpha, cv::ColorConversionCodes code=cv::COLOR_COLORCVT_MAX);
static std::map<std::pair<int, int>, std::function<float(cv::Vec3f)>> base_functions = {
{ std::make_pair<int, int>(0, 0), [](cv::Vec3f position)->float { (void) position; return 0.282095f; } },
{ std::make_pair<int, int>(1, 1), [](cv::Vec3f position)->float { return 0.488603f * position[0]; } },
{ std::make_pair<int, int>(1, 0), [](cv::Vec3f position)->float { return 0.488603f * position[2]; } },
{ std::make_pair<int, int>(1,-1), [](cv::Vec3f position)->float { return 0.488603f * position[1]; } },
{ std::make_pair<int, int>(2, 1), [](cv::Vec3f position)->float { return 1.092548f * position[0] * position[2]; } },
{ std::make_pair<int, int>(2,-1), [](cv::Vec3f position)->float { return 1.092548f * position[1] * position[2]; } },
{ std::make_pair<int, int>(2,-2), [](cv::Vec3f position)->float { return 1.092548f * position[0] * position[1]; } },
{ std::make_pair<int, int>(2, 0), [](cv::Vec3f position)->float { return 0.315392f * (3.f * position[2] * position[2] - 1.f); } },
{ std::make_pair<int, int>(2, 2), [](cv::Vec3f position)->float { return 0.546274f * (position[0] * position[0] - position[1] * position[1]); } },
};
template<typename InT, typename OutT>
OutT sphericalHarmonicMoment(cv::Mat &mat, int l, int m)
{
const int type = mat.type();
cv::Mat cdfs(cv::Mat const &mat, std::vector<int> channels, std::vector<const float *> ranges, float alpha, cv::ColorConversionCodes code=cv::COLOR_COLORCVT_MAX);
OutT mat_sh_moments(0.0);
auto baseFunc = base_functions[std::pair<int, int>(l, m)];
double d_theta = M_PI * (1. / double(mat.rows));
double d_phi = 2. * M_PI * (1. / double(mat.cols));
namespace SphericHarmonics
{
for (int i = 0; i < mat.rows; ++i)
static std::map<std::pair<int, int>, std::function<float(cv::Vec3f)>> base_functions = {
{ std::make_pair<int, int>(0, 0), [](cv::Vec3f position)->float { (void) position; return 0.282095f; } },
{ std::make_pair<int, int>(1, 1), [](cv::Vec3f position)->float { return 0.488603f * position[0]; } },
{ std::make_pair<int, int>(1, 0), [](cv::Vec3f position)->float { return 0.488603f * position[2]; } },
{ std::make_pair<int, int>(1,-1), [](cv::Vec3f position)->float { return 0.488603f * position[1]; } },
{ std::make_pair<int, int>(2, 1), [](cv::Vec3f position)->float { return 1.092548f * position[0] * position[2]; } },
{ std::make_pair<int, int>(2,-1), [](cv::Vec3f position)->float { return 1.092548f * position[1] * position[2]; } },
{ std::make_pair<int, int>(2,-2), [](cv::Vec3f position)->float { return 1.092548f * position[0] * position[1]; } },
{ std::make_pair<int, int>(2, 0), [](cv::Vec3f position)->float { return 0.315392f * (3.f * position[2] * position[2] - 1.f); } },
{ std::make_pair<int, int>(2, 2), [](cv::Vec3f position)->float { return 0.546274f * (position[0] * position[0] - position[1] * position[1]); } },
};
/**
* @brief Returns the spherical harmonic coefficient of the input matrix at order l and degree m
*
* Template type specifies the element type of the matrix
**/
template<typename T>
T SHCoefficient(const cv::Mat &mat, int l, int m)
{
for (int j = 0; j < mat.cols; ++j)
const int type = mat.type();
T sh_coeff(0.0);
auto baseFunc = base_functions[std::pair<int, int>(l, m)];
double d_theta = M_PI * (1. / double(mat.rows));
double d_phi = 2. * M_PI * (1. / double(mat.cols));
for (int i = 0; i < mat.rows; ++i)
{
double theta = M_PI * double(i) / double(mat.rows);
double phi = 2. * M_PI * double(j) / double(mat.cols);
for (int j = 0; j < mat.cols; ++j)
{
double theta = M_PI * double(i) / double(mat.rows);
double phi = 2. * M_PI * double(j) / double(mat.cols);
cv::Vec3f position(std::sin(theta) * std::cos(phi), std::sin(theta) * std::sin(phi), std::cos(theta));
cv::Vec3f position(std::sin(theta) * std::cos(phi), std::sin(theta) * std::sin(phi), std::cos(theta));
InT channel_sh_moments = mat.at<InT>(i, j) * baseFunc(position) * std::sin(theta) * d_theta * d_phi;
mat_sh_moments += OutT(channel_sh_moments.val);
T channel_sh_moments = mat.at<T>(i, j) * baseFunc(position) * std::sin(theta) * d_theta * d_phi;
sh_coeff += T(channel_sh_moments);
}
}
return sh_coeff;
}
return mat_sh_moments;
}
/**
* @brief Returns the spherical harmonic coefficients on the first orders whose maximal index is specified in parameter
*
* Template type specifies the element type of the matrix
**/
template<typename T>
std::map<std::pair<int, int>, T> SHCoefficients(const cv::Mat &mat, int max_order=2)
{
std::map<std::pair<int, int>, T> sh_coeffs;
template<typename InT, typename OutT>
std::map<std::pair<int, int>, OutT> sphericalHarmonicComponents(cv::Mat &mat, int max_order=2)
{
std::map<std::pair<int, int>, OutT> moments;
for (int l = 0; l <= max_order; ++l)
{
for (int m = -l; m <= l; ++m)
{
sh_coeffs[std::pair<int, int>(l, m)] = SHCoefficient<T>(mat, l, m);
}
}
return sh_coeffs;
}
for (int l = 0; l <= max_order; ++l)
/**
* @brief Reconstructs the input matrix from the given spherical harmonic moments
*
* The given matrix must be allocated
* Template type specifies the element type of the matrix
**/
template<typename T>
void matFromSHCoefficients(cv::Mat &mat, std::map<std::pair<int, int>, T> &sh_moments)
{
for (int m = -l; m <= l; ++m)
for (int i = 0; i < mat.rows; ++i)
{
moments[std::pair<int, int>(l, m)] = sphericalHarmonicMoment<InT, OutT>(mat, l, m);
for (int j = 0; j < mat.cols; ++j)
{
double theta = M_PI * double(i) / double(mat.rows);
double phi = 2. * M_PI * (double(j) / double(mat.cols));
cv::Vec3f position = cv::Vec3f(std::sin(theta) * std::cos(phi), std::sin(theta) * std::sin(phi), std::cos(theta));
T &element = mat.at<T>(i,j);
element = T(0.f);
for (auto moment : sh_moments)
{
T inner_prod = moment.second * base_functions[moment.first](position);
element += inner_prod;
}
}
}
}
return moments;
}
template<typename T>
void matFromSHComponents(cv::Mat &mat, std::map<std::pair<int, int>, T> &sh_moments)
{
for (int i = 0; i < mat.rows; ++i)
template<typename T>
std::map<std::pair<int, int>, T> resetSphericalHarmonicPowers(float nrows, float ncols, std::map<std::pair<int, int>, T> moments, std::map<std::pair<int, int>, T> original_moments, int max_order=2)
{
for (int j = 0; j < mat.cols; ++j)
std::map<std::pair<int, int>, T> new_moments;
std::vector<T> order_powers(max_order + 1, T(0.0f));
std::vector<T> order_original_powers(max_order + 1, T(0.0f));
// Compute original and current powers at each order
for (int l = 0; l <= max_order; ++l)
{
double theta = M_PI * double(i) / double(mat.rows);
double phi = 2. * M_PI * (double(j) / double(mat.cols));
for (int m = -l; m <= l; ++m)
{
auto moment_indices = std::pair<int, int>(l, m);
for (int i = 0; i < nrows; ++i)
{
for (int j = 0; j < ncols; ++j)
{
double theta = M_PI * double(i) / double(nrows);
double phi = 2. * M_PI * double(j) / double(ncols);
cv::Vec3f position(std::sin(theta) * std::cos(phi), std::sin(theta) * std::sin(phi), std::cos(theta));
order_original_powers[l] += original_moments[moment_indices] * base_functions[moment_indices](position);
order_powers[l] += moments[moment_indices] * base_functions[moment_indices](position);
}
}
}
}
cv::Vec3f position = cv::Vec3f(std::sin(theta) * std::cos(phi), std::sin(theta) * std::sin(phi), std::cos(theta));
T &element = mat.at<T>(i,j);
element = T(0.f, 0.f, 0.f, element[3]);
// Reset moment powers
for (int l = 0; l <= max_order; ++l)
{
T output_power;
cv::divide(order_original_powers[l], order_powers[l], output_power);
for (auto moment : sh_moments)
for (int m = -l; m <= l; ++m)
{
cv::Vec3f inner_prod = moment.second * base_functions[moment.first](position);
element[0] += inner_prod[0];
element[1] += inner_prod[1];
element[2] += inner_prod[2];
cv::multiply(moments[std::pair<int, int>(l, m)], output_power, new_moments[std::pair<int, int>(l, m)]);
}
}
return new_moments;
}
}
}// namespace SphericHarmonics
}
#endif // IMAGE_PROCESSING_H
#endif // PREPROCESSING_H
......@@ -22,7 +22,7 @@ using mrf::uint;
#include <chrono>
#include <sys/time.h>
#include "imnodes.h"
//#include "imnodes.h"
double getTime()
......@@ -162,14 +162,19 @@ void LightingGUI::initializeUI(int w, int h)
sphericalHarmonics->addFilter(sh_recomp);
sphericalHarmonics->setCurrentFilter(0);
_filter_managers.add(sphericalHarmonics);
FilterManager *equalHist = new FilterManager();
equalHist->init(envmap_width, envmap_height);
equalHist->addFilter(hist_equal);
equalHist->setCurrentFilter(0);
_filter_managers.add(equalHist);
_current_filter_manager = 0;
}
void LightingGUI::run()
{
_envmap_image = Preprocessing::matFromOptiXBuffer(_renderer->_envmap_buffer);
_sh_components = Preprocessing::sphericalHarmonicComponents<cv::Vec4f, cv::Vec3f>(_envmap_image);
_envmap_image = CvOptixInterop::matFromOptiXBuffer(_renderer->_envmap_buffer);
_sh_components = Preprocessing::SphericHarmonics::SHCoefficients<cv::Vec4f>(_envmap_image);
// Bind environment map buffer to a buffer object for OpenGL to use