WMMarchingCubes.cpp 41.8 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
//---------------------------------------------------------------------------
//
// Project: OpenWalnut ( http://www.openwalnut.org )
//
// Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
// For more information see http://www.openwalnut.org/copying
//
// This file is part of OpenWalnut.
//
// OpenWalnut is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// OpenWalnut is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
//
//---------------------------------------------------------------------------

#include <iostream>
26
#include <fstream>
Alexander Wiebel's avatar
Alexander Wiebel committed
27
#include <string>
28
#include <vector>
29

30 31 32
#include <cmath>


33
#include "WMMarchingCubes.h"
34
#include "iso_surface.xpm"
35
#include "marchingCubesCaseTables.h"
36
#include "../../common/datastructures/WTriangleMesh.h"
37
#include "../../common/WLimits.h"
38

39 40
#include <osg/Geode>
#include <osg/Geometry>
41 42 43
#include <osg/StateSet>
#include <osg/StateAttribute>
#include <osg/PolygonMode>
44 45
#include <osg/LightModel>
#include <osgDB/WriteFile>
46

Alexander Wiebel's avatar
Alexander Wiebel committed
47
#include "../../common/WProgress.h"
48
#include "../../common/WPreferences.h"
49
#include "../../common/math/WVector3D.h"
50
#include "../../common/math/WLinearAlgebraFunctions.h"
51
#include "../../dataHandler/WDataHandler.h"
52 53
#include "../../dataHandler/WSubject.h"
#include "../../dataHandler/WGridRegular3D.h"
54
#include "../../dataHandler/WDataTexture3D.h"
55
#include "../../kernel/WKernel.h"
56
#include "../../graphicsEngine/WShader.h"
57 58 59

#include "../data/WMData.h"

60

61
WMMarchingCubes::WMMarchingCubes():
62
    WModule(),
63
    m_recompute( boost::shared_ptr< WCondition >( new WCondition() ) ),
64 65 66 67 68 69
    m_nCellsX( 0 ),
    m_nCellsY( 0 ),
    m_nCellsZ( 0 ),
    m_fCellLengthX( 1 ),
    m_fCellLengthY( 1 ),
    m_fCellLengthZ( 1 ),
70
    m_tIsoLevel( 100 ),
71
    m_dataSet(),
72
    m_shaderUseLighting( false ),
73
    m_shaderUseTransparency( false ),
74 75
    m_moduleNode( new WGEGroupNode() ),
    m_surfaceGeode( 0 )
76 77 78 79 80
{
    // WARNING: initializing connectors inside the constructor will lead to an exception.
    // Implement WModule::initializeConnectors instead.
}

81
WMMarchingCubes::~WMMarchingCubes()
82 83 84 85 86
{
    // cleanup
    removeConnectors();
}

87 88 89 90 91
boost::shared_ptr< WModule > WMMarchingCubes::factory() const
{
    return boost::shared_ptr< WModule >( new WMMarchingCubes() );
}

92 93 94 95 96
const char** WMMarchingCubes::getXPMIcon() const
{
    return iso_surface_xpm;
}

97
const std::string WMMarchingCubes::getName() const
98
{
99
    return "Isosurface";
100 101
}

102
const std::string WMMarchingCubes::getDescription() const
103
{
104
    return "This module implement the marching cubes"
105
" algorithm with a consistent triangulation. It allows to compute isosurfaces"
106
" for a given isovalue on data given on a grid only consisting of cubes. It yields"
Alexander Wiebel's avatar
Alexander Wiebel committed
107
" the surface as triangle soup.";
108 109
}

