OpenWalnut  1.3.1
WJoinContourTree.cpp
1 //---------------------------------------------------------------------------
2 //
3 // Project: OpenWalnut ( http://www.openwalnut.org )
4 //
5 // Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
6 // For more information see http://www.openwalnut.org/copying
7 //
8 // This file is part of OpenWalnut.
9 //
10 // OpenWalnut is free software: you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // OpenWalnut is distributed in the hope that it will be useful,
16 // but WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public License
21 // along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
22 //
23 //---------------------------------------------------------------------------
24 
25 #include <algorithm>
26 #include <set>
27 #include <string>
28 #include <vector>
29 
30 #include "../../common/WStringUtils.h"
31 #include "../../common/WTransferable.h"
32 #include "../../common/datastructures/WUnionFind.h"
33 #include "../../common/exceptions/WNotImplemented.h"
34 #include "../WValueSet.h"
35 #include "WJoinContourTree.h"
36 
38  : WTransferable()
39 {
40 }
41 
42 WJoinContourTree::WJoinContourTree( boost::shared_ptr< WDataSetSingle > dataset )
43  : WTransferable(),
44  m_elementIndices( dataset->getValueSet()->size() ),
45  m_joinTree( dataset->getValueSet()->size() ),
46  m_lowestVoxel( dataset->getValueSet()->size() )
47 {
48  if( dataset->getValueSet()->order() != 0 || dataset->getValueSet()->dimension() != 1 )
49  {
50  throw WNotImplemented( std::string( "ATM there is only support for scalar fields" ) );
51  }
52  m_valueSet = boost::dynamic_pointer_cast< WValueSet< double > >( dataset->getValueSet() );
53  if( !m_valueSet )
54  {
55  throw WNotImplemented( std::string( "ATM there is only support for scalar fields with doubles as scalars" ) );
56  }
57  m_grid = boost::dynamic_pointer_cast< WGridRegular3D >( dataset->getGrid() );
58  if( !m_grid )
59  {
60  throw WNotImplemented( std::string( "Only WGridRegular3D is supported, despite that its not a simplicial mesh!" ) );
61  }
62  for( size_t i = 0; i < m_elementIndices.size(); ++i )
63  {
64  m_joinTree[i] = m_elementIndices[i] = i;
65  }
66 }
67 
69 {
71  std::sort( m_elementIndices.begin(), m_elementIndices.end(), comp );
72 }
73 
75 {
77 
78  WUnionFind uf( m_joinTree.size() );
79 
80  for( size_t i = 0; i < m_joinTree.size(); ++i )
81  {
82  m_lowestVoxel[ m_elementIndices[i] ] = m_elementIndices[i];
83  std::vector< size_t > neighbours = m_grid->getNeighbours( m_elementIndices[i] );
84  std::vector< size_t >::const_iterator n = neighbours.begin();
85  for( ; n != neighbours.end(); ++n )
86  {
87  if( uf.find( m_elementIndices[i] ) == uf.find( *n ) || m_valueSet->getScalar( *n ) <= m_valueSet->getScalar( m_elementIndices[i] ) )
88  {
89  continue;
90  }
91  else
92  {
93  uf.merge( m_elementIndices[i], *n );
94  m_joinTree[ m_lowestVoxel[ uf.find( *n ) ] ] = m_elementIndices[i];
95  m_lowestVoxel[ uf.find( *n ) ] = m_elementIndices[i];
96  }
97  }
98  }
99 }
100 
101 boost::shared_ptr< std::set< size_t > > WJoinContourTree::getVolumeVoxelsEnclosedByIsoSurface( const double isoValue ) const
102 {
103  boost::shared_ptr< std::vector< size_t > > result( new std::vector< size_t >( m_elementIndices ) );
104  WUnionFind uf( m_elementIndices.size() );
105 
106  // using string_utils::operator<<;
107  // std::cout << "m_element: " << m_elementIndices << std::endl;
108 
109  //std::stringstream ss;
110  // assume the m_elementIndices array is still sorted descending on its iso values in the valueset
111  for( size_t i = 0; i < m_elementIndices.size() && m_valueSet->getScalar( m_elementIndices[i] ) >= isoValue; ++i )
112  {
113  // std::cout << "processing element: " << i << std::endl;
114  // std::cout << "having index: " << m_elementIndices[i] << std::endl;
115  // using string_utils::operator<<;
116  // std::cout << "xyz: " << m_grid->getNeighbours( m_elementIndices[i] ) << std::endl;
117  // std::cout << "having isovalue: " << m_valueSet->getScalar( m_elementIndices[i] ) << std::endl;
118  // ss << " m_elementIndices[i]:isovalue=" << m_valueSet->getScalar( m_elementIndices[i] ) << ", ";
119  size_t target = m_joinTree[ m_elementIndices[i] ];
120  // std::cout << "having edge to: " << target << std::endl;
121  if( m_valueSet->getScalar( target ) >= isoValue )
122  {
123  uf.merge( target, m_elementIndices[i] );
124  }
125  }
126  //std::cout << ss.str() << std::endl;
127  return uf.getMaxSet();
128 }
129 
131  : m_valueSet( valueSet )
132 {
133 }
134 
136 {
137  return ( m_valueSet->getScalar( i ) > m_valueSet->getScalar( j ) );
138 }
139 
140 // make the class beeing a WTransferrable:
141 // ---------------------------------------
142 boost::shared_ptr< WPrototyped > WJoinContourTree::m_prototype = boost::shared_ptr< WPrototyped >();
143 
144 boost::shared_ptr< WPrototyped > WJoinContourTree::getPrototype()
145 {
146  if( !m_prototype )
147  {
148  m_prototype = boost::shared_ptr< WPrototyped >( new WJoinContourTree() );
149  }
150  return m_prototype;
151 }