Biomechanical Joint Model  Author: Anderson Maciel

---

testMetaball_inventor.cpp

Go to the documentation of this file.

#include <util.h>
#undef vector
#undef min
#undef max


#include <Inventor/nodes/SoMaterial.h>
#include <Inventor/nodes/SoSeparator.h>
#include <Inventor/Xt/SoXt.h>
#include <Inventor/nodes/SoLineSet.h>
#include <Inventor/Xt/viewers/SoXtExaminerViewer.h>
#include <Inventor/nodes/SoPointSet.h>
#include <Inventor/nodes/SoVertexProperty.h>
#include <Inventor/nodes/SoDrawStyle.h>
#include <Inventor/nodes/SoSphere.h>
#include <Inventor/actions/SoWriteAction.h>
#include <Inventor/nodes/SoIndexedFaceSet.h>
#include <Inventor/nodes/SoCoordinate3.h>
#include <Inventor/actions/SoSearchAction.h>
#include <Inventor/nodes/SoNormal.h>
#include <Inventor/nodes/SoIndexedLineSet.h>
#include <Inventor/nodes/SoEventCallback.h>
#include <Inventor/events/SoEvents.h>
#include <Inventor/events/SoMouseButtonEvent.h>

#include <LinAlg/linalg.h>
#include <MarchingCubes/marching_cubes.h>
#include <MarchingCubes/blob.h>
#include <MarchingCubes/mesh.h>

#include <vector>
#include <list>

using namespace LinAlg;
using namespace MarchingCubes;
using namespace std;
using namespace VL;


// global vars.
float skin_magic_sphere = 1.0;  // sphere threshold redius
float skin_magic_blob = 0.5;

const int MAX_STRING_SIZE = 256; // maximal lenght of the strings

void PrintHelp(char *Command){
  printf("Usage: %s -r res -m spheres -k pot. field param -out mesh \n",Command);
  printf("-r   : Mesh resouluton[0.1 high resolution mesh generated, 2.0 very corse mesh generated]\n");
  printf("-m   : File with sphere parameters[cx cy cz r] -> center and radious of the sphere.\n"\
         "     : At the begining of the file give the number of the spheres\n");
  printf("-k   : Parameter of the potential filed[small values(<1.0) -> implicit surface is big and mesh\n"\
         "       doesn't follow the form of the spheres, big values(>1.0) -> surface fallows the form of the spheres\n");
  printf("-a   : Parameter of the exponential potential filed[small values(<1.0) -> implicit surface is big and mesh\n"\
         "       doesn't follow the form of the spheres, big values(>1.0) -> surface fallows the form of the spheres\n");
  printf("out  : Output triangulation base name[generates name.pts(vertices), name.pts.tri(triangles), name.nor(normals)\n"); 

}

/**************************************************************************************************
  Mouse click callback function  

***************************************************************************************************/

void myMouseEvent(void *userData, SoEventCallback *event){

  
  const SoEvent *myEvent = event->getEvent();

  if(SO_MOUSE_PRESS_EVENT(myEvent,BUTTON1))
    cout<<"Left mouse button pressed"<<endl;
  if(SO_MOUSE_PRESS_EVENT(myEvent,BUTTON2))
    cout<<"Middle mouse button pressed"<<endl; 
  if(SO_MOUSE_PRESS_EVENT(myEvent,BUTTON3))
    cout<<"Right mouse button pressed"<<endl; 
}


/***************************************************************************************************
  Read ellipsoids from file which format is:
 
  Number of ellipsoids
  N
  Center  Radii     Principal axes(they with the center define transf. matrix from g.c.s -> l.c.s
  x y z   lx ly lz  (i1.x,i1.y,i1.z) (j1.x,j1.y,j1.z) (k1.x, k1.y, k1.z)

****************************************************************************************************/
void readPrimitives(char* fileInp, Metaballs &mb){
  FILE *pFile;
  VEC C(3), i1(3), j1(3), k1(3);
  double lx, ly, lz;
  HMAT M;
  int number;
  
  printf("Read metaballs\n");
  pFile = fopen ( fileInp , "r" );
  if(pFile==NULL){
    fprintf(stderr,"Error in reading the file %s\n",fileInp);
    exit (1);
  }
  
  // obtain the number of implicit surfaces
  fscanf(pFile,"%d",&number);
  printf("Numb. of blobs: %d\n",number);
   
  for (int i=0;i<number;i++){
    fscanf(pFile,"%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",\
            &C[0],&C[1],&C[2],&lx,&ly,&lz,&i1[0],&i1[1],&i1[2],&j1[0],&j1[1],&j1[2],&k1[0],&k1[1],&k1[2]);


    M(0,0)=i1[0]; M(0,1)=j1[0]; M(0,2)=k1[0]; M(0,3)= C[0];
    M(1,0)=i1[1]; M(1,1)=j1[1]; M(1,2)=k1[1]; M(1,3)= C[1]; 
    M(2,0)=i1[2]; M(2,1)=j1[2]; M(2,2)=k1[2]; M(2,3)= C[2];
 
    cout<<"M\n"<<M<<endl;
    mb.addBlob(Blob(M, lx, ly, lz));
  }
  getchar();
  fclose(pFile);
}

