Commit 9a76bebc authored by Christian Schulte zu Berge's avatar Christian Schulte zu Berge
Browse files

work on fiber visualization for Columbia module: introducing FiberData and StrainFiberTracker

parent d0eefc5f
......@@ -49,7 +49,7 @@ namespace campvis {
void ButtonPropertyWidget::onButtonClicked(bool) {
ButtonProperty* bp = static_cast<ButtonProperty*>(_property);
bp->s_clicked();
bp->click();
}
}
\ No newline at end of file
......@@ -41,4 +41,9 @@ namespace campvis {
ButtonProperty::~ButtonProperty() {
}
void ButtonProperty::click() {
s_clicked();
s_changed(this);
}
}
......@@ -51,6 +51,11 @@ namespace campvis {
**/
virtual ~ButtonProperty();
/**
* Simulates a click event.
*/
void click();
/**
* Signal emitted when button was clicked.
*/
......
......@@ -45,8 +45,12 @@ namespace tgt {
void OpenGLGarbageCollector::deleteGarbage() {
tbb::spin_mutex::scoped_lock lock(_mutex);
size_t backIndex = _currentFrontindex;
_currentFrontindex = (_currentFrontindex+1) % 2;
size_t backIndex;
{
backIndex = _currentFrontindex;
tbb::spin_mutex::scoped_lock lock(_addMutex);
_currentFrontindex = (_currentFrontindex+1) % 2;
}
// delete textures
if (! _texturesToDelete[backIndex].empty()) {
......
......@@ -2,12 +2,14 @@
# Source files:
FILE(GLOB ThisModSources RELATIVE ${ModulesDir}
modules/columbia/datastructures/*.cpp
modules/columbia/pipelines/*.cpp
modules/columbia/processors/*.cpp
)
# Header files (including GLSL files so that they'll appear in VS projects)
FILE(GLOB ThisModHeaders RELATIVE ${ModulesDir}
modules/columbia/datastructures/*.h
modules/columbia/glsl/*.frag
modules/columbia/glsl/*.vert
modules/columbia/pipelines/*.h
......
// ================================================================================================
//
// This file is part of the CAMPVis Software Framework.
//
// If not explicitly stated otherwise: Copyright (C) 2012, all rights reserved,
// Christian Schulte zu Berge <christian.szb@in.tum.de>
// Chair for Computer Aided Medical Procedures
// Technische Universitt Mnchen
// Boltzmannstr. 3, 85748 Garching b. Mnchen, Germany
// For a full list of authors and contributors, please refer to the file "AUTHORS.txt".
//
// The licensing of this softare is not yet resolved. Until then, redistribution in source or
// binary forms outside the CAMP chair is not permitted, unless explicitly stated in legal form.
// However, the names of the original authors and the above copyright notice must retain in its
// original state in any case.
//
// Legal disclaimer provided by the BSD license:
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
// IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
// AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// ================================================================================================
#include "fiberdata.h"
namespace campvis {
FiberData::FiberData()
: AbstractData()
{
}
FiberData::~FiberData() {
}
void FiberData::addFiber(const std::deque<tgt::vec3>& vertices) {
_vertices.insert(_vertices.end(), vertices.begin(), vertices.end());
_fibers.push_back(Fiber(_vertices.size() - vertices.size(), _vertices.size()));
}
void FiberData::addFiber(const std::vector<tgt::vec3>& vertices) {
_vertices.insert(_vertices.end(), vertices.begin(), vertices.end());
_fibers.push_back(Fiber(_vertices.size() - vertices.size(), _vertices.size()));
}
void FiberData::clear() {
_fibers.clear();
_vertices.clear();
}
void FiberData::updateLengths() const {
for (size_t i = 0; i < _fibers.size(); ++i) {
_fibers[i]._length = 0.0f;
for (size_t j = _fibers[i]._startIndex + 1; j < _fibers[i]._endIndex; ++j)
_fibers[i]._length += distance(_vertices[j-1], _vertices[j]);
}
}
size_t FiberData::numFibers() const {
return _fibers.size();
}
size_t FiberData::numSegments() const {
size_t sum = 0;
for (std::vector<Fiber>::const_iterator it = _fibers.begin(); it != _fibers.end(); ++it)
sum += (it->_endIndex - it->_startIndex);
return sum;
}
bool FiberData::empty() const {
return _fibers.empty();
}
FiberData* FiberData::clone() const {
FiberData* toReturn = new FiberData(*this);
return toReturn;
}
size_t FiberData::getLocalMemoryFootprint() const {
size_t sum = _vertices.size() * sizeof(tgt::vec3);
sum += _fibers.size() * sizeof(Fiber);
sum += sizeof(*this);
return sum;
}
size_t FiberData::getVideoMemoryFootprint() const {
return 0;
}
}
\ No newline at end of file
// ================================================================================================
//
// This file is part of the CAMPVis Software Framework.
//
// If not explicitly stated otherwise: Copyright (C) 2012, 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".
//
// The licensing of this softare is not yet resolved. Until then, redistribution in source or
// binary forms outside the CAMP chair is not permitted, unless explicitly stated in legal form.
// However, the names of the original authors and the above copyright notice must retain in its
// original state in any case.
//
// Legal disclaimer provided by the BSD license:
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
// IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
// AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// ================================================================================================
#ifndef FIBERDATA_H__
#define FIBERDATA_H__
#include "tgt/vector.h"
#include "core/datastructures/abstractdata.h"
#include <deque>
#include <vector>
namespace campvis {
/**
* Data object storing fiber data.
*/
class FiberData : public AbstractData {
public:
/**
* Struct storing meta information about a single fiber.
*/
struct Fiber {
size_t _startIndex; ///< Start index of the fiber
size_t _endIndex; ///< End index of the fiber (as in STL iterators: points to the element _behind_ the last vertex)
mutable float _length; ///< Length of the fiber (cached)
int _segmentId; ///< Label of the fiber
bool _visible; ///< Visibility flag of the fiber
bool _selected; ///< Selected flag of the fiber
Fiber(size_t startIndex, size_t endIndex)
: _startIndex(startIndex), _endIndex(endIndex), _length(0.f), _segmentId(0), _visible(true), _selected(false)
{};
};
/**
* Constructor.
*/
FiberData();
/**
* Destructor.
*/
virtual ~FiberData();
/**
* Generates a new fiber from the given vertices and adds it to this data structure.
* \param vertices Coordinates of the fiber points.
*/
void addFiber(const std::deque<tgt::vec3>& vertices);
/**
* Generates a new fiber from the given vertices and adds it to this data structure.
* \param vertices Coordinates of the fiber points.
*/
void addFiber(const std::vector<tgt::vec3>& vertices);
/**
* Clears this data structure.
*/
void clear();
/**
* Computes the lengths of each fiber in this data structure and stores it in the
* corresponding field.
* \note Since there is currently not automatism to do this, you're responsible to do
* this yourself when needed.
*/
void updateLengths() const;
/**
* Returns the number of fibers in this data structure.
* \return _fibers.size()
*/
size_t numFibers() const;
/**
* Returns the number of fiber segments in this data structure.
* \return (i.e. the number of vertices - number of fibers)
*/
size_t numSegments() const;
/**
* Returns whether this data structure is empty (i.e. has no fibers).
*/
bool empty() const;
/// \see AbstractData::clone()
virtual FiberData* clone() const;
/// \see AbstractData::getLocalMemoryFootprint()
virtual size_t getLocalMemoryFootprint() const;
/// \see AbstractData::getVideoMemoryFootprint()
virtual size_t getVideoMemoryFootprint() const;
protected:
std::vector<tgt::vec3> _vertices; ///< The fiber vertex (coordinates) data
std::vector<Fiber> _fibers; ///< The fiber meta data
};
}
#endif // FIBERDATA_H__
......@@ -48,6 +48,7 @@ namespace campvis {
, _sr(_effectiveRenderTargetSize)
, _src(_effectiveRenderTargetSize)
, _gr(_effectiveRenderTargetSize)
, _sft()
, _trackballEH(0)
{
addProperty(&_camera);
......@@ -58,10 +59,11 @@ namespace campvis {
addProcessor(&_imageReader);
addProcessor(&_vtkReader);
addProcessor(&_splitter);
addProcessor(&_vr);
addProcessor(&_src);
//addProcessor(&_vr);
//addProcessor(&_src);
addProcessor(&_sr);
addProcessor(&_gr);
addProcessor(&_sft);
}
Columbia1::~Columbia1() {
......@@ -80,10 +82,12 @@ namespace campvis {
_vr.p_outputImage.setValue("vr");
_sr.p_targetImageID.setValue("sr");
_src.p_targetImageID.setValue("src");
_renderTargetID.setValue("vr");
_renderTargetID.setValue("sr");
_imageReader.p_url.setValue("D:/Medical Data/Columbia/inputs/FullVolumeLV_3D_25Hz_[IM_0004]_NIF_diffused_crop_00.ltf");
_imageReader.p_url.setValue("D:/Medical Data/Columbia/outputs/FullVolumeLV_3D_25Hz_[IM_0004]_NIF_crop_flow_field_00_00.ltf");
_imageReader.p_size.setValue(tgt::ivec3(224, 176, 208));
_imageReader.p_numChannels.setValue(3);
_imageReader.p_baseType.selectById("float");
_imageReader.p_targetImageID.setValue("reader.output");
_imageReader.p_targetImageID.connect(&_splitter.p_inputID);
......@@ -94,7 +98,8 @@ namespace campvis {
_splitter.p_outputID.connect(&_vr.p_inputVolume);
_splitter.p_outputID.connect(&_src.p_sourceImageID);
_splitter.p_outputID.connect(&_sr.p_sourceImageID);
_splitter.p_outputID.connect(&_sft.p_strainId);
Geometry1DTransferFunction* dvrTF = new Geometry1DTransferFunction(128, tgt::vec2(0.f, 1.f));
dvrTF->addGeometry(TFGeometry1D::createQuad(tgt::vec2(.1f, .125f), tgt::col4(255, 0, 0, 32), tgt::col4(255, 0, 0, 32)));
dvrTF->addGeometry(TFGeometry1D::createQuad(tgt::vec2(.4f, .5f), tgt::col4(0, 255, 0, 128), tgt::col4(0, 255, 0, 128)));
......@@ -102,6 +107,8 @@ namespace campvis {
_gr.p_renderTargetID.setValue("gr");
_sft.p_outputID.setValue("fibers");
_trackballEH->setViewportSize(_effectiveRenderTargetSize.getValue());
_effectiveRenderTargetSize.s_changed.connect<Columbia1>(this, &Columbia1::onRenderTargetSizeChanged);
}
......
......@@ -39,6 +39,7 @@
#include "modules/io/processors/vtkimagereader.h"
#include "modules/columbia/processors/geometrystrainrenderer.h"
#include "modules/columbia/processors/imageseriessplitter.h"
#include "modules/columbia/processors/strainfibertracker.h"
#include "modules/columbia/processors/strainraycaster.h"
#include "modules/vis/processors/geometryrenderer.h"
#include "modules/vis/processors/sliceextractor.h"
......@@ -87,6 +88,8 @@ namespace campvis {
SliceExtractor _sr;
GeometryRenderer _gr;
StrainFiberTracker _sft;
TrackballNavigationEventHandler* _trackballEH;
};
......
// ================================================================================================
//
// This file is part of the CAMPVis Software Framework.
//
// If not explicitly stated otherwise: Copyright (C) 2012, all rights reserved,
// Christian Schulte zu Berge <christian.szb@in.tum.de>
// Chair for Computer Aided Medical Procedures
// Technische Universitt Mnchen
// Boltzmannstr. 3, 85748 Garching b. Mnchen, Germany
// For a full list of authors and contributors, please refer to the file "AUTHORS.txt".
//
// The licensing of this softare is not yet resolved. Until then, redistribution in source or
// binary forms outside the CAMP chair is not permitted, unless explicitly stated in legal form.
// However, the names of the original authors and the above copyright notice must retain in its
// original state in any case.
//
// Legal disclaimer provided by the BSD license:
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
// IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
// AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// ================================================================================================
#include "strainfibertracker.h"
#include "tbb/tbb.h"
#include "tbb/spin_mutex.h"
#include "core/datastructures/imagerepresentationlocal.h"
#include "modules/columbia/datastructures/fiberdata.h"
#include <deque>
namespace campvis {
const std::string StrainFiberTracker::loggerCat_ = "CAMPVis.modules.io.StrainFiberTracker";
class ApplyFiberTracking {
public:
ApplyFiberTracking(const ImageRepresentationLocal* input, const std::vector<tgt::vec3>& seeds, FiberData* output, int numSteps, float stepSize, float strainThreshold, float maximumAngle, tbb::spin_mutex& mutex)
: _input(input)
, _seeds(seeds)
, _output(output)
, _numSteps(numSteps)
, _stepSize(stepSize)
, _strainThreshold(strainThreshold * strainThreshold)
, _maxAngle(maximumAngle)
, _mutex(mutex)
{
_voxelSize = tgt::length(_input->getParent()->getMappingInformation().getVoxelSize());
}
/**
* Retrieves a vec3 from \a vol using trilinear interpolation.
* \param position voxel position
**/
inline tgt::vec3 getVec3FloatLinear(const tgt::vec3& position) const {
tgt::vec3 result;
result.x = _input->getElementNormalizedLinear(position, 0);
result.y = _input->getElementNormalizedLinear(position, 1);
result.z = _input->getElementNormalizedLinear(position, 2);
return result;
}
/**
* Checks whether the angle between \a a and \a b is lower that the given threshold.
*
* \param a direction of first tangent vector in world coordinates
* \param b direction of second tangent vector in world coordinates
*
* \return true, if angle is below the threshold.
**/
inline bool testTortuosity(const tgt::vec3& a, const tgt::vec3& b) const {
float angle = fabs(acos(tgt::dot(a, b) / (tgt::length(a) * tgt::length(b))));
return (angle < _maxAngle);
}
/**
* Checks whether \a position is within volume bounds.
*
* \param position position in voxel coordinates
*
* true, if \a position is within bounds of eigenvalues volume.
**/
inline bool testBounds(const tgt::vec3& position) const {
const tgt::svec3& dim = _input->getParent()->getSize();
tgt::svec3 pos(tgt::ceil(position));
return (tgt::hand(tgt::greaterThanEqual(pos, tgt::svec3::zero)) && tgt::hand(tgt::lessThanEqual(pos, dim)));
}
/**
* Performs fiber tracking of a single fiber to a single direction starting at \a worldPosition
* and stores path in \a result. \a result will NOT contain the start point \a worldPosition.
*
* \param worldPosition start position in world coordinates
* \param forwards flag whether to propagate forwards or backwards from \a worldPosition
* \param maxAngle maximum angle between two fiber points to continue tracking
* \param result fiber points will be stored in here - forward tracking will append, backward tracking push front
**/
void performSingleTracking(tgt::vec3 worldPosition, bool forwards, std::deque<tgt::vec3>& result) const {
const tgt::mat4& WtV = _input->getParent()->getMappingInformation().getWorldToVoxelMatrix();
tgt::vec3 voxelPosition = (WtV * tgt::vec4(worldPosition, 1.f)).xyz();
tgt::vec3 direction = getVec3FloatLinear(voxelPosition);
if (! forwards)
direction *= -1.f;
for (int i = 0; i < _numSteps; ++i) {
// apply second order runge-kutta integration (Heun method)
tgt::vec3 dir1 = getVec3FloatLinear((WtV * tgt::vec4(worldPosition, 1.f)).xyz()) * _stepSize * _voxelSize;
if (tgt::dot(direction, dir1) < 0)
dir1 *= -1.f;
tgt::vec3 dir2 = getVec3FloatLinear((WtV * tgt::vec4(worldPosition + dir1, 1.f)).xyz()) * _stepSize * _voxelSize;
if (tgt::dot(direction, dir2) < 0)
dir2 *= -1.f;
tgt::vec3 vProp = (dir1 + dir2) * .5f;
worldPosition += vProp;
voxelPosition = (WtV * tgt::vec4(worldPosition, 1.f)).xyz();
// check termination criteria
if (tgt::lengthSq(vProp) < _strainThreshold || !testBounds(voxelPosition) || !testTortuosity(direction, vProp))
break;
direction = vProp;
if (forwards)
result.push_back(worldPosition);
else
result.push_front(worldPosition);
}
}
void operator() (const tbb::blocked_range<size_t>& range) const {
for (size_t i = range.begin(); i != range.end(); ++i) {
const tgt::vec3& position = _seeds[i];
// perform fiber tracking in both directions
std::deque<tgt::vec3> vertices;
performSingleTracking(position, false, vertices);
vertices.push_back(position);
performSingleTracking(position, true, vertices);
if (vertices.size() > 1) {
tbb::spin_mutex::scoped_lock lock(_mutex);
_output->addFiber(vertices);
}
}
}
protected:
const ImageRepresentationLocal* _input;
const std::vector<tgt::vec3>& _seeds;
FiberData* _output;
int _numSteps;
float _stepSize;
float _voxelSize;
float _strainThreshold;
float _maxAngle;
tbb::spin_mutex& _mutex;
};
StrainFiberTracker::StrainFiberTracker()
: AbstractProcessor()
, p_strainId("StrainId", "Input Strain Data", "input", DataNameProperty::READ, AbstractProcessor::VALID)
, p_outputID("OutputId", "Output Fiber Data", "output", DataNameProperty::WRITE, AbstractProcessor::VALID)
, p_updateButton("UpdateButton", "Perform Tracking")
, p_numSteps("NumSteps", "Maximum Number of Steps", 256, 16, 1024, AbstractProcessor::VALID)
, p_stepSize("StepSize", "Base Step Size", 1.f, .01f, 10.f, AbstractProcessor::VALID)
, p_strainThreshold("StrainThreshold", "Local Strain Threshold", .5f, .1f, 1.f, AbstractProcessor::VALID)
, p_maximumAngle("MaximumAngle", "Maxium Angle", 25.f, 0.f, 100.f, AbstractProcessor::VALID)
{
addProperty(&p_strainId);
addProperty(&p_outputID);
addProperty(&p_updateButton);
addProperty(&p_numSteps);
addProperty(&p_stepSize);
addProperty(&p_strainThreshold);
addProperty(&p_maximumAngle);
}
StrainFiberTracker::~StrainFiberTracker() {
}
void StrainFiberTracker::process(DataContainer& data) {
ImageRepresentationLocal::ScopedRepresentation strainData(data, p_strainId.getValue());
if (strainData != 0) {
if (strainData.getImageData()->getNumChannels() == 3) {
LDEBUG("Generating seeds...");
std::vector<tgt::vec3> seeds = performUniformSeeding(*strainData);
LDEBUG("Generating fibers...");
FiberData* fibers = new FiberData();
tbb::spin_mutex mutex;
tbb::parallel_for(
tbb::blocked_range<size_t>(0, seeds.size()),
ApplyFiberTracking(strainData, seeds, fibers, p_numSteps.getValue(), p_stepSize.getValue(), p_strainThreshold.getValue(), p_maximumAngle.getValue(), mutex));
LDEBUG("done.");
data.addData(p_outputID.getValue(), fibers);
p_outputID.issueWrite();
}
else {
LERROR("Wrong number of channels.");
}
}
else {
LERROR("No input data.");
}
validate(INVALID_RESULT);
}
std::vector<tgt::vec3> StrainFiberTracker::performUniformSeeding(const ImageRepresentationLocal& strainData) const {
std::vector<tgt::vec3> seeds;
const tgt::mat4& VtW = strainData.getParent()->getMappingInformation().getVoxelToWorldMatrix();
float threshold = p_strainThreshold.getValue() * p_strainThreshold.getValue();
for (size_t z = 0; z < strainData.getSize().z; ++z) {
for (size_t y = 0; y < strainData.getSize().y; ++y) {