OpenWalnut  1.4.0
WFiberCluster.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 <list>
00026 #include <vector>
00027 #include <algorithm>
00028 
00029 #include <boost/shared_ptr.hpp>
00030 
00031 #include "../../common/math/WMath.h"
00032 #include "../../common/math/WPlane.h"
00033 #include "../../common/WLimits.h"
00034 #include "../../common/WTransferable.h"
00035 #include "../WDataSetFiberVector.h"
00036 #include "WFiberCluster.h"
00037 
00038 // The prototype as singleton. Created during first getPrototype() call
00039 boost::shared_ptr< WPrototyped > WFiberCluster::m_prototype = boost::shared_ptr< WPrototyped >();
00040 
00041 WFiberCluster::WFiberCluster()
00042     : WTransferable(),
00043     m_centerLineCreationLock( new boost::shared_mutex() ),
00044     m_longestLineCreationLock( new boost::shared_mutex() )
00045 {
00046 }
00047 
00048 WFiberCluster::WFiberCluster( size_t index )
00049     : WTransferable(),
00050     m_centerLineCreationLock( new boost::shared_mutex() ),
00051     m_longestLineCreationLock( new boost::shared_mutex() )
00052 {
00053     m_memberIndices.push_back( index );
00054 }
00055 
00056 WFiberCluster::WFiberCluster( const WFiberCluster::IndexList& indices, const WColor& color )
00057     : WTransferable(),
00058     m_memberIndices( indices ),
00059     m_color( color ),
00060     m_centerLineCreationLock( new boost::shared_mutex() ),
00061     m_longestLineCreationLock( new boost::shared_mutex() )
00062 {
00063     // init
00064 }
00065 
00066 WFiberCluster::WFiberCluster( WFiberCluster::IndexListConstIterator indicesBegin,
00067                               WFiberCluster::IndexListConstIterator indicesEnd, const WColor& color )
00068     : WTransferable(),
00069     m_color( color ),
00070     m_centerLineCreationLock( new boost::shared_mutex() ),
00071     m_longestLineCreationLock( new boost::shared_mutex() )
00072 {
00073     // now copy the index list
00074     std::copy( indicesBegin, indicesEnd, m_memberIndices.begin() );
00075 }
00076 
00077 WFiberCluster::WFiberCluster( const WFiberCluster& other )
00078     : WTransferable( other ),
00079     m_memberIndices( other.m_memberIndices ),
00080     m_fibs( other.m_fibs ),
00081     m_color( other.m_color ),
00082     m_centerLineCreationLock( new boost::shared_mutex() ),  // do not copy the mutex as both instances of WFiberCluster can be modifed at the
00083                                                             // same time
00084     m_longestLineCreationLock( new boost::shared_mutex() ),
00085     m_centerLine(),     // << we can't ensure that the centerline and longest line are initialized yet, but we want a deep copy
00086     m_longestLine()
00087 {
00088     // copy them only if they exist
00089     if( other.m_centerLine )
00090     {
00091         m_centerLine = boost::shared_ptr< WFiber >( new WFiber( *other.m_centerLine.get() ) );
00092     }
00093     if( other.m_longestLine )
00094     {
00095         m_longestLine = boost::shared_ptr< WFiber >( new WFiber( *other.m_longestLine.get() ) );
00096     }
00097 }
00098 
00099 WFiberCluster::~WFiberCluster()
00100 {
00101     delete m_centerLineCreationLock;
00102     delete m_longestLineCreationLock;
00103 }
00104 
00105 void WFiberCluster::merge( WFiberCluster& other ) // NOLINT
00106 {
00107     WFiberCluster::IndexList::const_iterator cit = other.m_memberIndices.begin();
00108     for( ; cit != other.m_memberIndices.end(); ++cit)
00109     {
00110         m_memberIndices.push_back( *cit );
00111     }
00112     // make sure that those indices aren't occuring anywhere else
00113     other.m_centerLine.reset();     // they are not valid anymore
00114     other.m_longestLine.reset();
00115     other.clear();
00116 }
00117 
00118 void WFiberCluster::merge( WFiberCluster::IndexListConstIterator indicesBegin,
00119                            WFiberCluster::IndexListConstIterator indicesEnd )
00120 {
00121     // now copy the index list
00122     m_memberIndices.insert( m_memberIndices.end(), indicesBegin, indicesEnd );
00123 }
00124 
00125 // NODOXYGEN
00126 // \cond Suppress_Doxygen
00127 void WFiberCluster::setDataSetReference( boost::shared_ptr< const WDataSetFiberVector > fibs )
00128 {
00129     m_fibs = fibs;
00130 }
00131 
00132 boost::shared_ptr< const WDataSetFiberVector > WFiberCluster::getDataSetReference() const
00133 {
00134     return m_fibs;
00135 }
00136 // TODO(math): The only reason why we store here a Reference to the fiber
00137 // dataset is, we need it in the WMVoxelizer module as well as the clustering
00138 // information. Since we don't have the possibility of multiple
00139 // InputConnectors we must agglomerate those into one object. Please remove this.
00140 // \endcond
00141 
00142 boost::shared_ptr< WPrototyped > WFiberCluster::getPrototype()
00143 {
00144     if( !m_prototype )
00145     {
00146         m_prototype = boost::shared_ptr< WPrototyped >( new WFiberCluster() );
00147     }
00148     return m_prototype;
00149 }
00150 
00151 void WFiberCluster::generateCenterLine() const
00152 {
00153     // ensure nobody changes the mutable m_centerline
00154     boost::unique_lock< boost::shared_mutex > lock = boost::unique_lock< boost::shared_mutex >( *m_centerLineCreationLock );
00155     // has the line been calculated while we waited?
00156     if( m_centerLine )
00157     {
00158         lock.unlock();
00159         return;
00160     }
00161 
00162     // make copies of the fibers
00163     boost::shared_ptr< WDataSetFiberVector > fibs( new WDataSetFiberVector() );
00164     size_t avgFiberSize = 0;
00165     for( WFiberCluster::IndexList::const_iterator cit = m_memberIndices.begin(); cit != m_memberIndices.end(); ++cit )
00166     {
00167         fibs->push_back( m_fibs->at( *cit ) );
00168         avgFiberSize += fibs->back().size();
00169     }
00170     avgFiberSize /= fibs->size();
00171 
00172     unifyDirection( fibs );
00173 
00174     for( WDataSetFiberVector::iterator cit = fibs->begin(); cit != fibs->end(); ++cit )
00175     {
00176         cit->resampleByNumberOfPoints( avgFiberSize );
00177     }
00178 
00179     m_centerLine = boost::shared_ptr< WFiber >( new WFiber() );
00180     m_centerLine->reserve( avgFiberSize );
00181     for( size_t i = 0; i < avgFiberSize; ++i )
00182     {
00183         WPosition avgPosition( 0, 0, 0 );
00184         for( WDataSetFiberVector::const_iterator cit = fibs->begin(); cit != fibs->end(); ++cit )
00185         {
00186             avgPosition += cit->at( i );
00187         }
00188         avgPosition /= static_cast< double >( fibs->size() );
00189         m_centerLine->push_back( avgPosition );
00190     }
00191 
00192     elongateCenterLine();
00193     lock.unlock();
00194 }
00195 
00196 void WFiberCluster::generateLongestLine() const
00197 {
00198     // ensure nobody changes the mutable m_longestline
00199     boost::unique_lock< boost::shared_mutex > lock = boost::unique_lock< boost::shared_mutex >( *m_longestLineCreationLock );
00200     // has the line been calculated while we waited?
00201     if( m_longestLine )
00202     {
00203         lock.unlock();
00204         return;
00205     }
00206 
00207     m_longestLine = boost::shared_ptr< WFiber >( new WFiber() );
00208 
00209     // empty datasets can be ignored
00210     if( m_fibs->size() == 0 )
00211     {
00212         return;
00213     }
00214 
00215     size_t longest = 0;
00216     size_t longestID = 0;
00217     for( size_t cit = 0; cit < m_fibs->size(); ++cit )
00218     {
00219         if( m_fibs->at( cit ).size() > longest )
00220         {
00221             longest = m_fibs->at( cit ).size();
00222             longestID = cit;
00223         }
00224     }
00225 
00226     for( WFiber::const_iterator cit = m_fibs->at( longestID ).begin(); cit != m_fibs->at( longestID ).end(); ++cit )
00227     {
00228         m_longestLine->push_back( *cit );
00229     }
00230 
00231     lock.unlock();
00232 }
00233 
00234 void WFiberCluster::elongateCenterLine() const
00235 {
00236     // first ending of the centerline
00237     assert( m_centerLine->size() > 1 );
00238     WFiber cL( *m_centerLine );
00239     WPlane p( cL[0] - cL[1], cL[0] + ( cL[0] - cL[1] ) );
00240     boost::shared_ptr< WPosition > cutPoint( new WPosition( 0, 0, 0 ) );
00241     bool intersectionFound = true;
00242 
00243     // in the beginning all fibers participate
00244     boost::shared_ptr< WDataSetFiberVector > fibs( new WDataSetFiberVector() );
00245     for( WFiberCluster::IndexList::const_iterator cit = m_memberIndices.begin(); cit != m_memberIndices.end(); ++cit )
00246     {
00247         fibs->push_back( m_fibs->at( *cit ) );
00248     }
00249 
00250     while( intersectionFound )
00251     {
00252         intersectionFound = false;
00253         size_t intersectingFibers = 0;
00254 //        WPosition avg( 0, 0, 0 );
00255         for( WDataSetFiberVector::iterator cit = fibs->begin(); cit != fibs->end(); )
00256         {
00257             if( intersectPlaneLineNearCP( p, *cit, cutPoint ) )
00258             {
00259                 if( length( *cutPoint - p.getPosition() ) < 20 )
00260                 {
00261 //                    avg += *cutPoint;
00262                     intersectingFibers++;
00263                     intersectionFound = true;
00264                     ++cit;
00265                 }
00266                 else
00267                 {
00268                     cit = fibs->erase( cit );
00269                 }
00270             }
00271             else
00272             {
00273                 cit = fibs->erase( cit );
00274             }
00275         }
00276         if( intersectingFibers > 10 )
00277         {
00278             cL.insert( cL.begin(), cL[0] + ( cL[0] - cL[1] ) );
00279             p.resetPosition( cL[0] + ( cL[0] -  cL[1] ) );
00280 //            avg[0] /= static_cast< double >( intersectingFibers );
00281 //            avg[1] /= static_cast< double >( intersectingFibers );
00282 //            avg[2] /= static_cast< double >( intersectingFibers );
00283 //            cL.insert( cL.begin(), 0.995 * ( cL[0] + ( cL[0] - cL[1] ) ) + 0.005 * avg );
00284 //            p.resetPosition( cL[0] + ( cL[0] -  cL[1] ) );
00285 //            p.setNormal( ( cL[0] -  cL[1] ) );
00286         }
00287         else // no intersections found => abort
00288         {
00289             break;
00290         }
00291     }
00292     // second ending of the centerline
00293     boost::shared_ptr< WDataSetFiberVector > fobs( new WDataSetFiberVector() );
00294     for( WFiberCluster::IndexList::const_iterator cit = m_memberIndices.begin(); cit != m_memberIndices.end(); ++cit )
00295     {
00296         fobs->push_back( m_fibs->at( *cit ) );
00297     }
00298 
00299     // try to discard other lines from other end
00300 
00301     WPlane q( cL.back() - cL[ cL.size() - 2 ], cL.back() + ( cL.back() - cL[ cL.size() - 2 ] ) );
00302     intersectionFound = true;
00303     while( intersectionFound )
00304     {
00305         intersectionFound = false;
00306         size_t intersectingFibers = 0;
00307 //        WPosition avg( 0, 0, 0 );
00308         for( WDataSetFiberVector::iterator cit = fobs->begin(); cit != fobs->end(); )
00309         {
00310             if( intersectPlaneLineNearCP( q, *cit, cutPoint ) )
00311             {
00312                 if( length( *cutPoint - q.getPosition() ) < 20 )
00313                 {
00314 //                    avg += *cutPoint;
00315                     intersectingFibers++;
00316                     intersectionFound = true;
00317                     ++cit;
00318                 }
00319                 else
00320                 {
00321                     cit = fobs->erase( cit );
00322                 }
00323             }
00324             else
00325             {
00326                 cit = fobs->erase( cit );
00327             }
00328         }
00329         if( intersectingFibers > 10 )
00330         {
00331             cL.push_back(  cL.back() + ( cL.back() - cL[ cL.size() - 2 ] ) );
00332             q.resetPosition(  cL.back() + ( cL.back() - cL[ cL.size() - 2 ] ) );
00333 //            avg[0] /= static_cast< double >( intersectingFibers );
00334 //            avg[1] /= static_cast< double >( intersectingFibers );
00335 //            avg[2] /= static_cast< double >( intersectingFibers );
00336 //            cL.push_back( 0.995 * ( cL.back() + ( cL.back() - cL[ cL.size() - 2 ] ) ) + 0.005 * avg );
00337 //            q.resetPosition( cL.back() + ( cL.back() - cL[ cL.size() - 2 ] ) );
00338 //            q.setNormal( cL.back() - cL[ cL.size() - 2 ] );
00339         }
00340         else // no intersections found => abort
00341         {
00342             break;
00343         }
00344     }
00345     *m_centerLine = cL;
00346 }
00347 
00348 void WFiberCluster::unifyDirection( boost::shared_ptr< WDataSetFiberVector > fibs ) const
00349 {
00350     if( fibs->size() < 2 )
00351     {
00352         return;
00353     }
00354 
00355     assert( !( fibs->at( 0 ).empty() ) && "WFiberCluster.unifyDirection: Empty fiber processed.. aborting" );
00356 
00357     // first fiber defines direction
00358     const WFiber& firstFib = fibs->front();
00359     const WPosition start = firstFib.front();
00360     const WPosition m1    = firstFib.at( firstFib.size() * 1.0 / 3.0 );
00361     const WPosition m2    = firstFib.at( firstFib.size() * 2.0 / 3.0 );
00362     const WPosition end   = firstFib.back();
00363 
00364     for( WDataSetFiberVector::iterator cit = fibs->begin() + 1; cit != fibs->end(); ++cit )
00365     {
00366         const WFiber& other = *cit;
00367         double        distance = length2( start - other.front() ) +
00368                                  length2( m1 - other.at( other.size() * 1.0 / 3.0 ) ) +
00369                                  length2( m2 - other.at( other.size() * 2.0 / 3.0 ) ) +
00370                                  length2( end - other.back() );
00371         double inverseDistance = length2( start - other.back() ) +
00372                                  length2( m1 - other.at( other.size() * 2.0 / 3.0 ) ) +
00373                                  length2( m2 - other.at( other.size() * 1.0 / 3.0 ) ) +
00374                                  length2( end - other.front() );
00375         distance        /= 4.0;
00376         inverseDistance /= 4.0;
00377         if( inverseDistance < distance )
00378         {
00379             cit->reverseOrder();
00380         }
00381     }
00382 }
00383 
00384 boost::shared_ptr< WFiber > WFiberCluster::getCenterLine() const
00385 {
00386     if( !m_centerLine )
00387     {
00388         generateCenterLine();
00389     }
00390     return m_centerLine;
00391 }
00392 
00393 boost::shared_ptr< WFiber > WFiberCluster::getLongestLine() const
00394 {
00395     if( !m_longestLine )
00396     {
00397         generateLongestLine();
00398     }
00399     return m_longestLine;
00400 }
00401 
00402 WBoundingBox WFiberCluster::getBoundingBox() const
00403 {
00404     WBoundingBox result;
00405     for( WFiberCluster::IndexList::const_iterator cit = m_memberIndices.begin(); cit != m_memberIndices.end(); ++cit )
00406     {
00407         result.expandBy( computeBoundingBox( m_fibs->at( *cit ) ) );
00408     }
00409     return result;
00410 }