// Read the metaballs(Raquel's predoc project) from file
void readSpheres(char* fileInp,  MetaballsSphere &mb, double k)
{
  // file with the number of implicit surfaces
  // for each implicit surface (x y z e k)-> params of the field function for the spherical metaballs

  FILE *pFile;
  //TsiSphere *vector;
  VEC C(3);
  double ro;
  int number;

  printf("Read spheres\n");
  printf("Reading the file %s\n",fileInp);
  pFile = fopen ( fileInp , "r" );
  if (pFile==NULL){
      printf("Error in reading the file\n");
      exit (1);
  }

  // obtain the number of implicit surfaces
  fscanf(pFile,"%d",&number);
  //printf("NUMBER %d\n",*number);

  //printf("number %d\n",*number);
  //vector=(TsiSphere *) malloc((*number)*sizeof(TsiSphere));
  for (int i=0;i<number;i++){
    fscanf(pFile,"%lf %lf %lf %lf\n",&C[0],&C[1],&C[2],&ro);
    printf("%lf %lf %lf %lf %lf\n",C[0],C[1],C[2],ro,k);
    mb.addSphere(Sphere(C, ro, k));
  }
  fclose(pFile);
  printf("Read spheres finished\n");
}

/**************************************************************************************************
  Save triangulated isosurface of the metaball field function int the Trig like format.
  This means make three files: *.pts with vertex coordinates
                               *.tri with facets
                               *.nx, *.ny, *.nz - binary files with vertex normals
**************************************************************************************************/

void SaveIsoSurf2Trig(const marchingCubes &surf, char *BaseName){
  vector<VEC> vertlist = surf.getVertices();
  vector<VEC> normlist = surf.getNormals();
  vector<Triangle> trilist = surf.getTriangles();
  char pts_fname[MAX_STRING_SIZE], tri_fname[MAX_STRING_SIZE], nor_fname[MAX_STRING_SIZE];
  FILE *fpts, *ftri, *fnor;


  // specific output file name creation

  strcpy(pts_fname,BaseName); strcat(pts_fname,".pts");
  strcpy(tri_fname,BaseName); strcat(tri_fname,".pts.tri");
  strcpy(nor_fname,BaseName); strcat(nor_fname,".nor");

  // oppening the output files for writting
  fpts = fopen (pts_fname,"wt");
  if(fpts==NULL){
    fprintf(stderr,"Error in reading the file %s\n",pts_fname);
    exit (1);
  }

  ftri = fopen (tri_fname,"wt");
  if(ftri==NULL){
    fprintf(stderr,"Error in reading the file %s\n",tri_fname);
    exit (1);
  }

  fnor = fopen (nor_fname,"wb");
  if(fnor==NULL){
    fprintf(stderr,"Error in reading the file %s\n",nor_fname);
    exit (1);
  }

  // copy vertices and normals ito the *.pts and *.nx, *.ny and *.nz files
  for (int i=0; i<vertlist.size(); i++){
    fprintf(fnor,"%lf %lf %lf\n",normlist[i][0],normlist[i][1],normlist[i][2]);
    fprintf(fpts,"%lf %lf %lf\n",vertlist[i][0],vertlist[i][1],vertlist[i][2]);
  }
  for (vector<Triangle>::const_iterator tri=trilist.begin(); tri!=trilist.end(); tri++){     
    fprintf(ftri,"%d %d %d\n",tri->indices[0], tri->indices[1], tri->indices[2]); 
  }
  fclose(fpts); fprintf(stdout,"%s saved ...\n",pts_fname);
  fclose(ftri); fprintf(stdout,"%s saved ...\n",tri_fname);
  fclose(fnor); fprintf(stdout,"%s saved ...\n",nor_fname);
}                                          
                                                                                                   
