OpenWalnut  1.4.0
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 #include <algorithm>
28 
29 #include <boost/shared_ptr.hpp>
30 
31 #include "../../common/math/WMath.h"
32 #include "../../common/math/WPlane.h"
33 #include "../../common/WLimits.h"
34 #include "../../common/WTransferable.h"
35 #include "../WDataSetFiberVector.h"
36 #include "WFiberCluster.h"
37 
38 // The prototype as singleton. Created during first getPrototype() call
39 boost::shared_ptr< WPrototyped > WFiberCluster::m_prototype = boost::shared_ptr< WPrototyped >();
40 
42  : WTransferable(),
43  m_centerLineCreationLock( new boost::shared_mutex() ),
44  m_longestLineCreationLock( new boost::shared_mutex() )
45 {
46 }
47 
49  : WTransferable(),
50  m_centerLineCreationLock( new boost::shared_mutex() ),
51  m_longestLineCreationLock( new boost::shared_mutex() )
52 {
53  m_memberIndices.push_back( index );
54 }
55 
56 WFiberCluster::WFiberCluster( const WFiberCluster::IndexList& indices, const WColor& color )
57  : WTransferable(),
58  m_memberIndices( indices ),
59  m_color( color ),
60  m_centerLineCreationLock( new boost::shared_mutex() ),
61  m_longestLineCreationLock( new boost::shared_mutex() )
62 {
63  // init
64 }
65 
67  WFiberCluster::IndexListConstIterator indicesEnd, const WColor& color )
68  : WTransferable(),
69  m_color( color ),
70  m_centerLineCreationLock( new boost::shared_mutex() ),
71  m_longestLineCreationLock( new boost::shared_mutex() )
72 {
73  // now copy the index list
74  std::copy( indicesBegin, indicesEnd, m_memberIndices.begin() );
75 }
76 
78  : WTransferable( other ),
79  m_memberIndices( other.m_memberIndices ),
80  m_fibs( other.m_fibs ),
81  m_color( other.m_color ),
82  m_centerLineCreationLock( new boost::shared_mutex() ), // do not copy the mutex as both instances of WFiberCluster can be modifed at the
83  // same time
84  m_longestLineCreationLock( new boost::shared_mutex() ),
85  m_centerLine(), // << we can't ensure that the centerline and longest line are initialized yet, but we want a deep copy
86  m_longestLine()
87 {
88  // copy them only if they exist
89  if( other.m_centerLine )
90  {
91  m_centerLine = boost::shared_ptr< WFiber >( new WFiber( *other.m_centerLine.get() ) );
92  }
93  if( other.m_longestLine )
94  {
95  m_longestLine = boost::shared_ptr< WFiber >( new WFiber( *other.m_longestLine.get() ) );
96  }
97 }
98 
100 {
103 }
104 
105 void WFiberCluster::merge( WFiberCluster& other ) // NOLINT
106 {
107  WFiberCluster::IndexList::const_iterator cit = other.m_memberIndices.begin();
108  for( ; cit != other.m_memberIndices.end(); ++cit)
109  {
110  m_memberIndices.push_back( *cit );
111  }
112  // make sure that those indices aren't occuring anywhere else
113  other.m_centerLine.reset(); // they are not valid anymore
114  other.m_longestLine.reset();
115  other.clear();
116 }
117 
120 {
121  // now copy the index list
122  m_memberIndices.insert( m_memberIndices.end(), indicesBegin, indicesEnd );
123 }
124 
125 // NODOXYGEN
126 // \cond Suppress_Doxygen
127 void WFiberCluster::setDataSetReference( boost::shared_ptr< const WDataSetFiberVector > fibs )
128 {
129  m_fibs = fibs;
130 }
131 
132 boost::shared_ptr< const WDataSetFiberVector > WFiberCluster::getDataSetReference() const
133 {
134  return m_fibs;
135 }
136 // TODO(math): The only reason why we store here a Reference to the fiber
137 // dataset is, we need it in the WMVoxelizer module as well as the clustering
138 // information. Since we don't have the possibility of multiple
139 // InputConnectors we must agglomerate those into one object. Please remove this.
140 // \endcond
141 
142 boost::shared_ptr< WPrototyped > WFiberCluster::getPrototype()
143 {
144  if( !m_prototype )
145  {
146  m_prototype = boost::shared_ptr< WPrototyped >( new WFiberCluster() );
147  }
148  return m_prototype;
149 }
150 
152 {
153  // ensure nobody changes the mutable m_centerline
154  boost::unique_lock< boost::shared_mutex > lock = boost::unique_lock< boost::shared_mutex >( *m_centerLineCreationLock );
155  // has the line been calculated while we waited?
156  if( m_centerLine )
157  {
158  lock.unlock();
159  return;
160  }
161 
162  // make copies of the fibers
163  boost::shared_ptr< WDataSetFiberVector > fibs( new WDataSetFiberVector() );
164  size_t avgFiberSize = 0;
165  for( WFiberCluster::IndexList::const_iterator cit = m_memberIndices.begin(); cit != m_memberIndices.end(); ++cit )
166  {
167  fibs->push_back( m_fibs->at( *cit ) );
168  avgFiberSize += fibs->back().size();
169  }
170  avgFiberSize /= fibs->size();
171 
172  unifyDirection( fibs );
173 
174  for( WDataSetFiberVector::iterator cit = fibs->begin(); cit != fibs->end(); ++cit )
175  {
176  cit->resampleByNumberOfPoints( avgFiberSize );
177  }
178 
179  m_centerLine = boost::shared_ptr< WFiber >( new WFiber() );
180  m_centerLine->reserve( avgFiberSize );
181  for( size_t i = 0; i < avgFiberSize; ++i )
182  {
183  WPosition avgPosition( 0, 0, 0 );
184  for( WDataSetFiberVector::const_iterator cit = fibs->begin(); cit != fibs->end(); ++cit )
185  {
186  avgPosition += cit->at( i );
187  }
188  avgPosition /= static_cast< double >( fibs->size() );
189  m_centerLine->push_back( avgPosition );
190  }
191 
193  lock.unlock();
194 }
195 
197 {
198  // ensure nobody changes the mutable m_longestline
199  boost::unique_lock< boost::shared_mutex > lock = boost::unique_lock< boost::shared_mutex >( *m_longestLineCreationLock );
200  // has the line been calculated while we waited?
201  if( m_longestLine )
202  {
203  lock.unlock();
204  return;
205  }
206 
207  m_longestLine = boost::shared_ptr< WFiber >( new WFiber() );
208 
209  // empty datasets can be ignored
210  if( m_fibs->size() == 0 )
211  {
212  return;
213  }
214 
215  size_t longest = 0;
216  size_t longestID = 0;
217  for( size_t cit = 0; cit < m_fibs->size(); ++cit )
218  {
219  if( m_fibs->at( cit ).size() > longest )
220  {
221  longest = m_fibs->at( cit ).size();
222  longestID = cit;
223  }
224  }
225 
226  for( WFiber::const_iterator cit = m_fibs->at( longestID ).begin(); cit != m_fibs->at( longestID ).end(); ++cit )
227  {
228  m_longestLine->push_back( *cit );
229  }
230 
231  lock.unlock();
232 }
233 
235 {
236  // first ending of the centerline
237  assert( m_centerLine->size() > 1 );
238  WFiber cL( *m_centerLine );
239  WPlane p( cL[0] - cL[1], cL[0] + ( cL[0] - cL[1] ) );
240  boost::shared_ptr< WPosition > cutPoint( new WPosition( 0, 0, 0 ) );
241  bool intersectionFound = true;
242 
243  // in the beginning all fibers participate
244  boost::shared_ptr< WDataSetFiberVector > fibs( new WDataSetFiberVector() );
245  for( WFiberCluster::IndexList::const_iterator cit = m_memberIndices.begin(); cit != m_memberIndices.end(); ++cit )
246  {
247  fibs->push_back( m_fibs->at( *cit ) );
248  }
249 
250  while( intersectionFound )
251  {
252  intersectionFound = false;
253  size_t intersectingFibers = 0;
254 // WPosition avg( 0, 0, 0 );
255  for( WDataSetFiberVector::iterator cit = fibs->begin(); cit != fibs->end(); )
256  {
257  if( intersectPlaneLineNearCP( p, *cit, cutPoint ) )
258  {
259  if( length( *cutPoint - p.getPosition() ) < 20 )
260  {
261 // avg += *cutPoint;
262  intersectingFibers++;
263  intersectionFound = true;
264  ++cit;
265  }
266  else
267  {
268  cit = fibs->erase( cit );
269  }
270  }
271  else
272  {
273  cit = fibs->erase( cit );
274  }
275  }
276  if( intersectingFibers > 10 )
277  {
278  cL.insert( cL.begin(), cL[0] + ( cL[0] - cL[1] ) );
279  p.resetPosition( cL[0] + ( cL[0] - cL[1] ) );
280 // avg[0] /= static_cast< double >( intersectingFibers );
281 // avg[1] /= static_cast< double >( intersectingFibers );
282 // avg[2] /= static_cast< double >( intersectingFibers );
283 // cL.insert( cL.begin(), 0.995 * ( cL[0] + ( cL[0] - cL[1] ) ) + 0.005 * avg );
284 // p.resetPosition( cL[0] + ( cL[0] - cL[1] ) );
285 // p.setNormal( ( cL[0] - cL[1] ) );
286  }
287  else // no intersections found => abort
288  {
289  break;
290  }
291  }
292  // second ending of the centerline
293  boost::shared_ptr< WDataSetFiberVector > fobs( new WDataSetFiberVector() );
294  for( WFiberCluster::IndexList::const_iterator cit = m_memberIndices.begin(); cit != m_memberIndices.end(); ++cit )
295  {
296  fobs->push_back( m_fibs->at( *cit ) );
297  }
298 
299  // try to discard other lines from other end
300 
301  WPlane q( cL.back() - cL[ cL.size() - 2 ], cL.back() + ( cL.back() - cL[ cL.size() - 2 ] ) );
302  intersectionFound = true;
303  while( intersectionFound )
304  {
305  intersectionFound = false;
306  size_t intersectingFibers = 0;
307 // WPosition avg( 0, 0, 0 );
308  for( WDataSetFiberVector::iterator cit = fobs->begin(); cit != fobs->end(); )
309  {
310  if( intersectPlaneLineNearCP( q, *cit, cutPoint ) )
311  {
312  if( length( *cutPoint - q.getPosition() ) < 20 )
313  {
314 // avg += *cutPoint;
315  intersectingFibers++;
316  intersectionFound = true;
317  ++cit;
318  }
319  else
320  {
321  cit = fobs->erase( cit );
322  }
323  }
324  else
325  {
326  cit = fobs->erase( cit );
327  }
328  }
329  if( intersectingFibers > 10 )
330  {
331  cL.push_back( cL.back() + ( cL.back() - cL[ cL.size() - 2 ] ) );
332  q.resetPosition( cL.back() + ( cL.back() - cL[ cL.size() - 2 ] ) );
333 // avg[0] /= static_cast< double >( intersectingFibers );
334 // avg[1] /= static_cast< double >( intersectingFibers );
335 // avg[2] /= static_cast< double >( intersectingFibers );
336 // cL.push_back( 0.995 * ( cL.back() + ( cL.back() - cL[ cL.size() - 2 ] ) ) + 0.005 * avg );
337 // q.resetPosition( cL.back() + ( cL.back() - cL[ cL.size() - 2 ] ) );
338 // q.setNormal( cL.back() - cL[ cL.size() - 2 ] );
339  }
340  else // no intersections found => abort
341  {
342  break;
343  }
344  }
345  *m_centerLine = cL;
346 }
347 
348 void WFiberCluster::unifyDirection( boost::shared_ptr< WDataSetFiberVector > fibs ) const
349 {
350  if( fibs->size() < 2 )
351  {
352  return;
353  }
354 
355  assert( !( fibs->at( 0 ).empty() ) && "WFiberCluster.unifyDirection: Empty fiber processed.. aborting" );
356 
357  // first fiber defines direction
358  const WFiber& firstFib = fibs->front();
359  const WPosition start = firstFib.front();
360  const WPosition m1 = firstFib.at( firstFib.size() * 1.0 / 3.0 );
361  const WPosition m2 = firstFib.at( firstFib.size() * 2.0 / 3.0 );
362  const WPosition end = firstFib.back();
363 
364  for( WDataSetFiberVector::iterator cit = fibs->begin() + 1; cit != fibs->end(); ++cit )
365  {
366  const WFiber& other = *cit;
367  double distance = length2( start - other.front() ) +
368  length2( m1 - other.at( other.size() * 1.0 / 3.0 ) ) +
369  length2( m2 - other.at( other.size() * 2.0 / 3.0 ) ) +
370  length2( end - other.back() );
371  double inverseDistance = length2( start - other.back() ) +
372  length2( m1 - other.at( other.size() * 2.0 / 3.0 ) ) +
373  length2( m2 - other.at( other.size() * 1.0 / 3.0 ) ) +
374  length2( end - other.front() );
375  distance /= 4.0;
376  inverseDistance /= 4.0;
377  if( inverseDistance < distance )
378  {
379  cit->reverseOrder();
380  }
381  }
382 }
383 
384 boost::shared_ptr< WFiber > WFiberCluster::getCenterLine() const
385 {
386  if( !m_centerLine )
387  {
389  }
390  return m_centerLine;
391 }
392 
393 boost::shared_ptr< WFiber > WFiberCluster::getLongestLine() const
394 {
395  if( !m_longestLine )
396  {
398  }
399  return m_longestLine;
400 }
401 
403 {
404  WBoundingBox result;
405  for( WFiberCluster::IndexList::const_iterator cit = m_memberIndices.begin(); cit != m_memberIndices.end(); ++cit )
406  {
407  result.expandBy( computeBoundingBox( m_fibs->at( *cit ) ) );
408  }
409  return result;
410 }
virtual ~WFiberCluster()
Destructs.
Represents a neural pathway.
Definition: WFiber.h:39
vector_type::const_iterator const_iterator
Compares to std::vector type.
Definition: WMixinVector.h:87
boost::shared_ptr< WFiber > m_centerLine
Average fiber for this cluster representing the main direction and curvature of this cluster...
const_reference at(size_type index) const
Wrapper around std::vector member function.
Definition: WMixinVector.h:413
void elongateCenterLine() const
The centerline may be shortened due to the averaging of outliers.
void expandBy(const WBoundingBoxImpl< VT > &bb)
Expands this bounding box to include the given bounding box.
Definition: WBoundingBox.h:253
Represents a cluster of indices of a WDataSetFiberVector.
Definition: WFiberCluster.h:47
boost::shared_ptr< WFiber > getCenterLine() const
Returns the center line of this cluster.
std::list< size_t > IndexList
This is the list of indices of fibers.
Definition: WFiberCluster.h:64
boost::shared_mutex * m_centerLineCreationLock
Lock the modification in the m_centerLine mutable.
vector_type::iterator iterator
Compares to std::vector type.
Definition: WMixinVector.h:92
This only is a 3d double vector.
Represents a simple set of WFibers.
WBoundingBox getBoundingBox() const
Recomputes on every call the axis aligned bounding box incorporating all tracts in this cluster...
Represents a plane with a normal vector and a position in space.
Definition: WPlane.h:41
const WPosition & getPosition() const
Returns a point in that plane.
Definition: WPlane.h:183
void resetPosition(WPosition newPos)
Reset the position of the plane, normal remains the same.
Definition: WPlane.cpp:62
void merge(WFiberCluster &other)
Merge the fibers of the other cluster with the fibers of this cluster.
Class building the interface for classes that might be transferred using WModuleConnector.
Definition: WTransferable.h:37
IndexList m_memberIndices
All indices in this set are members of this cluster.
void unifyDirection(boost::shared_ptr< WDataSetFiberVector > fibs) const
Alings all fibers within the given dataset to be in one main direction.
const_reference front() const
Wrapper around std::vector member function.
Definition: WMixinVector.h:557
void clear()
Make this cluster empty.
void generateCenterLine() const
Makes the hard work to compute the center line.
boost::shared_ptr< WFiber > m_longestLine
The longest fiber in the dataset.
static boost::shared_ptr< WPrototyped > getPrototype()
Returns a prototype instantiated with the true type of the deriving class.
boost::shared_mutex * m_longestLineCreationLock
Lock the modification in the m_longestLine mutable.
void generateLongestLine() const
Makes the hard work to find the longest line.
static boost::shared_ptr< WPrototyped > m_prototype
Prototype for this dataset.
boost::shared_ptr< WFiber > getLongestLine() const
Returns the center line of this cluster.
const_reference back() const
Wrapper around std::vector member function.
Definition: WMixinVector.h:537
boost::shared_ptr< const WDataSetFiberVector > m_fibs
Reference to the real fibers of the brain this cluster belongs to.
IndexList::const_iterator IndexListConstIterator
Const iterator on the index list.
Definition: WFiberCluster.h:69
size_type size() const
Wrapper around std::vector member function.
Definition: WMixinVector.h:267
WFiberCluster()
Constructs an empty cluster.