OpenWalnut  1.4.0
WMarchingLegoAlgorithm.cpp
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 #include <vector>
00026 
00027 #include "WMarchingLegoAlgorithm.h"
00028 
00029 WMarchingLegoAlgorithm::WMarchingLegoAlgorithm()
00030     : m_matrix( 4, 4 )
00031 {
00032 }
00033 
00034 WMarchingLegoAlgorithm::~WMarchingLegoAlgorithm()
00035 {
00036 }
00037 
00038 void WMarchingLegoAlgorithm::addSurface( size_t x, size_t y, size_t z, size_t surface )
00039 {
00040     WMLPointXYZId pt1;
00041     WMLPointXYZId pt2;
00042     WMLPointXYZId pt3;
00043     WMLPointXYZId pt4;
00044 
00045     pt1.newID = 0;
00046     pt2.newID = 0;
00047     pt3.newID = 0;
00048     pt4.newID = 0;
00049 
00050     switch( surface )
00051     {
00052         case 1:
00053         {
00054             pt1.x = x;
00055             pt1.y = y;
00056             pt1.z = z;
00057             unsigned int id1 = getVertexID( pt1.x , pt1.y, pt1.z );
00058             m_idToVertices.insert( ID2WMLPointXYZId::value_type( id1, pt1 ) );
00059 
00060             pt2.x = x;
00061             pt2.y = y + 1;
00062             pt2.z = z;
00063             unsigned int id2 = getVertexID( pt2.x , pt2.y, pt2.z );
00064             m_idToVertices.insert( ID2WMLPointXYZId::value_type( id2, pt2 ) );
00065 
00066             pt3.x = x;
00067             pt3.y = y + 1;
00068             pt3.z = z + 1;
00069             unsigned int id3 = getVertexID( pt3.x , pt3.y, pt3.z );
00070             m_idToVertices.insert( ID2WMLPointXYZId::value_type( id3, pt3 ) );
00071 
00072             pt4.x = x;
00073             pt4.y = y;
00074             pt4.z = z + 1;
00075             unsigned int id4 = getVertexID( pt4.x , pt4.y, pt4.z );
00076             m_idToVertices.insert( ID2WMLPointXYZId::value_type( id4, pt4 ) );
00077 
00078             WMLTriangle triangle1;
00079             triangle1.pointID[0] = id1;
00080             triangle1.pointID[1] = id2;
00081             triangle1.pointID[2] = id3;
00082             WMLTriangle triangle2;
00083             triangle2.pointID[0] = id3;
00084             triangle2.pointID[1] = id4;
00085             triangle2.pointID[2] = id1;
00086             m_trivecTriangles.push_back( triangle1 );
00087             m_trivecTriangles.push_back( triangle2 );
00088             break;
00089         }
00090         case 2:
00091         {
00092             pt1.x = x + 1;
00093             pt1.y = y;
00094             pt1.z = z;
00095             unsigned int id1 = getVertexID( pt1.x , pt1.y, pt1.z );
00096             m_idToVertices.insert( ID2WMLPointXYZId::value_type( id1, pt1 ) );
00097 
00098             pt2.x = x + 1;
00099             pt2.y = y;
00100             pt2.z = z + 1;
00101             unsigned int id2 = getVertexID( pt2.x , pt2.y, pt2.z );
00102             m_idToVertices.insert( ID2WMLPointXYZId::value_type( id2, pt2 ) );
00103 
00104             pt3.x = x + 1;
00105             pt3.y = y + 1;
00106             pt3.z = z + 1;
00107             unsigned int id3 = getVertexID( pt3.x , pt3.y, pt3.z );
00108             m_idToVertices.insert( ID2WMLPointXYZId::value_type( id3, pt3 ) );
00109 
00110             pt4.x = x + 1;
00111             pt4.y = y + 1;
00112             pt4.z = z;
00113             unsigned int id4 = getVertexID( pt4.x , pt4.y, pt4.z );
00114             m_idToVertices.insert( ID2WMLPointXYZId::value_type( id4, pt4 ) );
00115 
00116             WMLTriangle triangle1;
00117             triangle1.pointID[0] = id1;
00118             triangle1.pointID[1] = id2;
00119             triangle1.pointID[2] = id3;
00120             WMLTriangle triangle2;
00121             triangle2.pointID[0] = id3;
00122             triangle2.pointID[1] = id4;
00123             triangle2.pointID[2] = id1;
00124             m_trivecTriangles.push_back( triangle1 );
00125             m_trivecTriangles.push_back( triangle2 );
00126             break;
00127         }
00128         case 3:
00129         {
00130             pt1.x = x;
00131             pt1.y = y;
00132             pt1.z = z;
00133             unsigned int id1 = getVertexID( pt1.x , pt1.y, pt1.z );
00134             m_idToVertices.insert( ID2WMLPointXYZId::value_type( id1, pt1 ) );
00135 
00136             pt2.x = x;
00137             pt2.y = y;
00138             pt2.z = z + 1;
00139             unsigned int id2 = getVertexID( pt2.x , pt2.y, pt2.z );
00140             m_idToVertices.insert( ID2WMLPointXYZId::value_type( id2, pt2 ) );
00141 
00142             pt3.x = x + 1;
00143             pt3.y = y;
00144             pt3.z = z + 1;
00145             unsigned int id3 = getVertexID( pt3.x , pt3.y, pt3.z );
00146             m_idToVertices.insert( ID2WMLPointXYZId::value_type( id3, pt3 ) );
00147 
00148             pt4.x = x + 1;
00149             pt4.y = y;
00150             pt4.z = z;
00151             unsigned int id4 = getVertexID( pt4.x , pt4.y, pt4.z );
00152             m_idToVertices.insert( ID2WMLPointXYZId::value_type( id4, pt4 ) );
00153 
00154             WMLTriangle triangle1;
00155             triangle1.pointID[0] = id1;
00156             triangle1.pointID[1] = id2;
00157             triangle1.pointID[2] = id3;
00158             WMLTriangle triangle2;
00159             triangle2.pointID[0] = id3;
00160             triangle2.pointID[1] = id4;
00161             triangle2.pointID[2] = id1;
00162             m_trivecTriangles.push_back( triangle1 );
00163             m_trivecTriangles.push_back( triangle2 );
00164             break;
00165         }
00166         case 4:
00167         {
00168             pt1.x = x;
00169             pt1.y = y + 1;
00170             pt1.z = z;
00171             unsigned int id1 = getVertexID( pt1.x , pt1.y, pt1.z );
00172             m_idToVertices.insert( ID2WMLPointXYZId::value_type( id1, pt1 ) );
00173 
00174             pt2.x = x + 1;
00175             pt2.y = y + 1;
00176             pt2.z = z;
00177             unsigned int id2 = getVertexID( pt2.x , pt2.y, pt2.z );
00178             m_idToVertices.insert( ID2WMLPointXYZId::value_type( id2, pt2 ) );
00179 
00180             pt3.x = x + 1;
00181             pt3.y = y + 1;
00182             pt3.z = z + 1;
00183             unsigned int id3 = getVertexID( pt3.x , pt3.y, pt3.z );
00184             m_idToVertices.insert( ID2WMLPointXYZId::value_type( id3, pt3 ) );
00185 
00186             pt4.x = x;
00187             pt4.y = y + 1;
00188             pt4.z = z + 1;
00189             unsigned int id4 = getVertexID( pt4.x , pt4.y, pt4.z );
00190             m_idToVertices.insert( ID2WMLPointXYZId::value_type( id4, pt4 ) );
00191 
00192             WMLTriangle triangle1;
00193             triangle1.pointID[0] = id1;
00194             triangle1.pointID[1] = id2;
00195             triangle1.pointID[2] = id3;
00196             WMLTriangle triangle2;
00197             triangle2.pointID[0] = id3;
00198             triangle2.pointID[1] = id4;
00199             triangle2.pointID[2] = id1;
00200             m_trivecTriangles.push_back( triangle1 );
00201             m_trivecTriangles.push_back( triangle2 );
00202             break;
00203         }
00204         case 5:
00205         {
00206             pt1.x = x;
00207             pt1.y = y;
00208             pt1.z = z;
00209             unsigned int id1 = getVertexID( pt1.x , pt1.y, pt1.z );
00210             m_idToVertices.insert( ID2WMLPointXYZId::value_type( id1, pt1 ) );
00211 
00212             pt2.x = x + 1;
00213             pt2.y = y;
00214             pt2.z = z;
00215             unsigned int id2 = getVertexID( pt2.x , pt2.y, pt2.z );
00216             m_idToVertices.insert( ID2WMLPointXYZId::value_type( id2, pt2 ) );
00217 
00218             pt3.x = x + 1;
00219             pt3.y = y + 1;
00220             pt3.z = z;
00221             unsigned int id3 = getVertexID( pt3.x , pt3.y, pt3.z );
00222             m_idToVertices.insert( ID2WMLPointXYZId::value_type( id3, pt3 ) );
00223 
00224             pt4.x = x;
00225             pt4.y = y + 1;
00226             pt4.z = z;
00227             unsigned int id4 = getVertexID( pt4.x , pt4.y, pt4.z );
00228             m_idToVertices.insert( ID2WMLPointXYZId::value_type( id4, pt4 ) );
00229 
00230             WMLTriangle triangle1;
00231             triangle1.pointID[0] = id1;
00232             triangle1.pointID[1] = id2;
00233             triangle1.pointID[2] = id3;
00234             WMLTriangle triangle2;
00235             triangle2.pointID[0] = id3;
00236             triangle2.pointID[1] = id4;
00237             triangle2.pointID[2] = id1;
00238             m_trivecTriangles.push_back( triangle1 );
00239             m_trivecTriangles.push_back( triangle2 );
00240             break;
00241         }
00242         case 6:
00243         {
00244             pt1.x = x;
00245             pt1.y = y;
00246             pt1.z = z + 1;
00247             unsigned int id1 = getVertexID( pt1.x , pt1.y, pt1.z );
00248             m_idToVertices.insert( ID2WMLPointXYZId::value_type( id1, pt1 ) );
00249 
00250             pt2.x = x;
00251             pt2.y = y + 1;
00252             pt2.z = z + 1;
00253             unsigned int id2 = getVertexID( pt2.x , pt2.y, pt2.z );
00254             m_idToVertices.insert( ID2WMLPointXYZId::value_type( id2, pt2 ) );
00255 
00256             pt3.x = x + 1;
00257             pt3.y = y + 1;
00258             pt3.z = z + 1;
00259             unsigned int id3 = getVertexID( pt3.x , pt3.y, pt3.z );
00260             m_idToVertices.insert( ID2WMLPointXYZId::value_type( id3, pt3 ) );
00261 
00262             pt4.x = x + 1;
00263             pt4.y = y;
00264             pt4.z = z + 1;
00265             unsigned int id4 = getVertexID( pt4.x , pt4.y, pt4.z );
00266             m_idToVertices.insert( ID2WMLPointXYZId::value_type( id4, pt4 ) );
00267 
00268             WMLTriangle triangle1;
00269             triangle1.pointID[0] = id1;
00270             triangle1.pointID[1] = id2;
00271             triangle1.pointID[2] = id3;
00272             WMLTriangle triangle2;
00273             triangle2.pointID[0] = id3;
00274             triangle2.pointID[1] = id4;
00275             triangle2.pointID[2] = id1;
00276             m_trivecTriangles.push_back( triangle1 );
00277             m_trivecTriangles.push_back( triangle2 );
00278             break;
00279         }
00280         default:
00281             break;
00282     }
00283 }
00284 
00285 size_t WMarchingLegoAlgorithm::getVertexID( size_t nX, size_t nY, size_t nZ )
00286 {
00287     return nZ * ( m_nCellsY + 1 ) * ( m_nCellsX + 1) + nY * ( m_nCellsX + 1 ) + nX;
00288 }
00289 
00290 boost::shared_ptr<WTriangleMesh> WMarchingLegoAlgorithm::genSurfaceOneValue( size_t nbCoordsX, size_t nbCoordsY, size_t nbCoordsZ,
00291                                                                                                  const WMatrix< double >& mat,
00292                                                                                                  const std::vector< size_t >* vals,
00293                                                                                                  size_t isoValue,
00294                                                                                                  boost::shared_ptr< WProgressCombiner > mainProgress )
00295 {
00296     WAssert( vals, "No value set provided." );
00297 
00298     m_idToVertices.clear();
00299     m_trivecTriangles.clear();
00300 
00301     m_nCellsX = nbCoordsX - 1;
00302     m_nCellsY = nbCoordsY - 1;
00303     m_nCellsZ = nbCoordsZ - 1;
00304 
00305     m_matrix = mat;
00306 
00307     size_t nX = nbCoordsX;
00308     size_t nY = nbCoordsY;
00309 
00310     size_t nPointsInSlice = nX * nY;
00311 
00312     boost::shared_ptr< WProgress > progress;
00313     if( mainProgress )
00314     {
00315         progress = boost::shared_ptr< WProgress >( new WProgress( "Marching Legos", m_nCellsZ ) );
00316         mainProgress->addSubProgress( progress );
00317     }
00318 
00319     // Generate isosurface.
00320     for( size_t z = 0; z < m_nCellsZ; z++ )
00321     {
00322         if( progress )
00323         {
00324             ++*progress;
00325         }
00326         for( size_t y = 0; y < m_nCellsY; y++ )
00327         {
00328             for( size_t x = 0; x < m_nCellsX; x++ )
00329             {
00330                 if( ( *vals )[ z * nPointsInSlice + y * nX + x ] != isoValue )
00331                 {
00332                     continue;
00333                 }
00334 
00335                 if( x > 0 && ( ( *vals )[ z * nPointsInSlice + y * nX + x - 1 ] != isoValue ) )
00336                 {
00337                     addSurface( x, y, z, 1 );
00338                 }
00339                 if( x < m_nCellsX - 1 && ( ( *vals )[ z * nPointsInSlice + y * nX + x + 1 ] != isoValue ) )
00340                 {
00341                     addSurface( x, y, z, 2 );
00342                 }
00343 
00344                 if( y > 0 && ( ( *vals )[ z * nPointsInSlice + ( y - 1 ) * nX + x ] != isoValue ) )
00345                 {
00346                     addSurface( x, y, z, 3 );
00347                 }
00348 
00349                 if( y < m_nCellsY - 1 && ( ( *vals )[ z * nPointsInSlice + ( y + 1 ) * nX + x ] != isoValue ) )
00350                 {
00351                     addSurface( x, y, z, 4 );
00352                 }
00353 
00354                 if( z > 0 && ( ( *vals )[ ( z - 1 ) * nPointsInSlice + y * nX + x ] != isoValue ) )
00355                 {
00356                     addSurface( x, y, z, 5 );
00357                 }
00358 
00359                 if( z < m_nCellsZ - 1 && ( ( *vals )[ ( z + 1 ) * nPointsInSlice + y * nX + x ] != isoValue ) )
00360                 {
00361                     addSurface( x, y, z, 6 );
00362                 }
00363 
00364                 if( x == 0 )
00365                 {
00366                     addSurface( x, y, z, 1 );
00367                 }
00368                 if( x == m_nCellsX - 1 )
00369                 {
00370                     addSurface( x, y, z, 2 );
00371                 }
00372 
00373                 if( y == 0 )
00374                 {
00375                     addSurface( x, y, z, 3 );
00376                 }
00377 
00378                 if( y == m_nCellsY - 1 )
00379                 {
00380                     addSurface( x, y, z, 4 );
00381                 }
00382 
00383                 if( z == 0 )
00384                 {
00385                     addSurface( x, y, z, 5 );
00386                 }
00387 
00388                 if( z == m_nCellsZ - 1 )
00389                 {
00390                     addSurface( x, y, z, 6 );
00391                 }
00392             }
00393         }
00394     }
00395     unsigned int nextID = 0;
00396     boost::shared_ptr< WTriangleMesh > triMesh( new WTriangleMesh( m_idToVertices.size(), m_trivecTriangles.size() ) );
00397 
00398     // Rename vertices.
00399     ID2WMLPointXYZId::iterator mapIterator = m_idToVertices.begin();
00400     while( mapIterator != m_idToVertices.end() )
00401     {
00402         WPosition texCoord = WPosition( mapIterator->second.x / nbCoordsX,
00403                                                       mapIterator->second.y / nbCoordsY,
00404                                                       mapIterator->second.z / nbCoordsZ );
00405 
00406         // transform from grid coordinate system to world coordinates
00407         WPosition pos = WPosition( mapIterator->second.x, mapIterator->second.y, mapIterator->second.z );
00408 
00409         std::vector< double > resultPos4D( 4 );
00410         resultPos4D[0] = m_matrix( 0, 0 ) * pos[0] + m_matrix( 0, 1 ) * pos[1] + m_matrix( 0, 2 ) * pos[2] + m_matrix( 0, 3 ) * 1;
00411         resultPos4D[1] = m_matrix( 1, 0 ) * pos[0] + m_matrix( 1, 1 ) * pos[1] + m_matrix( 1, 2 ) * pos[2] + m_matrix( 1, 3 ) * 1;
00412         resultPos4D[2] = m_matrix( 2, 0 ) * pos[0] + m_matrix( 2, 1 ) * pos[1] + m_matrix( 2, 2 ) * pos[2] + m_matrix( 2, 3 ) * 1;
00413         resultPos4D[3] = m_matrix( 3, 0 ) * pos[0] + m_matrix( 3, 1 ) * pos[1] + m_matrix( 3, 2 ) * pos[2] + m_matrix( 3, 3 ) * 1;
00414 
00415         ( *mapIterator ).second.newID = nextID;
00416         triMesh->addVertex( resultPos4D[0] / resultPos4D[3],
00417                             resultPos4D[1] / resultPos4D[3],
00418                             resultPos4D[2] / resultPos4D[3] );
00419         triMesh->addTextureCoordinate( texCoord );
00420         nextID++;
00421         mapIterator++;
00422     }
00423 
00424     // Now rename triangles.
00425     WMLTriangleVECTOR::iterator vecIterator = m_trivecTriangles.begin();
00426     while( vecIterator != m_trivecTriangles.end() )
00427     {
00428         for( unsigned int i = 0; i < 3; i++ )
00429         {
00430             unsigned int newID = m_idToVertices[( *vecIterator ).pointID[i]].newID;
00431             ( *vecIterator ).pointID[i] = newID;
00432         }
00433         triMesh->addTriangle( ( *vecIterator ).pointID[0], ( *vecIterator ).pointID[1], ( *vecIterator ).pointID[2] );
00434         vecIterator++;
00435     }
00436 
00437     if( progress )
00438     {
00439         progress->finish();
00440     }
00441     return triMesh;
00442 }