← Back to branch summary

~ubuntu-branches/ubuntu/vivid/openwalnut/vivid

« back to all changes in this revision

Viewing changes to .pc/boost153/src/core/dataHandler/WThreadedTrackingFunction.cpp

  • Committer: Package Import Robot
  • Author(s): Dmitrijs Ledkovs
  • Date: 2013-05-24 03:12:03 UTC
  • Revision ID: package-import@ubuntu.com-20130524031203-l5g1lzm1vd83fupi
Tags: 1.3.1+hg5849-1ubuntu1
Cherrypick boost1.53 pointer cast fixes.

Show diffs side-by-side

added added

removed removed

Lines of Context:
.pc/boost153/src/core/dataHandler/WThreadedTrackingFunction.cpp
 
1
//---------------------------------------------------------------------------
 
2
//
 
3
// Project: OpenWalnut ( http://www.openwalnut.org )
 
4
//
 
5
// Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
 
6
// For more information see http://www.openwalnut.org/copying
 
7
//
 
8
// This file is part of OpenWalnut.
 
9
//
 
10
// OpenWalnut is free software: you can redistribute it and/or modify
 
11
// it under the terms of the GNU Lesser General Public License as published by
 
12
// the Free Software Foundation, either version 3 of the License, or
 
13
// (at your option) any later version.
 
14
//
 
15
// OpenWalnut is distributed in the hope that it will be useful,
 
16
// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
17
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
18
// GNU Lesser General Public License for more details.
 
19
//
 
20
// You should have received a copy of the GNU Lesser General Public License
 
21
// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
 
22
//
 
23
//---------------------------------------------------------------------------
 
24
 
 
25
#include <limits>
 
26
#include <string>
 
27
#include <vector>
 
28
 
 
29
#include "../common/WLimits.h"
 
30
#include "WThreadedTrackingFunction.h"
 
31
 
 
32
namespace wtracking
 
33
{
 
34
    bool WTrackingUtility::followToNextVoxel( DataSetPtr dataset, JobType& job, DirFunc const& dirFunc )
 
35
    {
 
36
        // TODO(reichenbach): Give those WAsserts correct strings incase they occur
 
37
        WAssert( dataset->getGrid(), "" );
 
38
        Grid3DPtr g = boost::shared_dynamic_cast< WGridRegular3D >( dataset->getGrid() );
 
39
        WAssert( g, "" );
 
40
        WAssert( g->encloses( job.first ), "" );
 
41
 
 
42
        // find matching direction
 
43
        WVector3d dir = dirFunc( dataset, job );
 
44
 
 
45
        if( fabs( length( dir ) - 1.0 ) > TRACKING_EPS )
 
46
        {
 
47
            return false;
 
48
        }
 
49
 
 
50
        // find t such that job.first() + t * dir is a point on the boundary of the current voxel
 
51
        double t = getDistanceToBoundary( g, job.first, dir );
 
52
        WAssert( !wlimits::isInf( t ) && !wlimits::isNaN( t ), "Warning in WTrackingUtility::followToNextVoxel NaN's or INF's occured" );
 
53
        WAssert( t > 0.0, "" );
 
54
        WAssert( onBoundary( g, job.first + dir * t ), "" );
 
55
 
 
56
        if( !g->encloses( job.first + dir * t ) )
 
57
        {
 
58
            return false;
 
59
        }
 
60
 
 
61
        // increase t, so that job.first() + t * dir does not lie on any boundary anymore
 
62
        // ( so we can find the corresponding voxel without ambiguities )
 
63
        // ( this also means the resulting new position will never be in the same
 
64
        //   voxel as the starting position )
 
65
        int i = 0;
 
66
        while( onBoundary( g, job.first + dir * t ) && i < 50 )
 
67
        {
 
68
            t += TRACKING_EPS;
 
69
            if( !g->encloses( job.first + dir * t ) )
 
70
            {
 
71
                return false;
 
72
            }
 
73
            ++i;
 
74
        }
 
75
        if( i == 50 )
 
76
        {
 
77
            return false;
 
78
        }
 
79
 
 
80
        // set the new position
 
81
        job.first += t * dir;
 
82
        job.second = dir;
 
83
        return true;
 
84
    }
 
85
 
 
86
    bool WTrackingUtility::onBoundary( Grid3DPtr grid, WVector3d const& pos )
 
87
    {
 
88
        WVector3d v = pos - grid->getOrigin();
 
89
        double p;
 
90
        double z;
 
91
        WVector3d d = grid->getDirectionX();
 
92
        d = normalize( d );
 
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 )
 
96
        {
 
97
            return true;
 
98
        }
 
99
        d = grid->getDirectionY();
 
100
        d = normalize( d );
 
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 )
 
