[CHANGE] - histogram euqalizer now makes use of value clamping

src/modules/histogramEqualization/WMHistogramEqualization.cpp

@@ -26,6 +26,8 @@
 
 #include "../../kernel/WKernel.h"
 #include "../../common/WPropertyHelper.h"
+#include "../../dataHandler/WDataHandler.h"
+#include "../../dataHandler/WDataTexture3D.h"
 
 #include "WMHistogramEqualization.h"
 #include "WMHistogramEqualization.xpm"
@@ -87,9 +89,18 @@ void WMHistogramEqualization::properties()
 {
     m_propCondition = boost::shared_ptr< WCondition >( new WCondition() );
 
-    // enable histogram equalization?
-    m_equalize = m_properties->addProperty( "Equalize Histogram", "If true, the dataset's cumulative histogram gets linearized.",
-                                            true, m_propCondition );
+    m_histogramResolution = m_properties->addProperty( "Histogram Resolution", "How many buckets should be used for the initial data histogram?",
+                                                       10000, m_propCondition );
+    m_histogramResolution->setMin( 10 );
+    m_histogramResolution->setMax( 1000000 );
+
+    m_clampPerc = m_properties->addProperty( "Min-Max clamping in %", "Percent that are clamped from the beginning and the end of the histogram.",
+                                             10, m_propCondition );
+    m_clampPerc->setMin( 0 );
+    m_clampPerc->setMax( 100 );
+
+    m_equalize  = m_properties->addProperty( "Equalize Histogram", "If true, the dataset's cumulative histogram gets linearized.",
+                                             true, m_propCondition );
 
     // call WModule's initialization
     WModule::properties();
@@ -122,20 +133,126 @@ void WMHistogramEqualization::moduleMain()
 
         // Remember the above criteria. We now need to check if the data is valid. After a connect-update, it might be NULL.
         boost::shared_ptr< WDataSetScalar > dataSet = m_input->getData();
+        boost::shared_ptr< WValueSetBase > valueSet = dataSet->getValueSet();
         dataUpdated = dataUpdated && dataSet;
 
-        // this ensures the below code only gets executed upon data or property change
-        if ( !dataUpdated && !m_equalize->changed() )
+        // The data is valid and we received an update. The data is not NULL but may be the same as in previous loops.
+        size_t histRes = m_histogramResolution->get( true );
+        debugLog() << "Calculating histogram with resolution " << histRes;
+
+        // Grab the histogram whose modus (interval with most of the action) is used as interval for histogram equalization
+        boost::shared_ptr< const WValueSetHistogram > hist = dataSet->getHistogram( histRes );
+
+        // find interval borders and remove first and last p%
+        double lower = hist->getMinimum();
+        double upper = hist->getMaximum();
+
+        size_t perc = m_clampPerc->get( true );
+
+        // should the histogram be clamped before processing?
+        if ( perc > 0 )
         {
-            continue;
+            debugLog() << "Clamping histogram";
+
+            size_t accumL = 0;  // accumulation of values from below
+            size_t accumU = 0;  // accumulation of values from upper side
+            size_t accumMax = hist->getTotalElementCount() * perc / 100;
+            bool foundL = false;
+            bool foundU = false;
+            size_t curI = 0;
+            size_t maxI = hist->size() - 1; // the largest index in hist
+            while ( curI <= maxI )
+            {
+                accumL += ( *hist )[ curI ];
+                accumU += ( *hist )[ maxI - curI ];
+
+                // lower border found?
+                if ( !foundL && ( accumL >= accumMax ) )
+                {
+                    foundL = true;
+                    lower = hist->getIntervalForIndex( curI ).first;
+                }
+                if ( !foundU && ( accumU >= accumMax ) )
+                {
+                    foundU = true;
+                    upper = hist->getIntervalForIndex( maxI - curI ).second;
+                }
+
+                curI++;
+            }
+
+            debugLog() << "Clamped " << perc << "% resulting in new interval ["<< lower << ", " << upper <<").";
+
+            // with this new interval, extract a new histogram and use it for equalization
+            hist = boost::shared_ptr< const WValueSetHistogram >( new WValueSetHistogram( valueSet, lower, upper, histRes ) );
         }
 
-        // The data is valid and we received an update. The data is not NULL but may be the same as in previous loops.
-        debugLog() << "Recalculating";
+        // the new data
+        std::vector< unsigned char > newData;
+        newData.resize( hist->getTotalElementCount(), 0 );
+
+        // equalize?
+        if ( m_equalize->get( true ) )
+        {
+            // calculate the cumulative distribution function
+            debugLog() << "Calculating cumulative distribution function";
+            std::vector< double > cdf;
+            cdf.resize( hist->size() );
+
+            // go through each CDF item and fill it, which is the sum of all previous items in hist
+            size_t accum = 0;
+            for( size_t i = 0; i < hist->size(); ++i )
+            {
+                // the CDF at i is the summed up histogram from 0 to i
+                // we additionally require the histogram to be normalized so divide by the total count
+                accum += ( *hist )[ i ];
+                cdf[ i ] = accum;
+            }
+
+            // finally, build the new dataset
+            debugLog() << "Calculating equalized value-set";
+            double cdfMin = cdf[ 0 ];
+            for ( size_t vi = 0; vi < valueSet->rawSize(); ++vi )
+            {
+                double cdfVI = cdf[ hist->getIndexForValue( valueSet->getScalarDouble( vi ) ) ];
+
+                newData[ vi ] = static_cast< unsigned char >( 255.0 * ( cdfVI - cdfMin ) / ( double( valueSet->rawSize() ) - cdfMin ) );
+            }
+        }
+        else
+        {
+            // finally, build the new dataset
+            debugLog() << "Calculating value-set";
+            size_t maxI = hist->size() - 1;
+            for ( size_t vi = 0; vi < valueSet->rawSize(); ++vi )
+            {
+                size_t idx = hist->getIndexForValue( valueSet->getScalarDouble( vi ) );
+                newData[ vi ] = static_cast< unsigned char >( idx / maxI * 255 );
+            }
+        }
+
+        // update output with a new dataset, reuse grid
+        debugLog() << "Updating output";
+
+        // de-register at datahandler
+        if ( m_lastOutputDataSet )
+        {
+            WDataHandler::deregisterDataSet( m_lastOutputDataSet );
+        }
+
+        // construct
+        m_lastOutputDataSet = boost::shared_ptr< WDataSetScalar >(
+            new WDataSetScalar( boost::shared_ptr< WValueSetBase >(
+                new WValueSet< unsigned char >( 0, 1, newData, W_DT_UNSIGNED_CHAR ) ), dataSet->getGrid() )
+        );
+
+        m_lastOutputDataSet->getTexture()->setSelectedColormap( 0 );
+
+        // register new
+        WDataHandler::registerDataSet( m_lastOutputDataSet );
+        m_output->updateData( m_lastOutputDataSet );
 
-        // Grab the histogram
-        boost::shared_ptr< const WValueSetHistogram > hist = dataSet->getHistogram();
-        debugLog() << *hist;
+        debugLog() << "Done";
     }
 }
 

src/modules/histogramEqualization/WMHistogramEqualization.h

@@ -109,7 +109,7 @@ private:
     /**
      * This is a pointer to the dataset the module is currently working on.
      */
-    boost::shared_ptr< WDataSetScalar > m_dataSet;
+    boost::shared_ptr< WDataSetScalar > m_lastOutputDataSet;
 
     /**
      * A condition used to notify about changes in several properties.
@@ -120,6 +120,16 @@ private:
      * If true, histogram equalization is turned on.
      */
     WPropBool m_equalize;
+
+    /**
+     * How many percent should be clamped from the histogram.
+     */
+    WPropInt m_clampPerc;
+
+    /**
+     * Resolution of the initial histogram.
+     */
+    WPropInt m_histogramResolution;
 };
 
 #endif  // WMHISTOGRAMEQUALIZATION_H