110
void WMMarchingCubes::moduleMain()
111
{
112
    // use the m_input "data changed" flag
113
    m_moduleState.setResetable( true, true );
114
    m_moduleState.add( m_input->getDataChangedCondition() );
115
    m_moduleState.add( m_recompute );
116

117 118 119
    // signal ready state
    ready();

120 121 122 123 124
    // now, to watch changing/new textures use WSubject's change condition
    boost::signals2::connection con = WDataHandler::getDefaultSubject()->getChangeCondition()->subscribeSignal(
            boost::bind( &WMMarchingCubes::notifyTextureChange, this )
    );

125 126
    // loop until the module container requests the module to quit
    while ( !m_shutdownFlag() )
127
    {
128 129 130 131
        // acquire data from the input connector
        m_dataSet = m_input->getData();
        if ( !m_dataSet.get() )
        {
132
            // OK, the output has not yet sent data
133
            // NOTE: see comment at the end of this while loop for m_moduleState
134
            debugLog() << "Waiting for data ...";
135 136 137
            m_moduleState.wait();
            continue;
        }
138

139 140
        // update ISO surface
        debugLog() << "Computing surface ...";
141

142
        generateSurfacePre( m_isoValueProp->get() );
143

144 145
        // TODO(wiebel): MC remove this from here
        //    renderMesh( load( "/tmp/isosurfaceTestMesh.vtk" ) );
146

147
        debugLog() << "Rendering surface ...";
148

149 150
        // settings for normal isosurface
        m_shaderUseLighting = true;
151
        m_shaderUseTransparency = true;
Alexander Wiebel's avatar
Alexander Wiebel committed
152

153 154
        renderSurface();
        debugLog() << "Done!";
Alexander Wiebel's avatar
Alexander Wiebel committed
155

156 157 158 159
        // this waits for m_moduleState to fire. By default, this is only the m_shutdownFlag condition.
        // NOTE: you can add your own conditions to m_moduleState using m_moduleState.add( ... )
        m_moduleState.wait();
    }
160 161
}

162
void WMMarchingCubes::connectors()
163 164
{
    // initialize connectors
165
    m_input = boost::shared_ptr< WModuleInputData < WDataSetSingle  > >(
166
        new WModuleInputData< WDataSetSingle >( shared_from_this(),
167 168 169 170 171 172
                                                               "in", "Dataset to compute isosurface for." )
        );

    // add it to the list of connectors. Please note, that a connector NOT added via addConnector will not work as expected.
    addConnector( m_input );

173 174 175 176 177
    m_output = boost::shared_ptr< WModuleOutputData< WTriangleMesh > >(
            new WModuleOutputData< WTriangleMesh >( shared_from_this(), "out", "The mesh representing the isosurface." ) );

    addConnector( m_output );

178 179 180 181
    // call WModules initialization
    WModule::connectors();
}

182
void WMMarchingCubes::properties()
183
{
184
    m_isoValueProp = m_properties2->addProperty( "Iso Value", "The surface will show the area that has this value.", 100., m_recompute );
185 186
    m_isoValueProp->setMin( wlimits::MIN_DOUBLE );
    m_isoValueProp->setMax( wlimits::MAX_DOUBLE );
187 188 189 190 191 192 193 194 195
    {
        // If set in config file use standard isovalue from config file
        double tmpIsoValue;
        if( WPreferences::getPreference( "modules.MC.isoValue", &tmpIsoValue ) )
        {
            m_isoValueProp->set( tmpIsoValue );
        }
    }

196 197 198
    m_opacityProp = m_properties2->addProperty( "Opacity %", "Opaqueness of surface.", 100 );
    m_opacityProp->setMin( 0 );
    m_opacityProp->setMax( 100 );
199

200
    m_useTextureProp = m_properties2->addProperty( "Use Texture", "Use texturing of the surface?", false );
201
}
202

203 204
void WMMarchingCubes::generateSurfacePre( double isoValue )
{
205
    debugLog() << "Isovalue: " << isoValue;
206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223
    switch( (*m_dataSet).getValueSet()->getDataType() )
    {
        case W_DT_UNSIGNED_CHAR:
        {
            boost::shared_ptr< WValueSet< unsigned char > > vals;
            vals =  boost::shared_dynamic_cast< WValueSet< unsigned char > >( ( *m_dataSet ).getValueSet() );
            assert( vals );
            generateSurface( ( *m_dataSet ).getGrid(), vals, isoValue );
            break;
        }
        case W_DT_INT16:
        {
            boost::shared_ptr< WValueSet< int16_t > > vals;
            vals =  boost::shared_dynamic_cast< WValueSet< int16_t > >( ( *m_dataSet ).getValueSet() );
            assert( vals );
            generateSurface( ( *m_dataSet ).getGrid(), vals, isoValue );
            break;
        }
224 225 226 227 228 229 230 231
        case W_DT_SIGNED_INT:
        {
            boost::shared_ptr< WValueSet< int32_t > > vals;
            vals =  boost::shared_dynamic_cast< WValueSet< int32_t > >( ( *m_dataSet ).getValueSet() );
            assert( vals );
            generateSurface( ( *m_dataSet ).getGrid(), vals, isoValue );
            break;
        }
232 233 234 235 236 237 238 239
        case W_DT_FLOAT:
        {
            boost::shared_ptr< WValueSet< float > > vals;
            vals =  boost::shared_dynamic_cast< WValueSet< float > >( ( *m_dataSet ).getValueSet() );
            assert( vals );
            generateSurface( ( *m_dataSet ).getGrid(), vals, isoValue );
            break;
        }
240 241 242 243 244 245 246 247
        case W_DT_DOUBLE:
        {
            boost::shared_ptr< WValueSet< double > > vals;
            vals =  boost::shared_dynamic_cast< WValueSet< double > >( ( *m_dataSet ).getValueSet() );
            assert( vals );
            generateSurface( ( *m_dataSet ).getGrid(), vals, isoValue );
            break;
        }
248 249 250 251 252 253
        default:
            assert( false && "Unknow data type in MarchingCubes module" );
    }
}


