OpenWalnut 1.3.1
WThreadedTrackingFunction.h
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 #ifndef WTHREADEDTRACKINGFUNCTION_H
00026 #define WTHREADEDTRACKINGFUNCTION_H
00027 
00028 #include <vector>
00029 #include <utility>
00030 
00031 #include <boost/array.hpp>
00032 
00033 #include "../common/math/linearAlgebra/WLinearAlgebra.h"
00034 #include "../common/WSharedObject.h"
00035 #include "../common/WThreadedJobs.h"
00036 
00037 #include "WDataSetSingle.h"
00038 
00039 class WThreadedTrackingFunctionTest; //! forward declaration
00040 
00041 
00042 /**
00043  * For tracking related functionality.
00044  */
00045 namespace wtracking // note that this is not final
00046 {
00047     // an epsilon value for various floating point comparisons
00048 #define TRACKING_EPS 0.0000001
00049 
00050     /**
00051      * \class WTrackingUtility
00052      *
00053      * A class that provides untility functions and typedefs for tracking algorithms.
00054      */
00055     class WTrackingUtility
00056     {
00057     public:
00058         //! define a job type for tracking algorithms
00059         typedef std::pair< WVector3d, WVector3d > JobType;
00060 
00061         //! the dataset type
00062         typedef WDataSetSingle DataSetType;
00063 
00064         //! a pointer to a dataset
00065         typedef boost::shared_ptr< DataSetType const > DataSetPtr;
00066 
00067         //! a function that calculates a direction to continue tracking
00068         typedef boost::function< WVector3d ( DataSetPtr, JobType const& ) > DirFunc;
00069 
00070         //! a pointer to a regular 3d grid
00071         // other grid types are not supported at the moment
00072         typedef boost::shared_ptr< WGridRegular3D > Grid3DPtr;
00073 
00074         /**
00075          * A function that follows a direction until leaving the current voxel.
00076          *
00077          * \param dataset A pointer to the input dataset.
00078          * \param job A pair of vectors, the position and the direction of the last integration.
00079          * \param dirFunc A function that computes the next direction.
00080          *
00081          * \return true, iff the calculated point is a valid position inside the grid
00082          */
00083         static bool followToNextVoxel( DataSetPtr dataset, JobType& job, DirFunc const& dirFunc );
00084 
00085         // one could add a runge-kutta-integrator too
00086 
00087         /**
00088          * Check if a point is on the boundary of the given grid, where boundary
00089          * means a distance less then TRACKING_EPS from any plane between
00090          * voxels. This does not check if the position is actually inside the grid.
00091          *
00092          * \param grid The grid.
00093          * \param pos The position to test.
00094          *
00095          * \return true, iff the position is on any voxel boundary
00096          */
00097         static bool onBoundary( Grid3DPtr grid, WVector3d const& pos );
00098 
00099         /**
00100          * Calculate the distance from a given position to the nearest voxel boundary
00101          * on the ray from the position along the given direction.
00102          *
00103          * \param grid The grid.
00104          * \param pos The starting position of the ray.
00105          * \param dir The normalized direction of the ray.
00106          *
00107          * \return The distance to the next voxel boundary.
00108          *
00109          * \note pos + getDistanceToBoundary( grid, pos, dir ) * dir will be a position on a voxel boundary
00110          */
00111         static double getDistanceToBoundary( Grid3DPtr grid, WVector3d const& pos, WVector3d const& dir );
00112     };
00113 
00114     //////////////////////////////////////////////////////////////////////////////////////////
00115 
00116     /**
00117      * \class WThreadedTrackingFunction
00118      *
00119      * Implements a generalized multithreaded tracking algorithm. A function that calculates the direction
00120      * and a function that calculates a new position have to be provided.
00121      *
00122      * Output values can be retrieved via two visitor functions that get called per fiber tracked and
00123      * per point calculated respectively.
00124      *
00125      * There are a certain number n of seeds per direction, this meens n*n*n seeds per voxel. For every
00126      * seed, m fibers get integrated. These two parameters are the seedPositions and seedsPerVoxel parameters
00127      * of the constructor, respectively.
00128      *
00129      * A 'cubic' region of the grid can be chosen for seeding. The v0 and v1 parameters of the constructor
00130      * are the starting/target voxel coords. Example:
00131      *
00132      * v0: 1, 1, 1
00133      * v1: 4, 5, 3
00134      *
00135      * In this case, only voxels between coords 1 to 3 in the x-direction, the voxels 1 to 4 in y- and the voxels 1 to 2
00136      * in z-direction are used for seeding.
00137      *
00138      * Note that voxels at the first (0) and last (grid->getNbCoords*()) position in any direction are
00139      * invalid seeding voxels as they are partially outside of the grid.
00140      */
00141     class WThreadedTrackingFunction : public WThreadedJobs< WTrackingUtility::DataSetType, WTrackingUtility::JobType >
00142     {
00143     //! make the test a friend
00144     friend class ::WThreadedTrackingFunctionTest;
00145 
00146     //! the job type
00147     typedef WTrackingUtility::JobType JobType;
00148 
00149     //! the dataset type
00150     typedef WTrackingUtility::DataSetType DataSetType;
00151 
00152     //! a pointer to a dataset
00153     typedef WTrackingUtility::DataSetPtr DataSetPtr;
00154 
00155     //! the grid type
00156     typedef WGridRegular3D GridType;
00157 
00158     //! a pointer to the grid
00159     typedef boost::shared_ptr< GridType > GridPtr;
00160 
00161     //! the direction calculation function
00162     typedef WTrackingUtility::DirFunc DirFunc;
00163 
00164     //! the path integration function
00165     typedef boost::function< bool ( DataSetPtr, JobType&, DirFunc const& ) > NextPositionFunc;
00166 
00167     //! a visitor function for fibers
00168     typedef boost::function< void ( std::vector< WVector3d > const& ) > FiberVisitorFunc;
00169 
00170     //! a visitor function type for points
00171     typedef boost::function< void ( WVector3d const& ) > PointVisitorFunc;
00172 
00173     //! the base class, a threaded job function
00174     typedef WThreadedJobs< DataSetType, JobType > Base;
00175 
00176     //! this type
00177     typedef WThreadedTrackingFunction This;
00178 
00179     public:
00180         /**
00181          * Constructor.
00182          *
00183          * \param dataset A pointer to a dataset.
00184          * \param dirFunc A direction calculation function.
00185          * \param nextFunc A position integration function.
00186          * \param fiberVst A visitor for fibers.
00187          * \param pointVst A visitor for points.
00188          * \param seedPositions The number of seed positions in every direction per voxel.
00189          * \param seedsPerPos The number of fibers startet from every seed position.
00190          * \param v0 A vector of starting voxel indices for every direction.
00191          * \param v1 A vector of target voxel indices for every direction.
00192          */
00193         WThreadedTrackingFunction( DataSetPtr dataset, DirFunc dirFunc, NextPositionFunc nextFunc,
00194                 FiberVisitorFunc fiberVst, PointVisitorFunc pointVst,
00195                 std::size_t seedPositions = 1, std::size_t seedsPerPos = 1,
00196                 std::vector< int > v0 = std::vector< int >(),
00197                 std::vector< int > v1 = std::vector< int >() );
00198 
00199         /**
00200          * Destructor.
00201          */
00202         virtual ~WThreadedTrackingFunction();
00203 
00204         /**
00205          * The job generator.
00206          *
00207          * \param job The next job (output).
00208          *
00209          * \return false, iff there are no more jobs.
00210          */
00211         virtual bool getJob( JobType& job ); // NOLINT
00212 
00213         /**
00214          * The calculation per job.
00215          *
00216          * \param input The input dataset.
00217          * \param job The job.
00218          */
00219         virtual void compute( DataSetPtr input, JobType const& job );
00220 
00221     private:
00222         /**
00223          * \class IndexType
00224          *
00225          * An index for seed positions.
00226          */
00227         class IndexType
00228         {
00229         friend class ::WThreadedTrackingFunctionTest;
00230         public:
00231             /**
00232              * Construct an invalid index.
00233              */
00234             IndexType();
00235 
00236             /**
00237              * Construct an index.
00238              *
00239              * \param grid The grid.
00240              * \param v0 A vector of starting voxel indices for every direction.
00241              * \param v1 A vector of target voxel indices for every direction.
00242              * \param seedPositions The number of seed positions in every direction per voxel.
00243              * \param seedsPerPosition The number of fibers startet from every seed position.
00244              */
00245             IndexType( GridPtr grid, std::vector< int > const& v0,
00246                     std::vector< int > const& v1, std::size_t seedPositions,
00247                     std::size_t seedsPerPosition );
00248 
00249             /**
00250              * Increase the index by one, effectively generating the next seed position.
00251              *
00252              * \return *this
00253              */
00254             IndexType& operator++ ();
00255 
00256             /**
00257              * Check if there aren't any more seed positions.
00258              *
00259              * \return true, iff there aren't any more seed positions.
00260              */
00261             bool done();
00262 
00263             /**
00264              * Create a job from this index.
00265              *
00266              * \return The job that is the current position.
00267              */
00268             JobType job();
00269 
00270         private:
00271             //! a pointer to the grid
00272             GridPtr m_grid;
00273 
00274             //! true, iff there are no more seeds
00275             bool m_done;
00276 
00277             //! the position in the seed space
00278             boost::array< std::size_t, 4 > m_pos;
00279 
00280             //! the minimum position in the seed space
00281             boost::array< std::size_t, 4 > m_min;
00282 
00283             //! the maximum position in the seed space
00284             boost::array< std::size_t, 4 > m_max;
00285 
00286             //! the relative (to the size of a voxel) distance between seeds
00287             double m_offset;
00288             };
00289 
00290             //! a pointer to the grid
00291             GridPtr m_grid;
00292 
00293             //! a function that returns the next direction
00294             DirFunc m_directionFunc;
00295 
00296             //! a function that calculates the next position
00297             NextPositionFunc m_nextPosFunc;
00298 
00299             //! the fiber visitor
00300             FiberVisitorFunc m_fiberVisitor;
00301 
00302             //! the point visitor
00303             PointVisitorFunc m_pointVisitor;
00304 
00305             //! the maximum number of points per forward/backward integration of a fiber
00306             std::size_t m_maxPoints;
00307 
00308             //! the current index/seed position
00309             WSharedObject< IndexType > m_currentIndex;
00310         };
00311 
00312 } /* namespace wtracking */
00313 
00314 #endif  // WTHREADEDTRACKINGFUNCTION_H