Skip to content
GitLab
Commit 634d2db1 authored by Christian Schulte zu Berge's avatar Christian Schulte zu Berge
Browse files

Further work on TensorAnalyzer: added computation of various anisotropy measures, yet untested.

parent 1d6510ae
Hide whitespace changes
Inline Side-by-side
modules/tensor/processors/tensoranalyzer.cpp
View file @ 634d2db1
......@@ -32,7 +32,7 @@
#include <tbb/atomic.h>
#include "core/datastructures/imagedata.h"
#include "core/tools/stringutils.h"
namespace campvis {
......@@ -43,8 +43,32 @@ namespace campvis {
GenericOption<TensorAnalyzer::DegeneratedEvHandling>("shift", "Shift", TensorAnalyzer::SHIFT)
};
static const GenericOption<std::string> measurementOptions[15] = {
GenericOption<std::string>("Disabled", "Disabled"),
GenericOption<std::string>("EigenValue1", "Eigenvalue 1"),
GenericOption<std::string>("EigenValue2", "Eigenvalue 2"),
GenericOption<std::string>("EigenValue3", "Eigenvalue 3"),
GenericOption<std::string>("MainEigenvector", "Main Eigenvector"),
GenericOption<std::string>("VolumeRatio", "Volume Ratio"),
GenericOption<std::string>("FractionalAnisotropy", "Fractional Anisotropy"),
GenericOption<std::string>("RelativeAnisotropy", "Relative Anisotropy"),
GenericOption<std::string>("MeanDiffusivity", "Mean Diffusivity"),
GenericOption<std::string>("Trace", "Trace"),
GenericOption<std::string>("AxialDiffusivity", "Axial Diffusivity"),
GenericOption<std::string>("RadialDiffusivity", "Radial Diffusivity"),
GenericOption<std::string>("LinearAnisotropy", "Linear Anisotropy"),
GenericOption<std::string>("PlanarAnisotropy", "Planar Anisotropy"),
GenericOption<std::string>("Isotropy", "Isotropy")
};
const std::string TensorAnalyzer::loggerCat_ = "CAMPVis.modules.classification.TensorAnalyzer";
TensorAnalyzer::OutputPropertyPair::OutputPropertyPair(size_t index)
: _imageId("OutputId" + StringUtils::toString(index), "Output " + StringUtils::toString(index) + " Image", "TensorAnalyzer.output" + StringUtils::toString(index), DataNameProperty::WRITE)
, _imageType("OutputType" + StringUtils::toString(index), "Output " + StringUtils::toString(index) + " Image Type", measurementOptions, 15)
{}
TensorAnalyzer::TensorAnalyzer()
: AbstractProcessor()
, p_inputImage("InputImage", "Input Tensor Image", "tensors", DataNameProperty::READ, AbstractProcessor::INVALID_RESULT | EIGENSYSTEM_INVALID)
......@@ -52,6 +76,7 @@ namespace campvis {
, p_evecsImage("EvecsImage", "Output Eigenvectors Image", "eigenvectors", DataNameProperty::WRITE)
, p_degeneratedHandling("DegeneratedHandling", "Handling of Degenerated Tensors", handlingModes, 4)
, p_maskMixedTensors("MaskMixedTensors", "Mask Mixed Tensors", true)
, p_addOutputButton("AddOutputButton", "Add Output", AbstractProcessor::VALID)
, _eigenvalues(0)
, _eigenvectors(0)
{
......@@ -60,6 +85,10 @@ namespace campvis {
addProperty(&p_evecsImage);
addProperty(&p_degeneratedHandling);
addProperty(&p_maskMixedTensors);
addProperty(&p_addOutputButton);
addOutput();
p_addOutputButton.s_clicked.connect(this, &TensorAnalyzer::addOutput);
}
TensorAnalyzer::~TensorAnalyzer() {
......@@ -72,7 +101,9 @@ namespace campvis {
}
if (_eigenvalues.getData() != 0 && _eigenvectors.getData() != 0) {
for (size_t i = 0; i < p_outputProperties.size(); ++i) {
computeOutput(data, i);
}
}
else {
LERROR("Could not compute Eigensystem");
......@@ -195,4 +226,253 @@ namespace campvis {
validate(EIGENSYSTEM_INVALID);
}
void TensorAnalyzer::computeOutput(DataContainer& data, size_t index) {
// sanity check
if (index >= p_outputProperties.size()) {
LERROR("Index out of bounds while computing output #" << index);
return;
}
// gather eigensystem
const GenericImageRepresentationLocal<float, 3>* evalRep = 0;
const GenericImageRepresentationLocal<float, 9>* evecRep = 0;
if (_eigenvalues.getData() != 0 && _eigenvectors.getData() != 0) {
evalRep = static_cast<const ImageData*>(_eigenvalues.getData())->getRepresentation< GenericImageRepresentationLocal<float, 3> >(false);
evecRep = static_cast<const ImageData*>(_eigenvectors.getData())->getRepresentation< GenericImageRepresentationLocal<float, 9> >(false);
}
if (evalRep == 0 || evecRep == 0) {
LERROR("Could not compute output, no eigensystem present.");
return;
}
OutputPropertyPair* opp = p_outputProperties[index];
const std::string& type = opp->_imageType.getOptionValue();
if (type == "Disabled") {
return;
}
else if (type == "EigenValue1") {
ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
for (size_t i = range.begin(); i != range.end(); ++i) {
const tgt::vec3& vals = evalRep->getElement(i);
float ev = vals.x;
if (ev == 0.0f || tgt::isNaN(ev))
output->setElement(i, 0.0f);
else
output->setElement(i, ev);
}
});
data.addData(opp->_imageId.getValue(), id);
}
else if (type == "EigenValue2") {
ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
for (size_t i = range.begin(); i != range.end(); ++i) {
const tgt::vec3& vals = evalRep->getElement(i);
float ev = vals.y;
if (ev == 0.0f || tgt::isNaN(ev))
output->setElement(i, 0.0f);
else
output->setElement(i, ev);
}
});
data.addData(opp->_imageId.getValue(), id);
}
else if (type == "EigenValue3") {
ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
for (size_t i = range.begin(); i != range.end(); ++i) {
const tgt::vec3& vals = evalRep->getElement(i);
float ev = vals.z;
if (ev == 0.0f || tgt::isNaN(ev))
output->setElement(i, 0.0f);
else
output->setElement(i, ev);
}
});
data.addData(opp->_imageId.getValue(), id);
}
else if (type == "MainEigenvector") {
ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 3);
GenericImageRepresentationLocal<float, 3>* output = GenericImageRepresentationLocal<float, 3>::create(id, 0);
tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
for (size_t i = range.begin(); i != range.end(); ++i) {
output->setElement(i, evecRep->getElement(i)[0]);
}
});
data.addData(opp->_imageId.getValue(), id);
}
// = Anisotropy Measures: =========================================================================
else if (type == "VolumeRatio") {
ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
for (size_t i = range.begin(); i != range.end(); ++i) {
const tgt::vec3& vals = evalRep->getElement(i);
if (vals == tgt::vec3::zero || tgt::isNaN(vals))
output->getElement(i) = 0;
else
output->getElement(i) = (vals.x * vals.y * vals.z) / pow((vals.x + vals.y + vals.z)/3.f, 3);
}
});
data.addData(opp->_imageId.getValue(), id);
}
else if (type == "FractionalAnisotropy") {
ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
for (size_t i = range.begin(); i != range.end(); ++i) {
const tgt::vec3& vals = evalRep->getElement(i);
const float root = sqrt(.5f);
if (vals == tgt::vec3::zero || tgt::isNaN(vals))
output->getElement(i) = 0;
else
output->getElement(i) = root * sqrt((vals.x-vals.y)*(vals.x-vals.y) + (vals.y-vals.z)*(vals.y-vals.z) + (vals.z-vals.x)*(vals.z-vals.x)) / sqrt(vals.x*vals.x + vals.y*vals.y + vals.z*vals.z);
}
});
data.addData(opp->_imageId.getValue(), id);
}
else if (type == "RelativeAnisotropy") {
ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
for (size_t i = range.begin(); i != range.end(); ++i) {
const tgt::vec3& vals = evalRep->getElement(i);
const float root = sqrt(.5f);
if (vals == tgt::vec3::zero || tgt::isNaN(vals))
output->getElement(i) = 0;
else
output->getElement(i) = root * sqrt((vals.x-vals.y)*(vals.x-vals.y) + (vals.y-vals.z)*(vals.y-vals.z) + (vals.z-vals.x)*(vals.z-vals.x)) / (vals.x + vals.y + vals.z);
}
});
data.addData(opp->_imageId.getValue(), id);
}
else if (type == "MeanDiffusivity") {
ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
for (size_t i = range.begin(); i != range.end(); ++i) {
const tgt::vec3& vals = evalRep->getElement(i);
if (vals == tgt::vec3::zero || tgt::isNaN(vals))
output->getElement(i) = 0;
else
output->getElement(i) = (vals.x + vals.y + vals.z) / 3;
}
});
data.addData(opp->_imageId.getValue(), id);
}
else if (type == "Trace") {
ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
for (size_t i = range.begin(); i != range.end(); ++i) {
const tgt::vec3& vals = evalRep->getElement(i);
if (vals == tgt::vec3::zero || tgt::isNaN(vals))
output->getElement(i) = 0;
else
output->getElement(i) = vals.x + vals.y + vals.z;
}
});
data.addData(opp->_imageId.getValue(), id);
}
else if (type == "AxialDiffusivity") {
ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
for (size_t i = range.begin(); i != range.end(); ++i) {
const tgt::vec3& vals = evalRep->getElement(i);
if (vals == tgt::vec3::zero || tgt::isNaN(vals))
output->getElement(i) = 0;
else
output->getElement(i) = vals.x;
}
});
data.addData(opp->_imageId.getValue(), id);
}
else if (type == "RadialDiffusivity") {
ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
for (size_t i = range.begin(); i != range.end(); ++i) {
const tgt::vec3& vals = evalRep->getElement(i);
if (vals == tgt::vec3::zero || tgt::isNaN(vals))
output->getElement(i) = 0;
else
output->getElement(i) = (vals.y + vals.z) / 2;
}
});
data.addData(opp->_imageId.getValue(), id);
}
else if (type == "LinearAnisotropy") {
ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
for (size_t i = range.begin(); i != range.end(); ++i) {
const tgt::vec3& vals = evalRep->getElement(i);
if (vals == tgt::vec3::zero || tgt::isNaN(vals))
output->getElement(i) = 0;
else
output->getElement(i) = (vals.x - vals.y) / (vals.x + vals.y + vals.z);
}
});
data.addData(opp->_imageId.getValue(), id);
}
else if (type == "PlanarAnisotropy") {
ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
for (size_t i = range.begin(); i != range.end(); ++i) {
const tgt::vec3& vals = evalRep->getElement(i);
if (vals == tgt::vec3::zero || tgt::isNaN(vals))
output->getElement(i) = 0;
else
output->getElement(i) = 2.f*(vals.y - vals.z) / (vals.x + vals.y + vals.z);
}
});
data.addData(opp->_imageId.getValue(), id);
}
else if (type == "Isotropy") {
ImageData* id = new ImageData(evalRep->getDimensionality(), evalRep->getSize(), 1);
GenericImageRepresentationLocal<float, 1>* output = GenericImageRepresentationLocal<float, 1>::create(id, 0);
tbb::parallel_for(tbb::blocked_range<size_t>(0, id->getNumElements()), [&] (const tbb::blocked_range<size_t>& range) {
for (size_t i = range.begin(); i != range.end(); ++i) {
const tgt::vec3& vals = evalRep->getElement(i);
if (vals == tgt::vec3::zero || tgt::isNaN(vals))
output->getElement(i) = 0;
else
output->getElement(i) = (3.f*vals.z) / (vals.x + vals.y + vals.z);
}
});
data.addData(opp->_imageId.getValue(), id);
}
}
void TensorAnalyzer::addOutput() {
OutputPropertyPair* opp = new OutputPropertyPair(p_outputProperties.size() + 1);
p_outputProperties.push_back(opp);
addProperty(&opp->_imageId);
addProperty(&opp->_imageType);
}
}
modules/tensor/processors/tensoranalyzer.h
View file @ 634d2db1
......@@ -29,6 +29,7 @@
#include "core/pipeline/visualizationprocessor.h"
#include "core/properties/datanameproperty.h"
#include "core/properties/buttonproperty.h"
#include "core/properties/genericproperty.h"
#include "core/properties/floatingpointproperty.h"
#include "core/properties/numericproperty.h"
......@@ -51,6 +52,14 @@ namespace campvis {
NONE, MASK, INVERT, SHIFT
};
// Pair of DataNameProperty for output image ID and OptionProperty for image type
struct OutputPropertyPair {
OutputPropertyPair(size_t index);;
DataNameProperty _imageId;
GenericOptionProperty<std::string> _imageType;
};
/**
* Constructs a new TensorAnalyzer Processor
**/
......@@ -79,6 +88,9 @@ namespace campvis {
GenericOptionProperty<DegeneratedEvHandling> p_degeneratedHandling; ///< Handling of degenerated tensors
BoolProperty p_maskMixedTensors;
ButtonProperty p_addOutputButton;
std::vector<OutputPropertyPair*> p_outputProperties;
protected:
/**
......@@ -88,6 +100,18 @@ namespace campvis {
*/
void computeEigensystem(DataContainer& data);
/**
* Computes the derived measurement for output number \a index.
* \param data DataContainer to store output image in.
* \param index Index of output to compute.
*/
void computeOutput(DataContainer& data, size_t index);
/**
* Adds another output for this processor (i.e. adds another OutputPropertyPair).
*/
void addOutput();
DataHandle _eigenvalues; ///< Current eigenvalues cached
DataHandle _eigenvectors; ///< Current eigenvectors cached
......
Supports Markdown
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

LRZ Homepage | Datenschutz | Dokumentation und Betriebsbedingungen | Impressum