254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272
void WMMarchingCubes::transformPositions( ID2WPointXYZId* positions )
{
    wmath::WMatrix< double > mat = m_grid->getTransformationMatrix();
    for( ID2WPointXYZId::iterator it = positions->begin(); it != positions->end(); ++it )
    {
        wmath::WPosition pos = wmath::WPosition( it->second.x, it->second.y, it->second.z );

        std::vector< double > resultPos4D( 4 );
        resultPos4D[0] = mat( 0, 0 ) * pos[0] + mat( 0, 1 ) * pos[1] + mat( 0, 2 ) * pos[2] + mat( 0, 3 ) * 1;
        resultPos4D[1] = mat( 1, 0 ) * pos[0] + mat( 1, 1 ) * pos[1] + mat( 1, 2 ) * pos[2] + mat( 1, 3 ) * 1;
        resultPos4D[2] = mat( 2, 0 ) * pos[0] + mat( 2, 1 ) * pos[1] + mat( 2, 2 ) * pos[2] + mat( 2, 3 ) * 1;
        resultPos4D[3] = mat( 3, 0 ) * pos[0] + mat( 3, 1 ) * pos[1] + mat( 3, 2 ) * pos[2] + mat( 3, 3 ) * 1;

        it->second.x = resultPos4D[0] / resultPos4D[3];
        it->second.y = resultPos4D[1] / resultPos4D[3];
        it->second.z = resultPos4D[2] / resultPos4D[3];
    }
}

273 274 275
template< typename T > void WMMarchingCubes::generateSurface( boost::shared_ptr< WGrid > inGrid,
                                                              boost::shared_ptr< WValueSet< T > > vals,
                                                              double isoValue )
