marchingcubes.cpp

ＲＢＦ平台

//_____________________________________________________________________________
// Calculating gradient

float MarchingCubes::get_x_grad( const int i, const int j, const int k ) const
//-----------------------------------------------------------------------------
{
  if( i > 0 )
  {
    if ( i < _size_x - 1 )
      return ( get_data( i+1, j, k ) - get_data( i-1, j, k ) ) / 2 ;
    else
      return get_data( i, j, k ) - get_data( i-1, j, k ) ;
  }
  else
    return get_data( i+1, j, k ) - get_data( i, j, k ) ;
}
//-----------------------------------------------------------------------------

float MarchingCubes::get_y_grad( const int i, const int j, const int k ) const
//-----------------------------------------------------------------------------
{
  if( j > 0 )
  {
    if ( j < _size_y - 1 )
      return ( get_data( i, j+1, k ) - get_data( i, j-1, k ) ) / 2 ;
    else
      return get_data( i, j, k ) - get_data( i, j-1, k ) ;
  }
  else
    return get_data( i, j+1, k ) - get_data( i, j, k ) ;
}
//-----------------------------------------------------------------------------

float MarchingCubes::get_z_grad( const int i, const int j, const int k ) const
//-----------------------------------------------------------------------------
{
  if( k > 0 )
  {
    if ( k < _size_z - 1 )
      return ( get_data( i, j, k+1 ) - get_data( i, j, k-1 ) ) / 2 ;
    else
      return get_data( i, j, k ) - get_data( i, j, k-1 ) ;
  }
  else
    return get_data( i, j, k+1 ) - get_data( i, j, k ) ;
}
//_____________________________________________________________________________


//_____________________________________________________________________________
// Adding vertices

void MarchingCubes::test_vertex_addition()
{
  if( _nverts >= _Nverts )
  {
    Vertex_mc *temp = _vertices ;
    _vertices = new Vertex_mc[ _Nverts*2 ] ;
    memcpy( _vertices, temp, _Nverts*sizeof(Vertex_mc) ) ;
    delete[] temp ;
    printf("%d allocated vertices\n", _Nverts) ;
    _Nverts *= 2 ;
  }
}


int MarchingCubes::add_x_vertex( )
//-----------------------------------------------------------------------------
{
  test_vertex_addition() ;
  Vertex_mc *vert = _vertices + _nverts++ ;

  float   u = ( _cube[0] ) / ( _cube[0] - _cube[1] ) ;

  vert->x      = (float)_i+u;
  vert->y      = (float) _j ;
  vert->z      = (float) _k ;

  vert->nx = (1-u)*get_x_grad(_i,_j,_k) + u*get_x_grad(_i+1,_j,_k) ;
  vert->ny = (1-u)*get_y_grad(_i,_j,_k) + u*get_y_grad(_i+1,_j,_k) ;
  vert->nz = (1-u)*get_z_grad(_i,_j,_k) + u*get_z_grad(_i+1,_j,_k) ;

  u = (float) sqrt( vert->nx * vert->nx + vert->ny * vert->ny +vert->nz * vert->nz ) ;
  if( u > 0 )
  {
    vert->nx /= u ;
    vert->ny /= u ;
    vert->nz /= u ;
  }

  return _nverts-1 ;
}
//-----------------------------------------------------------------------------

int MarchingCubes::add_y_vertex( )
//-----------------------------------------------------------------------------
{
  test_vertex_addition() ;
  Vertex_mc *vert = _vertices + _nverts++ ;

  float   u = ( _cube[0] ) / ( _cube[0] - _cube[3] ) ;

  vert->x      = (float) _i ;
  vert->y      = (float)_j+u;
  vert->z      = (float) _k ;

  vert->nx = (1-u)*get_x_grad(_i,_j,_k) + u*get_x_grad(_i,_j+1,_k) ;
  vert->ny = (1-u)*get_y_grad(_i,_j,_k) + u*get_y_grad(_i,_j+1,_k) ;
  vert->nz = (1-u)*get_z_grad(_i,_j,_k) + u*get_z_grad(_i,_j+1,_k) ;

  u = (float) sqrt( vert->nx * vert->nx + vert->ny * vert->ny +vert->nz * vert->nz ) ;
  if( u > 0 )
  {
    vert->nx /= u ;
    vert->ny /= u ;
    vert->nz /= u ;
  }

  return _nverts-1 ;
}
//-----------------------------------------------------------------------------

