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