marchingcubes.h

Go to the documentation of this file.

/*
 * Author: Tao Ju, 5/16/2007 <taoju@cs.wustl.edu>, code ported by Grant Tang
 * Copyright (c) 2000-2006 Baylor College of Medicine
 *
 * This software is issued under a joint BSD/GNU license. You may use the
 * source code in this file under either license. However, note that the
 * complete EMAN2 and SPARX software packages have some GPL dependencies,
 * so you are responsible for compliance with the licenses of these packages
 * if you opt to use BSD licensing. The warranty disclaimer below holds
 * in either instance.
 *
 * This complete copyright notice must be included in any revised version of the
 * source code. Additional authorship citations may be added, but existing
 * author citations must be preserved.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program 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 General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 *
 * */
 
#ifndef _MARCHING_CUBES_H_
#define _MARCHING_CUBES_H_
 
#include <vector>
using std::vector;
 
#include "vecmath.h"
#include "isosurface.h"
 
// Marching cubes debug will turn on debug and timing information
#define MARCHING_CUBES_DEBUG 0
 
#include <ostream>
using std::ostream;
 
#include <climits>
// for CHAR_BIT
 
#ifdef __APPLE__
        #include "OpenGL/gl.h"
        #include "OpenGL/glext.h"
#else // WIN32, LINUX
        #include "GL/gl.h"
        #include "GL/glext.h"
#endif  //__APPLE__
 
namespace EMAN
{
        template<typename type>
        class CustomVector
        {
                public:
                        explicit CustomVector(unsigned int starting_size=1024) : data(0), size(0), elements(0)
                        {
                                resize(starting_size);
                        }
 
                        ~CustomVector()
                        {
                                if(data) {free(data); data=0;}
                        }
 
                        inline void clear(unsigned int starting_size=1024)
                        {
                                if (data) {free(data); data = 0;}
                                size = 0;
                                elements = 0;
                                resize(starting_size);
                        }
 
                        inline void resize(const unsigned int n)
                        {
                                data = (type*)realloc(data, n*sizeof(type));
 
                                for(unsigned int i = size; i < n; ++i) data[i] = 0;
                                size = n;
                        }
 
                        inline void push_back(const type& t)
                        {
                                if ( elements == size ) resize(2*size);
                                data[elements] = t;
                                ++elements;
                        }
 
                        inline void push_back_3(const type* const p)
                        {
                                if ( elements+2 >= size ) resize(2*size);
                                memcpy( &data[elements], p, 3*sizeof(type));
                                elements = elements + 3;
                        }
 
                        inline void mult3(const type& v1, const type& v2, const type& v3)
                        {
                                for(unsigned int i = 0; (i + 2) < elements; i += 3 ){
                                        data[i] *= v1;
                                        data[i+1] *= v2;
                                        data[i+2] *= v3;
                                }
                        }
 
                        inline unsigned int elem() { return elements; }
 
                        inline type& operator[](const unsigned int idx)
                        {
                                unsigned int csize = size;
                                while (idx >= csize ) {
                                        csize *= 2;
                                }
                                if ( csize != size ) resize(csize);
                                return data[idx];
                        }
 
                        inline type* get_data() { return data; }
 
                private:
                        type* data;
                        unsigned int size;
                        unsigned int elements;
        };
 
        class ColorRGBGenerator{
        public:
                ColorRGBGenerator();
                ColorRGBGenerator(EMData* emdata);
                
                float* getRGBColor(int x, int y, int z);
                
                void set_data(EMData* data);
                
                void generateRadialColorMap();
                
                void set_cmap_data(EMData* data);
                
                inline void setOrigin(int orix, int oriy, int oriz)
                {
                        originx = orix;
                        originy = oriy;
                        originz = oriz;
                        needtorecolor = true;
                }
                inline void setScale(float i, float o)
                {
                        inner = i;
                        outer = o;
                        needtorecolor = true;
                }
                inline void setRGBmode(int mode)
                {
                        rgbmode = mode;
                        needtorecolor = true;
                }
                
                inline void setMinMax(float min, float max)
                {
                        minimum = min;
                        maximum = max;
                        needtorecolor = true;
                }
                
                inline int getRGBmode()
                {
                        return rgbmode;
                }
                