int MarchingCubes::add_z_vertex( )
//-----------------------------------------------------------------------------
{
  test_vertex_addition() ;
  Vertex_mc *vert = _vertices + _nverts++ ;

  float   u = ( _cube[0] ) / ( _cube[0] - _cube[4] ) ;

  vert->x      = (float) _i ;
  vert->y      = (float) _j ;
  vert->z      = (float)_k+u;

  vert->nx = (1-u)*get_x_grad(_i,_j,_k) + u*get_x_grad(_i,_j,_k+1) ;
  vert->ny = (1-u)*get_y_grad(_i,_j,_k) + u*get_y_grad(_i,_j,_k+1) ;
  vert->nz = (1-u)*get_z_grad(_i,_j,_k) + u*get_z_grad(_i,_j,_k+1) ;

  u = (float) sqrt( vert->nx * vert->nx + vert->ny * vert->ny +vert->nz * vert->nz ) ;
  if( u > 0 )
  {
    vert->nx /= u ;
    vert->ny /= u ;
    vert->nz /= u ;
  }

  return _nverts-1 ;
}


int MarchingCubes::add_c_vertex( )
//-----------------------------------------------------------------------------
{
  test_vertex_addition() ;
  Vertex_mc *vert = _vertices + _nverts++ ;

  float  u = 0 ;
  int   vid ;

  vert->x = vert->y = vert->z =  vert->nx = vert->ny = vert->nz = 0 ;

  // Computes the average of the intersection points of the cube
  vid = get_x_vert( _i , _j , _k ) ;
  if( vid != -1 ) { ++u ; const Vertex_mc &v = _vertices[vid] ; vert->x += v.x ;  vert->y += v.y ;  vert->z += v.z ;  vert->nx += v.nx ; vert->ny += v.ny ; vert->nz += v.nz ; }
  vid = get_y_vert(_i+1, _j , _k ) ;
  if( vid != -1 ) { ++u ; const Vertex_mc &v = _vertices[vid] ; vert->x += v.x ;  vert->y += v.y ;  vert->z += v.z ;  vert->nx += v.nx ; vert->ny += v.ny ; vert->nz += v.nz ; }
  vid = get_x_vert( _i ,_j+1, _k ) ;
  if( vid != -1 ) { ++u ; const Vertex_mc &v = _vertices[vid] ; vert->x += v.x ;  vert->y += v.y ;  vert->z += v.z ;  vert->nx += v.nx ; vert->ny += v.ny ; vert->nz += v.nz ; }
  vid = get_y_vert( _i , _j , _k ) ;
  if( vid != -1 ) { ++u ; const Vertex_mc &v = _vertices[vid] ; vert->x += v.x ;  vert->y += v.y ;  vert->z += v.z ;  vert->nx += v.nx ; vert->ny += v.ny ; vert->nz += v.nz ; }
  vid = get_x_vert( _i , _j ,_k+1) ;
  if( vid != -1 ) { ++u ; const Vertex_mc &v = _vertices[vid] ; vert->x += v.x ;  vert->y += v.y ;  vert->z += v.z ;  vert->nx += v.nx ; vert->ny += v.ny ; vert->nz += v.nz ; }
  vid = get_y_vert(_i+1, _j ,_k+1) ;
  if( vid != -1 ) { ++u ; const Vertex_mc &v = _vertices[vid] ; vert->x += v.x ;  vert->y += v.y ;  vert->z += v.z ;  vert->nx += v.nx ; vert->ny += v.ny ; vert->nz += v.nz ; }
  vid = get_x_vert( _i ,_j+1,_k+1) ;
  if( vid != -1 ) { ++u ; const Vertex_mc &v = _vertices[vid] ; vert->x += v.x ;  vert->y += v.y ;  vert->z += v.z ;  vert->nx += v.nx ; vert->ny += v.ny ; vert->nz += v.nz ; }
  vid = get_y_vert( _i , _j ,_k+1) ;
  if( vid != -1 ) { ++u ; const Vertex_mc &v = _vertices[vid] ; vert->x += v.x ;  vert->y += v.y ;  vert->z += v.z ;  vert->nx += v.nx ; vert->ny += v.ny ; vert->nz += v.nz ; }
  vid = get_z_vert( _i , _j , _k ) ;
  if( vid != -1 ) { ++u ; const Vertex_mc &v = _vertices[vid] ; vert->x += v.x ;  vert->y += v.y ;  vert->z += v.z ;  vert->nx += v.nx ; vert->ny += v.ny ; vert->nz += v.nz ; }
  vid = get_z_vert(_i+1, _j , _k ) ;
  if( vid != -1 ) { ++u ; const Vertex_mc &v = _vertices[vid] ; vert->x += v.x ;  vert->y += v.y ;  vert->z += v.z ;  vert->nx += v.nx ; vert->ny += v.ny ; vert->nz += v.nz ; }
  vid = get_z_vert(_i+1,_j+1, _k ) ;
  if( vid != -1 ) { ++u ; const Vertex_mc &v = _vertices[vid] ; vert->x += v.x ;  vert->y += v.y ;  vert->z += v.z ;  vert->nx += v.nx ; vert->ny += v.ny ; vert->nz += v.nz ; }
  vid = get_z_vert( _i ,_j+1, _k ) ;
  if( vid != -1 ) { ++u ; const Vertex_mc &v = _vertices[vid] ; vert->x += v.x ;  vert->y += v.y ;  vert->z += v.z ;  vert->nx += v.nx ; vert->ny += v.ny ; vert->nz += v.nz ; }

  vert->x  /= u ;
  vert->y  /= u ;
  vert->z  /= u ;

  u = (float) sqrt( vert->nx * vert->nx + vert->ny * vert->ny +vert->nz * vert->nz ) ;
  if( u > 0 )
  {
    vert->nx /= u ;
    vert->ny /= u ;
    vert->nz /= u ;
  }

  return _nverts-1 ;
}
//_____________________________________________________________________________