276
{
277
    assert( vals );
278

279 280 281
    m_idToVertices.clear();
    m_trivecTriangles.clear();

282
    boost::shared_ptr< WGridRegular3D > grid = boost::shared_dynamic_cast< WGridRegular3D >( inGrid );
283
    m_grid = grid;
284 285 286 287 288 289 290 291 292 293
    assert( grid );

    m_fCellLengthX = grid->getOffsetX();
    m_fCellLengthY = grid->getOffsetY();
    m_fCellLengthZ = grid->getOffsetZ();

    m_nCellsX = grid->getNbCoordsX() - 1;
    m_nCellsY = grid->getNbCoordsY() - 1;
    m_nCellsZ = grid->getNbCoordsZ() - 1;

294 295
    m_tIsoLevel = isoValue;

296 297 298
    unsigned int nX = m_nCellsX + 1;
    unsigned int nY = m_nCellsY + 1;

299 300 301

    unsigned int nPointsInSlice = nX * nY;

302
    boost::shared_ptr< WProgress > progress = boost::shared_ptr< WProgress >( new WProgress( "Marching Cubes", m_nCellsZ ) );
Alexander Wiebel's avatar
Alexander Wiebel committed
303
    m_progress->addSubProgress( progress );
304
    // Generate isosurface.
305
    for( unsigned int z = 0; z < m_nCellsZ; z++ )
306
    {
Alexander Wiebel's avatar
Alexander Wiebel committed
307
        ++*progress;
308
        for( unsigned int y = 0; y < m_nCellsY; y++ )
309
        {
310
            for( unsigned int x = 0; x < m_nCellsX; x++ )
311 312 313 314
            {
                // Calculate table lookup index from those
                // vertices which are below the isolevel.
                unsigned int tableIndex = 0;
315
                if( vals->getScalar( z * nPointsInSlice + y * nX + x ) < m_tIsoLevel )
316
                    tableIndex |= 1;
317
                if( vals->getScalar( z * nPointsInSlice + ( y + 1 ) * nX + x ) < m_tIsoLevel )
318
                    tableIndex |= 2;
319
                if( vals->getScalar( z * nPointsInSlice + ( y + 1 ) * nX + ( x + 1 ) ) < m_tIsoLevel )
320
                    tableIndex |= 4;
321
                if( vals->getScalar( z * nPointsInSlice + y * nX + ( x + 1 ) ) < m_tIsoLevel )
322
                    tableIndex |= 8;
323
                if( vals->getScalar( ( z + 1 ) * nPointsInSlice + y * nX + x ) < m_tIsoLevel )
324
                    tableIndex |= 16;
325
                if( vals->getScalar( ( z + 1 ) * nPointsInSlice + ( y + 1 ) * nX + x ) < m_tIsoLevel )
326
                    tableIndex |= 32;
327
                if( vals->getScalar( ( z + 1 ) * nPointsInSlice + ( y + 1 ) * nX + ( x + 1 ) ) < m_tIsoLevel )
328
                    tableIndex |= 64;
329
                if( vals->getScalar( ( z + 1 ) * nPointsInSlice + y * nX + ( x + 1 ) ) < m_tIsoLevel )
330 331 332 333
                    tableIndex |= 128;

                // Now create a triangulation of the isosurface in this
                // cell.
334
                if( wMarchingCubesCaseTables::edgeTable[tableIndex] != 0 )
335
                {
336
                    if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 8 )
337
                    {
338
                        WPointXYZId pt = calculateIntersection( vals, x, y, z, 3 );
339 340 341
                        unsigned int id = getEdgeID( x, y, z, 3 );
                        m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
                    }
342
                    if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 1 )
343
                    {
344
                        WPointXYZId pt = calculateIntersection( vals, x, y, z, 0 );
345 346 347
                        unsigned int id = getEdgeID( x, y, z, 0 );
                        m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
                    }
348
                    if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 256 )
349
                    {
350
                        WPointXYZId pt = calculateIntersection( vals, x, y, z, 8 );
351 352 353 354
                        unsigned int id = getEdgeID( x, y, z, 8 );
                        m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
                    }

355
                    if( x == m_nCellsX - 1 )
356
                    {
357
                        if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 4 )
358
                        {
359
                            WPointXYZId pt = calculateIntersection( vals, x, y, z, 2 );
360 361 362
                            unsigned int id = getEdgeID( x, y, z, 2 );
                            m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
                        }
363
                        if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 2048 )
364
                        {
365
                            WPointXYZId pt = calculateIntersection( vals, x, y, z, 11 );
366 367 368 369
                            unsigned int id = getEdgeID( x, y, z, 11 );
                            m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
                        }
                    }
370
                    if( y == m_nCellsY - 1 )
371
                    {
372
                        if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 2 )
373
                        {
374
                            WPointXYZId pt = calculateIntersection( vals, x, y, z, 1 );
375 376 377
                            unsigned int id = getEdgeID( x, y, z, 1 );
                            m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
                        }
378
                        if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 512 )
379
                        {
380
                            WPointXYZId pt = calculateIntersection( vals, x, y, z, 9 );
381 382 383 384
                            unsigned int id = getEdgeID( x, y, z, 9 );
                            m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
                        }
                    }
385
                    if( z == m_nCellsZ - 1 )
386
                    {
387
                        if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 16 )
388
                        {
389
                            WPointXYZId pt = calculateIntersection( vals, x, y, z, 4 );
390 391 392
                            unsigned int id = getEdgeID( x, y, z, 4 );
                            m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
                        }
393
                        if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 128 )
394
                        {
395
                            WPointXYZId pt = calculateIntersection( vals, x, y, z, 7 );
396 397 398 399
                            unsigned int id = getEdgeID( x, y, z, 7 );
                            m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
                        }
                    }
400
                    if( ( x == m_nCellsX - 1 ) && ( y == m_nCellsY - 1 ) )
401
                        if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 1024 )
402
                        {
403
                            WPointXYZId pt = calculateIntersection( vals, x, y, z, 10 );
404 405 406
                            unsigned int id = getEdgeID( x, y, z, 10 );
                            m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
                        }
407
                    if( ( x == m_nCellsX - 1 ) && ( z == m_nCellsZ - 1 ) )
408
                        if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 64 )
409
                        {
410
                            WPointXYZId pt = calculateIntersection( vals, x, y, z, 6 );
411 412 413
                            unsigned int id = getEdgeID( x, y, z, 6 );
                            m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
                        }
414
                    if( ( y == m_nCellsY - 1 ) && ( z == m_nCellsZ - 1 ) )
415
                        if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 32 )