                inline bool getNeedToRecolor()
                {
                        return needtorecolor;
                }
                
                inline void setNeedToRecolor(bool recolor)
                {
                        needtorecolor = recolor;
                }
        private:
                int rgbmode;
                int originx;
                int originy;
                int originz;
                float inner;
                float outer;
                float minimum;
                float maximum;
        
                bool needtorecolor;     // dirty bit to let the system know when we need to recolor
                
                float* colormap;        // pointer to a colormap
                EMData* em_data;        // pointer to EMdata
                EMData* cmap;           // pointer to colormap data
                
                float rgb[3];           //little array to hold RGB values;
        };
        
        class MarchingCubes : public Isosurface {
                friend class GLUtil;
        public:
                MarchingCubes();
 
                MarchingCubes(EMData * em);
                virtual ~MarchingCubes();
 
                void set_data(EMData* data);
 
                void set_surface_value(const float value);
 
                float get_surface_value() const { return _surf_value; }
 
                void set_sampling(const int rate) { drawing_level = rate; }
 
                int get_sampling() const { return drawing_level; }
 
                int get_sampling_range() { return minvals.size()-1; }
 
                void color_vertices();
                
                Dict get_isosurface();
 
                void surface_face_z();
                
                inline void setRGBorigin(int x, int y, int z)
                {
                        rgbgenerator.setOrigin(x, y, z);
                }
                
                inline void setRGBscale(float i, float o)
                {
                        rgbgenerator.setScale(i, o);
                }
                
                inline void setRGBmode(int mode)
                {
                        rgbgenerator.setRGBmode(mode);
                }
                
                inline int getRGBmode()
                {
                        return rgbgenerator.getRGBmode();
                }
                
                inline void setCmapData(EMData* data)
                {
                        rgbgenerator.set_cmap_data(data);
                }
                
                inline void setCmapMinMax(float min, float max)
                {
                        rgbgenerator.setMinMax(min, max);
                }
                
        private:
                map<int, int> point_map;
                unsigned long _isodl;
                GLuint buffer[4];
 
                bool calculate_min_max_vals();
                
        
                void clear_min_max_vals();
 
                vector<EMData*> minvals, maxvals;
 
                int drawing_level;
 
                void draw_cube(const int x, const int y, const int z, const int cur_level );
 
 
                void marching_cube(int fX, int fY, int fZ, const int cur_level);
 
                void calculate_surface();
                
                float get_offset(float fValue1, float fValue2, float fValueDesired);
 
                int get_edge_num(int x, int y, int z, int edge);
 
                void get_normal(Vector3 &normal, int fX, int fY, int fZ);
 
                CustomVector<float> pp;
                CustomVector<float> cc;
                CustomVector<int> vv;
                CustomVector<float> nn;
                CustomVector<unsigned int> ff;
                
                ColorRGBGenerator rgbgenerator;
                
                bool needtobind;        // A dirty bit to signal when the the MC algorithm or color has chaged and hence a need to update GPU buffers
        };
 
//      class CudaMarchingCubes : public Isosurface {
//              public:
//
//                      /** Most commonly used constructor
//                      * calls set_data(em)
//                      * @param em the EMData object to generate triangles and normals for
//                      */
//                      CudaMarchingCubes(EMData * em);
//                      virtual ~CudaMarchingCubes();
//                      /**
//                       * Set Isosurface value
//                       */
//                      virtual void set_surface_value(const float value);
//
//                      /** Sets Voxel data for Isosurface implementation
//                      * Calls calculate_min_max_vals which generates the tree of data
//                      * @param data the emdata object to be rendered in 3D
//                      * @exception ImageDimensionException if the image z dimension is 1
//                      */
//                      void set_data(EMData* data);
//
//                      virtual float get_surface_value() const { return _surf_value; }
//
//                      /**
//                       * Set Grid Size
//                       */
//                      virtual void set_sampling(const int size) {} ;
//
//                      virtual int get_sampling() const { return 0; };
//
//                      /** Get the number of feasible samplings
//                      *
//                       */
//                      virtual int get_sampling_range() {return 0; }
//
//                      virtual Dict get_isosurface()  { return Dict;}
//
//      }
        