/*************************************************************************************
  Save caracteristic points defining metaballs into the file in world coordinates.
  These points are on the ends of the half radius of the metaballs(in our case
  exponential metaballs-elipsoids with exponential field function)
**************************************************************************************/
void SaveMb2Pts(char *in_fname, char *out_fname){
  VEC X(3), Lxp(3), Lxn(3), Lyp(3), Lyn(3), Lzp(3), Lzn(3), C(3);
  FILE *fin, *fout;
  int number;
  HMAT M;
  VEC  i1(3), j1(3), k1(3);
  double lx, ly, lz;
  
  // oppening the output files for writting
  fin = fopen(in_fname,"rt");
  if(fin==NULL){
    fprintf(stderr,"Error in reading the file %s\n",in_fname);
    exit (1);
  }
  
  // oppening the output files for writting
  fout = fopen(out_fname,"wt");
  if(fout==NULL){
    fprintf(stderr,"Error in reading the file %s\n",out_fname);
    exit (1);
  }

  // obtain the number of implicit surfaces
  fscanf(fin,"%d",&number);
  printf("Numb. of blobs: %d\n",number);
   
  for (int i=0;i<number;i++){
    // read metaball params from file
    fscanf(fin,"%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",\
            &C[0],&C[1],&C[2],&lx,&ly,&lz,&i1[0],&i1[1],&i1[2],&j1[0],&j1[1],&j1[2],&k1[0],&k1[1],&k1[2]);

    M(0,0)=i1[0]; M(0,1)=j1[0]; M(0,2)=k1[0]; M(0,3)= C[0];
    M(1,0)=i1[1]; M(1,1)=j1[1]; M(1,2)=k1[1]; M(1,3)= C[1]; 
    M(2,0)=i1[2]; M(2,1)=j1[2]; M(2,2)=k1[2]; M(2,3)= C[2];

    // set 8 points in the l.c.s defining ellipsoid and its center
    Lxp[0]=lx/2.0;  Lxp[1]=0.0;    Lxp[2]=0.0;
    Lyp[0]=0.0;     Lyp[1]=ly/2.0; Lyp[2]=0.0;
    Lzp[0]=0.0;     Lzp[1]=0.0;    Lzp[2]=lz/2.0;
    Lxn = -Lxp;
    Lyn = -Lyp;
    Lzn = -Lzp;
  
    // find their coordinates in g.c.s and write them to file
    fprintf(fout,"%lf %lf %lf\n",C[0],C[1],C[2]);
    X = M.fullMult(Lxp);
    fprintf(fout,"%lf %lf %lf\n",X[0],X[1],X[2]);
    X = M.fullMult(Lxn);
    fprintf(fout,"%lf %lf %lf\n",X[0],X[1],X[2]);
    X = M.fullMult(Lyp);
    fprintf(fout,"%lf %lf %lf\n",X[0],X[1],X[2]);
    X = M.fullMult(Lyn);
    fprintf(fout,"%lf %lf %lf\n",X[0],X[1],X[2]);
    X = M.fullMult(Lzp);
    fprintf(fout,"%lf %lf %lf\n",X[0],X[1],X[2]);
    X = M.fullMult(Lzn);
    fprintf(fout,"%lf %lf %lf\n",X[0],X[1],X[2]);
 
  }
}
                                   
                                                                                                   