//_____________________________________________________________________________
//_____________________________________________________________________________




//_____________________________________________________________________________
// Grid exportation
void MarchingCubes::writeISO(const char *fn )
//-----------------------------------------------------------------------------
{
  unsigned char buf[sizeof(float)] ;

  FILE *fp = fopen( fn, "wb" ) ;

  // header
  * (int*) buf = _size_x ;
  fwrite(buf, sizeof(float), 1, fp);
  * (int*) buf = _size_y ;
  fwrite(buf, sizeof(float), 1, fp);
  * (int*) buf = _size_z ;
  fwrite(buf, sizeof(float), 1, fp);

  for( int i = 0 ; i < _size_x ; i++ )
  {
    for( int j = 0 ; j < _size_y ; j++ )
    {
      for( int k = 0 ; k < _size_z ; k++ )
      {
        * (float*) buf = get_data( i,j,k ) ;
        fwrite(buf, sizeof(float), 1, fp);
      }
    }
  }

  fclose(fp) ;
}
//_____________________________________________________________________________



//_____________________________________________________________________________
// PLY exportation
typedef struct PlyFace
{
  unsigned char nverts;    /* number of Vertex indices in list */
  int *verts;              /* Vertex index list */
} PlyFace;

PlyProperty vert_props[]  = { /* list of property information for a PlyVertex */
  {"x", Float32, Float32, offsetof( Vertex_mc,x ), 0, 0, 0, 0},
  {"y", Float32, Float32, offsetof( Vertex_mc,y ), 0, 0, 0, 0},
  {"z", Float32, Float32, offsetof( Vertex_mc,z ), 0, 0, 0, 0},
  {"nx", Float32, Float32, offsetof( Vertex_mc,nx ), 0, 0, 0, 0},
  {"ny", Float32, Float32, offsetof( Vertex_mc,ny ), 0, 0, 0, 0},
  {"nz", Float32, Float32, offsetof( Vertex_mc,nz ), 0, 0, 0, 0}
};

PlyProperty face_props[]  = { /* list of property information for a PlyFace */
  {"vertex_indices", Int32, Int32, offsetof( PlyFace,verts ),
    1, Uint8, Uint8, offsetof( PlyFace,nverts )},
};


