29 #include "../common/WLimits.h"
30 #include "WThreadedTrackingFunction.h"
37 WAssert( dataset->getGrid(),
"" );
40 WAssert( g->encloses( job.first ),
"" );
45 if( fabs( length( dir ) - 1.0 ) > TRACKING_EPS )
53 WAssert( t > 0.0,
"" );
54 WAssert(
onBoundary( g, job.first + dir * t ),
"" );
56 if( !g->encloses( job.first + dir * t ) )
66 while(
onBoundary( g, job.first + dir * t ) && i < 50 )
69 if( !g->encloses( job.first + dir * t ) )
93 double c = dot( d, v ) / grid->getOffsetX() + 0.5;
94 p = c >= 0.0 ? modf( c, &z ) : 1.0 + modf( c, &z );
95 if( fabs( p ) < TRACKING_EPS || fabs( p - 1 ) < TRACKING_EPS )
99 d = grid->getDirectionY();
101 c = dot( d, v ) / grid->getOffsetY() + 0.5;
102 p = c >= 0.0 ? modf( c, &z ) : 1.0 + modf( c, &z );
103 if( fabs( p ) < TRACKING_EPS || fabs( p - 1 ) < TRACKING_EPS )
107 d = grid->getDirectionZ();
109 c = dot( d, v ) / grid->getOffsetZ() + 0.5;
110 p = c >= 0.0 ? modf( c, &z ) : 1.0 + modf( c, &z );
111 if( fabs( p ) < TRACKING_EPS || fabs( p - 1 ) < TRACKING_EPS )
124 double dist = std::numeric_limits< double >::infinity();
129 WVector3d o( grid->getOffsetX(), grid->getOffsetY(), grid->getOffsetZ() );
130 WVector3d s[ 3 ] = { grid->getDirectionX() / o[ 0 ], grid->getDirectionY() / o[ 1 ], grid->getDirectionZ() / o[ 2 ] };
131 double c = dot( s[ 0 ], v ) / o[ 0 ] + 0.5;
132 p[ 0 ] = c >= 0.0 ? modf( c, &z ) : 1.0 + modf( c, &z );
133 c = dot( s[ 1 ], v ) / o[ 1 ] + 0.5;
134 p[ 1 ] = c >= 0.0 ? modf( c, &z ) : 1.0 + modf( c, &z );
135 c = dot( s[ 2 ], v ) / o[ 2 ] + 0.5;
136 p[ 2 ] = c >= 0.0 ? modf( c, &z ) : 1.0 + modf( c, &z );
137 WVector3d d( dot( dir, s[ 0 ] ), dot( dir, s[ 1 ] ), dot( dir, s[ 2 ] ) );
140 for(
int i = 0; i < 3; ++i )
142 WAssert( -TRACKING_EPS < p[ i ] && p[ i ] < o[ i ] + TRACKING_EPS,
"" );
143 for(
double j = 0.0; j <= 1.0; j += 1.0 )
145 if( fabs( d[ i ] ) > TRACKING_EPS )
147 t = ( j - p[ i ] ) * o[ i ] / d[ i ];
148 if( t > 0.0 && t < dist )
156 WAssert( dist >= TRACKING_EPS,
"" );
165 std::size_t seedPositions, std::size_t seedsPerPos,
166 std::vector< int > v0, std::vector< int > v1 )
168 m_grid( boost::shared_dynamic_cast<
GridType >( dataset->getGrid() ) ),
169 m_directionFunc( dirFunc ),
170 m_nextPosFunc( nextFunc ),
171 m_fiberVisitor( fiberVst ),
172 m_pointVisitor( pointVst ),
179 throw WException( std::string(
"Cannot find WGridRegular3D. Are you sure the dataset has the correct grid type?" ) );
182 m_maxPoints =
static_cast< std::size_t
>( 5 * pow( static_cast< double >(
m_grid->size() ), 1.0 / 3.0 ) );
194 if( t->get().done() )
200 job = t->get().job();
212 std::vector< WVector3d > fiber;
213 if( fabs( length( e ) - 1.0 ) > TRACKING_EPS )
222 fiber.push_back( j.first );
231 if( fabs( length( j.second ) - 1.0 ) > TRACKING_EPS )
239 fiber.push_back( j.first );
247 std::reverse( fiber.begin(), fiber.end() );
252 if( fabs( length( j.second ) - 1.0 ) > TRACKING_EPS )
260 fiber.push_back( j.first );
281 std::vector< int >
const& v1, std::size_t seedPositions,
282 std::size_t seedsPerPosition )
292 WAssert( v0[ 0 ] >= 1 && static_cast< std::size_t >( v0[ 0 ] ) <
m_grid->getNbCoordsX() - 1,
"" );
293 WAssert( v0[ 1 ] >= 1 && static_cast< std::size_t >( v0[ 1 ] ) <
m_grid->getNbCoordsY() - 1,
"" );
294 WAssert( v0[ 2 ] >= 1 && static_cast< std::size_t >( v0[ 2 ] ) <
m_grid->getNbCoordsZ() - 1,
"" );
295 for(
int i = 0; i < 3; ++i )
297 m_min[ i ] = v0[ i ] * seedPositions;
302 m_max[ 0 ] = (
m_grid->getNbCoordsX() - 1 ) * seedPositions;
303 m_max[ 1 ] = (
m_grid->getNbCoordsY() - 1 ) * seedPositions;
304 m_max[ 2 ] = (
m_grid->getNbCoordsZ() - 1 ) * seedPositions;
308 WAssert( v1[ 0 ] > v0[ 0 ] && static_cast< std::size_t >( v1[ 0 ] ) < grid->getNbCoordsX(),
"" );
309 WAssert( v1[ 1 ] > v0[ 1 ] && static_cast< std::size_t >( v1[ 1 ] ) < grid->getNbCoordsY(),
"" );
310 WAssert( v1[ 2 ] > v0[ 2 ] && static_cast< std::size_t >( v1[ 2 ] ) < grid->getNbCoordsZ(),
"" );
311 for(
int i = 0; i < 3; ++i )
313 m_max[ i ] = v1[ i ] * seedPositions;
318 m_max[ 3 ] = seedsPerPosition;
329 for(
int i = 3; i > -1; --i )
332 if( m_pos[ i ] == m_max[ i ] )
334 m_pos[ i ] = m_min[ i ];
355 job.first =
m_grid->getOrigin() +
m_grid->getDirectionX() * ( m_offset * ( 0.5 + m_pos[ 0 ] ) - 0.5 )
356 +
m_grid->getDirectionY() * ( m_offset * ( 0.5 + m_pos[ 1 ] ) - 0.5 )
357 +
m_grid->getDirectionZ() * ( m_offset * ( 0.5 + m_pos[ 2 ] ) - 0.5 );