source: proiecte/pmake3d/make3d_original/Make3dSingleImageStanford_version0.1/third_party/vrippack-0.31/src/plytools/plyvolume.cc @ 37

/*

Compute the volume of a mesh

Sean Anderson, February 1999

*/

#include <assert.h>
#include <stdio.h>
#include <ply.h>
#include <stdlib.h>
#include <strings.h>
#include <iostream>
#include <limits.h>

#include <Linear.h>

#include <vector.h>

//#include "plyio.h"


typedef Vec3f PlyVertex;


struct PlyFace 
{
   uchar nverts;
   int * verts;
};

struct PlyTriStrip 
{
   int nverts;
   int * verts;
};

#define MAX_VERT_PROPS 3

static PlyProperty vert_prop_x =  
{"x", PLY_FLOAT, PLY_FLOAT, 0, 0, PLY_START_TYPE, PLY_START_TYPE, 0};
static PlyProperty vert_prop_y =  
{"y", PLY_FLOAT, PLY_FLOAT, 0, 0, PLY_START_TYPE, PLY_START_TYPE, 0};
static PlyProperty vert_prop_z =  
{"z", PLY_FLOAT, PLY_FLOAT, 0, 0, PLY_START_TYPE, PLY_START_TYPE, 0};


static PlyProperty vert_props[MAX_VERT_PROPS];


static PlyProperty face_props[] = 
{ 
   {"vertex_indices", PLY_INT, PLY_INT, 0, 1, PLY_UCHAR, PLY_UCHAR, 0},
};

static PlyProperty tristrip_props[] = 
{
   {"vertex_indices", PLY_INT, PLY_INT, 0, 1, PLY_INT, PLY_INT},
};

#define voffset(field) ((int) &(((PlyVertex *) 0)->field))
#define foffset(field) ((int) &(((PlyFace *) 0)->field))



class Plane
{
public:
   Plane()
   {
      n = Vec3f(0, 0, 1);
      d = 0;
   }
   Vec3f n;
   float d;
};

inline void
computeFaceContribution(const vector<PlyVertex> & verts, 
                        const PlyFace & plyFace, 
                        const Plane & plane, 
                        double * subvolume, 
                        double * subarea)
{
   const int * v = plyFace.verts;
   const int k = plyFace.nverts;
   Vec3f area;      
   float vol;
   Vec3f avePt;
   
   if (k == 3) // Make the common case fast:  triangle
   {
      area = 0.5 * 
         (verts[v[1]] - verts[v[0]]).cross(verts[v[2]] - verts[v[0]]);
      avePt = (1. / 3) * (verts[v[0]] + verts[v[1]] + verts[v[2]]);      
      vol = (plane.n.dot(avePt) - plane.d) * plane.n.dot(area);
   }
   else
   {
      area.clear();
      vol = 0;
      const Vec3f v0 = verts[v[0]];
      for (int i = 0; i < k - 1; i++)
      {
         Vec3f a = (verts[v[i]] - v0).cross(verts[v[i + 1]] - v0);
         area += a;
         avePt = (1. / 3) * (v0 + verts[v[i]] + verts[v[i + 1]]);      
         vol += (plane.n.dot(avePt) - plane.d) * plane.n.dot(a);
      }
      vol *= 0.5;
      area *= 0.5;
   }

   if (subvolume)
   {
      *subvolume += vol;
   }
   if (subarea)
   {
      *subarea += area.length();
   } 
}
            

//
// Read in data from the ply file and compute the volume and area of faces
// with plane.  The volume and area contributions are ADDED to what was 
// passed in.
//