        class U3DWriter{
        public:
                typedef unsigned int U32;
                typedef long unsigned int U64;
                typedef double F64;
                typedef float F32;
                typedef short int I16;
                typedef short unsigned int U16;
                typedef unsigned char U8;
 
 
 
                U3DWriter();
                ~U3DWriter();
 
                int write(const string& filename);
 
                ostream& write(ostream&);
 
        private:
                unsigned int size_of_in_bytes();
 
                void test_type_sizes();
 
                ostream& write_header(ostream&);
                ostream& write_clod_mesh_generator_node(ostream& os);
 
 
                template<typename type>
                ostream& write(ostream& os, const type& T) {
                        os.write( (const char*)(&T), sizeof(type) );
                        return os;
                }
 
 
//              U8 reverse_bits_3d( U8 val ) {
//                      U8 ret = 0;
//                      U8 n_bits = sizeof ( val ) * CHAR_BIT;
//
//                      for ( unsigned i = 0; i < n_bits; ++i ) {
//                              ret = ( ret << 1 ) | ( val & 1 );
//                              val >>= 1;
//                      }
//
//                      return ret;
//              }
 
//              U16 reverse_bits_u3d( U16 val ) {
//                      bool is_big_endian = ByteOrder::is_host_big_endian();
//
//                      U16 ret = 0;
//
//                      U8* pu8 = (U8*)&val;
//                      U8* pu8ret = (U8*)&ret;
//                      U8 high_order = 0;
//                      U8 low_order = 0;
//                      if (is_big_endian) {
//                              high_order = pu8[0];
//                              low_order = pu8[1];
//                      }
//                      else {
//                              high_order = pu8[1];
//                              low_order = pu8[0];
//                      }
//
//                      pu8ret[0] = low_order;
//                      pu8ret[1] = higher_order;
//
//                      high_order = reverse_bits_u3d(high_order);
//                      low_order = reverse_bits_u3d(low_order);
//
//                      return ret;
//
//              }
//
//              U32 reverse_bits_u3d( U32 val ) {
//                      bool is_big_endian = ByteOrder::is_host_big_endian();
//
//                      U32 ret = 0;
//
//                      U16 * pu16 = (U16*)&val;
//                      U16* pu16ret = (U16*)&ret;
//                      U16 high_order = 0;
//                      U16 low_order = 0;
//                      if (is_big_endian) {
//                              high_order = pu16[0];
//                              low_order = pu16[1];
//                      }
//                      else {
//                              high_order = pu16[1];
//                              low_order = pu16[0];
//                      }
//
//                      high_order = reverse_bits_u3d(high_order);
//                      low_order = reverse_bits_u3d(low_order);
//
//                      pu16ret[0] = low_order;
//                      pu16ret[1] = higher_order;
//
//                      return ret;
//
//              }
//
 
                U32 DIFFUSE_COLOR_COUNT;
                U32 SPECULAR_COLOR_COUNT;
 
 
                CustomVector<F32> pp;
                CustomVector<F32> nn;
                CustomVector<unsigned int> ff;
        };
 
        // Template specialization. Have to be careful when dealing with strings
        template<> ostream& U3DWriter::write(ostream& os, const string& );
 
//      class U3DBitStreamWriter {
//      public:
//              typedef unsigned int U32;
//              typedef long unsigned int U64;
//              typedef double F64;
//              typedef float F32;
//              typedef short int I16;
//              typedef short unsigned int U16;
//              typedef unsigned char U8;
//
//
//              U3D_BitStreamWriter();
//              ~U3D_BitStreamWriter();
//
//              // We could do this kind of thing with templates,  but
//              // we would have to specialize each function anyhow,
//              // so it makes sense just to do it the C way.
//              void writeU8( U8 value );
//              void writeU16( U16 value );
//              void writeU32( U32 value );
//              void writeU64( U64 value );
//              void writeI32( I32 value );
//              void writeF32( F32 value );
//              void writeF64( F64 value );
//
//              void WriteCompressedU32( U32 context, U32 value );
//              void WriteCompressedU16( U32 context, U32 value );
//              void WriteCompressedU8( U32 context, U32 value );
//
//      };
 
 
 
}
 
#endif