104
        {
 
105
            return true;
 
106
        }
 
107
        d = grid->getDirectionZ();
 
108
        d = normalize( d );
 
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 )
 
112
        {
 
113
            return true;
 
114
        }
 
115
        return false;
 
116
    }
 
117
 
 
118
    double WTrackingUtility::getDistanceToBoundary( Grid3DPtr grid, WVector3d const& pos, WVector3d const& dir )
 
119
    {
 
120
        // the input pos is at least TRACKING_EPS away from the boundary
 
121
        // so distance calculations don't get screwed
 
122
        WAssert( !onBoundary( grid, pos ), "" );
 
123
 
 
124
        double dist = std::numeric_limits< double >::infinity();
 
125
 
 
126
        WVector3d v = pos - grid->getOrigin();
 
127
        WVector3d p;
 
128
        double z;
 
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 ] ) );
 
138
 
 
139
        double t;
 
140
        for( int i = 0; i < 3; ++i )
 
141
        {
 
142
            WAssert( -TRACKING_EPS < p[ i ] && p[ i ] < o[ i ] + TRACKING_EPS, "" );
 
143
            for( double j = 0.0; j <= 1.0; j += 1.0 )
 
144
            {
 
145
                if( fabs( d[ i ] ) > TRACKING_EPS )
 
146
                {
 
147
                    t = ( j - p[ i ] ) * o[ i ] / d[ i ];
 
148
                    if( t > 0.0 && t < dist )
 
149
                    {
 
150
                        dist = t;
 
151
                    }
 
152
                }
 
153
            }
 
154
        }
 
155
        // distance absolute must be at least TRACKING_EPS!
 
156
        WAssert( dist >= TRACKING_EPS, "" );
 
157
        return dist;
 
158
    }
 
159
 
 
160
    //////////////////////////////////////////////////////////////////////////////////////////
 
161
 
 
162
    WThreadedTrackingFunction::WThreadedTrackingFunction( DataSetPtr dataset, DirFunc dirFunc,
 
163
            NextPositionFunc nextFunc,
 
164
            FiberVisitorFunc fiberVst, PointVisitorFunc pointVst,
 
165
            std::size_t seedPositions, std::size_t seedsPerPos,
 
166
            std::vector< int > v0, std::vector< int > v1 )
 
167
        : Base( dataset ),
 
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 ),
 
173
        m_maxPoints(),
 
174
        m_currentIndex()
 
175
    {
 
176
        // dataset != 0 is tested by the base constructor
 
177
        if( !m_grid )
 
178
        {
 
179
            throw WException( std::string( "Cannot find WGridRegular3D. Are you sure the dataset has the correct grid type?" ) );
 
180
        }
 
181
 
 
182
        m_maxPoints = static_cast< std::size_t >( 5 * pow( static_cast< double >( m_grid->size() ), 1.0 / 3.0 ) );
 
183
 
 
184
        m_currentIndex.getWriteTicket()->get() = IndexType( m_grid, v0, v1, seedPositions, seedsPerPos );
 
185
    }
 
186
 
 
187
    WThreadedTrackingFunction::~WThreadedTrackingFunction()
 
188
    {
 
189
    }
 
190
 
 
191
    bool WThreadedTrackingFunction::getJob( JobType& job )  // NOLINT
 
192
    {
 
193
        WSharedObject< IndexType >::WriteTicket t = m_currentIndex.getWriteTicket();
 
194
        if( t->get().done() )
 
195
        {
 
196
            return false;
 
197
        }
 
198
        else
 
199
        {
 
200
            job = t->get().job();
 
201
            ++t->get();
 
202
            return true;
 
203
        }
 
204
    }
 
205
 
 
206
    void WThreadedTrackingFunction::compute( DataSetPtr input, JobType const& job )
 
