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