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

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

Simplified equation system creation code for confidence maps

parent c470e61c
......@@ -101,40 +101,30 @@ struct ComputeLaplacianData
int centralDiagonal;
int offsets[9];
float gammaList[9];
};
static inline float _getAttenuation(int y, int height, float alpha) {
return (1 - exp(-alpha * ((float)y / (float)(height-1))));
}
std::vector<float> attenuationLUT;
};
static inline float _getGradient(const ComputeLaplacianData &data, int idx, int offset)
static inline float _getWeight(const ComputeLaplacianData &data, int x, int y, int diagonal, bool isUpsideDown)
{
const unsigned char *image = data.image;
int y1 = idx / data.width;
int y2 = (idx+offset) / data.width;
float a1 = _getAttenuation(y1, data.height, data.alpha);
float a2 = _getAttenuation(y2, data.height, data.alpha);
int idx1 = y * data.width + x;
int idx2 = idx1 + data.offsets[diagonal];
return abs(image[idx]*a1/255.0f - image[idx+offset]*a2/255.0f);
}
float attenuation1 = data.attenuationLUT[idx1 / data.width];
float attenuation2 = data.attenuationLUT[idx2 / data.width];
static inline float _calculateWeight(float gradient, float beta, float gamma, float scaling)
{
return exp(-beta * (gradient*scaling + gamma));
}
float gradient = abs(image[idx1]*attenuation1/255.0f - image[idx2]*attenuation2/255.0f) * data.gradientScaling;
static inline float _getWeight(const ComputeLaplacianData &data, int x, int y, int diagonal)
{
float gradient = _getGradient(data, y * data.width + x, data.offsets[diagonal]);
float weight = _calculateWeight(gradient, data.beta, data.gammaList[diagonal], 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)
float gradientScaling,
bool isUpsideDown)
{
// Gather all of the options together
ComputeLaplacianData data;
......@@ -155,38 +145,46 @@ void CUSP_CM_buildEquationSystem(CuspGPUData &gpudata, const unsigned char* imag
gpudata.L_h.diagonal_offsets[i] = data.offsets[i];
}
// Reset B
for (int x = 0; x < width * height; ++x) {
gpudata.b_h[x] = 0.0f;
}
// Fill in first row of b, which sholud be one
for (int x = 0; x < width; ++x) {
gpudata.b_h[x] = 1.0f;
// 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);
}
// Fill in rows of A corresponding to first and last row of the image
for (int x = 0; x < width; ++x) {
for (int d = 0; d < 9; ++d) {
gpudata.L_h.values(x, d) = (d == data.centralDiagonal ? 1.0f : 0.0f);
gpudata.L_h.values(width*(height-1) + x, d) = (d == data.centralDiagonal ? 1.0f : 0.0f);
}
// 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 = 1; y < height-1; ++y) {
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 == 1) filter &= 47; // 000 101 111
if (y == height-2) filter &= 488; // 111 101 000
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;
......@@ -194,17 +192,8 @@ void CUSP_CM_buildEquationSystem(CuspGPUData &gpudata, const unsigned char* imag
float value = 0.0f;
if (((256>>d) & filter) != 0) {
value = _getWeight(data, x, y, d);
value = _getWeight(data, x, y, d, isUpsideDown);
gpudata.L_h.values(idx, d) = -value;
} else if(y == 1 || y == height - 2) {
unsigned short filter2 = 495; // 111 101 11
if (x == 0) filter2 &= 203; // 011 001 011
if (x == width-1) filter2 &= 422; // 110 100 110
if (((256>>d) & filter2) != 0) {
value = _getWeight(data, x, y, d);
if (y == 1) gpudata.b_h[idx] += value;
}
}
valueSum += value;
......
#pragma once
#include <vector>
#include <GL/gl.h>
class CuspGPUData;
......@@ -9,9 +10,13 @@ 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);
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);
\ No newline at end of file
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
......@@ -29,8 +29,11 @@
#include "core/datastructures/imagedata.h"
#include "core/datastructures/genericimagerepresentationlocal.h"
#include "core/datastructures/imagerepresentationgl.h"
#include "core/tools/stringutils.h"
#include "modules/cudaconfidencemaps/core/cm_cusp_cuda_exports.h"
namespace campvis {
const std::string CudaConfidenceMapsSolver::loggerCat_ = "CAMPVis.modules.cudaconfidencemaps.solver";
......@@ -58,27 +61,16 @@ namespace campvis {
addProperty(p_paramGamma);
}
CudaConfidenceMapsSolver::~CudaConfidenceMapsSolver() {
/*if (_receiverRunning)
stopReceiver();
if (_receiverThread)
delete _receiverThread;*/
}
CudaConfidenceMapsSolver::~CudaConfidenceMapsSolver() { }
void CudaConfidenceMapsSolver::init() {/*
p_connect.s_clicked.connect(this, &CudaConfidenceMapsSolver::connect);
p_disconnect.s_clicked.connect(this, &CudaConfidenceMapsSolver::disconnect);
p_disconnect.setVisible(false);*/
void CudaConfidenceMapsSolver::init() {
//CUSP_CM_initCuda();
}
void CudaConfidenceMapsSolver::deinit() {/*
stopReceiver();
p_connect.s_clicked.disconnect(this);
p_disconnect.s_clicked.disconnect(this);*/
}
void CudaConfidenceMapsSolver::deinit() { }
void CudaConfidenceMapsSolver::updateResult(DataContainer& data) {
ImageRepresentationLocal::ScopedRepresentation img(data, p_inputImage.getValue());
if (img != 0) {
int iterations = p_iterations.getValue();
......@@ -86,124 +78,27 @@ namespace campvis {
float alpha = p_paramAlpha.getValue();
float beta = p_paramBeta.getValue();
float gamma = p_paramGamma.getValue();
bool isFlipped = true;
cgt::ivec3 size = img->getSize();
size_t elementCount = cgt::hmul(size);
unsigned char *flippedData = new unsigned char[elementCount];
for (size_t i = 0; i < elementCount; ++i) {
flippedData[elementCount-i-1] = ((unsigned char*)img->getWeaklyTypedPointer()._pointer)[i];
}
_solver.createSystem(flippedData, size.x, size.y,
gradientScaling, alpha, beta, gamma);
delete[] flippedData;
auto image = (unsigned char*)img->getWeaklyTypedPointer()._pointer;
_solver.createSystem(image, size.x, size.y, gradientScaling, alpha, beta, gamma);
_solver.solve(iterations);
const float *solution = _solver.getSolution(size.x, size.y);
float *solutionCopy = new float[elementCount];
memcpy(solutionCopy, solution, sizeof(float) * elementCount);
WeaklyTypedPointer wtpData(WeaklyTypedPointer::FLOAT, 1, solutionCopy);
ImageData *id = new ImageData(img->getParent()->getDimensionality(), size, img->getParent()->getNumChannels());
float *imgData = new float[elementCount];
for (size_t i = 0; i < elementCount; ++i) {
imgData[elementCount-i-1] = solution[i];
}
WeaklyTypedPointer wtpData(WeaklyTypedPointer::FLOAT, 1, imgData);
ImageRepresentationLocal::create(id, wtpData);
id->setMappingInformation(img->getParent()->getMappingInformation());
data.addData(p_outputConfidenceMap.getValue(), id);
}
/*
if(p_receiveTransforms.getValue())
{
_transformMutex.lock();
for(auto it = _receivedTransforms.begin(), end = _receivedTransforms.end(); it != end; ++it) {
TransformData * td = new TransformData(it->second);
data.addData(p_targetTransformPrefix.getValue() + it->first, td);
#ifdef IGTL_CLIENT_DEBUGGING
LDEBUG("Transform data put into container. ");
#endif
}
_receivedTransforms.clear();
_transformMutex.unlock();
}
if(p_receiveImages.getValue())
{
_imageMutex.lock();
for(auto it = _receivedImages.begin(), end = _receivedImages.end(); it != end; ++it)
{
igtl::ImageMessage::Pointer imageMessage = it->second;
WeaklyTypedPointer wtp;
wtp._pointer = new uint8_t[imageMessage->GetImageSize()];
#ifdef IGTL_CLIENT_DEBUGGING
LDEBUG("Image has " << imageMessage->GetNumComponents() << " components and is of size " << imageMessage->GetImageSize());
#endif
memcpy(wtp._pointer, imageMessage->GetScalarPointer(), imageMessage->GetImageSize());
wtp._numChannels = imageMessage->GetNumComponents();
switch (imageMessage->GetScalarType()) {
case igtl::ImageMessage::TYPE_INT8:
wtp._baseType = WeaklyTypedPointer::INT8; break;
case igtl::ImageMessage::TYPE_UINT8:
wtp._baseType = WeaklyTypedPointer::UINT8; break;
case igtl::ImageMessage::TYPE_INT16:
wtp._baseType = WeaklyTypedPointer::INT16; break;
case igtl::ImageMessage::TYPE_UINT16:
wtp._baseType = WeaklyTypedPointer::UINT16; break;
case igtl::ImageMessage::TYPE_INT32:
wtp._baseType = WeaklyTypedPointer::INT32; break;
case igtl::ImageMessage::TYPE_UINT32:
wtp._baseType = WeaklyTypedPointer::UINT32; break;
case igtl::ImageMessage::TYPE_FLOAT32:
wtp._baseType = WeaklyTypedPointer::FLOAT; break;
default:
LERROR("Error while receiving IGTL IMAGE message: unsupported type: " << imageMessage->GetScalarType());
return;
}
tgt::vec3 imageOffset(0.f);
tgt::vec3 voxelSize(1.f);
tgt::ivec3 size_i(1);
imageMessage->GetSpacing(voxelSize.elem);
imageMessage->GetDimensions(size_i.elem);
tgt::svec3 size(size_i);
imageMessage->GetOrigin(imageOffset.elem);
size_t dimensionality = (size_i[2] == 1) ? ((size_i[1] == 1) ? 1 : 2) : 3;
ImageData* image = new ImageData(dimensionality, size, wtp._numChannels);
ImageRepresentationLocal::create(image, wtp);
image->setMappingInformation(ImageMappingInformation(size, p_imageOffset.getValue(), voxelSize * p_voxelSize.getValue()));
data.addData(p_targetImagePrefix.getValue() + it->first, image);
}
_receivedImages.clear();
_imageMutex.unlock();
}
if(p_receivePositions.getValue())
{
_positionMutex.lock();
for(auto it = _receivedPositions.begin(), end = _receivedPositions.end(); it != end; ++it)
{
PositionData * pd = new PositionData(it->second._position, it->second._quaternion);
data.addData(p_targetPositionPrefix.getValue() + it->first, pd);
}
_receivedPositions.clear();
_positionMutex.unlock();
}
validate(INVALID_RESULT);*/
}
void CudaConfidenceMapsSolver::updateProperties(DataContainer& dataContainer) {
/*p_targetImagePrefix.setVisible(p_receiveImages.getValue());
p_imageOffset.setVisible(p_receiveImages.getValue());
p_voxelSize.setVisible(p_receiveImages.getValue());
p_targetTransformPrefix.setVisible(p_receiveImages.getValue() || p_receiveTransforms.getValue());
p_targetPositionPrefix.setVisible(p_receivePositions.getValue());
validate(INVALID_PROPERTIES);*/
}
void CudaConfidenceMapsSolver::updateProperties(DataContainer& dataContainer) { }
}
\ No newline at end of file
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