void MarchingCubes::writePLY(const char *fn, bool bin )
//-----------------------------------------------------------------------------
{
  Ply_tk plytk;

  PlyFile    *ply;
  FILE       *fp = fopen( fn, "w" );

  int          i ;
  PlyFace     face ;
  int         verts[3] ;
  char       *elem_names[]  = { "vertex", "face" };
  printf("Marching Cubes::writePLY(%s)...", fn ) ;
  ply = plytk.write_ply ( fp, 2, elem_names, bin? PLY_BINARY_LE : PLY_ASCII );

  /* describe what properties go into the PlyVertex elements */
  plytk.describe_element_ply ( ply, "vertex", _nverts );
  plytk.describe_property_ply ( ply, &vert_props[0] );
  plytk.describe_property_ply ( ply, &vert_props[1] );
  plytk.describe_property_ply ( ply, &vert_props[2] );
  plytk.describe_property_ply ( ply, &vert_props[3] );
  plytk.describe_property_ply ( ply, &vert_props[4] );
  plytk.describe_property_ply ( ply, &vert_props[5] );

  /* describe PlyFace properties (just list of PlyVertex indices) */
  plytk.describe_element_ply ( ply, "face", _ntrigs );
  plytk.describe_property_ply ( ply, &face_props[0] );

  plytk.header_complete_ply ( ply );

  /* set up and write the PlyVertex elements */
  plytk.put_element_setup_ply ( ply, "vertex" );
  for ( i = 0; i < _nverts; i++ )
    plytk.put_element_ply ( ply, ( void * ) &(_vertices[i]) );
  printf("   %d vertices written\n", _nverts ) ;

  /* set up and write the PlyFace elements */
  plytk.put_element_setup_ply ( ply, "face" );
  face.nverts = 3 ;
  face.verts  = verts ;
  for ( i = 0; i < _ntrigs; i++ )
  {
    face.verts[0] = _triangles[i].v1 ;
    face.verts[1] = _triangles[i].v2 ;
    face.verts[2] = _triangles[i].v3 ;
    plytk.put_element_ply ( ply, ( void * ) &face );
  }
  printf("   %d triangles written\n", _ntrigs ) ;

  plytk.close_ply ( ply );
  plytk.free_ply ( ply );
  fclose( fp ) ;
}
//_____________________________________________________________________________

//_____________________________________________________________________________
// Open Inventor / VRML 1.0 ascii exportation
void MarchingCubes::writeIV(const char *fn )
//-----------------------------------------------------------------------------
{
  FILE *fp = fopen( fn, "w" ) ;
  int   i ;

  printf("Marching Cubes::exportIV(%s)...", fn) ;

  fprintf( fp, "#Inventor V2.1 ascii \n\nSeparator { \n    ShapeHints {\n        vertexOrdering  COUNTERCLOCKWISE\n        shapeType       UNKNOWN_SHAPE_TYPE\n        creaseAngle     0.0\n    }\n Coordinate3 { \n point [  \n" ) ;
  for ( i = 0; i < _nverts; i++ )
    fprintf( fp, " %f %f %f,\n", _vertices[i].x, _vertices[i].y, _vertices[i].z ) ;
  printf("   %d vertices written\n", _nverts ) ;

  fprintf( fp, "\n ] \n} \nNormal { \nvector [ \n" ) ;
  for ( i = 0; i < _nverts; i++ )
    fprintf( fp, " %f %f %f,\n", _vertices[i].nx, _vertices[i].ny, _vertices[i].nz ) ;

  fprintf( fp, "\n ] \n} \nIndexedFaceSet { \ncoordIndex [ \n" ) ;
  for ( i = 0; i < _ntrigs; i++ )
    fprintf( fp, "%d, %d, %d, -1,\n", _triangles[i].v1, _triangles[i].v2, _triangles[i].v3 ) ;

  fprintf( fp, " ] \n } \n } \n" ) ;
  fclose( fp ) ;
  printf("   %d triangles written\n", _ntrigs ) ;
}
//_____________________________________________________________________________

void MarchingCubes::convert2trimesh(TriMesh * trimesh)
{
	if(!(trimesh->vertices.empty()))
		trimesh->vertices.clear();
	if(!(trimesh->faces.empty()))
		trimesh->faces.clear();
	if(!(trimesh->normals.empty()))
		trimesh->normals.clear();
	if(!(trimesh->tstrips.empty()))
		trimesh->tstrips.clear();
		
	int nf = ntrigs();
	int nv = nverts();
	Vertex_mc * varray = vertices();
	Triangle * triarray = triangles();
	trimesh->vertices.resize(nv);

	for(int i = 0; i < nf; ++i)
	{
		Triangle tridx = triarray[i];
		TriMesh::Face fsidx(tridx.v1,tridx.v2,tridx.v3);
		trimesh->faces.push_back(fsidx);

		//vertex1
		Vertex_mc vtxidx = varray[tridx.v1];
		point vtx1(vtxidx.x, vtxidx.y, vtxidx.z);
		trimesh->vertices[tridx.v1] = vtx1;
			
		//vertex2
		vtxidx = varray[tridx.v2];
		point vtx2(vtxidx.x, vtxidx.y, vtxidx.z);
		trimesh->vertices[tridx.v2] = vtx2;
		
		//vertex3
		vtxidx = varray[tridx.v3];
		point vtx3(vtxidx.x, vtxidx.y, vtxidx.z);
		trimesh->vertices[tridx.v3] = vtx3;
	}

	clean_temps();
	clean_all();
}