OpenWalnut  1.4.0
WJoinContourTree.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 <set>
00027 #include <string>
00028 #include <vector>
00029 
00030 #include "../../common/WStringUtils.h"
00031 #include "../../common/WTransferable.h"
00032 #include "../../common/datastructures/WUnionFind.h"
00033 #include "../../common/exceptions/WNotImplemented.h"
00034 #include "../WValueSet.h"
00035 #include "WJoinContourTree.h"
00036 
00037 WJoinContourTree::WJoinContourTree()
00038     : WTransferable()
00039 {
00040 }
00041 
00042 WJoinContourTree::WJoinContourTree( boost::shared_ptr< WDataSetSingle > dataset )
00043     : WTransferable(),
00044       m_elementIndices( dataset->getValueSet()->size() ),
00045       m_joinTree( dataset->getValueSet()->size() ),
00046       m_lowestVoxel( dataset->getValueSet()->size() )
00047 {
00048     if( dataset->getValueSet()->order() != 0 || dataset->getValueSet()->dimension() != 1 )
00049     {
00050         throw WNotImplemented( std::string( "ATM there is only support for scalar fields" ) );
00051     }
00052     m_valueSet = boost::dynamic_pointer_cast< WValueSet< double > >( dataset->getValueSet() );
00053     if( !m_valueSet )
00054     {
00055         throw WNotImplemented( std::string( "ATM there is only support for scalar fields with doubles as scalars" ) );
00056     }
00057     m_grid = boost::dynamic_pointer_cast< WGridRegular3D >( dataset->getGrid() );
00058     if( !m_grid )
00059     {
00060         throw WNotImplemented( std::string( "Only WGridRegular3D is supported, despite that its not a simplicial mesh!" ) );
00061     }
00062     for( size_t i = 0; i < m_elementIndices.size(); ++i )
00063     {
00064         m_joinTree[i] = m_elementIndices[i] = i;
00065     }
00066 }
00067 
00068 void WJoinContourTree::sortIndexArray()
00069 {
00070     IndirectCompare comp( m_valueSet );
00071     std::sort( m_elementIndices.begin(), m_elementIndices.end(), comp );
00072 }
00073 
00074 void WJoinContourTree::buildJoinTree()
00075 {
00076     sortIndexArray();
00077 
00078     WUnionFind uf( m_joinTree.size() );
00079 
00080     for( size_t i = 0; i < m_joinTree.size(); ++i )
00081     {
00082         m_lowestVoxel[ m_elementIndices[i] ] = m_elementIndices[i];
00083         std::vector< size_t > neighbours = m_grid->getNeighbours( m_elementIndices[i] );
00084         std::vector< size_t >::const_iterator n = neighbours.begin();
00085         for( ; n != neighbours.end(); ++n )
00086         {
00087             if( uf.find( m_elementIndices[i] ) == uf.find( *n ) || m_valueSet->getScalar( *n ) <= m_valueSet->getScalar( m_elementIndices[i] ) )
00088             {
00089                 continue;
00090             }
00091             else
00092             {
00093                 uf.merge( m_elementIndices[i], *n );
00094                 m_joinTree[ m_lowestVoxel[ uf.find( *n ) ] ] = m_elementIndices[i];
00095                 m_lowestVoxel[ uf.find( *n ) ] = m_elementIndices[i];
00096             }
00097         }
00098     }
00099 }
00100 
00101 boost::shared_ptr< std::set< size_t > > WJoinContourTree::getVolumeVoxelsEnclosedByIsoSurface( const double isoValue ) const
00102 {
00103     boost::shared_ptr< std::vector< size_t > > result( new std::vector< size_t >( m_elementIndices ) );
00104     WUnionFind uf( m_elementIndices.size() );
00105 
00106     // using string_utils::operator<<;
00107     // std::cout << "m_element: " << m_elementIndices << std::endl;
00108 
00109     //std::stringstream ss;
00110     // assume the m_elementIndices array is still sorted descending on its iso values in the valueset
00111     for( size_t i = 0; i < m_elementIndices.size() && m_valueSet->getScalar( m_elementIndices[i] ) >= isoValue; ++i )
00112     {
00113         // std::cout << "processing element: " << i << std::endl;
00114         // std::cout << "having index: " << m_elementIndices[i] << std::endl;
00115         // using string_utils::operator<<;
00116         // std::cout << "xyz: " << m_grid->getNeighbours( m_elementIndices[i] ) << std::endl;
00117         // std::cout << "having isovalue: " << m_valueSet->getScalar( m_elementIndices[i] ) << std::endl;
00118         // ss << " m_elementIndices[i]:isovalue=" << m_valueSet->getScalar( m_elementIndices[i] ) << ", ";
00119         size_t target = m_joinTree[ m_elementIndices[i] ];
00120         // std::cout << "having edge to: " << target << std::endl;
00121         if( m_valueSet->getScalar( target ) >= isoValue )
00122         {
00123             uf.merge( target, m_elementIndices[i] );
00124         }
00125     }
00126     //std::cout << ss.str() << std::endl;
00127     return uf.getMaxSet();
00128 }
00129 
00130 WJoinContourTree::IndirectCompare::IndirectCompare( boost::shared_ptr< WValueSet< double > > valueSet )
00131     : m_valueSet( valueSet )
00132 {
00133 }
00134 
00135 bool WJoinContourTree::IndirectCompare::operator()( size_t i, size_t j )
00136 {
00137     return ( m_valueSet->getScalar( i ) > m_valueSet->getScalar( j ) );
00138 }
00139 
00140 // make the class beeing a WTransferrable:
00141 // ---------------------------------------
00142 boost::shared_ptr< WPrototyped > WJoinContourTree::m_prototype = boost::shared_ptr< WPrototyped >();
00143 
00144 boost::shared_ptr< WPrototyped > WJoinContourTree::getPrototype()
00145 {
00146     if( !m_prototype )
00147     {
00148         m_prototype = boost::shared_ptr< WPrototyped >( new WJoinContourTree() );
00149     }
00150     return m_prototype;
00151 }