OpenWalnut  1.2.5
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
WFiberCluster.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 <list>
26 #include <vector>
27 
28 #include <boost/shared_ptr.hpp>
29 
30 #include "../../common/math/WMath.h"
31 #include "../../common/math/WPlane.h"
32 #include "../../common/WLimits.h"
33 #include "../../common/WTransferable.h"
34 #include "../WDataSetFiberVector.h"
35 #include "WFiberCluster.h"
36 
37 // TODO(math): The only reason why we store here a Reference to the fiber
38 // dataset is, we need it in the WMVoxelizer module as well as the clustering
39 // information. Since we don't have the possibility of multiple
40 // InputConnectors we must agglomerate those into one object. Please remove this.
41 // initializes the variable and provides a linker reference
42 // \cond Suppress_Doxygen
43 boost::shared_ptr< WPrototyped > WFiberCluster::m_prototype = boost::shared_ptr< WPrototyped >();
44 // \endcond
45 
47  : WTransferable(),
48  m_centerLineCreationLock( new boost::shared_mutex() ),
49  m_longestLineCreationLock( new boost::shared_mutex() )
50 {
51 }
52 
54  : WTransferable(),
55  m_centerLineCreationLock( new boost::shared_mutex() ),
56  m_longestLineCreationLock( new boost::shared_mutex() )
57 {
58  m_memberIndices.push_back( index );
59 }
60 
62  : WTransferable( other ),
63  m_memberIndices( other.m_memberIndices ),
64  m_fibs( other.m_fibs ),
65  m_color( other.m_color ),
66  m_centerLineCreationLock( new boost::shared_mutex() ), // do not copy the mutex as both instances of WFiberCluster can be modifed at the
67  // same time
68  m_longestLineCreationLock( new boost::shared_mutex() ),
69  m_centerLine(), // << we can't ensure that the centerline and longest line are initialized yet, but we want a deep copy
70  m_longestLine()
71 {
72  // copy them only if they exist
73  if( other.m_centerLine )
74  {
75  m_centerLine = boost::shared_ptr< WFiber >( new WFiber( *other.m_centerLine.get() ) );
76  }
77  if( other.m_longestLine )
78  {
79  m_longestLine = boost::shared_ptr< WFiber >( new WFiber( *other.m_longestLine.get() ) );
80  }
81 }
82 
84 {
87 }
88 
89 void WFiberCluster::merge( WFiberCluster& other ) // NOLINT
90 {
91  std::list< size_t >::const_iterator cit = other.m_memberIndices.begin();
92  for( ; cit != other.m_memberIndices.end(); ++cit)
93  {
94  m_memberIndices.push_back( *cit );
95  }
96  // make sure that those indices aren't occuring anywhere else
97  other.m_centerLine.reset(); // they are not valid anymore
98  other.m_longestLine.reset();
99  other.clear();
100 }
101 
102 // NODOXYGEN
103 // \cond Suppress_Doxygen
104 void WFiberCluster::setDataSetReference( boost::shared_ptr< const WDataSetFiberVector > fibs )
105 {
106  m_fibs = fibs;
107 }
108 
109 boost::shared_ptr< const WDataSetFiberVector > WFiberCluster::getDataSetReference() const
110 {
111  return m_fibs;
112 }
113 
114 // TODO(math): The only reason why we store here a Reference to the fiber
115 // dataset is, we need it in the WMVoxelizer module as well as the clustering
116 // information. Since we don't have the possibility of multiple
117 // InputConnectors we must agglomerate those into one object. Please remove this.
118 boost::shared_ptr< WPrototyped > WFiberCluster::getPrototype()
119 {
120  if( !m_prototype )
121  {
122  m_prototype = boost::shared_ptr< WPrototyped >( new WFiberCluster() );
123  }
124  return m_prototype;
125 }
126 // \endcond
127 
129 {
130  // ensure nobody changes the mutable m_centerline
131  boost::unique_lock< boost::shared_mutex > lock = boost::unique_lock< boost::shared_mutex >( *m_centerLineCreationLock );
132  // has the line been calculated while we waited?
133  if( m_centerLine )
134  {
135  lock.unlock();
136  return;
137  }
138 
139  // make copies of the fibers
140  boost::shared_ptr< WDataSetFiberVector > fibs( new WDataSetFiberVector() );
141  size_t avgFiberSize = 0;
142  for( std::list< size_t >::const_iterator cit = m_memberIndices.begin(); cit != m_memberIndices.end(); ++cit )
143  {
144  fibs->push_back( m_fibs->at( *cit ) );
145  avgFiberSize += fibs->back().size();
146  }
147  avgFiberSize /= fibs->size();
148 
149  unifyDirection( fibs );
150 
151  for( WDataSetFiberVector::iterator cit = fibs->begin(); cit != fibs->end(); ++cit )
152  {
153  cit->resampleByNumberOfPoints( avgFiberSize );
154  }
155 
156  m_centerLine = boost::shared_ptr< WFiber >( new WFiber() );
157  m_centerLine->reserve( avgFiberSize );
158  for( size_t i = 0; i < avgFiberSize; ++i )
159  {
160  WPosition avgPosition( 0, 0, 0 );
161  for( WDataSetFiberVector::const_iterator cit = fibs->begin(); cit != fibs->end(); ++cit )
162  {
163  avgPosition += cit->at( i );
164  }
165  avgPosition /= static_cast< double >( fibs->size() );
166  m_centerLine->push_back( avgPosition );
167  }
168 
170  lock.unlock();
171 }
172 
174 {
175  // ensure nobody changes the mutable m_longestline
176  boost::unique_lock< boost::shared_mutex > lock = boost::unique_lock< boost::shared_mutex >( *m_longestLineCreationLock );
177  // has the line been calculated while we waited?
178  if( m_longestLine )
179  {
180  lock.unlock();
181  return;
182  }
183 
184  m_longestLine = boost::shared_ptr< WFiber >( new WFiber() );
185 
186  // empty datasets can be ignored
187  if( m_fibs->size() == 0 )
188  {
189  return;
190  }
191 
192  size_t longest = 0;
193  size_t longestID = 0;
194  for( size_t cit = 0; cit < m_fibs->size(); ++cit )
195  {
196  if( m_fibs->at( cit ).size() > longest )
197  {
198  longest = m_fibs->at( cit ).size();
199  longestID = cit;
200  }
201  }
202 
203  for( WFiber::const_iterator cit = m_fibs->at( longestID ).begin(); cit != m_fibs->at( longestID ).end(); ++cit )
204  {
205  m_longestLine->push_back( *cit );
206  }
207 
208  lock.unlock();
209 }
210 
212 {
213  // first ending of the centerline
214  assert( m_centerLine->size() > 1 );
215  WFiber cL( *m_centerLine );
216  WPlane p( cL[0] - cL[1], cL[0] + ( cL[0] - cL[1] ) );
217  boost::shared_ptr< WPosition > cutPoint( new WPosition( 0, 0, 0 ) );
218  bool intersectionFound = true;
219 
220  // in the beginning all fibers participate
221  boost::shared_ptr< WDataSetFiberVector > fibs( new WDataSetFiberVector() );
222  for( std::list< size_t >::const_iterator cit = m_memberIndices.begin(); cit != m_memberIndices.end(); ++cit )
223  {
224  fibs->push_back( m_fibs->at( *cit ) );
225  }
226 
227  while( intersectionFound )
228  {
229  intersectionFound = false;
230  size_t intersectingFibers = 0;
231 // WPosition avg( 0, 0, 0 );
232  for( WDataSetFiberVector::iterator cit = fibs->begin(); cit != fibs->end(); )
233  {
234  if( intersectPlaneLineNearCP( p, *cit, cutPoint ) )
235  {
236  if( length( *cutPoint - p.getPosition() ) < 20 )
237  {
238 // avg += *cutPoint;
239  intersectingFibers++;
240  intersectionFound = true;
241  ++cit;
242  }
243  else
244  {
245  cit = fibs->erase( cit );
246  }
247  }
248  else
249  {
250  cit = fibs->erase( cit );
251  }
252  }
253  if( intersectingFibers > 10 )
254  {
255  cL.insert( cL.begin(), cL[0] + ( cL[0] - cL[1] ) );
256  p.resetPosition( cL[0] + ( cL[0] - cL[1] ) );
257 // avg[0] /= static_cast< double >( intersectingFibers );
258 // avg[1] /= static_cast< double >( intersectingFibers );
259 // avg[2] /= static_cast< double >( intersectingFibers );
260 // cL.insert( cL.begin(), 0.995 * ( cL[0] + ( cL[0] - cL[1] ) ) + 0.005 * avg );
261 // p.resetPosition( cL[0] + ( cL[0] - cL[1] ) );
262 // p.setNormal( ( cL[0] - cL[1] ) );
263  }
264  else // no intersections found => abort
265  {
266  break;
267  }
268  }
269  // second ending of the centerline
270  boost::shared_ptr< WDataSetFiberVector > fobs( new WDataSetFiberVector() );
271  for( std::list< size_t >::const_iterator cit = m_memberIndices.begin(); cit != m_memberIndices.end(); ++cit )
272  {
273  fobs->push_back( m_fibs->at( *cit ) );
274  }
275 
276  // try to discard other lines from other end
277 
278  WPlane q( cL.back() - cL[ cL.size() - 2 ], cL.back() + ( cL.back() - cL[ cL.size() - 2 ] ) );
279  intersectionFound = true;
280  while( intersectionFound )
281  {
282  intersectionFound = false;
283  size_t intersectingFibers = 0;
284 // WPosition avg( 0, 0, 0 );
285  for( WDataSetFiberVector::iterator cit = fobs->begin(); cit != fobs->end(); )
286  {
287  if( intersectPlaneLineNearCP( q, *cit, cutPoint ) )
288  {
289  if( length( *cutPoint - q.getPosition() ) < 20 )
290  {
291 // avg += *cutPoint;
292  intersectingFibers++;
293  intersectionFound = true;
294  ++cit;
295  }
296  else
297  {
298  cit = fobs->erase( cit );
299  }
300  }
301  else
302  {
303  cit = fobs->erase( cit );
304  }
305  }
306  if( intersectingFibers > 10 )
307  {
308  cL.push_back( cL.back() + ( cL.back() - cL[ cL.size() - 2 ] ) );
309  q.resetPosition( cL.back() + ( cL.back() - cL[ cL.size() - 2 ] ) );
310 // avg[0] /= static_cast< double >( intersectingFibers );
311 // avg[1] /= static_cast< double >( intersectingFibers );
312 // avg[2] /= static_cast< double >( intersectingFibers );
313 // cL.push_back( 0.995 * ( cL.back() + ( cL.back() - cL[ cL.size() - 2 ] ) ) + 0.005 * avg );
314 // q.resetPosition( cL.back() + ( cL.back() - cL[ cL.size() - 2 ] ) );
315 // q.setNormal( cL.back() - cL[ cL.size() - 2 ] );
316  }
317  else // no intersections found => abort
318  {
319  break;
320  }
321  }
322  *m_centerLine = cL;
323 }
324 
325 void WFiberCluster::unifyDirection( boost::shared_ptr< WDataSetFiberVector > fibs ) const
326 {
327  if( fibs->size() < 2 )
328  {
329  return;
330  }
331 
332  assert( !( fibs->at( 0 ).empty() ) && "WFiberCluster.unifyDirection: Empty fiber processed.. aborting" );
333 
334  // first fiber defines direction
335  const WFiber& firstFib = fibs->front();
336  const WPosition start = firstFib.front();
337  const WPosition m1 = firstFib.at( firstFib.size() * 1.0 / 3.0 );
338  const WPosition m2 = firstFib.at( firstFib.size() * 2.0 / 3.0 );
339  const WPosition end = firstFib.back();
340 
341  for( WDataSetFiberVector::iterator cit = fibs->begin() + 1; cit != fibs->end(); ++cit )
342  {
343  const WFiber& other = *cit;
344  double distance = length2( start - other.front() ) +
345  length2( m1 - other.at( other.size() * 1.0 / 3.0 ) ) +
346  length2( m2 - other.at( other.size() * 2.0 / 3.0 ) ) +
347  length2( end - other.back() );
348  double inverseDistance = length2( start - other.back() ) +
349  length2( m1 - other.at( other.size() * 2.0 / 3.0 ) ) +
350  length2( m2 - other.at( other.size() * 1.0 / 3.0 ) ) +
351  length2( end - other.front() );
352  distance /= 4.0;
353  inverseDistance /= 4.0;
354  if( inverseDistance < distance )
355  {
356  cit->reverseOrder();
357  }
358  }
359 }
360 
361 boost::shared_ptr< WFiber > WFiberCluster::getCenterLine() const
362 {
363  if( !m_centerLine )
364  {
366  }
367  return m_centerLine;
368 }
369 
370 boost::shared_ptr< WFiber > WFiberCluster::getLongestLine() const
371 {
372  if( !m_longestLine )
373  {
375  }
376  return m_longestLine;
377 }
378 
380 {
381  WBoundingBox result;
382  for( std::list< size_t >::const_iterator cit = m_memberIndices.begin(); cit != m_memberIndices.end(); ++cit )
383  {
384  result.expandBy( computeBoundingBox( m_fibs->at( *cit ) ) );
385  }
386  return result;
387 }