416
                        {
417
                            WPointXYZId pt = calculateIntersection( vals, x, y, z, 5 );
418 419 420 421
                            unsigned int id = getEdgeID( x, y, z, 5 );
                            m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
                        }

422
                    for( int i = 0; wMarchingCubesCaseTables::triTable[tableIndex][i] != -1; i += 3 )
423 424 425
                    {
                        WMCTriangle triangle;
                        unsigned int pointID0, pointID1, pointID2;
426 427 428
                        pointID0 = getEdgeID( x, y, z, wMarchingCubesCaseTables::triTable[tableIndex][i] );
                        pointID1 = getEdgeID( x, y, z, wMarchingCubesCaseTables::triTable[tableIndex][i + 1] );
                        pointID2 = getEdgeID( x, y, z, wMarchingCubesCaseTables::triTable[tableIndex][i + 2] );
429 430 431 432 433 434 435 436 437
                        triangle.pointID[0] = pointID0;
                        triangle.pointID[1] = pointID1;
                        triangle.pointID[2] = pointID2;
                        m_trivecTriangles.push_back( triangle );
                    }
                }
            }
        }
    }
438 439

    transformPositions( &m_idToVertices );
Alexander Wiebel's avatar
Alexander Wiebel committed
440
    progress->finish();
441 442
}

443 444 445
template< typename T > WPointXYZId WMMarchingCubes::calculateIntersection( boost::shared_ptr< WValueSet< T > > vals,
                                                                           unsigned int nX, unsigned int nY, unsigned int nZ,
                                                                           unsigned int nEdgeNo )
446
{
447 448 449 450 451 452 453 454 455 456 457 458 459 460 461
    double x1;
    double y1;
    double z1;

    double x2;
    double y2;
    double z2;

    unsigned int v1x = nX;
    unsigned int v1y = nY;
    unsigned int v1z = nZ;

    unsigned int v2x = nX;
    unsigned int v2y = nY;
    unsigned int v2z = nZ;
462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533

    switch ( nEdgeNo )
    {
        case 0:
            v2y += 1;
            break;
        case 1:
            v1y += 1;
            v2x += 1;
            v2y += 1;
            break;
        case 2:
            v1x += 1;
            v1y += 1;
            v2x += 1;
            break;
        case 3:
            v1x += 1;
            break;
        case 4:
            v1z += 1;
            v2y += 1;
            v2z += 1;
            break;
        case 5:
            v1y += 1;
            v1z += 1;
            v2x += 1;
            v2y += 1;
            v2z += 1;
            break;
        case 6:
            v1x += 1;
            v1y += 1;
            v1z += 1;
            v2x += 1;
            v2z += 1;
            break;
        case 7:
            v1x += 1;
            v1z += 1;
            v2z += 1;
            break;
        case 8:
            v2z += 1;
            break;
        case 9:
            v1y += 1;
            v2y += 1;
            v2z += 1;
            break;
        case 10:
            v1x += 1;
            v1y += 1;
            v2x += 1;
            v2y += 1;
            v2z += 1;
            break;
        case 11:
            v1x += 1;
            v2x += 1;
            v2z += 1;
            break;
    }

    x1 = v1x * m_fCellLengthX;
    y1 = v1y * m_fCellLengthY;
    z1 = v1z * m_fCellLengthZ;
    x2 = v2x * m_fCellLengthX;
    y2 = v2y * m_fCellLengthY;
    z2 = v2z * m_fCellLengthZ;

534
    unsigned int nPointsInSlice = ( m_nCellsX + 1 ) * ( m_nCellsY + 1 );
535 536
    double val1 = vals->getScalar( v1z * nPointsInSlice + v1y * ( m_nCellsX + 1 ) + v1x );
    double val2 = vals->getScalar( v2z * nPointsInSlice + v2y * ( m_nCellsX + 1 ) + v2x );
537

538 539 540 541 542
    WPointXYZId intersection = interpolate( x1, y1, z1, x2, y2, z2, val1, val2 );
    intersection.newID = 0;
    return intersection;
}

543
WPointXYZId WMMarchingCubes::interpolate( double fX1, double fY1, double fZ1, double fX2, double fY2, double fZ2,
544 545 546 547 548
                                             double tVal1, double tVal2 )
{
    WPointXYZId interpolation;
    double mu;

Alexander Wiebel's avatar
Alexander Wiebel committed
549
    mu = static_cast<double>( ( m_tIsoLevel - tVal1 ) ) / ( tVal2 - tVal1 );
550 551 552 553
    interpolation.x = fX1 + mu * ( fX2 - fX1 );
    interpolation.y = fY1 + mu * ( fY2 - fY1 );
    interpolation.z = fZ1 + mu * ( fZ2 - fZ1 );
    interpolation.newID = 0;
554

555 556 557
    return interpolation;
}