207
    {
 
208
        WVector3d e = m_directionFunc( input, job );
 
209
        JobType j = job;
 
210
        j.second = e;
 
211
 
 
212
        std::vector< WVector3d > fiber;
 
213
        if( fabs( length( e ) - 1.0 ) > TRACKING_EPS )
 
214
        {
 
215
            if( m_fiberVisitor )
 
216
            {
 
217
                m_fiberVisitor( fiber );
 
218
            }
 
219
            return;
 
220
        }
 
221
 
 
222
        fiber.push_back( j.first );
 
223
        if( m_pointVisitor )
 
224
        {
 
225
            m_pointVisitor( j.first );
 
226
        }
 
227
 
 
228
        // forward integration
 
229
        for( std::size_t k = 0; k < m_maxPoints; ++k )
 
230
        {
 
231
            if( fabs( length( j.second ) - 1.0 ) > TRACKING_EPS )
 
232
            {
 
233
                break;
 
234
            }
 
235
            if( !m_nextPosFunc( input, j, m_directionFunc ) )
 
236
            {
 
237
                break;
 
238
            }
 
239
            fiber.push_back( j.first );
 
240
            if( m_pointVisitor )
 
241
            {
 
242
                m_pointVisitor( j.first );
 
243
            }
 
244
        }
 
245
 
 
246
        // backward integration
 
247
        std::reverse( fiber.begin(), fiber.end() );
 
248
        j = job;
 
249
        j.second = e * -1.0;
 
250
        for( std::size_t k = 0; k < m_maxPoints; ++k )
 
251
        {
 
252
            if( fabs( length( j.second ) - 1.0 ) > TRACKING_EPS )
 
253
            {
 
254
                break;
 
255
            }
 
256
            if( !m_nextPosFunc( input, j, m_directionFunc ) )
 
257
            {
 
258
                break;
 
259
            }
 
260
            fiber.push_back( j.first );
 
261
            if( m_pointVisitor )
 
262
            {
 
263
                m_pointVisitor( j.first );
 
264
            }
 
265
        }
 
266
 
 
267
        // output result
 
268
        if( m_fiberVisitor )
 
269
        {
 
270
            m_fiberVisitor( fiber );
 
271
        }
 
272
    }
 
273
 
 
274
    WThreadedTrackingFunction::IndexType::IndexType()
 
275
        :   m_grid(),
 
276
        m_done( true )
 
277
    {
 
278
    }
 
279
 
 
280
    WThreadedTrackingFunction::IndexType::IndexType( GridPtr grid, std::vector< int > const& v0,
 
281
            std::vector< int > const& v1, std::size_t seedPositions,
 
282
            std::size_t seedsPerPosition )
 
283
        : m_grid( grid ),
 
284
        m_done( false )
 
285
    {
 
286
        if( v0.size() != 3 )
 
287
        {
 
288
            m_min[ 0 ] = m_min[ 1 ] = m_min[ 2 ] = seedPositions;
 
289
        }
 
290
        else
 
291
        {
 
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 )
 
296
            {
 
297
                m_min[ i ] = v0[ i ] * seedPositions;
 
298
            }
 
299
        }
 
300
        if( v1.size() != 3 )
 
301
        {
 
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;
 
305
        }
 
306
        else
 
307
        {
 
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 )
 
312
            {
 
313
                m_max[ i ] = v1[ i ] * seedPositions;
 
314
            }
 
315
        }
 
316
        m_offset = 1.0 / seedPositions;
 
317
        m_min[ 3 ] = 0;
 
318
        m_max[ 3 ] = seedsPerPosition;
 
319
        m_pos[ 0 ] = m_min[ 0 ];
 
320
        m_pos[ 1 ] = m_min[ 1 ];
 
321
        m_pos[ 2 ] = m_min[ 2 ];
 
322
        m_pos[ 3 ] = 0;
 
323
    }
 
324
 
 
325
    WThreadedTrackingFunction::IndexType& WThreadedTrackingFunction::IndexType::operator++ ()
 
326
    {
 
327
        if( !m_done )
 
328
        {
 
329
            for( int i = 3; i > -1; --i )
 
330
            {
 
331
                ++m_pos[ i ];
 
332
                if( m_pos[ i ] == m_max[ i ] )
 
333
                {
 
334
                    m_pos[ i ] = m_min[ i ];
 
335
                    m_done = i == 0;
 
336
                }
 
337
                else
 
338
                {
 
339
                    break;
 
340
                };
 
341
            }
 
342
        }
 
343
        return *this;
 
344
    }
 
345
 
 
346
    bool WThreadedTrackingFunction::IndexType::done()
 
347
    {
 
348
        return m_done;
 
349
    }
 
350
 
 
351
    WThreadedTrackingFunction::JobType WThreadedTrackingFunction::IndexType::job()
 
352
    {
 
353
        JobType job;
 
354
        job.second = WVector3d( 0.0, 0.0, 0.0 );
 
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 );
 
358
        return job;
 
359
    }
 
360
 
 
361
}  // namespace wtracking