void displayBBox(SoSeparator *group, const BBox &bb, float lineWidth, uint32_t color){
   // Line Set
   SoDrawStyle *drawStyle= new SoDrawStyle;
   drawStyle->style.setValue(SoDrawStyle::LINES);
   drawStyle->lineWidth.setValue(lineWidth);
   group->addChild(drawStyle);
   
   SoIndexedLineSet *ls = new SoIndexedLineSet;
   SoVertexProperty *vp = new SoVertexProperty;
   vp->materialBinding = SoVertexProperty::OVERALL;  
   ls->vertexProperty.setValue(vp);
   group->addChild(ls);
   vp->orderedRGBA.set1Value(0, color);

   int idx_vtx=0, idx=0;

   VEC IX(bb.getDim()), IY(bb.getDim()), IZ(bb.getDim());
   IX[1]=IX[2]=IY[0]=IY[2]=IZ[0]=IZ[1]=0;
   VEC O(bb.getMin());
             
   // construct points
   for (int p=0; p<8; p++)
   {
      int x= ((p&4)>>2);
      int y= ((p&2)>>1);
      int z= (p&1);
                
      VEC P = O + x*IX +  y*IY + z*IZ;
      vp->vertex.set1Value(idx_vtx+p, P);
   }
             
   // construct lines
   ls->coordIndex.set1Value(idx++, idx_vtx+0);
   ls->coordIndex.set1Value(idx++, idx_vtx+4);
   ls->coordIndex.set1Value(idx++, idx_vtx+6);
   ls->coordIndex.set1Value(idx++, idx_vtx+2);
   ls->coordIndex.set1Value(idx++, idx_vtx+0);
   ls->coordIndex.set1Value(idx++, SO_END_LINE_INDEX);

   ls->coordIndex.set1Value(idx++, idx_vtx+0);
   ls->coordIndex.set1Value(idx++, idx_vtx+1);
   ls->coordIndex.set1Value(idx++, idx_vtx+3);
   ls->coordIndex.set1Value(idx++, idx_vtx+2);
   ls->coordIndex.set1Value(idx++, idx_vtx+0);
   ls->coordIndex.set1Value(idx++, SO_END_LINE_INDEX);

   ls->coordIndex.set1Value(idx++, idx_vtx+1);
   ls->coordIndex.set1Value(idx++, idx_vtx+5);
   ls->coordIndex.set1Value(idx++, idx_vtx+7);
   ls->coordIndex.set1Value(idx++, idx_vtx+3);
   ls->coordIndex.set1Value(idx++, idx_vtx+1);
   ls->coordIndex.set1Value(idx++, SO_END_LINE_INDEX);

   ls->coordIndex.set1Value(idx++, idx_vtx+3);
   ls->coordIndex.set1Value(idx++, idx_vtx+2);
   ls->coordIndex.set1Value(idx++, idx_vtx+6);
   ls->coordIndex.set1Value(idx++, idx_vtx+7);
   ls->coordIndex.set1Value(idx++, idx_vtx+3);
   ls->coordIndex.set1Value(idx++, SO_END_LINE_INDEX);

   ls->coordIndex.set1Value(idx++, idx_vtx+0);
   ls->coordIndex.set1Value(idx++, idx_vtx+1);
   ls->coordIndex.set1Value(idx++, idx_vtx+5);
   ls->coordIndex.set1Value(idx++, idx_vtx+4);
   ls->coordIndex.set1Value(idx++, idx_vtx+0);
   ls->coordIndex.set1Value(idx++, SO_END_LINE_INDEX);

   ls->coordIndex.set1Value(idx++, idx_vtx+5);
   ls->coordIndex.set1Value(idx++, idx_vtx+4);
   ls->coordIndex.set1Value(idx++, idx_vtx+6);
   ls->coordIndex.set1Value(idx++, idx_vtx+7);
   ls->coordIndex.set1Value(idx++, idx_vtx+5);
   ls->coordIndex.set1Value(idx++, SO_END_LINE_INDEX);

   idx_vtx += 8;
          
   vp->vertex.setNum(idx_vtx);
   ls->coordIndex.setNum(idx);
}

/*
void
displayBBoxHierarchy(SoSeparator *group, const BBoxHierarchy *bbh)
{   
   int level = bbh->level;
   float scale = 1 - (0.02*level);
   BBox bb(bbh->bb);
   bb.scale(scale);
        
   displayBBox(group, bb, 3./(float)(1.+level), (uint32_t)((0xFF00 << 4*level)+ 0xFF));

   for (int i=0; i< bbh->bboxlist.size(); i++)
      displayBBoxHierarchy(group, bbh->bboxlist[i] );
}
*/