558
int WMMarchingCubes::getEdgeID( unsigned int nX, unsigned int nY, unsigned int nZ, unsigned int nEdgeNo )
559 560 561 562
{
    switch ( nEdgeNo )
    {
        case 0:
563
            return 3 * getVertexID( nX, nY, nZ ) + 1;
564
        case 1:
565
            return 3 * getVertexID( nX, nY + 1, nZ );
566
        case 2:
567
            return 3 * getVertexID( nX + 1, nY, nZ ) + 1;
568
        case 3:
569
            return 3 * getVertexID( nX, nY, nZ );
570
        case 4:
571
            return 3 * getVertexID( nX, nY, nZ + 1 ) + 1;
572
        case 5:
573
            return 3 * getVertexID( nX, nY + 1, nZ + 1 );
574
        case 6:
575
            return 3 * getVertexID( nX + 1, nY, nZ + 1 ) + 1;
576
        case 7:
577
            return 3 * getVertexID( nX, nY, nZ + 1 );
578
        case 8:
579
            return 3 * getVertexID( nX, nY, nZ ) + 2;
580
        case 9:
581
            return 3 * getVertexID( nX, nY + 1, nZ ) + 2;
582
        case 10:
583
            return 3 * getVertexID( nX + 1, nY + 1, nZ ) + 2;
584
        case 11:
585
            return 3 * getVertexID( nX + 1, nY, nZ ) + 2;
586 587 588 589 590 591
        default:
            // Invalid edge no.
            return -1;
    }
}

592
unsigned int WMMarchingCubes::getVertexID( unsigned int nX, unsigned int nY, unsigned int nZ )
593
{
594
    return nZ * ( m_nCellsY + 1 ) * ( m_nCellsX + 1) + nY * ( m_nCellsX + 1 ) + nX;
595
}
596

597
void WMMarchingCubes::renderSurface()
598 599
{
    unsigned int nextID = 0;
600 601 602 603
    boost::shared_ptr< WTriangleMesh > triMesh( new WTriangleMesh() );
    triMesh->clearMesh();
    triMesh->resizeVertices( m_idToVertices.size() );
    triMesh->resizeTriangles( m_trivecTriangles.size() );
604

605
    // TODO(wiebel): MC what is this for?
606 607 608 609 610
    float xOff = 0.5f;
    float yOff = 0.5f;
    float zOff = 0.5f;

    // Rename vertices.
611
    ID2WPointXYZId::iterator mapIterator = m_idToVertices.begin();
612 613 614
    while ( mapIterator != m_idToVertices.end() )
    {
        ( *mapIterator ).second.newID = nextID;
615 616 617
        triMesh->fastAddVert( wmath::WPosition( ( *mapIterator ).second.x + xOff,
                                                ( *mapIterator ).second.y + yOff,
                                                ( *mapIterator ).second.z + zOff ) );
618 619 620 621 622
        nextID++;
        mapIterator++;
    }

    // Now rename triangles.
623
    WMCTriangleVECTOR::iterator vecIterator = m_trivecTriangles.begin();
624 625
    while ( vecIterator != m_trivecTriangles.end() )
    {
626
        for( unsigned int i = 0; i < 3; i++ )
627 628 629 630
        {
            unsigned int newID = m_idToVertices[( *vecIterator ).pointID[i]].newID;
            ( *vecIterator ).pointID[i] = newID;
        }
631
        triMesh->fastAddTriangle( ( *vecIterator ).pointID[0], ( *vecIterator ).pointID[1],
632 633 634 635
                                 ( *vecIterator ).pointID[2] );
        vecIterator++;
    }

636 637 638
    renderMesh( triMesh );
    m_triMesh = triMesh;
    m_output->updateData( m_triMesh );
639 640 641 642 643 644 645 646 647

    // TODO(wiebel): MC make the filename set automatically
//     bool saved = save( "/tmp/isosurfaceTestMesh.vtk", triMesh );
//     if( saved )
//     {
//         WLogger::getLogger()->addLogMessage( "File opened", "Marching Cubes", LL_DEBUG );
//     }
}

648

