Introducing various processors to preprocessing module all using the mighty power of OpenGL:

* GlGaussianFilter performs a Gaussian blur
* GlIntensityQuantizer quantizes image intensities into bins using a transfer function
* GlSignalToNoiseRatioFilter computes a very simple variance-based SnR value
* GlVesselnessFilter computes a Frangi-similar vesselness measure
parent 44ac5827
......@@ -14,7 +14,8 @@ namespace tgt {
public:
enum TargetType {
ARRAY_BUFFER = GL_ARRAY_BUFFER,
ELEMENT_ARRAY_BUFFER = GL_ELEMENT_ARRAY_BUFFER
ELEMENT_ARRAY_BUFFER = GL_ELEMENT_ARRAY_BUFFER,
TEXTURE_BUFFER = GL_TEXTURE_BUFFER
};
enum UsageType {
......
// ================================================================================================
//
// This file is part of the CAMPVis Software Framework.
//
// If not explicitly stated otherwise: Copyright (C) 2012-2013, all rights reserved,
// Christian Schulte zu Berge <christian.szb@in.tum.de>
// Chair for Computer Aided Medical Procedures
// Technische Universität München
// Boltzmannstr. 3, 85748 Garching b. München, Germany
//
// For a full list of authors and contributors, please refer to the file "AUTHORS.txt".
//
// Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file
// except in compliance with the License. You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software distributed under the
// License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND,
// either express or implied. See the License for the specific language governing permissions
// and limitations under the License.
//
// ================================================================================================
in vec3 ex_TexCoord;
out vec4 out_Color;
#include "tools/texture3d.frag"
uniform sampler3D _texture;
uniform TextureParameters3D _textureParams;
uniform float _zTexCoord;
uniform ivec3 _direction;
uniform int _halfKernelSize;
uniform samplerBuffer _kernel;
void main() {
ivec3 texel = ivec3(vec3(ex_TexCoord.xy, _zTexCoord) * _textureParams._size);
vec4 result = vec4(0.0, 0.0, 0.0, 0.0);
float norm = 0.0;
for (int i = -_halfKernelSize; i <= _halfKernelSize; ++i) {
// TODO: why the fuck does abs(i) not work here?!?
int absi = (i < 0) ? -i : i;
ivec3 lookupTexel = texel + (_direction * i);
if (all(greaterThanEqual(lookupTexel, ivec3(0, 0, 0))) && all(lessThan(lookupTexel, _textureParams._size))) {
vec4 curValue = texelFetch(_texture, lookupTexel, 0);
float mult = texelFetch(_kernel, absi).r;
result += curValue * mult;
norm += mult;
}
}
out_Color = result / norm;
}
// ================================================================================================
//
// This file is part of the CAMPVis Software Framework.
//
// If not explicitly stated otherwise: Copyright (C) 2012-2013, all rights reserved,
// Christian Schulte zu Berge <christian.szb@in.tum.de>
// Chair for Computer Aided Medical Procedures
// Technische Universität München
// Boltzmannstr. 3, 85748 Garching b. München, Germany
//
// For a full list of authors and contributors, please refer to the file "AUTHORS.txt".
//
// Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file
// except in compliance with the License. You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software distributed under the
// License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND,
// either express or implied. See the License for the specific language governing permissions
// and limitations under the License.
//
// ================================================================================================
in vec3 ex_TexCoord;
out vec4 out_Color;
#include "tools/gradient.frag"
#include "tools/texture3d.frag"
#include "tools/transferfunction.frag"
uniform sampler3D _texture;
uniform TextureParameters3D _textureParams;
uniform sampler1D _transferFunction;
uniform TFParameters1D _transferFunctionParams;
uniform float _zTexCoord;
uniform int _numberOfBins;
void main() {
vec4 intensity = texture(_texture, vec3(ex_TexCoord.xy, _zTexCoord));
intensity.x = lookupTF(_transferFunction, _transferFunctionParams, intensity.x).a;
intensity.y = lookupTF(_transferFunction, _transferFunctionParams, intensity.y).a;
intensity.z = lookupTF(_transferFunction, _transferFunctionParams, intensity.z).a;
intensity.w = lookupTF(_transferFunction, _transferFunctionParams, intensity.w).a;
out_Color = floor(intensity * float(_numberOfBins + 1)) / float(_numberOfBins);
}
// ================================================================================================
//
// This file is part of the CAMPVis Software Framework.
//
// If not explicitly stated otherwise: Copyright (C) 2012-2013, all rights reserved,
// Christian Schulte zu Berge <christian.szb@in.tum.de>
// Chair for Computer Aided Medical Procedures
// Technische Universität München
// Boltzmannstr. 3, 85748 Garching b. München, Germany
//
// For a full list of authors and contributors, please refer to the file "AUTHORS.txt".
//
// Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file
// except in compliance with the License. You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software distributed under the
// License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND,
// either express or implied. See the License for the specific language governing permissions
// and limitations under the License.
//
// ================================================================================================
in vec3 ex_TexCoord;
out float out_result;
#include "tools/gradient.frag"
#include "tools/texture3d.frag"
uniform sampler3D _texture;
uniform TextureParameters3D _textureParams;
uniform float _zTexCoord;
void sampleOffset(inout float sum, inout float sumsq, inout int count, in float referenceValue, const in ivec3 offset) {
float value = abs(referenceValue - textureOffset(_texture, vec3(ex_TexCoord.xy, _zTexCoord), offset).r);
sum += value;
sumsq += value * value;
++count;
}
void main() {
float refVal = texture(_texture, vec3(ex_TexCoord.xy, _zTexCoord)).r;
float sum = refVal;
float sumsq = refVal * refVal;
int count = 1;
sampleOffset(sum, sumsq, count, refVal, ivec3(-1, 0, 0));
sampleOffset(sum, sumsq, count, refVal, ivec3( 1, 0, 0));
sampleOffset(sum, sumsq, count, refVal, ivec3( 0, -1, 0));
sampleOffset(sum, sumsq, count, refVal, ivec3( 0, 1, 0));
sampleOffset(sum, sumsq, count, refVal, ivec3( 0, 0, -1));
sampleOffset(sum, sumsq, count, refVal, ivec3( 0, 0, 1));
sampleOffset(sum, sumsq, count, refVal, ivec3(-2, -2, -2));
sampleOffset(sum, sumsq, count, refVal, ivec3( 2, -2, -2));
sampleOffset(sum, sumsq, count, refVal, ivec3(-2, 2, -2));
sampleOffset(sum, sumsq, count, refVal, ivec3( 2, 2, -2));
sampleOffset(sum, sumsq, count, refVal, ivec3(-2, -2, 2));
sampleOffset(sum, sumsq, count, refVal, ivec3( 2, -2, 2));
sampleOffset(sum, sumsq, count, refVal, ivec3(-2, 2, 2));
sampleOffset(sum, sumsq, count, refVal, ivec3( 2, 2, 2));
sampleOffset(sum, sumsq, count, refVal, ivec3(-3, 0, 0));
sampleOffset(sum, sumsq, count, refVal, ivec3( 3, 0, 0));
sampleOffset(sum, sumsq, count, refVal, ivec3( 0, -3, 0));
sampleOffset(sum, sumsq, count, refVal, ivec3( 0, 3, 0));
sampleOffset(sum, sumsq, count, refVal, ivec3( 0, 0, -3));
sampleOffset(sum, sumsq, count, refVal, ivec3( 0, 0, 3));
float mean = sum / float(count);
float sdev = sqrt( (1.0 / float(count - 1)) * (sumsq - (sum * sum / float(count))));
out_result = (sdev != 0.0) ? mean / sdev : 0.0;
}
// ================================================================================================
//
// This file is part of the CAMPVis Software Framework.
//
// If not explicitly stated otherwise: Copyright (C) 2012-2013, all rights reserved,
// Christian Schulte zu Berge <christian.szb@in.tum.de>
// Chair for Computer Aided Medical Procedures
// Technische Universität München
// Boltzmannstr. 3, 85748 Garching b. München, Germany
//
// For a full list of authors and contributors, please refer to the file "AUTHORS.txt".
//
// Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file
// except in compliance with the License. You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software distributed under the
// License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND,
// either express or implied. See the License for the specific language governing permissions
// and limitations under the License.
//
// ================================================================================================
in vec3 ex_TexCoord;
out float out_vesselness;
//out vec3 out_evals;
#include "tools/gradient.frag"
#include "tools/texture3d.frag"
uniform sampler3D _texture;
uniform TextureParameters3D _textureParams;
uniform sampler3D _labels;
uniform TextureParameters3D _labelsTexParams;
uniform int _labelBit = 3;
uniform float _zTexCoord;
uniform vec2 _lod = vec2(1.0, 3.0);
uniform float _alpha = 0.5;
uniform float _beta = 0.5;
uniform float _gamma = 0.001;
uniform float _theta = 0.5;
vec4 computeVesselness(in vec3 texCoords, in float lod) {
if (true || texture(_texture, texCoords).rgb != vec3(0.0, 0.0, 0.0)) {
mat3 hessian = computeHessianCentralDifferencesLod(_texture, _textureParams, texCoords, lod);
vec3 evals = computeEigenvalues(hessian);
if (abs(evals.x) < abs(evals.y)) {
if (abs(evals.x) < abs(evals.z)) {
evals.xyz = (abs(evals.y) < abs(evals.z)) ? evals.xyz : evals.xzy;
}
else {
evals.xyz = evals.zxy;
}
}
else {
if (abs(evals.x) < abs(evals.z)) {
evals.xyz = evals.yxz;
}
else {
evals.xyz = (abs(evals.y) < abs(evals.z)) ? evals.yzx : evals.zyx;
}
}
float l1 = evals.x;
float l2 = evals.y;
float l3 = evals.z;
float rb = abs(l1) / sqrt(abs(l2 * l3));
float ra = abs(l2) / abs(l3);
float s = sqrt(l1 * l1 + l2 * l2 + l3 * l3);
vec3 v3 = (hessian - mat3(l1)) * (hessian - mat3(l2))[0];
float fbdir = abs(dot(normalize(v3), vec3(0, 1, 0)));
float vesselness = 0.0;
if (l2 < 0.0 && l3 < 0.0) {
vesselness = (1.0 - exp(- (ra*ra)/(2.0 * _alpha * _alpha))) * exp(- (rb*rb)/(2.0 * _beta * _beta)) * (1.0 - exp(- (s*s)/(2.0 * _gamma * _gamma))) * (1.0 - exp(- (fbdir*fbdir)/(2.0 * _theta * _theta)));
}
//out_evals = vec3(fbdir);//max(out_evals, evals);
return vec4(vesselness, ra, rb, s);
}
else {
return vec4(0.0, 0.0, 0.0, 0.0);
}
}
void main() {
vec3 texCoords = vec3(ex_TexCoord.xy, _zTexCoord);
//out_evals = vec3(1337.0);
out_vesselness = 0.0;
uint l = uint(texture(_labels, texCoords).r * 255);
if (true || bitfieldExtract(l, 4, 1) != 0U) {
for (float lod = _lod.x; lod < _lod.y; lod += 0.5) {
float level_vesselness = computeVesselness(texCoords, lod);
out_vesselness = max(out_vesselness, level_vesselness.x);
}
}
}
// ================================================================================================
//
// This file is part of the CAMPVis Software Framework.
//
// If not explicitly stated otherwise: Copyright (C) 2012-2013, all rights reserved,
// Christian Schulte zu Berge <christian.szb@in.tum.de>
// Chair for Computer Aided Medical Procedures
// Technische Universität München
// Boltzmannstr. 3, 85748 Garching b. München, Germany
//
// For a full list of authors and contributors, please refer to the file "AUTHORS.txt".
//
// Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file
// except in compliance with the License. You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software distributed under the
// License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND,
// either express or implied. See the License for the specific language governing permissions
// and limitations under the License.
//
// ================================================================================================
#include "glgaussianfilter.h"
#include "tgt/buffer.h"
#include "tgt/logmanager.h"
#include "tgt/shadermanager.h"
#include "tgt/textureunit.h"
#include "tgt/texture.h"
#include "core/datastructures/imagedata.h"
#include "core/datastructures/imagerepresentationgl.h"
#include "core/datastructures/renderdata.h"
#include "core/tools/quadrenderer.h"
#include "core/tools/stringutils.h"
namespace campvis {
const std::string GlGaussianFilter::loggerCat_ = "CAMPVis.modules.classification.GlGaussianFilter";
#define MAX_SIGMA 20.f
#define MAX_HALF_KERNEL_SIZE static_cast<int>(MAX_SIGMA * 2.5f) + 2
GlGaussianFilter::GlGaussianFilter(IVec2Property* viewportSizeProp)
: VisualizationProcessor(viewportSizeProp)
, p_inputImage("InputImage", "Input Image", "", DataNameProperty::READ)
, p_outputImage("OutputImage", "Output Image", "GlGaussianFilter.out", DataNameProperty::WRITE)
, p_sigma("Sigma", "Sigma (relates to kernel size)", 2.5f, 1.f, MAX_SIGMA, .1f, 1)
, _shader(0)
, _kernelBuffer(0)
, _kernelBufferTexture(0)
{
addProperty(&p_inputImage);
addProperty(&p_outputImage);
addProperty(&p_sigma);
}
GlGaussianFilter::~GlGaussianFilter() {
}
void GlGaussianFilter::init() {
VisualizationProcessor::init();
_shader = ShdrMgr.load("core/glsl/passthrough.vert", "modules/preprocessing/glsl/glgaussianfilter.frag", "");
_shader->setAttributeLocation(0, "in_Position");
_shader->setAttributeLocation(1, "in_TexCoord");
// create kernel buffer
tgt::TextureUnit inputUnit;
inputUnit.activate();
_kernelBuffer = new tgt::BufferObject(tgt::BufferObject::TEXTURE_BUFFER, tgt::BufferObject::USAGE_STATIC_DRAW);
glGenTextures(1, &_kernelBufferTexture);
LGL_ERROR;
}
void GlGaussianFilter::deinit() {
ShdrMgr.dispose(_shader);
delete _kernelBuffer;
glDeleteTextures(1, &_kernelBufferTexture);
VisualizationProcessor::deinit();
}
void GlGaussianFilter::updateResult(DataContainer& data) {
ImageRepresentationGL::ScopedRepresentation img(data, p_inputImage.getValue());
if (img != 0) {
tgt::ivec3 size = img->getSize();
int halfKernelSize = static_cast<int>(2.5 * p_sigma.getValue());
tgtAssert(halfKernelSize < MAX_HALF_KERNEL_SIZE, "halfKernelSize too big -> kernel uniform buffer will be out of bounds!")
tgt::TextureUnit inputUnit, kernelUnit;
inputUnit.activate();
// create texture for result
tgt::Texture* resultTextures[2];
for (size_t i = 0; i < 2; ++i) {
resultTextures[i] = new tgt::Texture(0, size, img->getTexture()->getFormat(), img->getTexture()->getInternalFormat(), img->getTexture()->getDataType(), tgt::Texture::LINEAR);
resultTextures[i]->uploadTexture();
}
// create and upload kernel buffer
GLfloat kernel[MAX_HALF_KERNEL_SIZE];
for (int i = 0; i <= halfKernelSize; ++i) {
kernel[i] = exp(- static_cast<GLfloat>(i*i) / (2.f * p_sigma.getValue() * p_sigma.getValue()));
}
_kernelBuffer->data(kernel, (halfKernelSize + 1) * sizeof(GLfloat), tgt::BufferObject::FLOAT, 1);
// activate shader
_shader->activate();
_shader->setUniform("_halfKernelSize", halfKernelSize);
// bind kernel buffer texture
kernelUnit.activate();
glBindTexture(GL_TEXTURE_BUFFER, _kernelBufferTexture);
glTexBuffer(GL_TEXTURE_BUFFER, GL_R32F, _kernelBuffer->getId());
_shader->setUniform("_kernel", kernelUnit.getUnitNumber());
LGL_ERROR;
// activate FBO and attach texture
_fbo->activate();
glViewport(0, 0, static_cast<GLsizei>(size.x), static_cast<GLsizei>(size.y));
// start 3 passes of convolution: in X, Y and Z direction:
{
// X pass
_shader->setUniform("_direction", tgt::ivec3(1, 0, 0));
img->bind(_shader, inputUnit);
// render quad to compute difference measure by shader
for (int z = 0; z < size.z; ++z) {
float zTexCoord = static_cast<float>(z)/static_cast<float>(size.z) + .5f/static_cast<float>(size.z);
_shader->setUniform("_zTexCoord", zTexCoord);
_fbo->attachTexture(resultTextures[0], GL_COLOR_ATTACHMENT0, 0, z);
QuadRdr.renderQuad();
}
}
{
// Y pass
_shader->setUniform("_direction", tgt::ivec3(0, 1, 0));
inputUnit.activate();
resultTextures[0]->bind();
// render quad to compute difference measure by shader
for (int z = 0; z < size.z; ++z) {
float zTexCoord = static_cast<float>(z)/static_cast<float>(size.z) + .5f/static_cast<float>(size.z);
_shader->setUniform("_zTexCoord", zTexCoord);
_fbo->attachTexture(resultTextures[1], GL_COLOR_ATTACHMENT0, 0, z);
QuadRdr.renderQuad();
}
}
{
// Z pass
_shader->setUniform("_direction", tgt::ivec3(0, 0, 1));
inputUnit.activate();
resultTextures[1]->bind();
// render quad to compute difference measure by shader
for (int z = 0; z < size.z; ++z) {
float zTexCoord = static_cast<float>(z)/static_cast<float>(size.z) + .5f/static_cast<float>(size.z);
_shader->setUniform("_zTexCoord", zTexCoord);
_fbo->attachTexture(resultTextures[0], GL_COLOR_ATTACHMENT0, 0, z);
QuadRdr.renderQuad();
}
}
_fbo->detachAll();
_fbo->deactivate();
_shader->deactivate();
// put resulting image into DataContainer
ImageData* id = new ImageData(3, size, img->getParent()->getNumChannels());
ImageRepresentationGL::create(id, resultTextures[0]);
id->setMappingInformation(img->getParent()->getMappingInformation());
data.addData(p_outputImage.getValue(), id);
delete resultTextures[1];
tgt::TextureUnit::setZeroUnit();
LGL_ERROR;
}
else {
LERROR("No suitable input image found.");
}
validate(INVALID_RESULT);
}
}
// ================================================================================================
//
// This file is part of the CAMPVis Software Framework.
//
// If not explicitly stated otherwise: Copyright (C) 2012-2013, all rights reserved,
// Christian Schulte zu Berge <christian.szb@in.tum.de>
// Chair for Computer Aided Medical Procedures
// Technische Universität München
// Boltzmannstr. 3, 85748 Garching b. München, Germany
//
// For a full list of authors and contributors, please refer to the file "AUTHORS.txt".
//
// Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file
// except in compliance with the License. You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software distributed under the
// License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND,
// either express or implied. See the License for the specific language governing permissions
// and limitations under the License.
//
// ================================================================================================
#ifndef GLGAUSSIANFILTER_H__
#define GLGAUSSIANFILTER_H__
#include <string>
#include "core/pipeline/abstractprocessordecorator.h"
#include "core/pipeline/visualizationprocessor.h"
#include "core/properties/datanameproperty.h"
#include "core/properties/floatingpointproperty.h"
#include "core/properties/optionproperty.h"
namespace tgt {
class BufferObject;
class Shader;
}
namespace campvis {
/**
* Performs a gaussian filtering on the input image using OpenGL.
*/
class GlGaussianFilter : public VisualizationProcessor {
public:
/**
* Constructs a new GlGaussianFilter Processor
**/
GlGaussianFilter(IVec2Property* viewportSizeProp);
/**
* Destructor
**/
virtual ~GlGaussianFilter();
/// \see AbstractProcessor::init
virtual void init();
/// \see AbstractProcessor::deinit
virtual void deinit();
/// \see AbstractProcessor::getName()
virtual const std::string getName() const { return "GlGaussianFilter"; };
/// \see AbstractProcessor::getDescription()
virtual const std::string getDescription() const { return "Performs a gaussian filtering on the input image using OpenGL."; };
/// \see AbstractProcessor::getAuthor()
virtual const std::string getAuthor() const { return "Christian Schulte zu Berge <christian.szb@in.tum.de>"; };
/// \see AbstractProcessor::getProcessorState()
virtual ProcessorState getProcessorState() const { return AbstractProcessor::TESTING; };
DataNameProperty p_inputImage; ///< ID for input volume
DataNameProperty p_outputImage; ///< ID for output gradient volume
FloatProperty p_sigma; ///< Sigma for specifying kernel size