* added ConcurrentGenericHistogramND

 * some work on LHHistogram processor

git-svn-id: https://camplinux.in.tum.de/svn/campvis/trunk@307 bb408c1c-ae56-11e1-83d9-df6b3e0c105e

src/tumvis/core/tools/concurrenthistogram.h

+
+#include "tgt/tgt_math.h"
+#include "tbb/include/tbb/atomic.h"
+
+namespace TUMVis {
+
+    template<typename T, size_t ND>
+    class ConcurrentGenericHistogramND {
+    public:
+        ConcurrentGenericHistogramND(T mins[ND], T maxs[ND], size_t numBuckets[ND]);
+
+        virtual ~ConcurrentGenericHistogramND();
+
+        void addSample(T sample1[ND]);
+
+        const tbb::atomic<size_t>* getBuckets() const { return _buckets; };
+
+        size_t getNumSamples() const { return _numSamples; };
+
+    protected:
+        size_t getBucketNumber(size_t dimension, T sample) const;
+
+        size_t getArrayIndex(size_t* bucketNumbers) const;
+
+        T _min[ND];
+        T _max[ND];
+        size_t _numBuckets[ND];
+        size_t _arraySize;
+        tbb::atomic<size_t>* _buckets;
+        tbb::atomic<size_t> _numSamples;
+    };
+
+// ================================================================================================
+
+    template<typename T, size_t ND>
+    TUMVis::ConcurrentGenericHistogramND<T, ND>::ConcurrentGenericHistogramND(T mins[ND], T maxs[ND], size_t numBuckets[ND]) 
+        : _arraySize(1)
+        , _buckets(0)
+    {
+        for (size_t i = 0; i < ND; ++i) {
+            _min[i] = mins[i];
+            _max[i] = maxs[i];
+            tgtAssert(_min[i] < _max[i], "Min must be smaller than Max!");
+            _numBuckets[i] = numBuckets[i];
+            _arraySize *= _numBuckets[i];
+        }
+
+        _buckets = new tbb::atomic<size_t>[_arraySize];
+        for (size_t i = 0; i < _arraySize; ++i)
+            _buckets[i] = 0;
+    }
+
+    template<typename T, size_t ND>
+    TUMVis::ConcurrentGenericHistogramND<T, ND>::~ConcurrentGenericHistogramND() {
+        delete [] _buckets;
+    }
+
+    template<typename T, size_t ND>
+    void TUMVis::ConcurrentGenericHistogramND<T, ND>::addSample(T sample1[ND]) {
+        size_t bucketNumbers[ND];
+        for(int i = 0; i < ND; ++i)
+            bucketNumbers[i] = getBucketNumber(i, sample1[i]);
+
+        size_t index = getArrayIndex(bucketNumbers);
+        ++(_buckets[index]);
+        ++_numSamples;
+    }
+
+
+    template<typename T, size_t ND>
+    size_t TUMVis::ConcurrentGenericHistogramND<T, ND>::getBucketNumber(size_t dimension, T sample) const {
+        double ratio = static_cast<double>(sample - _min[dimension]) / static_cast<double>(_max[dimension] - _min[dimension]);
+        return static_cast<size_t>(tgt::clamp(static_cast<int>(ratio * _numBuckets[dimension]), static_cast<int>(0), static_cast<int>(_numBuckets[dimension] - 1)));
+    }
+
+    template<typename T, size_t ND>
+    size_t TUMVis::ConcurrentGenericHistogramND<T, ND>::getArrayIndex(size_t* bucketNumbers) const {
+        size_t index = 0;
+        size_t multiplier = 1;
+        for (size_t i = 0; i < ND; ++i) {
+            index += multiplier * bucketNumbers[i];
+            multiplier *= _numBuckets[i];
+        }
+        return index;
+    }
+
+
+}

src/tumvis/core/tools/weaklytypedpointer.h

