2.12.2021, 9:00 - 11:00: Due to updates GitLab may be unavailable for some minutes between 09:00 and 11:00.

Commit d309ea25 authored by Declara Denis's avatar Declara Denis Committed by Christian Schulte zu Berge
Browse files

Cleaned up solver code

parent abca9934
#include "cm_cusp.h"
#include "cm_cusp_cuda_exports.h"
CuspConfidenceMapSolver::CuspConfidenceMapSolver(int width, int height)
: width(width), height(height)
, gpuData(CUSP_CM_createGpuData(width, height)) { }
CuspConfidenceMapSolver::~CuspConfidenceMapSolver()
{
CUSP_CM_destroyGpuData(gpuData);
}
void CuspConfidenceMapSolver::createSystem(const unsigned char* image, int width, int height,
float gradientScaling,
float alpha, float beta, float gamma)
{
if (width != this->width || height != this->height) {
CUSP_CM_resizeGpuData(*gpuData, width, height);
}
this->width = width;
this->height = height;
CUSP_CM_buildEquationSystem(*gpuData, image, width, height, alpha, beta, gamma, gradientScaling);
}
void CuspConfidenceMapSolver::setInitialSolution(const std::vector<float> &val)
{
CUSP_CM_setInitialSolution(*gpuData, val);
}
void CuspConfidenceMapSolver::solve(int iterations)
{
CUSP_CM_uploadSystem(*gpuData);
CUSP_CM_solveSystem(*gpuData, iterations, 1e-20);
CUSP_CM_downloadSolution(*gpuData);
}
const float* CuspConfidenceMapSolver::getSolution(int &width, int &height) const
{
const float *ptr = CUSP_CM_getSolutionPtr(*gpuData);
width = width; height = height; // FIXME: height and width are not actually checked against the size of the buffer...
return ptr;
}
#include <vector>
#include <cmath>
#include <cusp/dia_matrix.h>
#include <cusp/krylov/cg.h>
#include <cusp/precond/diagonal.h>
#include "cm_cusp_cuda_exports.h"
typedef cusp::dia_matrix<int, float, cusp::host_memory> CuspHostDiaMatrix;
typedef cusp::dia_matrix<int, float, cusp::device_memory> CuspDeviceDiaMatrix;
typedef cusp::array1d<float, cusp::host_memory> CuspHostVector;
typedef cusp::array1d<float, cusp::device_memory> CuspDeviceVector;
class CuspGPUData {
public:
CuspGPUData(int width, int height)
: L_h(width*height, width*height, width*height*9, 9),
b_h(width*height), x_h(width*height),
L_d(width*height, width*height, width*height*9, 9),
b_d(width*height), x_d(width*height) {};
void resize(int width, int height)
{
this->L_h.resize(width*height, width*height, width*height*9, 9);
this->b_h.resize(width*height);
this->x_h.resize(width*height);
}
CuspHostDiaMatrix L_h;
CuspHostVector b_h;
CuspHostVector x_h;
CuspDeviceDiaMatrix L_d;
CuspDeviceVector b_d;
CuspDeviceVector x_d;
};
CuspGPUData* CUSP_CM_createGpuData(int width, int height)
{
CuspGPUData *data = new CuspGPUData(width, height);
return data;
}
void CUSP_CM_destroyGpuData(CuspGPUData *data)
{
delete data;
}
void CUSP_CM_resizeGpuData(CuspGPUData& data, int width, int height)
{
data.resize(width, height);
}
void CUSP_CM_setInitialSolution(CuspGPUData& data, const std::vector<float>& values)
{
for (int i = 0; i < data.x_h.size() && i < values.size(); ++i) {
data.x_h[i] = values[i];
}
}
void CUSP_CM_uploadSystem(CuspGPUData &data)
{
data.x_d = data.x_h;
data.b_d = data.b_h;
data.L_d = data.L_h;
}
void CUSP_CM_downloadSolution(CuspGPUData &data)
{
data.x_h = data.x_d;
// FIXME: Should not be needed... (Clamp solution)
for (int i = 0; i < data.x_h.size(); ++i) {
data.x_h[i] = min(1.0f, data.x_h[i]);
}
}
void CUSP_CM_solveSystem(CuspGPUData& data, int iterations, float precision)
{
cusp::default_monitor<float> monitor(data.b_d, iterations, precision);
//cusp::precond::bridson_ainv<float, cusp::device_memory> M(A, 1);
cusp::precond::diagonal<float, cusp::device_memory> M(data.L_d);
//cusp::krylov::bicgstab(A, x, b, monitor);
cusp::krylov::cg(data.L_d, data.x_d, data.b_d, monitor, M);
}
const float* CUSP_CM_getSolutionPtr(const CuspGPUData& data)
{
//const_cast<CuspGPUData&>(data).x_h = const_cast<CuspGPUData&>(data).x_d;
return &data.x_h[0];
}
struct ComputeLaplacianData
{
float alpha, beta, gamma;
float gradientScaling;
const unsigned char *image;
int width, height;
int centralDiagonal;
int offsets[9];
float gammaList[9];
std::vector<float> attenuationLUT;
};
static inline float _getWeight(const ComputeLaplacianData &data, int x, int y, int diagonal, bool isUpsideDown)
{
const unsigned char *image = data.image;
int idx1 = y * data.width + x;
int idx2 = idx1 + data.offsets[diagonal];
float attenuation1 = data.attenuationLUT[idx1 / data.width];
float attenuation2 = data.attenuationLUT[idx2 / data.width];
float gradient = abs(image[idx1]*attenuation1/255.0f - image[idx2]*attenuation2/255.0f) * data.gradientScaling;
float weight = exp(-data.beta * (gradient + data.gammaList[diagonal]));
return weight + 1e-4;
}
void CUSP_CM_buildEquationSystem(CuspGPUData &gpudata, const unsigned char* image, int width, int height,
float alpha, float beta, float gamma,
float gradientScaling,
bool isUpsideDown)
{
// Gather all of the options together
ComputeLaplacianData data;
data.alpha = alpha; data.beta = beta; data.gamma = gamma;
data.gradientScaling = gradientScaling;
data.image = image;
data.width = width; data.height = height;
data.centralDiagonal = 4;
int offsets[9] = {-width-1, -width, -width+1, -1, 0, 1, width-1, width, width+1};
float gammaList[9] = {sqrt(2.0f)*gamma, 0.0f, sqrt(2.0f)*gamma, gamma, 0.0f, gamma, sqrt(2.0f)*gamma, 0.0f, sqrt(2.0f)*gamma};
for (int i = 0; i < 9; ++i) {
data.offsets[i] = offsets[i];
data.gammaList[i] = gammaList[i];
}
// Prepare equation system data structure
for (int i = 0; i < 9; ++i) {
gpudata.L_h.diagonal_offsets[i] = data.offsets[i];
}
// Precompute attenuation tables
data.attenuationLUT = std::vector<float>(height);
for (int i = 0; i < height; ++i) {
float y = (float)i / (float)(height-1);
if (isUpsideDown) y = 1 - y;
data.attenuationLUT[i] = 1 - exp(-alpha * y);
}
// Initialize B
for (int x = 0; x < width * height; ++x) {
if (x < width)
gpudata.b_h[x] = isUpsideDown ? 0.0f : 1.0f;
else if (x >= (width*(height-1)))
gpudata.b_h[x] = isUpsideDown ? 1.0f : 0.0f;
else
gpudata.b_h[x] = 0.0f;
}
// Fill in the rest of the matrix
for (int y = 0; y < height; ++y) {
for (int x = 0; x < width; ++x) {
int idx = y * width + x;
// Filter off out-of-bounds edges
unsigned short filter = 495; // 111 101 111
// 8 - neighbourhood filter
if (x == 0) filter &= 203; // 011 001 011
if (x == width-1) filter &= 422; // 110 100 110
if (y == 0) filter &= 47; // 000 101 111
if (y == height-1) filter &= 488; // 111 101 000
// 4 - neighbourhood filter
//if (x == 0) filter &= 138; // 010 001 010
//if (x == width-1) filter &= 162; // 010 100 010
//if (y == 0) filter &= 42; // 000 101 010
//if (y == height-1) filter &= 168; // 010 101 000
float valueSum = 0.0f;
if (y == 0 || y == height - 1) valueSum = 1.0f;
for (int d = 0; d < 9; ++d) {
gpudata.L_h.values(idx, d) = 0;
float value = 0.0f;
if (((256>>d) & filter) != 0) {
value = _getWeight(data, x, y, d, isUpsideDown);
gpudata.L_h.values(idx, d) = -value;
}
valueSum += value;
}
gpudata.L_h.values(idx, data.centralDiagonal) = valueSum;
}
}
}
\ No newline at end of file
#pragma once
#include <vector>
class CuspGPUData;
class CuspConfidenceMapSolver
{
public:
CuspConfidenceMapSolver(int width, int height);
virtual ~CuspConfidenceMapSolver();
virtual void createSystem(const unsigned char* image, int width, int height,
float gradientScaling, float alpha, float beta, float gamma);
virtual void setInitialSolution(const std::vector<float> &val);
virtual void solve(int iterations);
virtual const float* getSolution(int &width, int &height) const;
// Input image data
int width, height;
private:
// Matrices and Vectors
CuspGPUData* gpuData;
};
\ No newline at end of file
#pragma once
#include <vector>
#include <GL/gl.h>
class CuspGPUData;
CuspGPUData* CUSP_CM_createGpuData(int width, int height);
void CUSP_CM_destroyGpuData(CuspGPUData *data);
void CUSP_CM_resizeGpuData(CuspGPUData& data, int width, int height);
void CUSP_CM_buildEquationSystem(CuspGPUData& data, const unsigned char* image, int width, int height,
float alpha, float beta, float gamma,
float gradientScaling,
bool isUpsideDown = true);
void CUSP_CM_setInitialSolution(CuspGPUData& data, const std::vector<float>& values);
void CUSP_CM_uploadSystem(CuspGPUData &data);
void CUSP_CM_downloadSolution(CuspGPUData &data);
void CUSP_CM_solveSystem(CuspGPUData& data, int iterations, float precision);
const float* CUSP_CM_getSolutionPtr(const CuspGPUData& data);
void CUSP_CM_copyTex(GLuint tex_in, GLuint tex_out, int width, int height);
void CUSP_CM_initCuda();
\ No newline at end of file
#include "cudaconfidencemaps_cuda.h"
#include <cusp/dia_matrix.h>
#include <cusp/krylov/cg.h>
#include <cusp/precond/diagonal.h>
typedef cusp::dia_matrix<int, float, cusp::host_memory> CuspHostDiaMatrix;
typedef cusp::dia_matrix<int, float, cusp::device_memory> CuspDeviceDiaMatrix;
typedef cusp::array1d<float, cusp::host_memory> CuspHostVector;
typedef cusp::array1d<float, cusp::device_memory> CuspDeviceVector;
namespace campvis {
namespace cuda {
struct CudaConfidenceMapsSystemGPUData {
// Data on the host (only solution vector x and matrix L)
CuspHostDiaMatrix L_h;
CuspHostVector x_h;
// Data on the device, a.k.a. GPU (Lx = b)
CuspDeviceDiaMatrix L_d;
CuspDeviceVector b_d;
CuspDeviceVector x_d;
// Information about the system
bool isUpsideDown;
int imageWidth;
int imageHeight;
float solutionResidualNorm;
};
struct ComputeLaplacianData
{
float alpha, beta, gamma;
float gradientScaling;
const unsigned char *image;
int width, height;
int centralDiagonal;
int offsets[9];
float gammaList[9];
std::vector<float> attenuationLUT;
};
CudaConfidenceMapsSystemSolver::CudaConfidenceMapsSystemSolver()
: _gpuData(new CudaConfidenceMapsSystemGPUData())
{
_gpuData->imageWidth = 0;
_gpuData->imageHeight = 0;
_gpuData->solutionResidualNorm = 0.0f;
}
static inline float _getWeight(const ComputeLaplacianData &data, int x, int y, int diagonal)
{
const unsigned char *image = data.image;
int idx1 = y * data.width + x;
int idx2 = idx1 + data.offsets[diagonal];
float attenuation1 = data.attenuationLUT[idx1 / data.width];
float attenuation2 = data.attenuationLUT[idx2 / data.width];
float gradient = abs(image[idx1]*attenuation1/255.0f - image[idx2]*attenuation2/255.0f) * data.gradientScaling;
float weight = exp(-data.beta * (gradient + data.gammaList[diagonal]));
return weight + 1e-4;
}
void CudaConfidenceMapsSystemSolver::uploadImage(const unsigned char* imageData, int imageWidth, int imageHeight,
float gradientScaling, float alpha, float beta, float gamma,
bool isUpsideDown) {
resizeDataStructures(imageWidth, imageHeight, isUpsideDown);
// Gather all the parameters needed to create the system in one place
ComputeLaplacianData data;
data.alpha = alpha;
data.beta = beta;
data.gamma = gamma;
data.gradientScaling = gradientScaling;
data.image = imageData;
data.width = imageWidth;
data.height = imageHeight;
data.centralDiagonal = 4;
int offsets[9] = {-imageWidth-1, -imageWidth, -imageWidth+1, -1, 0, 1, imageWidth-1, imageWidth, imageWidth+1};
float gammaList[9] = {sqrt(2.0f)*gamma, 0.0f, sqrt(2.0f)*gamma, gamma, 0.0f, gamma, sqrt(2.0f)*gamma, 0.0f, sqrt(2.0f)*gamma};
for (int i = 0; i < 9; ++i) {
data.offsets[i] = offsets[i];
data.gammaList[i] = gammaList[i];
_gpuData->L_h.diagonal_offsets[i] = offsets[i];
}
// Precompute attenuation tables
data.attenuationLUT = std::vector<float>(imageHeight);
for (int i = 0; i < imageHeight; ++i) {
float y = (float)i / (float)(imageHeight-1);
if (isUpsideDown) y = 1 - y;
data.attenuationLUT[i] = 1 - exp(-alpha * y);
}
// Fill in the rest of the matrix
for (int y = 0; y < imageHeight; ++y) {
for (int x = 0; x < imageWidth; ++x) {
int idx = y * imageWidth + x;
// Filter off out-of-bounds edges
unsigned short filter = 495; // 111 101 111
// 8 - neighbourhood filter
if (x == 0) filter &= 203; // 011 001 011
if (x == imageWidth-1) filter &= 422; // 110 100 110
if (y == 0) filter &= 47; // 000 101 111
if (y == imageHeight-1) filter &= 488; // 111 101 000
float valueSum = 0.0f;
if (y == 0 || y == imageHeight - 1) valueSum = 1.0f;
for (int d = 0; d < 9; ++d) {
_gpuData->L_h.values(idx, d) = 0;
float value = 0.0f;
if (((256>>d) & filter) != 0) {
value = _getWeight(data, x, y, d);
_gpuData->L_h.values(idx, d) = -value;
}
valueSum += value;
}
_gpuData->L_h.values(idx, data.centralDiagonal) = valueSum;
}
}
// Upload system
_gpuData->L_d = _gpuData->L_h;
}
void CudaConfidenceMapsSystemSolver::resetSolution() {
int width = _gpuData->imageWidth;
int height = _gpuData->imageHeight;
bool isUpsideDown = _gpuData->isUpsideDown;
for (int y = 0; y < height; ++y) {
for (int x = 0; x < width; ++x) {
float value;
if (!isUpsideDown)
value = y / (height - 1.0f);
else
value = 1.0f - y / (height - 1.0f);
_gpuData->x_h[y * width + x] = value;
}
}
// Uplaod solution vector on the GPU
_gpuData->x_d = _gpuData->x_h;
}
void CudaConfidenceMapsSystemSolver::solve(int maximumIterations, float errorTolerance) {
cusp::default_monitor<float> monitor(_gpuData->b_d, maximumIterations, errorTolerance);
cusp::precond::diagonal<float, cusp::device_memory> M(_gpuData->L_d);
cusp::krylov::cg(_gpuData->L_d, _gpuData->x_d, _gpuData->b_d, monitor, M);
_gpuData->solutionResidualNorm = monitor.residual_norm();
}
const float* CudaConfidenceMapsSystemSolver::getSolution(int& width, int& height) {
// Downlaod data
_gpuData->x_h = _gpuData->x_d;
return &_gpuData->x_h[0];
}
float CudaConfidenceMapsSystemSolver::getSolutionResidualNorm() const {
return _gpuData->solutionResidualNorm;
}
void CudaConfidenceMapsSystemSolver::resizeDataStructures(int imageWidth, int imageHeight, bool isUpsideDown) {
// If the problem size changed, reset the solution vector, as well as resizing the vectors and matrices
if (_gpuData->imageWidth != imageWidth || _gpuData->imageHeight != imageHeight || _gpuData->isUpsideDown != isUpsideDown) {
// Resize the system vectors and matrices to accomodate the different image sze
_gpuData->imageWidth = imageWidth;
_gpuData->imageHeight = imageHeight;
int numElements = imageWidth * imageHeight;
_gpuData->x_h.resize(numElements);
_gpuData->b_d.resize(numElements);
_gpuData->x_d.resize(numElements);
_gpuData->L_d.resize(numElements, numElements, numElements * 9, 9);
_gpuData->L_h.resize(numElements, numElements, numElements * 9, 9);
// Set the b vector to 0 (except for the row corresponding to the seed points set to 1)
if (isUpsideDown) {
CuspDeviceVector::view lastRow(_gpuData->b_d.begin() + (numElements - imageWidth), _gpuData->b_d.end());
CuspDeviceVector::view rest(_gpuData->b_d.begin(), _gpuData->b_d.begin() + (numElements - imageWidth));
cusp::blas::fill(rest, 0.0f);
cusp::blas::fill(lastRow, 1.0f);
}
else {
CuspDeviceVector::view firstRow(_gpuData->b_d.begin(), _gpuData->b_d.begin() + imageWidth);
CuspDeviceVector::view rest(_gpuData->b_d.begin() + imageWidth, _gpuData->b_d.end());
cusp::blas::fill(firstRow, 1.0f);
cusp::blas::fill(rest, 0.0f);
}
resetSolution();
}
}
} // cuda
} // campvis
\ No newline at end of file
#ifndef CUDADCONFIDENCEMAPS_CUDA_H__
#define CUDADCONFIDENCEMAPS_CUDA_H__
namespace campvis {
namespace cuda {
class CudaConfidenceMapsSystemGPUData;
class CudaConfidenceMapsSystemSolver {
public:
CudaConfidenceMapsSystemSolver();
/**
* Uploads an image to the solver. Additionally the matrices and vectors that are needed to
* solve the system are created. After calling this function it is possible to call \see solve() in order
* to start computing confidence maps.
* \param imageData pointer to a buffer containing the grayscale image 8-bit image data
* \param imageWidth width of the image
* \param imageHeight height of the image
* \param gradientScaling multiplication applied to the computed image gradients
* \param alpha controls the depth attenuation correction
* \param beta controls the non-linear mapping of gradients to weights
* \param gamma controls how much diagonal connections are penalized
* \param isUpsideDown if set to true, the image is interpreted as being upside down (as common in OpenGL)
*/
void uploadImage(const unsigned char* imageData, int imageWidth, int imageHeight,
float gradientScaling, float alpha, float beta, float gamma, bool isUpsideDown=true);
/**
* Resets the current solution vector to a linear fallof from white (top of the image) to black
* (bottom of the image). Note the calling \see uploadImage() for the first time or with a different
* image size, also results in the solution being reset.
*/
void resetSolution();
/**
* After calling \see uploadImage(), this functions launches a solver on the GPU that will solve
* the diffusion problem.
* \param maximumIterations maximum number of iterations the solver will preform
* \param errorTolerance if the solution error sinks below this value, the solver stops early
*/
void solve(int maximumIterations, float errorTolerance);
/**
* Returns a host buffer of the last solution computed by the solver. The pointer is guaranteed to
* remain valid until the next call of \see uploadImage() or the object is destructed.
* \param width outputs the width of the image
* \param height outputs the height of the image
* \return a pointer to a float buffer containing the confidence map
*/
const float* getSolution(int& width, int& height);
/**
* Returns the residual norm of the solution as a measure of error;
*/
float getSolutionResidualNorm() const;
private:
void resizeDataStructures(int imageWidth, int imageHeight, bool isUpsideDown);
CudaConfidenceMapsSystemGPUData *_gpuData;
};
} // cuda
} // campvis
#endif // CUDADCONFIDENCEMAPS_CUDA_H__
\ No newline at end of file