bool
computeSignedVolume(FILE * fp, const Plane & plane, 
                    double * volume, double * area)
{
   int     j;
   int     nelems;
   char ** elist;
   int     file_type;
   float   version;
   char *  elem_name;
   int     num_vert_props;
   int     num_elems;
   
   vector<PlyVertex> verts;
   
   face_props[0].offset = foffset(verts);
   face_props[0].count_offset = foffset(nverts);  /* count offset */
    
   PlyFile * ply = 
      ply_read(fp, &nelems, &elist);

   
   if (!ply)
   {
      return true;
   }

   int nvp = 0;
   
   if (ply_is_valid_property(ply, "vertex", vert_prop_x.name)) 
   {
      vert_props[nvp] = vert_prop_x;
      vert_props[nvp].offset = voffset(x); 
      nvp++;
   }
   
   if (ply_is_valid_property(ply, "vertex", vert_prop_y.name)) 
   {
      vert_props[nvp] = vert_prop_y;
      vert_props[nvp].offset = voffset(y); 
      nvp++;
   }
   
   if (ply_is_valid_property(ply, "vertex", vert_prop_z.name)) 
   {
      vert_props[nvp] = vert_prop_z;
      vert_props[nvp].offset = voffset(z); 
      nvp++;
   }
   

   num_vert_props = nvp;

   PlyVertex plyVert;
   PlyFace plyFace;
   int faceIndicesMem[256];
   plyFace.verts = faceIndicesMem;

   for (int i = 0; i < nelems; i++) 
   {
      /* get the description of the first element */
      elem_name = elist[i];
      int nprops;
      ply_get_element_description(ply, elem_name, &num_elems, &nprops);
      
      /* if we're on vertex elements, read them in */
      if (equal_strings("vertex", elem_name)) 
      {
         /* set up for getting vertex elements */
         ply_get_element_setup(ply, elem_name, num_vert_props, vert_props);
         
         cerr << "Reading vertices . . .\n";

         int percent = -1;
         
         /* grab all the vertex elements */
         for (j = 0; j < num_elems; j++) 
         {
            int newpercent = 100 * j / num_elems;
            if (newpercent != percent)
            {
               cerr << "\r" << (percent = newpercent) << "% complete" ;
               cerr.flush();
            }
            
            ply_get_element_noalloc(ply, (void *) &plyVert);
            
            verts.push_back(plyVert);
         }
         
         cerr << "\rDone.        \n";
      } 
      else if (equal_strings("face", elem_name)) 
      {      
         ply_get_element_setup(ply, elem_name, 1, face_props);
         
         if (num_elems == 0)
         {
            continue;
         }

         cerr << "Reading faces . . .\n";
         
         if (verts.size() == 0)
         {
            cerr << "No vertices read yet; can't read faces\n";
            continue;
         }
         
         int percent = -1;
         
         for (j = 0; j < num_elems; j++) 
         {              
            int newpercent = 100 * j / num_elems;
            if (newpercent != percent)
            {
               cerr << "\r" << (percent = newpercent) << "% complete";
               cerr.flush();
            }
            
            ply_get_element_noalloc(ply, (void *) &plyFace);
            assert(plyFace.verts == faceIndicesMem);
            
            computeFaceContribution(verts, plyFace, plane, volume, area);
             
            //faces->push_back(plyFace);
         }
         
         cerr << "\rDone.        \n";
      }
      else if (equal_strings("tristrips", elem_name)) 
      {      
         ply_get_element_setup(ply, elem_name, 1, tristrip_props);
         
         if (num_elems == 0)
         {
            continue;
         }
                 
         cerr << "Reading tristrips . . .\n";
         
         if (verts.size() == 0)
         {
            cerr << "No vertices read yet; can't read tristrip\n";
            continue;
         }
         
         PlyTriStrip plyTriStrip;
         
         for (j = 0; j < num_elems; j++) 
         {
            ply_get_element(ply, (void *) &plyTriStrip);
            
            int percent = -1;
         
            bool clockwise = false;
            for (int k = 2; k < plyTriStrip.nverts; )
            {
               int newpercent = 100 * j / num_elems;
               if (newpercent != percent)
               {
                  cerr << "\r" << (percent = newpercent) << "% complete";
                  cerr.flush();
               }
            
               plyFace.nverts = 3;
               //plyFace.verts = new int[3];
               plyFace.verts[0]             = plyTriStrip.verts[k - 2];
               plyFace.verts[1 + clockwise] = plyTriStrip.verts[k - 1];
               plyFace.verts[2 - clockwise] = plyTriStrip.verts[k - 0];

               computeFaceContribution(verts, plyFace, plane, 
                                       volume, area);
               
               //faces->push_back(plyFace);            
               
               clockwise = !clockwise;

               k++;

               if (k < plyTriStrip.nverts && plyTriStrip.verts[k] == -1)
               {
                  k += 3;
                  clockwise = false;
               }
            }
            
            free(plyTriStrip.verts);
         }
         
         cerr << "\rDone.        \n";
      }
   }

   ply_close(ply);

   return false;
}

void
printUsage(const char * name)
{
   cerr << "Usage:\n"
        << "  " << name 
        << " -h | [-p nx ny nz d] (- | file1.ply [file2.ply [...]])\n"
        << "where:\n"
        << "  -h         Displays this help\n"
        << "  -p         Specifies a plane to use as the base; \n"
           "             the plane has normal (nx, ny, nz) and is distance \n"
           "             d from the origin in the normal's direction\n"
        << "  -          Reads one ply file from standard input\n"
        << "  file.ply   Is a ply file whose volume is added to the total.\n";
}

int 
main(int argc, char * argv[])
{
   double volume = 0, area = 0;
   Plane plane;
   
   for (int c = 1; c < argc; c++)
   {
      if (!strcmp("-h", argv[c]))
      {
         printUsage(argv[0]);
         return 1;
      }
      else if (!strcmp("-p", argv[c]))
      {
         if (c + 4 >= argc)
         {
            cerr << "Insufficient parameters for plane description\n";
            return 1;
         }
         else
         {
            plane.n.x = atof(argv[c + 1]);
            plane.n.y = atof(argv[c + 2]);
            plane.n.z = atof(argv[c + 3]);
            plane.d  = atof(argv[c + 4]);
            c += 4;
         }
      }
      else if (!strcmp("-", argv[c]))
      {
         cerr << "stdin:\n";
         if (computeSignedVolume(stdin, plane, &volume, &area))
         {
            cerr << "Couldn't read ply file from stdin\n";
         }
      }
      else
      {
         cerr << argv[c] << ":\n";
         FILE * fp = fopen(argv[c], "r");
         if (!fp)
         {
            cerr << "Couldn't open file '" << argv[c] << "'\n";
         }
         else
         {
            if (computeSignedVolume(fp, plane, &volume, &area))
            {
               cerr << "Couldn't read ply file '" << argv[c] << "'\n";
            }
            fclose(fp);
         }
      }
   }
   
   //cout.flags(ios::fixed);
   cout.precision(10);
   cout << "Volume:  " << volume << endl
        << "Area:  " << area << endl;
   
   return 0;
}