Commit b0737b5f authored by Declara Denis's avatar Declara Denis Committed by Christian Schulte zu Berge

Added 4-neighbourhood option for US-CM

CudaConfidenceMaps now support a 4 neighbourhood version. This results
in a smaller number of weights and should theoretically be faster.
parent 4247cef0
......@@ -60,6 +60,7 @@ namespace cuda {
// Information about the system as well as statistics of the solution
bool isUpsideDown;
bool useAlphaBetaFiltering;
bool use8Neighbourhood; ///< If set to true the full 8-neighbourhood is used, otherwise only the 4-neighbourhood
int imageWidth;
int imageHeight;
int iterationCount;
......@@ -77,6 +78,7 @@ namespace cuda {
_gpuData->systemCreationTime = 0.0f;
_gpuData->systemSolveTime = 0.0f;
_gpuData->isUpsideDown = true;
_gpuData->use8Neighbourhood = true;
_gpuData->useAlphaBetaFiltering = false;
_gpuData->abFilterAlpha = 0.125;
......@@ -91,8 +93,8 @@ namespace cuda {
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);
bool use8Neighbourhood, bool isUpsideDown) {
resizeDataStructures(imageWidth, imageHeight, isUpsideDown, use8Neighbourhood);
// Measure execution time and record it in the _gpuData datastructure
CUDAClock clock; clock.start();
......@@ -202,21 +204,24 @@ namespace cuda {
}
void CudaConfidenceMapsSystemSolver::resizeDataStructures(int imageWidth, int imageHeight, bool isUpsideDown) {
void CudaConfidenceMapsSystemSolver::resizeDataStructures(int imageWidth, int imageHeight, bool isUpsideDown, bool use8Neighbourhood) {
// If the problem size changed, reset the solution vector, as well as all
// the vectors and matrices
if (_gpuData->imageWidth != imageWidth || _gpuData->imageHeight != imageHeight || _gpuData->isUpsideDown != isUpsideDown) {
if (_gpuData->imageWidth != imageWidth || _gpuData->imageHeight != imageHeight ||
_gpuData->isUpsideDown != isUpsideDown || _gpuData->use8Neighbourhood != use8Neighbourhood) {
// Resize the system vectors and matrices to accomodate the different image size
_gpuData->imageWidth = imageWidth;
_gpuData->imageHeight = imageHeight;
_gpuData->isUpsideDown = isUpsideDown;
_gpuData->use8Neighbourhood = use8Neighbourhood;
int numElements = imageWidth * imageHeight;
int numDiagonals = use8Neighbourhood ? 9 : 5;
_gpuData->x_h.resize(numElements);
_gpuData->b_d.resize(numElements);
_gpuData->x_d.resize(numElements);
_gpuData->imageBuffer_d.resize(numElements);
_gpuData->L_d.resize(numElements, numElements, numElements * 9, 9);
_gpuData->L_h.resize(numElements, numElements, numElements * 9, 9);
_gpuData->L_d.resize(numElements, numElements, numElements * numDiagonals, numDiagonals);
_gpuData->L_h.resize(numElements, numElements, numElements * numDiagonals, numDiagonals);
_gpuData->abFilterR_d.resize(numElements);
_gpuData->abFilterV_d.resize(numElements);
......@@ -306,13 +311,82 @@ namespace cuda {
L[4 * pitch + pidx] = weightSum;
}
// TODO: Unify system creation kernel code for 4 and 8 neighbourhood. DRYness!
static __global__ void k_buildSystem4Neighbourhood(float* L, int pitch, const unsigned char* image, int width, int height,
float gradientScaling, float alpha, float beta, float gamma, bool isUpsideDown)
{
const int x = blockDim.x * blockIdx.x + threadIdx.x;
const int y = blockDim.y * blockIdx.y + threadIdx.y;
const int pidx = y * width + x;
if (x >= width || y >= height) return;
const int offsets[] = {-width, -1, 0, 1, width};
const float gammas[] = {0.0f, gamma, 0.0f, gamma, 0.0f};
// Precompute the three attenuation values for the curernt pixel (the row above, current row, and row below)
float attenuations[3];
if (isUpsideDown) {
attenuations[0] = 1.0f - exp(-alpha * (1.0f - (y - 1.0f) / (height - 1.0f)));
attenuations[1] = 1.0f - exp(-alpha * (1.0f - (y ) / (height - 1.0f)));
attenuations[2] = 1.0f - exp(-alpha * (1.0f - (y + 1.0f) / (height - 1.0f)));
} else {
attenuations[0] = 1.0f - exp(-alpha * (y - 1.0f) / (height - 1.0f));
attenuations[1] = 1.0f - exp(-alpha * (y ) / (height - 1.0f));
attenuations[2] = 1.0f - exp(-alpha * (y + 1.0f) / (height - 1.0f));
}
// Filter off out-of-bounds edges
unsigned char filter = 27; // 1 101 1
// 4 - neighbourhood filter
if (x == 0) filter &= 19; // 1 001 1
if (x == width-1) filter &= 25; // 1 100 1
if (y == 0) filter &= 11; // 0 101 1
if (y == height-1) filter &= 26; // 1 101 0
// get central pixel
float centralValue = image[pidx] * attenuations[1];
// If the pixel is at the top or at the bottom, add a value of 1 to the diagonal, to
// account for the edge to the seed points
float weightSum = 0.0f;
if (y == 0 || y == height - 1)
weightSum = 1.0f;
for (int d = 0; d < 5; ++d) {
float weight = 0.0f;
if (((16>>d) & filter) != 0) {
int pidx_2 = pidx + offsets[d];
float v = image[pidx_2] * attenuations[(d+2)/3];
weight = d_getWeight(centralValue, v, gradientScaling, beta, gammas[d]);
}
// The matrix stores the data, so that values on the same diagonal are sequential.
// This means that all the values from [0, pitch) are on the first diagonal, [pitch, 2*pitch)
// are on the second diagonal and so on...
L[d * pitch + pidx] = -weight;
weightSum += weight;
}
// Store sum of weights in the central diagonal
L[2 * pitch + pidx] = weightSum;
}
void CudaConfidenceMapsSystemSolver::createSystemGPU(const unsigned char* imageData, int imageWidth, int imageHeight,
float gradientScaling, float alpha, float beta, float gamma,
bool isUpsideDown) {
// Initialize the DIA matrix diagonal offsets
int offsets[9] = {-imageWidth-1, -imageWidth, -imageWidth+1, -1, 0, 1, imageWidth-1, imageWidth, imageWidth+1};
for (int i = 0; i < 9; ++i) {
_gpuData->L_d.diagonal_offsets[i] = offsets[i];
if (_gpuData->use8Neighbourhood) {
int offsets[9] = {-imageWidth-1, -imageWidth, -imageWidth+1, -1, 0, 1, imageWidth-1, imageWidth, imageWidth+1};
for (int i = 0; i < 9; ++i) {
_gpuData->L_d.diagonal_offsets[i] = offsets[i];
}
} else {
int offsets[5] = {-imageWidth, -1, 0, 1, imageWidth};
for (int i = 0; i < 5; ++i) {
_gpuData->L_d.diagonal_offsets[i] = offsets[i];
}
}
int numElements = imageWidth * imageHeight;
......@@ -321,10 +395,18 @@ namespace cuda {
// Since the image will be needed by the CUDA kernel, it needs to be copied on the GPU first
cudaMemcpy(thrust::raw_pointer_cast(&_gpuData->imageBuffer_d[0]), imageData, numElements, cudaMemcpyHostToDevice);
k_buildSystem<<<dimGrid, dimBlock>>>(thrust::raw_pointer_cast(&_gpuData->L_d.values.values[0]), _gpuData->L_d.values.pitch,
thrust::raw_pointer_cast(&_gpuData->imageBuffer_d[0]),
imageWidth, imageHeight,
gradientScaling, alpha, beta, gamma, _gpuData->isUpsideDown);
if (_gpuData->use8Neighbourhood) {
k_buildSystem<<<dimGrid, dimBlock>>>(thrust::raw_pointer_cast(&_gpuData->L_d.values.values[0]), _gpuData->L_d.values.pitch,
thrust::raw_pointer_cast(&_gpuData->imageBuffer_d[0]),
imageWidth, imageHeight,
gradientScaling, alpha, beta, gamma, _gpuData->isUpsideDown);
}
else {
k_buildSystem4Neighbourhood<<<dimGrid, dimBlock>>>(thrust::raw_pointer_cast(&_gpuData->L_d.values.values[0]), _gpuData->L_d.values.pitch,
thrust::raw_pointer_cast(&_gpuData->imageBuffer_d[0]),
imageWidth, imageHeight,
gradientScaling, alpha, beta, gamma, _gpuData->isUpsideDown);
}
}
} // cuda
......
......@@ -23,10 +23,12 @@ namespace cuda {
* \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 use8Neighbourhood wether to use a graph connecting all 8 neighbours of a pixel or not. The original
* confidence maps problem formulation uses 8 neighbours.
* \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);
float gradientScaling, float alpha, float beta, float gamma, bool use8Neighbourhood=false, bool isUpsideDown=true);
/**
* Resets the current solution vector to a linear fallof from white (top of the image) to black
......@@ -94,9 +96,7 @@ namespace cuda {
float getSystemSolveTime() const;
private:
void resizeDataStructures(int imageWidth, int imageHeight, bool isUpsideDown);
void createSystemCPU(const unsigned char* imageData, int imageWidth, int imageHeight,
float gradientScaling, float alpha, float beta, float gamma, bool isUpsideDown=true);
void resizeDataStructures(int imageWidth, int imageHeight, bool isUpsideDown, bool use8Neighbourhood);
void createSystemGPU(const unsigned char* imageData, int imageWidth, int imageHeight,
float gradientScaling, float alpha, float beta, float gamma, bool isUpsideDown=true);
......
......@@ -40,6 +40,7 @@ namespace campvis {
, p_inputImage("InputImage", "Input Image", "", DataNameProperty::READ)
, p_outputConfidenceMap("OutputConfidenceMap", "Output Confidence Map", "us.confidence", DataNameProperty::WRITE)
, p_resetResult("ResetSolution", "Reset solution vector")
, p_use8Neighbourhood("Use8Neighbourhood", "Use 8 Neighbourhood (otherwise 4)", true)
, p_iterations("IterationCount", "Conjugate Gradient Iterations", 200, 1, 500)
, p_gradientScaling("GradientScaling", "Scaling factor for gradients", 2.0f, 0.001, 10)
, p_paramAlpha("Alpha", "Alpha (TGC)", 2.0f, 0.001, 10)
......@@ -55,6 +56,7 @@ namespace campvis {
addProperty(p_outputConfidenceMap);
addProperty(p_resetResult);
addProperty(p_use8Neighbourhood);
addProperty(p_iterations);
addProperty(p_gradientScaling);
addProperty(p_paramAlpha);
......@@ -81,6 +83,7 @@ namespace campvis {
ImageRepresentationLocal::ScopedRepresentation img(data, p_inputImage.getValue());
if (img != 0) {
bool use8Neighbourhood = p_use8Neighbourhood.getValue();
int iterations = p_iterations.getValue();
float gradientScaling = p_gradientScaling.getValue();
float alpha = p_paramAlpha.getValue();
......@@ -95,7 +98,7 @@ namespace campvis {
size_t elementCount = cgt::hmul(size);
auto image = (unsigned char*)img->getWeaklyTypedPointer()._pointer;
_solver.uploadImage(image, size.x, size.y, gradientScaling, alpha, beta, gamma);
_solver.uploadImage(image, size.x, size.y, gradientScaling, alpha, beta, gamma, use8Neighbourhood);
_solver.solve(iterations, 1e-10);
const float *solution = _solver.getSolution(size.x, size.y);
......@@ -109,12 +112,6 @@ namespace campvis {
ImageRepresentationGL::create(id, resultTexture);
id->setMappingInformation(img->getParent()->getMappingInformation());
data.addData(p_outputConfidenceMap.getValue(), id);
/*float *solutionCopy = new float[elementCount];
memcpy(solutionCopy, solution, sizeof(float) * elementCount);
WeaklyTypedPointer wtpData(WeaklyTypedPointer::FLOAT, 1, solutionCopy);
ImageRepresentationLocal::create(id, wtpData);*/
}
}
......
......@@ -88,6 +88,8 @@ namespace campvis {
ButtonProperty p_resetResult;
BoolProperty p_use8Neighbourhood; ///< Wether to use 8- or 4-neighbourhood
NumericProperty<int> p_iterations; ///< Number of CG-Iterations to do
FloatProperty p_gradientScaling;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment