EMAN2
marchingcubes.h
Go to the documentation of this file.
00001 
00005 /*
00006  * Author: Tao Ju, 5/16/2007 <taoju@cs.wustl.edu>, code ported by Grant Tang
00007  * Copyright (c) 2000-2006 Baylor College of Medicine
00008  *
00009  * This software is issued under a joint BSD/GNU license. You may use the
00010  * source code in this file under either license. However, note that the
00011  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00012  * so you are responsible for compliance with the licenses of these packages
00013  * if you opt to use BSD licensing. The warranty disclaimer below holds
00014  * in either instance.
00015  *
00016  * This complete copyright notice must be included in any revised version of the
00017  * source code. Additional authorship citations may be added, but existing
00018  * author citations must be preserved.
00019  *
00020  * This program is free software; you can redistribute it and/or modify
00021  * it under the terms of the GNU General Public License as published by
00022  * the Free Software Foundation; either version 2 of the License, or
00023  * (at your option) any later version.
00024  *
00025  * This program is distributed in the hope that it will be useful,
00026  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00027  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00028  * GNU General Public License for more details.
00029  *
00030  * You should have received a copy of the GNU General Public License
00031  * along with this program; if not, write to the Free Software
00032  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00033  *
00034  * */
00035 
00036 #ifndef _MARCHING_CUBES_H_
00037 #define _MARCHING_CUBES_H_
00038 
00039 #include <vector>
00040 using std::vector;
00041 
00042 #include "vecmath.h"
00043 #include "isosurface.h"
00044 
00045 // Marching cubes debug will turn on debug and timing information
00046 #define MARCHING_CUBES_DEBUG 0
00047 
00048 #include <ostream>
00049 using std::ostream;
00050 
00051 #include <climits>
00052 // for CHAR_BIT
00053 
00054 #ifdef __APPLE__
00055         #include "OpenGL/gl.h"
00056         #include "OpenGL/glext.h"
00057 #else // WIN32, LINUX
00058         #include "GL/gl.h"
00059         #include "GL/glext.h"
00060 #endif  //__APPLE__
00061 
00062 namespace EMAN
00063 {
00077         template<typename type>
00078         class CustomVector
00079         {
00080                 public:
00084                         explicit CustomVector(unsigned int starting_size=1024) : data(0), size(0), elements(0)
00085                         {
00086                                 resize(starting_size);
00087                         }
00088 
00089                         ~CustomVector()
00090                         {
00091                                 if(data) {free(data); data=0;}
00092                         }
00093 
00098                         inline void clear(unsigned int starting_size=1024)
00099                         {
00100                                 if (data) {free(data); data = 0;}
00101                                 size = 0;
00102                                 elements = 0;
00103                                 resize(starting_size);
00104                         }
00105 
00112                         inline void resize(const unsigned int n)
00113                         {
00114                                 data = (type*)realloc(data, n*sizeof(type));
00115 
00116                                 for(unsigned int i = size; i < n; ++i) data[i] = 0;
00117                                 size = n;
00118                         }
00119 
00124                         inline void push_back(const type& t)
00125                         {
00126                                 if ( elements == size ) resize(2*size);
00127                                 data[elements] = t;
00128                                 ++elements;
00129                         }
00130 
00135                         inline void push_back_3(const type* const p)
00136                         {
00137                                 if ( elements+2 >= size ) resize(2*size);
00138                                 memcpy( &data[elements], p, 3*sizeof(type));
00139                                 elements = elements + 3;
00140                         }
00141 
00144                         inline void mult3(const type& v1, const type& v2, const type& v3)
00145                         {
00146                                 for(unsigned int i = 0; (i + 2) < elements; i += 3 ){
00147                                         data[i] *= v1;
00148                                         data[i+1] *= v2;
00149                                         data[i+2] *= v3;
00150                                 }
00151                         }
00152 
00159                         inline unsigned int elem() { return elements; }
00160 
00165                         inline type& operator[](const unsigned int idx)
00166                         {
00167                                 unsigned int csize = size;
00168                                 while (idx >= csize ) {
00169                                         csize *= 2;
00170                                 }
00171                                 if ( csize != size ) resize(csize);
00172                                 return data[idx];
00173                         }
00174 
00178                         inline type* get_data() { return data; }
00179 
00180                 private:
00182                         type* data;
00184                         unsigned int size;
00186                         unsigned int elements;
00187         };
00188 
00192         class ColorRGBGenerator{
00193         public:
00194                 ColorRGBGenerator();
00196                 ColorRGBGenerator(EMData* emdata);
00197                 
00199                 float* getRGBColor(int x, int y, int z);
00200                 
00202                 void set_data(EMData* data);
00203                 
00205                 void generateRadialColorMap();
00206                 
00208                 void set_cmap_data(EMData* data);
00209                 
00211                 inline void setOrigin(int orix, int oriy, int oriz)
00212                 {
00213                         originx = orix;
00214                         originy = oriy;
00215                         originz = oriz;
00216                         needtorecolor = true;
00217                 }
00219                 inline void setScale(float i, float o)
00220                 {
00221                         inner = i;
00222                         outer = o;
00223                         needtorecolor = true;
00224                 }
00226                 inline void setRGBmode(int mode)
00227                 {
00228                         rgbmode = mode;
00229                         needtorecolor = true;
00230                 }
00231                 
00233                 inline void setMinMax(float min, float max)
00234                 {
00235                         minimum = min;
00236                         maximum = max;
00237                         needtorecolor = true;
00238                 }
00239                 
00241                 inline int getRGBmode()
00242                 {
00243                         return rgbmode;
00244                 }
00245                 
00247                 inline bool getNeedToRecolor()
00248                 {
00249                         return needtorecolor;
00250                 }
00251                 
00252                 inline void setNeedToRecolor(bool recolor)
00253                 {
00254                         needtorecolor = recolor;
00255                 }
00256         private:
00257                 int rgbmode;
00258                 int originx;
00259                 int originy;
00260                 int originz;
00261                 float inner;
00262                 float outer;
00263                 float minimum;
00264                 float maximum;
00265         
00266                 bool needtorecolor;     // dirty bit to let the sytem know when we need to recolor
00267                 
00268                 float* colormap;        // pointer to a colormap
00269                 EMData* em_data;        // pointer to EMdata
00270                 EMData* cmap;           // pointer to colormap data
00271                 
00272                 float rgb[3];           //little array to hold RGB values;
00273         };
00274         
00275         class MarchingCubes : public Isosurface {
00276                 friend class GLUtil;
00277         public:
00280                 MarchingCubes();
00281 
00286                 MarchingCubes(EMData * em);
00287                 virtual ~MarchingCubes();
00288 
00294                 void set_data(EMData* data);
00295 
00299                 void set_surface_value(const float value);
00300 
00304                 float get_surface_value() const { return _surf_value; }
00305 
00314                 void set_sampling(const int rate) { drawing_level = rate; }
00315 
00319                 int get_sampling() const { return drawing_level; }
00320 
00323                 int get_sampling_range() { return minvals.size()-1; }
00324 
00327                 void color_vertices();
00328                 
00333                 Dict get_isosurface();
00334 
00335                 void surface_face_z();
00336                 
00339                 inline void setRGBorigin(int x, int y, int z)
00340                 {
00341                         rgbgenerator.setOrigin(x, y, z);
00342                 }
00343                 
00344                 inline void setRGBscale(float i, float o)
00345                 {
00346                         rgbgenerator.setScale(i, o);
00347                 }
00348                 
00349                 inline void setRGBmode(int mode)
00350                 {
00351                         rgbgenerator.setRGBmode(mode);
00352                 }
00353                 
00355                 inline int getRGBmode()
00356                 {
00357                         return rgbgenerator.getRGBmode();
00358                 }
00359                 
00361                 inline void setCmapData(EMData* data)
00362                 {
00363                         rgbgenerator.set_cmap_data(data);
00364                 }
00365                 
00367                 inline void setCmapMinMax(float min, float max)
00368                 {
00369                         rgbgenerator.setMinMax(min, max);
00370                 }
00371                 
00372         private:
00373                 map<int, int> point_map;
00374                 unsigned long _isodl;
00375                 GLuint buffer[4];
00376 
00382                 bool calculate_min_max_vals();
00383                 
00384         
00388                 void clear_min_max_vals();
00389 
00391                 vector<EMData*> minvals, maxvals;
00392 
00394                 int drawing_level;
00395 
00405                 void draw_cube(const int x, const int y, const int z, const int cur_level );
00406 
00407 
00415                 void marching_cube(int fX, int fY, int fZ, const int cur_level);
00416 
00420                 void calculate_surface();
00421                 
00430                 float get_offset(float fValue1, float fValue2, float fValueDesired);
00431 
00435                 int get_edge_num(int x, int y, int z, int edge);
00436 
00447                 void get_normal(Vector3 &normal, int fX, int fY, int fZ);
00448 
00450                 CustomVector<float> pp;
00451                 CustomVector<float> cc;
00452                 CustomVector<int> vv;
00453                 CustomVector<float> nn;
00454                 CustomVector<unsigned int> ff;
00455                 
00457                 ColorRGBGenerator rgbgenerator;
00458                 
00459                 bool needtobind;        // A dirty bit to signal when the the MC algorithm or color has chaged and hence a need to update GPU buffers
00460         };
00461 
00465 //      class CudaMarchingCubes : public Isosurface {
00466 //              public:
00467 //
00468 //                      /** Most commonly used constructor
00469 //                      * calls set_data(em)
00470 //                      * @param em the EMData object to generate triangles and normals for
00471 //                      */
00472 //                      CudaMarchingCubes(EMData * em);
00473 //                      virtual ~CudaMarchingCubes();
00474 //                      /**
00475 //                       * Set Isosurface value
00476 //                       */
00477 //                      virtual void set_surface_value(const float value);
00478 //
00479 //                      /** Sets Voxel data for Isosurface implementation
00480 //                      * Calls calculate_min_max_vals which generates the tree of data
00481 //                      * @param data the emdata object to be rendered in 3D
00482 //                      * @exception ImageDimensionException if the image z dimension is 1
00483 //                      */
00484 //                      void set_data(EMData* data);
00485 //
00486 //                      virtual float get_surface_value() const { return _surf_value; }
00487 //
00488 //                      /**
00489 //                       * Set Grid Size
00490 //                       */
00491 //                      virtual void set_sampling(const int size) {} ;
00492 //
00493 //                      virtual int get_sampling() const { return 0; };
00494 //
00495 //                      /** Get the number of feasible samplings
00496 //                      *
00497 //                       */
00498 //                      virtual int get_sampling_range() {return 0; }
00499 //
00500 //                      virtual Dict get_isosurface()  { return Dict;}
00501 //
00502 //      }
00503         
00506         class U3DWriter{
00507         public:
00508                 typedef unsigned int U32;
00509                 typedef long unsigned int U64;
00510                 typedef double F64;
00511                 typedef float F32;
00512                 typedef short int I16;
00513                 typedef short unsigned int U16;
00514                 typedef unsigned char U8;
00515 
00516 
00517 
00518                 U3DWriter();
00519                 ~U3DWriter();
00520 
00521                 int write(const string& filename);
00522 
00523                 ostream& write(ostream&);
00524 
00525         private:
00526                 unsigned int size_of_in_bytes();
00527 
00528                 void test_type_sizes();
00529 
00530                 ostream& write_header(ostream&);
00531                 ostream& write_clod_mesh_generator_node(ostream& os);
00532 
00533 
00534                 template<typename type>
00535                 ostream& write(ostream& os, const type& T) {
00536                         os.write( (const char*)(&T), sizeof(type) );
00537                         return os;
00538                 }
00539 
00540 
00541 //              U8 reverse_bits_3d( U8 val ) {
00542 //                      U8 ret = 0;
00543 //                      U8 n_bits = sizeof ( val ) * CHAR_BIT;
00544 //
00545 //                      for ( unsigned i = 0; i < n_bits; ++i ) {
00546 //                              ret = ( ret << 1 ) | ( val & 1 );
00547 //                              val >>= 1;
00548 //                      }
00549 //
00550 //                      return ret;
00551 //              }
00552 
00553 //              U16 reverse_bits_u3d( U16 val ) {
00554 //                      bool is_big_endian = ByteOrder::is_host_big_endian();
00555 //
00556 //                      U16 ret = 0;
00557 //
00558 //                      U8* pu8 = (U8*)&val;
00559 //                      U8* pu8ret = (U8*)&ret;
00560 //                      U8 high_order = 0;
00561 //                      U8 low_order = 0;
00562 //                      if (is_big_endian) {
00563 //                              high_order = pu8[0];
00564 //                              low_order = pu8[1];
00565 //                      }
00566 //                      else {
00567 //                              high_order = pu8[1];
00568 //                              low_order = pu8[0];
00569 //                      }
00570 //
00571 //                      pu8ret[0] = low_order;
00572 //                      pu8ret[1] = higher_order;
00573 //
00574 //                      high_order = reverse_bits_u3d(high_order);
00575 //                      low_order = reverse_bits_u3d(low_order);
00576 //
00577 //                      return ret;
00578 //
00579 //              }
00580 //
00581 //              U32 reverse_bits_u3d( U32 val ) {
00582 //                      bool is_big_endian = ByteOrder::is_host_big_endian();
00583 //
00584 //                      U32 ret = 0;
00585 //
00586 //                      U16 * pu16 = (U16*)&val;
00587 //                      U16* pu16ret = (U16*)&ret;
00588 //                      U16 high_order = 0;
00589 //                      U16 low_order = 0;
00590 //                      if (is_big_endian) {
00591 //                              high_order = pu16[0];
00592 //                              low_order = pu16[1];
00593 //                      }
00594 //                      else {
00595 //                              high_order = pu16[1];
00596 //                              low_order = pu16[0];
00597 //                      }
00598 //
00599 //                      high_order = reverse_bits_u3d(high_order);
00600 //                      low_order = reverse_bits_u3d(low_order);
00601 //
00602 //                      pu16ret[0] = low_order;
00603 //                      pu16ret[1] = higher_order;
00604 //
00605 //                      return ret;
00606 //
00607 //              }
00608 //
00609 
00610                 U32 DIFFUSE_COLOR_COUNT;
00611                 U32 SPECULAR_COLOR_COUNT;
00612 
00613 
00614                 CustomVector<F32> pp;
00615                 CustomVector<F32> nn;
00616                 CustomVector<unsigned int> ff;
00617         };
00618 
00619         // Template specialization. Have to be careful when dealing with strings
00620         template<> ostream& U3DWriter::write(ostream& os, const string& );
00621 
00622 //      class U3DBitStreamWriter {
00623 //      public:
00624 //              typedef unsigned int U32;
00625 //              typedef long unsigned int U64;
00626 //              typedef double F64;
00627 //              typedef float F32;
00628 //              typedef short int I16;
00629 //              typedef short unsigned int U16;
00630 //              typedef unsigned char U8;
00631 //
00632 //
00633 //              U3D_BitStreamWriter();
00634 //              ~U3D_BitStreamWriter();
00635 //
00636 //              // We could do this kind of thing with templates,  but
00637 //              // we would have to specialize each function anyhow,
00638 //              // so it makes sense just to do it the C way.
00639 //              void writeU8( U8 value );
00640 //              void writeU16( U16 value );
00641 //              void writeU32( U32 value );
00642 //              void writeU64( U64 value );
00643 //              void writeI32( I32 value );
00644 //              void writeF32( F32 value );
00645 //              void writeF64( F64 value );
00646 //
00647 //              void WriteCompressedU32( U32 context, U32 value );
00648 //              void WriteCompressedU16( U32 context, U32 value );
00649 //              void WriteCompressedU8( U32 context, U32 value );
00650 //
00651 //      };
00652 
00653 
00654 
00655 }
00656 
00657 #endif