649
void WMMarchingCubes::renderMesh( boost::shared_ptr< WTriangleMesh > mesh )
650
{
651 652
//    WKernel::getRunningKernel()->getGraphicsEngine()->getScene()
    m_moduleNode->remove( m_surfaceGeode );
653
    osg::Geometry* surfaceGeometry = new osg::Geometry();
654
    m_surfaceGeode = osg::ref_ptr< osg::Geode >( new osg::Geode );
655

656
    osg::Vec3Array* vertices = new osg::Vec3Array;
657
    for( size_t i = 0; i < mesh->getNumVertices(); ++i )
658 659
    {
        wmath::WPosition vertPos;
660
        vertPos = mesh->getVertex( i );
661 662 663 664 665 666 667
        vertices->push_back( osg::Vec3( vertPos[0], vertPos[1], vertPos[2] ) );
    }
    surfaceGeometry->setVertexArray( vertices );

    osg::DrawElementsUInt* surfaceElement;

    surfaceElement = new osg::DrawElementsUInt( osg::PrimitiveSet::TRIANGLES, 0 );
668
    for( unsigned int triId = 0; triId < mesh->getNumTriangles(); ++triId )
669 670 671
    {
        for( unsigned int vertId = 0; vertId < 3; ++vertId )
        {
672
            surfaceElement->push_back( mesh->getTriangleVertexId( triId, vertId ) );
673 674 675 676
        }
    }
    surfaceGeometry->addPrimitiveSet( surfaceElement );

677
    // ------------------------------------------------
678 679
    // normals
    osg::ref_ptr< osg::Vec3Array> normals( new osg::Vec3Array() );
680 681 682

    mesh->computeVertNormals(); // time consuming
    for( unsigned int vertId = 0; vertId < mesh->getNumVertices(); ++vertId )
683
    {
684
        wmath::WVector3D tmpNormal = mesh->getVertexNormal( vertId );
685
        normals->push_back( osg::Vec3( tmpNormal[0], tmpNormal[1], tmpNormal[2] ) );
686 687 688 689
    }
    surfaceGeometry->setNormalArray( normals.get() );
    surfaceGeometry->setNormalBinding( osg::Geometry::BIND_PER_VERTEX );

690 691
    m_surfaceGeode->addDrawable( surfaceGeometry );
    osg::StateSet* state = m_surfaceGeode->getOrCreateStateSet();
692

693 694
    // ------------------------------------------------
    // colors
695
    osg::Vec4Array* colors = new osg::Vec4Array;
696

697
    colors->push_back( osg::Vec4( .9f, .9f, 0.9f, 1.0f ) );
698 699 700
    surfaceGeometry->setColorArray( colors );
    surfaceGeometry->setColorBinding( osg::Geometry::BIND_OVERALL );

701 702 703
    osg::ref_ptr<osg::LightModel> lightModel = new osg::LightModel();
    lightModel->setTwoSided( true );
    state->setAttributeAndModes( lightModel.get(), osg::StateAttribute::ON );
704
    state->setMode(  GL_BLEND, osg::StateAttribute::ON  );
705

706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774
    // ------------------------------------------------
    // Shader stuff

    // TODO(wiebel): fix texture coords.
    double xext = m_grid->getOffsetX() * m_grid->getNbCoordsX();
    double yext = m_grid->getOffsetY() * m_grid->getNbCoordsY();
    double zext = m_grid->getOffsetZ() * m_grid->getNbCoordsZ();

    osg::Vec3Array* texCoords = new osg::Vec3Array;
    for( size_t i = 0; i < mesh->getNumVertices(); ++i )
    {
        wmath::WPosition vertPos;
        vertPos = mesh->getVertex( i );
        texCoords->push_back( osg::Vec3( vertPos[0]/xext, vertPos[1]/yext, vertPos[2]/zext ) );
    }
    surfaceGeometry->setTexCoordArray( 0, texCoords );

    m_typeUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "type0", 0 ) ) );
    m_typeUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "type1", 0 ) ) );
    m_typeUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "type2", 0 ) ) );
    m_typeUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "type3", 0 ) ) );
    m_typeUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "type4", 0 ) ) );
    m_typeUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "type5", 0 ) ) );
    m_typeUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "type6", 0 ) ) );
    m_typeUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "type7", 0 ) ) );
    m_typeUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "type8", 0 ) ) );
    m_typeUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "type9", 0 ) ) );

    m_alphaUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "alpha0", 1.0f ) ) );
    m_alphaUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "alpha1", 1.0f ) ) );
    m_alphaUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "alpha2", 1.0f ) ) );
    m_alphaUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "alpha3", 1.0f ) ) );
    m_alphaUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "alpha4", 1.0f ) ) );
    m_alphaUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "alpha5", 1.0f ) ) );
    m_alphaUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "alpha6", 1.0f ) ) );
    m_alphaUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "alpha7", 1.0f ) ) );
    m_alphaUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "alpha8", 1.0f ) ) );
    m_alphaUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "alpha9", 1.0f ) ) );

    m_thresholdUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "threshold0", 0.0f ) ) );
    m_thresholdUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "threshold1", 0.0f ) ) );
    m_thresholdUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "threshold2", 0.0f ) ) );
    m_thresholdUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "threshold3", 0.0f ) ) );
    m_thresholdUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "threshold4", 0.0f ) ) );
    m_thresholdUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "threshold5", 0.0f ) ) );
    m_thresholdUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "threshold6", 0.0f ) ) );
    m_thresholdUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "threshold7", 0.0f ) ) );
    m_thresholdUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "threshold8", 0.0f ) ) );
    m_thresholdUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "threshold9", 0.0f ) ) );

    m_samplerUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "tex0", 0 ) ) );
    m_samplerUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "tex1", 1 ) ) );
    m_samplerUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "tex2", 2 ) ) );
    m_samplerUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "tex3", 3 ) ) );
    m_samplerUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "tex4", 4 ) ) );
    m_samplerUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "tex5", 5 ) ) );
    m_samplerUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "tex6", 6 ) ) );
    m_samplerUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "tex7", 7 ) ) );
    m_samplerUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "tex8", 8 ) ) );
    m_samplerUniforms.push_back( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "tex9", 9 ) ) );

    for ( int i = 0; i < 10; ++i )
    {
        state->addUniform( m_typeUniforms[i] );
        state->addUniform( m_thresholdUniforms[i] );
        state->addUniform( m_alphaUniforms[i] );
        state->addUniform( m_samplerUniforms[i] );
    }