@@ -93,7 +93,7 @@ namespace TUMVis {
          * Constructs a new weakly typed pointer.
          * \param pt            Base data type of the pointer.
          * \param numChannels   Number of channels, must be in [1, 4]!
-         * \param ptr           Pointer to the data, WeaklyTypedPointer will take ownership of it.
+         * \param ptr           Pointer to the data, WeaklyTypedPointer will _not_ take ownership of it.
          */
         WeaklyTypedPointer(BaseType pt, size_t numChannels, void* ptr);
 

src/tumvis/ext/tgt/texture.cpp

@@ -147,6 +147,7 @@ int Texture::calcBpp(GLint internalformat) {
         case GL_GREEN:
         case GL_BLUE:
         case GL_ALPHA:
+        case GL_ALPHA8:
         case GL_INTENSITY:
         case GL_LUMINANCE:
         case GL_DEPTH_COMPONENT:

src/tumvis/modules/classification/lhhistogram.cpp

@@ -36,21 +36,44 @@
 #include "tbb/include/tbb/tbb.h"
 
 #include "core/datastructures/genericimagedatalocal.h"
+#include "core/datastructures/imagedatagl.h"
+#include "core/tools/concurrenthistogram.h"
 
 namespace TUMVis {
     class LHGenerator {
     public:
-        LHGenerator(const ImageDataLocal* intensities, const GenericImageDataLocal<float, 4>* gradients, ImageDataLocal* fh, ImageDataLocal* fl, float epsilon)
+        LHGenerator(const ImageDataLocal* intensities, const GenericImageDataLocal<float, 4>* gradients, ImageDataLocal* fl, ImageDataLocal* fh, float epsilon)
             : _intensities(intensities)
             , _gradients(gradients)
-            , _fh(fh)
             , _fl(fl)
+            , _fh(fh)
             , _epsilon(epsilon)
         {
             tgtAssert(_intensities->getDimensionality() == _gradients->getDimensionality(), "Dimensionality of intensities volumes must match!");
             tgtAssert(_intensities->getSize() == _gradients->getSize(), "Size of intensities volumes must match!");
         }
 
+        void operator() (const tbb::blocked_range<size_t>& range) const {
+            for (size_t i = range.begin(); i != range.end(); ++i) {
+                tgt::svec3 pos = _intensities->indexToPosition(i);
+                const tgt::svec3& size = _intensities->getSize();
+
+                const tgt::vec4& gradient = _gradients->getElement(i);
+                float fl = _intensities->getElementNormalized(pos, 0);
+                float fh = fl;
+
+                float forwardIntensity = integrateHeun(tgt::vec3(static_cast<float>(pos.x), static_cast<float>(pos.y), static_cast<float>(pos.z)), gradient);
+                float backwardIntensity = integrateHeun(tgt::vec3(static_cast<float>(pos.x), static_cast<float>(pos.y), static_cast<float>(pos.z)), gradient * -1.f);
+
+                fh = std::max(forwardIntensity, backwardIntensity);
+                fl = std::min(forwardIntensity, backwardIntensity);
+
+                _fl->setElementNormalized(pos, 0, fl);
+                _fh->setElementNormalized(pos, 0, fh);
+            }
+        }
+
+    protected:
         tgt::vec4 getGradientLinear(const tgt::vec3& position) const {
             tgt::vec4 result;
             result.x = _gradients->getElementNormalizedLinear(position, 0);
@@ -64,45 +87,53 @@ namespace TUMVis {
             tgt::vec4 gradient1 = direction;
             tgt::vec3 stepSize(1.f);
             tgt::vec3 size(_intensities->getSize());
+            size_t numSteps = 0;
 
             while (gradient1.w > _epsilon) {
                 tgt::vec4 gradient2 = getGradientLinear(position + tgt::normalize(gradient1.xyz()) * stepSize);
                 position += tgt::normalize((gradient1 + gradient2).xyz()) * stepSize;
                 gradient1 = getGradientLinear(position);
+                ++numSteps;
 
-                if (tgt::hor(tgt::lessThan(position, tgt::vec3::zero)) || tgt::hor(tgt::greaterThan(position, size)))
+                if (numSteps > 128 || tgt::hor(tgt::lessThan(position, tgt::vec3::zero)) || tgt::hor(tgt::greaterThan(position, size)))
                     break;
             }
 
             return _intensities->getElementNormalizedLinear(position, 0);;
         }
 
-        void operator() (const tbb::blocked_range<size_t>& range) const {
-            for (size_t i = range.begin(); i != range.end(); ++i) {
-                tgt::svec3 pos = _intensities->indexToPosition(i);
-                const tgt::svec3& size = _intensities->getSize();
+        const ImageDataLocal* _intensities;
+        const GenericImageDataLocal<float, 4>* _gradients;
+        ImageDataLocal* _fl;
+        ImageDataLocal* _fh;
+        float _epsilon;
+    };
 
-                const tgt::vec4& gradient = _gradients->getElement(i);
-                float fl = _intensities->getElementNormalized(pos, 0);
-                float fh = fl;
+// ================================================================================================
 
-                float forwardIntensity = integrateHeun(tgt::vec3(static_cast<float>(pos.x), static_cast<float>(pos.y), static_cast<float>(pos.z)), gradient);
-                float backwardIntensity = integrateHeun(tgt::vec3(static_cast<float>(pos.x), static_cast<float>(pos.y), static_cast<float>(pos.z)), gradient * -1.f);
+    class LHHistogramGenerator {
+    public:
+        LHHistogramGenerator(const ImageDataLocal* fl, const ImageDataLocal* fh, ConcurrentGenericHistogramND<float, 2>* histogram)
+            : _fl(fl)
+            , _fh(fh)
+            , _histogram(histogram)
+        {
+            tgtAssert(_fh->getDimensionality() == _fl->getDimensionality(), "Dimensionality of input volumes must match!");
+            tgtAssert(_fh->getSize() == _fl->getSize(), "Size of input volumes must match!");
+        }
 
-                fh = std::max(forwardIntensity, backwardIntensity);
-                fl = std::min(forwardIntensity, backwardIntensity);
 
-                _fl->setElementNormalized(pos, 0, fl);
-                _fh->setElementNormalized(pos, 0, fh);
-            }
+    void operator() (const tbb::blocked_range<size_t>& range) const {
+        for (size_t i = range.begin(); i != range.end(); ++i) {
+            tgt::svec3 pos = _fl->indexToPosition(i);
+            float values[2] = { _fl->getElementNormalized(pos, 0), _fh->getElementNormalized(pos, 0) };
+            _histogram->addSample(values);
         }
-
+    }
     protected:
-        const ImageDataLocal* _intensities;
-        const GenericImageDataLocal<float, 4>* _gradients;
-        ImageDataLocal* _fh;
-        ImageDataLocal* _fl;
-        float _epsilon;
+        const ImageDataLocal* _fl;
+        const ImageDataLocal* _fh;
+        ConcurrentGenericHistogramND<float, 2>* _histogram;
     };
 
 // ================================================================================================
@@ -113,13 +144,13 @@ namespace TUMVis {
         : AbstractProcessor()
         , _inputVolume("InputVolume", "Input Volume ID", "volume", DataNameProperty::READ)
         , _inputGradients("InputGradients", "Input Gradient Volume ID", "gradients", DataNameProperty::READ)
-        , _outputFH("OutputFH", "FH Output Volume", "fh", DataNameProperty::WRITE)
         , _outputFL("OutputFL", "FL Output Volume", "fl", DataNameProperty::WRITE)
+        , _outputFH("OutputFH", "FH Output Volume", "fh", DataNameProperty::WRITE)
     {
         addProperty(&_inputVolume);
         addProperty(&_inputGradients);
-        addProperty(&_outputFH);
         addProperty(&_outputFL);
+        addProperty(&_outputFH);
     }
 
     LHHistogram::~LHHistogram() {
@@ -131,10 +162,25 @@ namespace TUMVis {
         DataContainer::ScopedTypedData< GenericImageDataLocal<float, 4> > gradients(data, _inputGradients.getValue());
 
         if (intensities != 0 && gradients != 0) {
-            ImageDataLocal* fh = intensities->clone();
             ImageDataLocal* fl = intensities->clone();
-            tbb::parallel_for(tbb::blocked_range<size_t>(0, intensities->getNumElements()), LHGenerator(intensities, gradients, fh, fl, .01f));
-
+            ImageDataLocal* fh = intensities->clone();
+            tbb::parallel_for(tbb::blocked_range<size_t>(0, intensities->getNumElements()), LHGenerator(intensities, gradients, fl, fh, .005f));
+
+            float mins[2] = { 0.f, 0.f };
+            float maxs[2] = { .01f, .01f };
+            size_t numBuckets[2] = { 256, 256 };
+            ConcurrentGenericHistogramND<float, 2> lhHistogram(mins, maxs, numBuckets);
+            tbb::parallel_for(tbb::blocked_range<size_t>(0, intensities->getNumElements()), LHHistogramGenerator(fl, fh, &lhHistogram));
+
+            // TODO: ugly hack...
+            int16_t* tmp = new int16_t[256*256];
+            for (size_t i = 0; i < 256*256; ++i)
+                tmp[i] = static_cast<int16_t>(lhHistogram.getBuckets()[i]);
+            WeaklyTypedPointer wtp(WeaklyTypedPointer::INT16, 1, tmp);
+            ImageDataGL* tex = new ImageDataGL(2, tgt::svec3(256, 256, 1), wtp);
+            delete [] tmp;
+
+            data.addData("foo", tex);
             data.addData(_outputFH.getValue(), fh);
             data.addData(_outputFL.getValue(), fl);
             _outputFH.issueWrite();

src/tumvis/modules/classification/lhhistogram.h

@@ -64,8 +64,8 @@ namespace TUMVis {
         DataNameProperty _inputVolume;      ///< ID for input volume
         DataNameProperty _inputGradients;   ///< ID for input gradient volume
 
-        DataNameProperty _outputFH;         ///< ID for output FH volume
         DataNameProperty _outputFL;         ///< ID for output FL volume
+        DataNameProperty _outputFH;         ///< ID for output FH volume
 
     protected:
 

src/tumvis/modules/pipelines/slicevis.cpp

@@ -54,7 +54,7 @@ namespace TUMVis {
     void SliceVis::init() {
         VisualizationPipeline::init();
 
-        _imageReader._url.setValue("D:\\Medical Data\\smallHeart.mhd");
+        _imageReader._url.setValue("D:\\Medical Data\\Dentalscan\\dental.mhd");
         _imageReader._targetImageID.setValue("reader.output");
 
         _gvg._inputVolume.setValue("se.input");
@@ -95,7 +95,7 @@ namespace TUMVis {
             executeProcessor(&_gvg);
         }
         if (! _lhh.getInvalidationLevel().isValid()) {
-            executeProcessor(&_lhh);
+            lockGLContextAndExecuteProcessor(&_lhh);
         }
         if (! _sliceExtractor.getInvalidationLevel().isValid()) {
             lockGLContextAndExecuteProcessor(&_sliceExtractor);