MarchingCubeUtility.cpp

Go to the documentation of this file.
00001 /******************************************************************************
00002 *       SOFA, Simulation Open-Framework Architecture, version 1.0 beta 3      *
00003 *                (c) 2006-2008 MGH, INRIA, USTL, UJF, CNRS                    *
00004 *                                                                             *
00005 * This library is free software; you can redistribute it and/or modify it     *
00006 * under the terms of the GNU Lesser General Public License as published by    *
00007 * the Free Software Foundation; either version 2.1 of the License, or (at     *
00008 * your option) any later version.                                             *
00009 *                                                                             *
00010 * This library is distributed in the hope that it will be useful, but WITHOUT *
00011 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or       *
00012 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License *
00013 * for more details.                                                           *
00014 *                                                                             *
00015 * You should have received a copy of the GNU Lesser General Public License    *
00016 * along with this library; if not, write to the Free Software Foundation,     *
00017 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301 USA.          *
00018 *******************************************************************************
00019 *                              SOFA :: Framework                              *
00020 *                                                                             *
00021 * Authors: M. Adam, J. Allard, B. Andre, P-J. Bensoussan, S. Cotin, C. Duriez,*
00022 * H. Delingette, F. Falipou, F. Faure, S. Fonteneau, L. Heigeas, C. Mendoza,  *
00023 * M. Nesme, P. Neumann, J-P. de la Plata Alcade, F. Poyer and F. Roy          *
00024 *                                                                             *
00025 * Contact information: contact@sofa-framework.org                             *
00026 ******************************************************************************/
00027 /**** Geometry Shader Marching Cubes
00028     * Copyright Cyril Crassin, tt 2007.
00029     * This code is partially based on the example of 
00030     * Paul Bourke "Polygonising a scalar field" located at :
00031     * http://local.wasp.uwa.edu.au/~pbourke/geometry/polygonise/
00032 ****/
00033 
00034 #include <sofa/helper/MarchingCubeUtility.h>
00035         
00036 #define PRECISION 16384.0
00037 
00038 #include <string.h>
00039 
00040 namespace sofa
00041 {
00042 
00043 namespace helper
00044 {
00045 
00046     const int MarchingCubeEdgeTable[256] = {
00047     0x0  , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
00048     0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
00049     0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
00050     0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
00051     0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c,
00052     0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
00053     0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac,
00054     0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
00055     0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c,
00056     0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
00057     0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc,
00058     0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
00059     0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c,
00060     0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
00061     0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc ,
00062     0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
00063     0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc,
00064     0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
00065     0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c,
00066     0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
00067     0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc,
00068     0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
00069     0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c,
00070     0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460,
00071     0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
00072     0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0,
00073     0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c,
00074     0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230,
00075     0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,
00076     0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190,
00077     0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
00078     0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0   };
00079 
00080     const int MarchingCubeTriTable[256][16] =
00081     {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00082     {0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00083     {0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00084     {1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00085     {1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00086     {0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00087     {9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00088     {2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1},
00089     {3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00090     {0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00091     {1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00092     {1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1},
00093     {3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00094     {0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1},
00095     {3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1},
00096     {9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00097     {4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00098     {4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00099     {0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00100     {4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1},
00101     {1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00102     {3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1},
00103     {9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
00104     {2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1},
00105     {8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00106     {11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1},
00107     {9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
00108     {4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1},
00109     {3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1},
00110     {1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1},
00111     {4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1},
00112     {4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1},
00113     {9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00114     {9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00115     {0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00116     {8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1},
00117     {1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00118     {3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
00119     {5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1},
00120     {2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1},
00121     {9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00122     {0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
00123     {0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
00124     {2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1},
00125     {10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1},
00126     {4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1},
00127     {5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1},
00128     {5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1},
00129     {9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00130     {9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1},
00131     {0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1},
00132     {1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00133     {9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1},
00134     {10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1},
00135     {8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1},
00136     {2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1},
00137     {7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1},
00138     {9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1},
00139     {2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1},
00140     {11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1},
00141     {9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1},
00142     {5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1},
00143     {11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1},
00144     {11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00145     {10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00146     {0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00147     {9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00148     {1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
00149     {1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00150     {1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1},
00151     {9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1},
00152     {5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1},
00153     {2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00154     {11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
00155     {0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
00156     {5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1},
00157     {6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1},
00158     {0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1},
00159     {3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1},
00160     {6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1},
00161     {5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00162     {4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1},
00163     {1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
00164     {10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1},
00165     {6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1},
00166     {1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1},
00167     {8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1},
00168     {7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1},
00169     {3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
00170     {5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1},
00171     {0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1},
00172     {9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1},
00173     {8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1},
00174     {5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1},
00175     {0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1},
00176     {6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1},
00177     {10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00178     {4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1},
00179     {10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1},
00180     {8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1},
00181     {1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1},
00182     {3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1},
00183     {0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00184     {8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1},
00185     {10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1},
00186     {0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1},
00187     {3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1},
00188     {6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1},
00189     {9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1},
00190     {8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1},
00191     {3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1},
00192     {6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00193     {7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1},
00194     {0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1},
00195     {10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1},
00196     {10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1},
00197     {1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1},
00198     {2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1},
00199     {7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1},
00200     {7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00201     {2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1},
00202     {2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1},
00203     {1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1},
00204     {11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1},
00205     {8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1},
00206     {0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00207     {7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1},
00208     {7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00209     {7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00210     {3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00211     {0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00212     {8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
00213     {10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00214     {1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
00215     {2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
00216     {6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1},
00217     {7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00218     {7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1},
00219     {2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1},
00220     {1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1},
00221     {10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1},
00222     {10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1},
00223     {0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1},
00224     {7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1},
00225     {6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00226     {3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1},
00227     {8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1},
00228     {9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1},
00229     {6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1},
00230     {1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1},
00231     {4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1},
00232     {10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1},
00233     {8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1},
00234     {0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00235     {1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1},
00236     {1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1},
00237     {8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1},
00238     {10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1},
00239     {4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1},
00240     {10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00241     {4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00242     {0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
00243     {5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
00244     {11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1},
00245     {9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
00246     {6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1},
00247     {7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1},
00248     {3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1},
00249     {7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1},
00250     {9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1},
00251     {3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1},
00252     {6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1},
00253     {9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1},
00254     {1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1},
00255     {4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1},
00256     {7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1},
00257     {6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1},
00258     {3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1},
00259     {0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1},
00260     {6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1},
00261     {1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1},
00262     {0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1},
00263     {11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1},
00264     {6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1},
00265     {5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1},
00266     {9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1},
00267     {1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1},
00268     {1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00269     {1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1},
00270     {10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1},
00271     {0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00272     {10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00273     {11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00274     {11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1},
00275     {5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1},
00276     {10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1},
00277     {11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1},
00278     {0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1},
00279     {9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1},
00280     {7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1},
00281     {2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1},
00282     {8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1},
00283     {9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1},
00284     {9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1},
00285     {1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00286     {0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1},
00287     {9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1},
00288     {9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00289     {5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1},
00290     {5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1},
00291     {0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1},
00292     {10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1},
00293     {2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1},
00294     {0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1},
00295     {0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1},
00296     {9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00297     {2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1},
00298     {5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1},
00299     {3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1},
00300     {5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1},
00301     {8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1},
00302     {0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00303     {8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1},
00304     {9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00305     {4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1},
00306     {0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1},
00307     {1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1},
00308     {3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1},
00309     {4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1},
00310     {9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1},
00311     {11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1},
00312     {11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1},
00313     {2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1},
00314     {9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1},
00315     {3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1},
00316     {1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00317     {4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1},
00318     {4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1},
00319     {4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00320     {4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00321     {9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00322     {3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1},
00323     {0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1},
00324     {3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00325     {1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1},
00326     {3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1},
00327     {0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00328     {3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00329     {2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1},
00330     {9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00331     {2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1},
00332     {1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00333     {1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00334     {0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00335     {0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00336     {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}};
00337 
00338 
00339     MarchingCubeUtility::MarchingCubeUtility()
00340     : cubeStep(1), convolutionSize(1), 
00341     dataResolution(0,0,0), dataVoxelSize(1.0f,1.0f,1.0f)
00342     {}
00343 
00344     /*
00345     Linearly interpolate the position where an isosurface cuts
00346     an edge between two vertices, each with their own scalar value
00347     */
00348     void MarchingCubeUtility::vertexInterp(const float isolevel, 
00349                                            const Vector3 &p1, const Vector3 &p2, 
00350                                            const float valp1, const float valp2, Vector3 &p) const
00351     {
00352         float mu = (isolevel - valp1) / (valp2 - valp1);
00353         p = p1 + (p2 - p1) * mu;
00354         p = p.linearProduct(dataVoxelSize*dataResolution[0]*0.5f);
00355         p[0] = (int)(p[0]*PRECISION)/PRECISION;
00356         p[1] = (int)(p[1]*PRECISION)/PRECISION;
00357         p[2] = (int)(p[2]*PRECISION)/PRECISION;
00358     }
00359 
00360 
00361     bool MarchingCubeUtility::testGrid(const float v, const float isolevel) const
00362     {
00363         return (v<isolevel);
00364     }
00365 
00366 
00367     /*
00368     Given a grid cell and an isolevel, calculate the triangular
00369     facets required to represent the isosurface through the cell.
00370     Return the number of triangular facets, the array "triangles"
00371     will be loaded up with the vertices at most 5 triangular facets.
00372     0 will be returned if the grid cell is either totally above
00373     of totally below the isolevel.
00374     */
00375     int MarchingCubeUtility::polygonise(const GridCell &grid, float isolevel, sofa::helper::vector< PointID > &triangles,
00376                                         std::map< Vector3, PointID> &map_vertices,
00377                                         std::map< PointID, Vector3> &map_indices, 
00378                                         unsigned int &ID) const
00379     {
00380 
00381         int i,ntriang;
00382         int cubeindex;
00383         Vector3 vertindex[12];
00384         /*
00385         Determine the index into the edge table which
00386         tells us which vertices are inside of the surface
00387         */
00388         cubeindex = 0;
00389         if (testGrid(grid.val[0], isolevel)) cubeindex |= 1;
00390         if (testGrid(grid.val[1], isolevel)) cubeindex |= 2;
00391         if (testGrid(grid.val[2], isolevel)) cubeindex |= 4;
00392         if (testGrid(grid.val[3], isolevel)) cubeindex |= 8;
00393         if (testGrid(grid.val[4], isolevel)) cubeindex |= 16;
00394         if (testGrid(grid.val[5], isolevel)) cubeindex |= 32;
00395         if (testGrid(grid.val[6], isolevel)) cubeindex |= 64;
00396         if (testGrid(grid.val[7], isolevel)) cubeindex |= 128;
00397 
00398         /* Cube is entirely in/out of the surface */
00399         if (MarchingCubeEdgeTable[cubeindex] == 0) return(0);
00400 
00401         /* Find the vertices where the surface intersects the cube */
00402         if (MarchingCubeEdgeTable[cubeindex] & 1)
00403             vertexInterp(isolevel,grid.pos[0],grid.pos[1],grid.val[0],grid.val[1], vertindex[0]);
00404         if (MarchingCubeEdgeTable[cubeindex] & 2)                            
00405             vertexInterp(isolevel,grid.pos[1],grid.pos[2],grid.val[1],grid.val[2], vertindex[1]);
00406         if (MarchingCubeEdgeTable[cubeindex] & 4)                            
00407             vertexInterp(isolevel,grid.pos[2],grid.pos[3],grid.val[2],grid.val[3], vertindex[2]);
00408         if (MarchingCubeEdgeTable[cubeindex] & 8)                            
00409             vertexInterp(isolevel,grid.pos[3],grid.pos[0],grid.val[3],grid.val[0], vertindex[3]);
00410         if (MarchingCubeEdgeTable[cubeindex] & 16)                           
00411             vertexInterp(isolevel,grid.pos[4],grid.pos[5],grid.val[4],grid.val[5], vertindex[4]);
00412         if (MarchingCubeEdgeTable[cubeindex] & 32)                           
00413             vertexInterp(isolevel,grid.pos[5],grid.pos[6],grid.val[5],grid.val[6], vertindex[5]);
00414         if (MarchingCubeEdgeTable[cubeindex] & 64)                           
00415             vertexInterp(isolevel,grid.pos[6],grid.pos[7],grid.val[6],grid.val[7], vertindex[6]);
00416         if (MarchingCubeEdgeTable[cubeindex] & 128)                          
00417             vertexInterp(isolevel,grid.pos[7],grid.pos[4],grid.val[7],grid.val[4], vertindex[7]);
00418         if (MarchingCubeEdgeTable[cubeindex] & 256)                          
00419             vertexInterp(isolevel,grid.pos[0],grid.pos[4],grid.val[0],grid.val[4], vertindex[8]);
00420         if (MarchingCubeEdgeTable[cubeindex] & 512)                          
00421             vertexInterp(isolevel,grid.pos[1],grid.pos[5],grid.val[1],grid.val[5], vertindex[9]);
00422         if (MarchingCubeEdgeTable[cubeindex] & 1024)                             
00423             vertexInterp(isolevel,grid.pos[2],grid.pos[6],grid.val[2],grid.val[6], vertindex[10]);
00424         if (MarchingCubeEdgeTable[cubeindex] & 2048)                             
00425             vertexInterp(isolevel,grid.pos[3],grid.pos[7],grid.val[3],grid.val[7], vertindex[11]);
00426 
00427         /* Create the triangle */
00428         ntriang = 0;
00429         std::map< Vector3, PointID>::iterator iter;
00430         Vector3 current_P;
00431         PointID current_ID;
00432         for (i=0;MarchingCubeTriTable[cubeindex][i]!=-1;i+=3) 
00433         {
00434             Vec<3,unsigned int> current_triangle;
00435             for (PointID j=0;j<3;++j)
00436             { 
00437                 current_P = vertindex[MarchingCubeTriTable[cubeindex][i+j]];
00438                 //Search if the current Vertex P is already stored with an ID
00439                 iter = map_vertices.find(current_P);
00440                 if (iter != map_vertices.end()) current_ID = iter->second;
00441                 else
00442                 {
00443                     //Add new Vertex in map 
00444                     current_ID = ++ID;
00445                     map_vertices.insert(std::make_pair(current_P, current_ID));   
00446                     map_indices.insert (std::make_pair(current_ID,current_P));
00447                 }
00448                 current_triangle[j]=current_ID;
00449             }
00450 
00451             if (current_triangle[0] == current_triangle[1] || current_triangle[0] == current_triangle[2] || current_triangle[2] == current_triangle[1]) continue;
00452             triangles.push_back(current_triangle[0]);
00453             triangles.push_back(current_triangle[1]);
00454             triangles.push_back(current_triangle[2]);
00455 
00456             ntriang+=3;
00457         }
00458         return(ntriang);
00459     }
00460 
00461 
00462 
00463     void MarchingCubeUtility::run( const unsigned char *_data, const float isolevel,  
00464                                   sofa::helper::vector< PointID > &mesh,
00465                                   sofa::helper::vector< Vector3 >  &vertices) const 
00466     {
00467         std::map<PointID, Vector3>  map_indices;
00468 
00469         vector< float > data(dataResolution[0]*dataResolution[1]*dataResolution[2]);
00470         if (data.size() == 0) 
00471             return;
00472 
00473         for (unsigned int i=0; i<data.size(); ++i) 
00474             data[i] = (float) getVoxel(i,_data);           
00475 
00476         if (convolutionSize != 0) 
00477             smoothData(&data[0]);
00478 
00479         unsigned int ID = 0;
00480         std::map< Vector3, PointID> map_vertices;
00481 
00482         Vec3i gridSize = Vec3i(dataResolution[0]/cubeStep, dataResolution[1]/cubeStep, dataResolution[2]/cubeStep);
00483 
00484         Vector3 gridStep = Vector3(2.0f/((float)gridSize[0]), 2.0f/((float)gridSize[1]), 2.0f/((float)gridSize[2]));
00485 
00486         Vec3i dataGridStep(dataResolution[0]/gridSize[0],dataResolution[1]/gridSize[1],dataResolution[2]/gridSize[2]);
00487 
00488         PointID counter_triangles=0;
00489         for(int k=0; k<gridSize[2]-1; k++)
00490             for(int j=0; j<gridSize[1]-1; j++)
00491                 for(int i=0; i<gridSize[0]-1; i++)
00492                 {
00493                     GridCell cell;
00494                     Vector3 vcurf((float)i, (float)j, (float)k);
00495                     Vec3i vcuri(i, j, k);
00496 
00497                     cell.pos[0]=vcurf.linearProduct(gridStep)-Vector3(1.0f,1.0f,1.0f);
00498                     Vec3i valPos0=vcuri.linearProduct(dataGridStep);
00499                     cell.val[0]=data[valPos0[0] + valPos0[1]*dataResolution[0] + valPos0[2]*dataResolution[0]*dataResolution[1]];
00500 
00501                     Vec3i valPos;
00502 
00503                     cell.pos[1]=cell.pos[0]+Vector3(gridStep[0], 0, 0);
00504                     if(i==gridSize[0]-1)
00505                         valPos=valPos0;
00506                     else
00507                         valPos=valPos0+Vec3i(dataGridStep[0], 0, 0);
00508                     cell.val[1]=data[valPos[0] + valPos[1]*dataResolution[0] + valPos[2]*dataResolution[0]*dataResolution[1]];
00509 
00510                     cell.pos[2]=cell.pos[0]+Vector3(gridStep[0], gridStep[1], 0);
00511                     valPos=valPos0+Vec3i(i==gridSize[0]-1 ? 0 : dataGridStep[0], j==gridSize[1]-1 ? 0 : dataGridStep[1], 0);
00512                     cell.val[2]=data[valPos[0] + valPos[1]*dataResolution[0] + valPos[2]*dataResolution[0]*dataResolution[1]];
00513 
00514                     cell.pos[3]=cell.pos[0]+Vector3(0, gridStep[1], 0);
00515                     valPos=valPos0+Vec3i(0, j==gridSize[1]-1 ? 0 : dataGridStep[1], 0);
00516                     cell.val[3]=data[valPos[0] + valPos[1]*dataResolution[0] + valPos[2]*dataResolution[0]*dataResolution[1]];
00517 
00518                     cell.pos[4]=cell.pos[0]+Vector3(0, 0, gridStep[2]);
00519                     valPos=valPos0+Vec3i(0, 0, k==gridSize[2]-1 ? 0 : dataGridStep[2]);
00520                     cell.val[4]=data[valPos[0] + valPos[1]*dataResolution[0] + valPos[2]*dataResolution[0]*dataResolution[1]];
00521 
00522                     cell.pos[5]=cell.pos[0]+Vector3(gridStep[0], 0, gridStep[2]);
00523                     valPos=valPos0+Vec3i(i==gridSize[0]-1 ? 0 : dataGridStep[0], 0, k==gridSize[2]-1 ? 0 : dataGridStep[2]);
00524                     cell.val[5]=data[valPos[0] + valPos[1]*dataResolution[0] + valPos[2]*dataResolution[0]*dataResolution[1]];
00525 
00526                     cell.pos[6]=cell.pos[0]+Vector3(gridStep[0], gridStep[1], gridStep[2]);
00527                     valPos=valPos0+Vec3i(i==gridSize[0]-1 ? 0 : dataGridStep[0], j==gridSize[1]-1 ? 0 : dataGridStep[1], k==gridSize[2]-1 ? 0 : dataGridStep[2]);
00528                     cell.val[6]=data[valPos[0] + valPos[1]*dataResolution[0] + valPos[2]*dataResolution[0]*dataResolution[1]];
00529 
00530                     cell.pos[7]=cell.pos[0]+Vector3(0, gridStep[1], gridStep[2]);
00531                     valPos=valPos0+Vec3i(0, j==gridSize[1]-1 ? 0 : dataGridStep[1], k==gridSize[2]-1 ? 0 : dataGridStep[2]);
00532                     cell.val[7]=data[valPos[0] + valPos[1]*dataResolution[0] + valPos[2]*dataResolution[0]*dataResolution[1]];
00533 
00534 
00535                     int numvert = polygonise(cell, isolevel, mesh, map_vertices, map_indices, ID);
00536                     counter_triangles += numvert/3;
00537                 }
00538 
00539         //WARNING mesh from Marching Cube has indices starting from ID 1, a sofa mesh begins with ID 0
00540         vertices.resize(map_indices.size());
00541         for(unsigned int i=0; i<map_indices.size(); ++i)
00542         {
00543             vertices[i] = map_indices[i+1]; 
00544         }
00545 
00546         for(unsigned int i=0; i<mesh.size(); ++i)
00547         {
00548             mesh[i] -= 1;
00549         }
00550     }
00551 
00552     void MarchingCubeUtility::run( const unsigned char *data,  const float isolevel, 
00553                                                 sofa::helper::io::Mesh &m) const
00554     {           
00555         using sofa::helper::vector;
00556         using sofa::defaulttype::Vector3;
00557 
00558         std::cout << "Creating Mesh using Marching Cubes\n";        
00559         vector<Vector3> &vertices                 = m.getVertices();
00560         vector< vector < vector <int> > > &facets = m.getFacets();
00561 
00562         vector< PointID >               triangles;
00563 
00564         //Do the Marching Cube
00565         run(data, isolevel, triangles, vertices);
00566 
00567         const unsigned int numTriangles = triangles.size()/3;
00568         facets.resize(numTriangles, vector< vector < int > >(3, vector<int>(3, 0)));
00569         for (unsigned int i=0; i<triangles.size(); /*i+=3*/)
00570         {
00571             vector< vector< int > > &vertNormTexIndices = facets[i/3];
00572             vector<int> &vIndices = vertNormTexIndices[0];
00573 
00574             vIndices[0] = triangles[i++]; 
00575             vIndices[1] = triangles[i++];
00576             vIndices[2] = triangles[i++];
00577         }   
00578     }
00579 
00580     void MarchingCubeUtility::smoothData( float *data) const
00581     {
00582         std::cout << "Smoothing Data using " << convolutionSize << "x"<< convolutionSize << "x"<< convolutionSize << " as gaussian convolution kernel\n";
00583         vector< float >convolutionKernel;
00584         createConvolutionKernel(convolutionKernel);
00585 
00586         float *original_data=new float[dataResolution[0]*dataResolution[1]*dataResolution[2]];
00587         memcpy(original_data, data, sizeof(float)*dataResolution[0]*dataResolution[1]*dataResolution[2]);
00588         unsigned int limit = convolutionSize/2;  
00589         unsigned int z,y,x;
00590         for (z=limit;z<dataResolution[2]-limit;++z)
00591         {
00592             for (y=limit;y<dataResolution[1]-limit;++y)
00593             {
00594                 for (x=limit;x<dataResolution[0]-limit;++x)
00595                 {
00596                     applyConvolution(x,y,z, original_data, data, convolutionKernel);
00597                 }
00598             }
00599         }
00600 
00601         delete [] original_data;
00602     }
00603 
00604     void  MarchingCubeUtility::applyConvolution(unsigned int x, unsigned int y, unsigned int z, 
00605                                                 const float *original_data, 
00606                                                 float *data, 
00607                                                 const vector< float >  &convolutionKernel) const
00608     {
00609 
00610         const unsigned int index = dataResolution[0]*(z*dataResolution[1] +y) +x;
00611         const unsigned int step=dataResolution[0]*dataResolution[1];
00612         data[index] = 0;
00613         const unsigned int c=convolutionSize/2;
00614         unsigned int _z,_y,_x;
00615         unsigned int i=0;
00616         for (  _z=0;_z<convolutionSize;++_z)
00617         {
00618             for (  _y=0;_y<convolutionSize;++_y)
00619             {
00620                 for (  _x=0;_x<convolutionSize;++_x)
00621                 {
00622                     data[index] += convolutionKernel[i++]*original_data[index+step*(_z-c)+dataResolution[0]*(_y-c)+(_x-c)];
00623                 }
00624             }
00625         }
00626     }
00627 
00628 
00629     void MarchingCubeUtility::createConvolutionKernel(vector< float >  &convolutionKernel) const
00630     {
00631         int c = (1 + convolutionSize/2);
00632         float total = 0.0;
00633 
00634         convolutionKernel.resize(convolutionSize*convolutionSize*convolutionSize);
00635         const float var = (float)convolutionSize/*/2.0f*/;
00636         unsigned int i=0;
00637         for (unsigned int z=0;z<convolutionSize;++z)
00638         {
00639             for (unsigned int y=0;y<convolutionSize;++y)
00640             {
00641                 for (unsigned int x=0;x<convolutionSize;++x)
00642                 {           
00643                     convolutionKernel[i] = (float)(exp( -(pow((float)(x+1-c),2) + pow((float)(y+1-c),2) + pow((float)(z+1-c),2))/(2.0f*var)));  
00644                     total += convolutionKernel[i++];
00645                 }
00646             }
00647         }
00648 
00649         i=0;
00650         for (unsigned int z=0;z<convolutionSize;++z)      
00651             for (unsigned int y=0;y<convolutionSize;++y)        
00652                 for (unsigned int x=0;x<convolutionSize;++x)          
00653                     convolutionKernel[i++] /= total;
00654 
00655 
00656     }
00657 
00658 
00659 }  
00660 
00661 }