775
    state->addUniform( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "useTexture", m_properties->getValue< bool >( "Use Texture" ) ) ) );
776
    state->addUniform( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "useLighting", m_shaderUseLighting ) ) );
777 778
    if( m_shaderUseTransparency )
    {
779
        state->addUniform( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "opacity", m_opacityProp->get( true ) ) ) );
780 781 782 783 784
    }
    else
    {
        state->addUniform( osg::ref_ptr<osg::Uniform>( new osg::Uniform( "opacity", 100 ) ) );
    }
785

786
    // NOTE: the following code should not be necessary. The update callback does this job just before the mesh is rendered
787
    // initially. Just set the texture changed flag to true. If this however might be needed use WSubject::getDataTextures.
788 789 790 791 792 793 794 795
    m_textureChanged = true;
    //std::vector< boost::shared_ptr< WDataSet > > dsl = WKernel::getRunningKernel()->getGui()->getDataSetList( 0, true );*/
    //if ( dsl.size() > 0 )
    //{
        //for ( int i = 0; i < 10; ++i )
        //{
            //m_typeUniforms[i]->set( 0 );
        //}
796

797 798 799 800
        //int c = 0;
        //for ( size_t i = 0; i < dsl.size(); ++i )
        //{
            //osg::ref_ptr<osg::Texture3D> texture3D = dsl[i]->getTexture()->getTexture();
801

802
            //state->setTextureAttributeAndModes( c, texture3D, osg::StateAttribute::ON );
803

804 805
            //float t = dsl[i]->getTexture()->getThreshold()/ 100.0;
            //float a = dsl[i]->getTexture()->getAlpha();
806

807 808 809
            //m_typeUniforms[c]->set( boost::shared_dynamic_cast<WDataSetSingle>( dsl[i] )->getValueSet()->getDataType() );
            //m_thresholdUniforms[c]->set( t );
            //m_alphaUniforms[c]->set( a );
810

811 812 813
            //++c;
        //}
    //}
814

815
    m_shader = osg::ref_ptr< WShader > ( new WShader( "surface" ) );
816
    m_shader->apply( m_surfaceGeode );
817

818 819
    m_moduleNode->insert( m_surfaceGeode );
    WKernel::getRunningKernel()->getGraphicsEngine()->getScene()->insert( m_moduleNode );
820

821
    m_moduleNode->addUpdateCallback( new SurfaceNodeCallback( boost::shared_dynamic_cast< WMMarchingCubes >( shared_from_this() ) ) );
822

823
    //osgDB::writeNodeFile( *m_surfaceGeode, "/tmp/saved.osg" ); //for debugging
824 825
}

826 827 828 829 830 831 832 833 834 835 836 837 838
// TODO(wiebel): move this somewhere, where more classes can use it