SoSeparator *drawVoxelisation(const marchingCubes &surf)
{
 
  SoSeparator *group = new SoSeparator;
  SoMaterial *material = new SoMaterial; 
  material->diffuseColor.setValue(0.8, 0.8, 0.8);
  group->addChild(material);
  SoDrawStyle *style = new SoDrawStyle;
  style->style = SoDrawStyle::FILLED;
  group->addChild(style);

  SoIndexedFaceSet *set = new SoIndexedFaceSet;
  group->addChild(set);
  SoVertexProperty *vp = new SoVertexProperty;
  vp->normalBinding = SoVertexProperty::PER_VERTEX_INDEXED;
  set->vertexProperty = vp;

  vector<VEC> vertlist = surf.getVertices();
  vector<VEC> normlist = surf.getNormals();
  vector<Triangle> trilist = surf.getTriangles();

  // copy vertices and normals into Vertex Property node
  for (int i=0; i<vertlist.size(); i++)
  {
     vp->normal.set1Value(i, normlist[i]);
     vp->vertex.set1Value(i, vertlist[i]);
  }
  vp->vertex.setNum(vertlist.size());
  vp->normal.setNum(normlist.size());
   
  int num=0;
  for (vector<Triangle>::const_iterator tri=trilist.begin(); tri!=trilist.end(); tri++)
  {     
     for (int k=0; k<3; k++)   
     {
        set->coordIndex.set1Value(num++, tri->indices[k]); 
     }
     set->coordIndex.set1Value(num++, SO_END_FACE_INDEX);
  }
  set->coordIndex.setNum(num);

  cout << "# tris = " << trilist.size() << " # points = " << vertlist.size() << endl;
  return group;
}

SoSeparator * readFile(const char *filename) 
{
   // Open the input file       
   SoInput mySceneInput;
   if (!mySceneInput.openFile(filename)) {
      fprintf(stderr, "Cannot open file %s\n", filename);
      return NULL;
   } 
   // Read the whole file into the database
   SoSeparator *myGraph = SoDB::readAll(&mySceneInput);
   if (myGraph == NULL) {
      fprintf(stderr, "Problem reading file\n");
      return NULL;
   } 
   mySceneInput.closeFile();
   return myGraph;
}

void writeFile(SoNode *node, const char *filename)
{
  SoWriteAction wa;
  wa.getOutput()->openFile(filename);
  wa.getOutput()->setBinary(0);
  wa.apply(node);
}

void analyseScene (SoSeparator *root, SoCoordinate3 **coords, SoNormal **norms, SoIndexedFaceSet **ifs)
{
   SoSearchAction mySearchAction;
   mySearchAction.setType(SoCoordinate3::getClassTypeId());
   mySearchAction.setInterest(SoSearchAction::FIRST);
   mySearchAction.apply(root);
   assert( mySearchAction.getPath()!=NULL );
   *coords = (SoCoordinate3*) mySearchAction.getPath()->getTail();
   (*coords)->ref();

   mySearchAction.setType(SoNormal::getClassTypeId());
   mySearchAction.setInterest(SoSearchAction::FIRST);
   mySearchAction.apply(root);
   assert( mySearchAction.getPath()!=NULL );
   *norms = (SoNormal*) mySearchAction.getPath()->getTail();
   (*norms)->ref();

   mySearchAction.setType(SoIndexedFaceSet::getClassTypeId());
   mySearchAction.setInterest(SoSearchAction::FIRST);
   mySearchAction.apply(root);
   assert( mySearchAction.getPath()!=NULL );
   *ifs = (SoIndexedFaceSet*) mySearchAction.getPath()->getTail();
   (*ifs)->ref();
}

#include <stdio.h>
void saveToPLY(const marchingCubes &surf, const char *filename)
{
   FILE *fp = fopen(filename, "w");
   vector<VEC> vertlist = surf.getVertices();
   vector<VEC> normlist = surf.getNormals();
   vector<Triangle> trilist = surf.getTriangles();

   fprintf(fp, "ply\nformat ascii 1.0\n");
   fprintf(fp, "element vertex %d\n", vertlist.size());
   fprintf(fp, "property float x\n");
   fprintf(fp, "property float y\n");
   fprintf(fp, "property float z\n");
   fprintf(fp, "element face %d\n", trilist.size());
   fprintf(fp, "property list uchar int vertex_indices\n");
   fprintf(fp, "end_header\n");

   for (int i=0; i<vertlist.size(); i++)
      fprintf(fp,"%.5f %.5f %.5f \n", vertlist[i][0], vertlist[i][1], vertlist[i][2]);

   for (int i=0; i<trilist.size(); i++)
      fprintf(fp,"3 %d %d %d \n", 
              trilist[i].indices[0], 
              trilist[i].indices[1], 
              trilist[i].indices[2]); 

   fclose(fp);
}

