OpenWalnut  1.4.0
WTransferFunction.cpp
00001 //---------------------------------------------------------------------------
00002 //
00003 // Project: OpenWalnut ( http://www.openwalnut.org )
00004 //
00005 // Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
00006 // For more information see http://www.openwalnut.org/copying
00007 //
00008 // This file is part of OpenWalnut.
00009 //
00010 // OpenWalnut is free software: you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as published by
00012 // the Free Software Foundation, either version 3 of the License, or
00013 // (at your option) any later version.
00014 //
00015 // OpenWalnut is distributed in the hope that it will be useful,
00016 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00018 // GNU Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public License
00021 // along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
00022 //
00023 //---------------------------------------------------------------------------
00024 
00025 #include <algorithm>
00026 #include <cmath>
00027 #include <iostream>
00028 #include <vector>
00029 
00030 #include "WAssert.h"
00031 #include "WTransferFunction.h"
00032 
00033 bool WTransferFunction::operator==( const WTransferFunction &rhs ) const
00034 {
00035     if( m_histogram.size() != rhs.m_histogram.size() )
00036     {
00037         return false;
00038     }
00039     {
00040         std::vector< double >::const_iterator ait1 = m_histogram.begin();
00041         std::vector< double >::const_iterator ait2 = rhs.m_histogram.begin();
00042         for( ;
00043                 ait1 != m_histogram.end();
00044                 ++ait1, ++ait2 )
00045         {
00046             if( *ait1 != *ait2 )
00047             {
00048                 return false;
00049             }
00050         }
00051     }
00052 
00053 
00054     if( m_colors.size() != rhs.m_colors.size() || m_alphas.size() != rhs.m_alphas.size() )
00055     {
00056         return false;
00057     }
00058 
00059     if( m_isomin != rhs.m_isomin && m_isomax != rhs.m_isomax )
00060     {
00061         return false;
00062     }
00063 
00064     std::vector< ColorEntry >::const_iterator it1 = m_colors.begin();
00065     std::vector< ColorEntry >::const_iterator it2 = rhs.m_colors.begin();
00066     for( ;
00067             it1 != m_colors.end();
00068             ++it1, ++it2 )
00069     {
00070         if( !( *it1 == *it2 ) )
00071         {
00072             return false;
00073         }
00074     }
00075 
00076     std::vector< AlphaEntry >::const_iterator ait1 = m_alphas.begin();
00077     std::vector< AlphaEntry >::const_iterator ait2 = rhs.m_alphas.begin();
00078     for( ;
00079             ait1 != m_alphas.end();
00080             ++ait1, ++ait2 )
00081     {
00082         if( !( *ait1 == *ait2 ) )
00083         {
00084             return false;
00085         }
00086     }
00087 
00088     return true;
00089 }
00090 
00091 bool WTransferFunction::operator!=( const WTransferFunction &rhs ) const
00092 {
00093     return !( ( *this ) == rhs );
00094 }
00095 
00096 namespace
00097 {
00098     //! linear blend between two colors in rgb space if ta = 1.-tb
00099     WColor blend( const WColor&a, double ta, const WColor &b, double tb )
00100     {
00101         return WColor(
00102                 ta*a[ 0 ]+tb*b[ 0 ],
00103                 ta*a[ 1 ]+tb*b[ 1 ],
00104                 ta*a[ 2 ]+tb*b[ 2 ], 1. );
00105     }
00106 
00107     //! linear blend between two variables
00108     double ablend( const double a, const double ta, const double b, const double tb )
00109     {
00110         return a*ta + b*tb;
00111     }
00112 } // namespace
00113 
00114 
00115 void WTransferFunction::sample1DTransferFunction( unsigned char*array, int width, double min, double max ) const
00116 {
00117     if( m_colors.size() < 1 ) return;
00118     if( m_alphas.size() < 1 ) return;
00119 
00120     std::vector< ColorEntry >::const_iterator c1 = m_colors.begin();
00121     std::vector< ColorEntry >::const_iterator c2 = c1+1;
00122 
00123     std::vector< AlphaEntry >::const_iterator a1 = m_alphas.begin();
00124     std::vector< AlphaEntry >::const_iterator a2 = a1+1;
00125 
00126     for( int i = 0; i < width; ++i )
00127     {
00128         WColor color;
00129         double iso = ( double )i/( double )width * ( max-min ) + min;
00130 
00131         if( iso <= m_isomin )
00132         {
00133             color = m_colors.begin()->color;
00134             color[ 3 ] = m_alphas.begin()->alpha;
00135         }
00136         else if( iso >= m_isomax )
00137         {
00138             color = m_colors.back().color;
00139             color[ 3 ] = m_alphas.back().alpha;
00140         }
00141         else
00142         {
00143             while( c2 != m_colors.end() && iso > c2->iso )
00144             {
00145                 WAssert( c2 != m_colors.end(), "Corruption of internal data structure." );
00146                 c1++;
00147                 c2++;
00148                 WAssert( c1 != m_colors.end(), "Corruption of internal data structure." );
00149             }
00150 
00151             while( a2 != m_alphas.end() && iso > a2->iso )
00152             {
00153                 WAssert( a2 != m_alphas.end(), "Corruption of internal data structure." );
00154                 a1++;
00155                 a2++;
00156                 WAssert( a1 != m_alphas.end(), "Corruption of internal data structure." );
00157             }
00158 
00159             if( c2 == m_colors.end() )
00160             {
00161                 color = c1->color;
00162             }
00163             else
00164             {
00165                 double colorParameter = ( iso - c1->iso )/( c2->iso - c1->iso );
00166                 color = blend( c1->color, 1.-colorParameter, c2->color, colorParameter );
00167             }
00168             if( a2 == m_alphas.end() )
00169             {
00170                 color[ 3 ] = a1->alpha;
00171             }
00172             else
00173             {
00174                 double alphaParameter = ( iso - a1->iso )/( a2->iso - a1->iso );
00175                 color[ 3 ] = ablend( a1->alpha, 1.-alphaParameter, a2->alpha, alphaParameter );
00176             }
00177         }
00178         for( int j = 0; j < 4; ++j )
00179         {
00180             array[ 4*i + j ] = color[ j ]*255.;
00181         }
00182     }
00183 }
00184 
00185 
00186 void WTransferFunction::addColor( double iso, const WColor& color )
00187 {
00188     if( m_colors.size() == 0 )
00189     {
00190         m_colors.push_back( ColorEntry( iso, color ) );
00191     }
00192     else
00193     {
00194         std::vector<ColorEntry>::iterator e = find_if( m_colors.begin(), m_colors.end(), LessPred<ColorEntry>( iso ) );
00195         m_colors.insert( e, ColorEntry( iso, color ) );
00196     }
00197 
00198     if( m_alphas.size() >= 1 )
00199     {
00200         m_isomin = std::min( m_colors.front().iso, m_alphas.front().iso );
00201         m_isomax = std::max( m_colors.back().iso, m_alphas.back().iso );
00202     }
00203     else
00204     {
00205         m_isomin = m_colors.front().iso;
00206         m_isomax = m_colors.back().iso;
00207     }
00208 }
00209 
00210 void WTransferFunction::addAlpha( double iso, double alpha )
00211 {
00212     if( m_alphas.size() == 0 )
00213     {
00214         m_alphas.push_back( AlphaEntry( iso, alpha ) );
00215     }
00216     else
00217     {
00218         std::vector<AlphaEntry>::iterator e = find_if( m_alphas.begin(), m_alphas.end(), LessPred<AlphaEntry>( iso ) );
00219         m_alphas.insert( e, AlphaEntry( iso, alpha ) );
00220     }
00221 
00222     if( m_colors.size() >= 1 )
00223     {
00224         m_isomin = std::min( m_colors.front().iso, m_alphas.front().iso );
00225         m_isomax = std::max( m_colors.back().iso, m_alphas.back().iso );
00226     }
00227     else
00228     {
00229         m_isomin = m_alphas.front().iso;
00230         m_isomax = m_alphas.back().iso;
00231     }
00232 }
00233 
00234 WTransferFunction WTransferFunction::createFromRGBA( unsigned char const * const rgba, size_t size )
00235 {
00236     // we create a linear match to the transfer funciton given by rgba by scanning the
00237     // alpha and color values in two passes.
00238     // each pass starts at the left of the picture and moves right looking for a good match of
00239     // the line between start point and end point to the given function in between. The error is
00240     // normalized to a per-sample basis. If the maximum error is below MIN_ERROR_THRESHOLD, the
00241     // line is accepted and the next segment is analyzed.
00242     const double MIN_ERROR_THRESHOLD = 5.0;
00243     WTransferFunction rgbatf;
00244     std::vector < float > values( size );
00245 
00246     // copy channel
00247     for( size_t i = 0; i < size; ++i )
00248     {
00249         values[ i ] =  static_cast<double>( rgba[ i*4+3 ] );
00250     }
00251 
00252     // add first and last alpha
00253     rgbatf.addAlpha( 0.0, values[ 0 ]/255. );
00254     rgbatf.addAlpha( 1.0, values[ size-1 ]/255. );
00255 
00256     std::vector < float > errors( size );
00257 
00258     size_t seed = 0;
00259     while( seed < size-1 )
00260     {
00261         // start at first pixel and fit a line to the data
00262         size_t to = seed+1;
00263         while( to < size )
00264         {
00265             double error = 0.0;
00266             double incline =  ( values[ to ] - values[ seed ] )/( to-seed );
00267             for( size_t j = seed+1; j < to; ++j )
00268             {
00269                 error +=  std::sqrt( std::pow( values[ j ] - values[ seed ] - incline * ( j-seed ), 2 ) );
00270             }
00271             errors[ to ] =  error/( to-seed ); // compute square error per pixel length of line
00272             ++to;
00273         }
00274         size_t minElement = size-1;
00275         double minerror = errors[ minElement ];
00276         for( to = size-1; to > seed; --to )
00277         {
00278             if( errors[ to ] <  minerror )
00279             {
00280                 minElement = to;
00281                 minerror = errors[ to ];
00282                 if( minerror < MIN_ERROR_THRESHOLD )
00283                 {
00284                     break;
00285                 }
00286             }
00287         }
00288         if( minElement < size-1 )
00289         {
00290             rgbatf.addAlpha( ( double )minElement/( double )( size-1 ), values[ minElement ]/255. );
00291         }
00292         seed = minElement;
00293     }
00294 
00295 
00296     // same for color
00297     // add first and last color
00298     rgbatf.addColor( 0.0, WColor( rgba[ 0*4+0 ]/255.f, rgba[ 0*4+1 ]/255.f, rgba[ 0*4+2 ]/255.f, 0.f ) );
00299     rgbatf.addColor( 1.0, WColor( rgba[ ( size-1 )*4+0 ]/255.f, rgba[ ( size-1 )*4+1 ]/255.f, rgba[ ( size-1 )*4+2 ]/255.f, 0.f ) );
00300 
00301     // first try of code: use combined RGB errors
00302 
00303     seed = 0;
00304     while( seed < size-1 )
00305     {
00306         // start at first pixel and fit a line to the data
00307         size_t to = seed+1;
00308         while( to < size )
00309         {
00310             double error = 0.0;
00311             double inclineR =  ( rgba[ to*4+0 ] - rgba[ seed*4+0 ] )/( to-seed );
00312             double inclineG =  ( rgba[ to*4+1 ] - rgba[ seed*4+1 ] )/( to-seed );
00313             double inclineB =  ( rgba[ to*4+2 ] - rgba[ seed*4+2 ] )/( to-seed );
00314 
00315             for( size_t j = seed; j < to; ++j )
00316             {
00317                 error +=  std::sqrt(
00318                         std::pow( rgba[ 4*j+0 ] - rgba[ 4*seed+0 ] - inclineR * ( j-seed ), 2 ) +
00319                         std::pow( rgba[ 4*j+1 ] - rgba[ 4*seed+1 ] - inclineG * ( j-seed ), 2 ) +
00320                         std::pow( rgba[ 4*j+2 ] - rgba[ 4*seed+2 ] - inclineB * ( j-seed ), 2 )
00321                         );
00322             }
00323             errors[ to ] =  error/( to-seed ); // compute square error per pixel length of line
00324             ++to;
00325         }
00326 
00327         size_t minElement = size-1;
00328         double minerror = errors[ size-1 ];
00329         // traverse from back
00330         for( to = size-2; to > seed; --to )
00331         {
00332             if( errors[ to ] <  minerror )
00333             {
00334                 minElement = to;
00335                 minerror = errors[ to ];
00336             }
00337             if( minerror < MIN_ERROR_THRESHOLD*2.0 ) //! the threshold here is larger than for alpha, becuase we compare all colors at once
00338             {
00339                 break;
00340             }
00341         }
00342         if( minElement < size-1 )
00343         {
00344             rgbatf.addColor( ( double )minElement/( double )( size-1 ),
00345                 WColor( rgba[ minElement*4+0 ]/255.f, rgba[ minElement*4+1 ]/255.f, rgba[ minElement*4+2 ]/255.f, 0.f ) );
00346         }
00347         seed = minElement;
00348     }
00349 
00350     std::cout <<  "New Transfer Function: " << rgbatf << "." <<  std::endl;
00351     return rgbatf;
00352 }
00353 
00354 std::ostream& operator << ( std::ostream& out, const WTransferFunction& tf )
00355 {
00356     size_t numColors = tf.numColors();
00357     for( size_t i = 0; i < numColors; ++i )
00358     {
00359         double iso = tf.getColorIsovalue( i );
00360         WColor c = tf.getColor( i );
00361         out << "c:" << iso << ":" << c[ 0 ] << ":" << c[ 1 ] << ":" << c[ 2 ] << ";";
00362     }
00363     size_t numAlphas = tf.numAlphas();
00364     for( size_t i = 0; i < numAlphas; ++i )
00365     {
00366         double iso = tf.getAlphaIsovalue( i );
00367         double alpha = tf.getAlpha( i );
00368         out << "a:" << iso << ":" << alpha;
00369         if( i != numAlphas-1 )
00370         {
00371             out << ";";
00372         }
00373     }
00374     return out;
00375 }