Mesh *buildMeshFromIv(const SoCoordinate3 *coords, const SoNormal *norms, const SoIndexedFaceSet *ifs)
{
  int i;
  // construct coords
  VECArray vertlist(coords->point.getNum());
  for (i=0; i<vertlist.size(); i++)
     vertlist[i] = VEC(coords->point[i].getValue(), 3);

  // construct normals
  VECArray normlist(norms->vector.getNum());
  for (i=0; i<normlist.size(); i++)
     normlist[i] = VEC(norms->vector[i].getValue(), 3);
  
  vector<int> triangles;
  for (int k, i=0; i<ifs->coordIndex.getNum();i=k)
  {
     k = i+1;
     while (ifs->coordIndex[k] != SO_END_FACE_INDEX) k++;

     if (k-i != 3) 
        cout << "found face other than triangle\n";
     else
     {          
        triangles.push_back(ifs->coordIndex[i]);
        triangles.push_back(ifs->coordIndex[i+1]);
        triangles.push_back(ifs->coordIndex[i+2]);
        while (ifs->coordIndex[k] == SO_END_FACE_INDEX) k++;
     }
  }

  vector<int> normIndex;
  for (int k, i=0; i<ifs->normalIndex.getNum();i=k)
  {
     k = i+1;
     while (ifs->normalIndex[k] != SO_END_FACE_INDEX) k++;

     if (k-i != 3) 
        cout << "found face other than triangle\n";
     else
     {          
        normIndex.push_back(ifs->normalIndex[i]);
        normIndex.push_back(ifs->normalIndex[i+1]);
        normIndex.push_back(ifs->normalIndex[i+2]);
        while (ifs->normalIndex[k] == SO_END_FACE_INDEX) k++;
     }
  }
  
  return new Mesh(vertlist, normlist, triangles, normIndex);
}


void main(int argc, char **argv)
{
  double res, coeff, stiff;
  char MbFile[256], OutFile[256];
  Widget myWindow = SoXt::init(argv[0]);
  if (argc<3) { PrintHelp(argv[0]); exit(-1); }
  ReadOptionDouble("-r",&res);
  ReadOptionDouble("-k",&stiff);
  ReadOptionDouble("-a",&coeff);
  ReadOptionString("-m",MbFile);
  ReadOptionString("-out",OutFile);
  int sP = IsOption("-e");
  MetaballsSphere_exp mb_e(skin_magic_blob,coeff);
  MetaballsSphere mb_s(skin_magic_sphere);

  if (myWindow == NULL) exit(1); 

  // construct a simple scene graph
  SoSeparator *root = new SoSeparator;
  root->ref();
  vector<Primitive*> primitives;

  if(sP){
    readSpheres(MbFile,mb_e,stiff);
    primitives.push_back(&mb_e);
    fprintf(stderr,"Added to primitives[-e]\n"); 
  }
  else{
    readSpheres(MbFile,mb_s,stiff);
    primitives.push_back(&mb_s);
    fprintf(stderr,"Added to primitives[-s]\n");
  }

     /*
  // }else if(IsOption("-e")){
    Metaballs_exp mb_e(skin_magic_blob); 
    readPrimitives(MbFile,mb_e);
    primitives.push_back(&mb_e);
    
  }else{
    printf("Error: Input file with metaball primitives is missing\n");
    exit(-1);
  }
    */
  marchingCubes surface(primitives, res);
  fprintf(stderr,"Marching cubes performed\n");
  /*
  SoSeparator *voxelisation = drawVoxelisation(surface);
  root->addChild(voxelisation);

  // Adding event handler for mouse clicked
  SoEventCallback *mouseEvent = new  SoEventCallback;
  mouseEvent->addEventCallback(SoMouseButtonEvent::getClassTypeId(),myMouseEvent,voxelisation);
  voxelisation->addChild(mouseEvent);


//   root->addChild(bone); 
//   displayBBoxHierarchy(root, &(mesh->bbox));

  writeFile(root, "./metaballs.iv");
  */
  if(strcmp(OutFile,""))
    SaveIsoSurf2Trig(surface, OutFile);
  //SaveMb2Pts(argv[2],argv[3]); 

  // Set up viewer:
  /*
  SoXtExaminerViewer *myViewer;
  myViewer = new SoXtExaminerViewer(myWindow);
  myViewer->setSceneGraph(root);
  myViewer->setTitle("Examiner Viewer");
  myViewer->show();

  SoXt::show(myWindow);
  SoXt::mainLoop();
  */
}

---