Initial Commit.

CMakeLists.txt

@@ -0,0 +1,10 @@
+
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8)
+
+project(voxelizer)
+
+file(GLOB_RECURSE HEADER_CODE ${voxelizer_SOURCE_DIR}/include/*.h)
+file(GLOB_RECURSE SRC_CODE ${voxelizer_SOURCE_DIR}/source/*.cpp)
+
+ADD_EXECUTABLE(voxelizer ${SRC_CODE} ${HEADER_CODE})
+

README.md

@@ -0,0 +1,6 @@
+**Procedural City Generation**
+
+A C++ engine for procedural city generation by manipulating .obj models and outputing an .obj model of a city.
+
+Use cmake file included to generate project.
+Current command line arguments required: (project directory)/data/(desired .obj model to use) output.obj

data/CSG.obj

@@ -0,0 +1,41 @@
+# File produced by Open Asset Import Library (http://www.assimp.sf.net)
+# (assimp v3.1.187496374)
+
+mtllib CSG.obj.mtl
+
+# 8 vertex positions
+v  0 1 3
+v  1 0 3
+v  1 1 3
+v  0 0 3
+v  0 0 0
+v  1 1 0
+v  1 0 0
+v  0 1 0
+
+# 0 UV coordinates
+
+# 6 vertex normals
+vn 0 0 1
+vn 0 0 -1
+vn 0 -1 0
+vn 1 0 0
+vn 0 1 0
+vn -1 0 0
+
+# Mesh 'OpenSCAD_Model' with 12 faces
+g OpenSCAD_Model
+usemtl DefaultMaterial
+f  1//1 2//1 3//1
+f  2//1 1//1 4//1
+f  5//2 6//2 7//2
+f  6//2 5//2 8//2
+f  5//3 2//3 4//3
+f  2//3 5//3 7//3
+f  2//4 6//4 3//4
+f  6//4 2//4 7//4
+f  6//5 1//5 3//5
+f  1//5 6//5 8//5
+f  5//6 1//6 8//6
+f  1//6 5//6 4//6
+

data/test.obj

@@ -0,0 +1,61 @@
+# File produced by Open Asset Import Library (http://www.assimp.sf.net)
+# (assimp v3.1.187496374)
+
+mtllib test.obj.mtl
+
+# 16 vertex positions
+v  0 0 0
+v  0 1 3
+v  0 1 0
+v  0 0 3
+v  1 0 3
+v  1 1 3
+v  1 1 0
+v  1 0 0
+v  2 0 0
+v  2 1 3
+v  2 1 0
+v  2 0 3
+v  3 0 3
+v  3 1 3
+v  3 1 0
+v  3 0 0
+
+# 0 UV coordinates
+
+# 6 vertex normals
+vn -1 0 0
+vn 0 0 1
+vn 1 0 0
+vn 0 1 0
+vn 0 0 -1
+vn 0 -1 0
+
+# Mesh 'OpenSCAD_Model' with 24 faces
+g OpenSCAD_Model
+usemtl DefaultMaterial
+f  1//1 2//1 3//1
+f  2//1 1//1 4//1
+f  2//2 5//2 6//2
+f  5//2 2//2 4//2
+f  5//3 7//3 6//3
+f  7//3 5//3 8//3
+f  7//4 2//4 6//4
+f  2//4 7//4 3//4
+f  1//5 7//5 8//5
+f  7//5 1//5 3//5
+f  1//6 5//6 4//6
+f  5//6 1//6 8//6
+f  9//1 10//1 11//1
+f  10//1 9//1 12//1
+f  10//2 13//2 14//2
+f  13//2 10//2 12//2
+f  13//3 15//3 14//3
+f  15//3 13//3 16//3
+f  15//4 10//4 14//4
+f  10//4 15//4 11//4
+f  9//5 15//5 16//5
+f  15//5 9//5 11//5
+f  9//6 13//6 12//6
+f  13//6 9//6 16//6
+

data/test.stl

@@ -0,0 +1,170 @@
+solid OpenSCAD_Model
+  facet normal -1 0 0
+    outer loop
+      vertex 0 0 0
+      vertex 0 1 3
+      vertex 0 1 0
+    endloop
+  endfacet
+  facet normal -1 -0 0
+    outer loop
+      vertex 0 1 3
+      vertex 0 0 0
+      vertex 0 0 3
+    endloop
+  endfacet
+  facet normal -0 0 1
+    outer loop
+      vertex 0 1 3
+      vertex 1 0 3
+      vertex 1 1 3
+    endloop
+  endfacet
+  facet normal 0 0 1
+    outer loop
+      vertex 1 0 3
+      vertex 0 1 3
+      vertex 0 0 3
+    endloop
+  endfacet
+  facet normal 1 -0 0
+    outer loop
+      vertex 1 0 3
+      vertex 1 1 0
+      vertex 1 1 3
+    endloop
+  endfacet
+  facet normal 1 0 0
+    outer loop
+      vertex 1 1 0
+      vertex 1 0 3
+      vertex 1 0 0
+    endloop
+  endfacet
+  facet normal 0 1 -0
+    outer loop
+      vertex 1 1 0
+      vertex 0 1 3
+      vertex 1 1 3
+    endloop
+  endfacet
+  facet normal 0 1 0
+    outer loop
+      vertex 0 1 3
+      vertex 1 1 0
+      vertex 0 1 0
+    endloop
+  endfacet
+  facet normal 0 0 -1
+    outer loop
+      vertex 0 0 0
+      vertex 1 1 0
+      vertex 1 0 0
+    endloop
+  endfacet
+  facet normal -0 0 -1
+    outer loop
+      vertex 1 1 0
+      vertex 0 0 0
+      vertex 0 1 0
+    endloop
+  endfacet
+  facet normal 0 -1 0
+    outer loop
+      vertex 0 0 0
+      vertex 1 0 3
+      vertex 0 0 3
+    endloop
+  endfacet
+  facet normal 0 -1 -0
+    outer loop
+      vertex 1 0 3
+      vertex 0 0 0
+      vertex 1 0 0
+    endloop
+  endfacet
+  facet normal -1 0 0
+    outer loop
+      vertex 2 0 0
+      vertex 2 1 3
+      vertex 2 1 0
+    endloop
+  endfacet
+  facet normal -1 -0 0
+    outer loop
+      vertex 2 1 3
+      vertex 2 0 0
+      vertex 2 0 3
+    endloop
+  endfacet
+  facet normal -0 0 1
+    outer loop
+      vertex 2 1 3
+      vertex 3 0 3
+      vertex 3 1 3
+    endloop
+  endfacet
+  facet normal 0 0 1
+    outer loop
+      vertex 3 0 3
+      vertex 2 1 3
+      vertex 2 0 3
+    endloop
+  endfacet
+  facet normal 1 -0 0
+    outer loop
+      vertex 3 0 3
+      vertex 3 1 0
+      vertex 3 1 3
+    endloop
+  endfacet
+  facet normal 1 0 0
+    outer loop
+      vertex 3 1 0
+      vertex 3 0 3
+      vertex 3 0 0
+    endloop
+  endfacet
+  facet normal 0 1 -0
+    outer loop
+      vertex 3 1 0
+      vertex 2 1 3
+      vertex 3 1 3
+    endloop
+  endfacet
+  facet normal 0 1 0
+    outer loop
+      vertex 2 1 3
+      vertex 3 1 0
+      vertex 2 1 0
+    endloop
+  endfacet
+  facet normal 0 0 -1
+    outer loop
+      vertex 2 0 0
+      vertex 3 1 0
+      vertex 3 0 0
+    endloop
+  endfacet
+  facet normal -0 0 -1
+    outer loop
+      vertex 3 1 0
+      vertex 2 0 0
+      vertex 2 1 0
+    endloop
+  endfacet
+  facet normal 0 -1 0
+    outer loop
+      vertex 2 0 0
+      vertex 3 0 3
+      vertex 2 0 3
+    endloop
+  endfacet
+  facet normal 0 -1 -0
+    outer loop
+      vertex 3 0 3
+      vertex 2 0 0
+      vertex 3 0 0
+    endloop
+  endfacet
+endsolid OpenSCAD_Model

include/CompFab.h

@@ -0,0 +1,139 @@
+//
+//  CompFab.h
+//  voxelizer
+//
+//  Created by David Levin on 2/3/14.
+//  Modified by Shinjiro Sueda on 8/20/16.
+//
+//
+
+#ifndef voxelizer_CompFab_h
+#define voxelizer_CompFab_h
+
+#define EPSILON 1e-9
+
+#include <cmath>
+
+namespace CompFab
+{
+    
+    //Data Types
+    typedef struct Vec3Struct
+    {
+        
+        Vec3Struct();
+        Vec3Struct(double x, double y, double z);
+        
+        union
+        {
+            double m_pos[3];
+            struct { double m_x,m_y,m_z; };
+        };
+        
+        inline double & operator[](unsigned int index) { return m_pos[index]; }
+        inline const double & operator[](unsigned int index) const { return m_pos[index]; }
+        inline void operator+=(const Vec3Struct &a)
+        {
+            m_x += a.m_x;
+            m_y += a.m_y;
+            m_z += a.m_z;
+        }
+        
+        void normalize();
+        
+    }Vec3;
+    
+    //Data Types
+    typedef struct Vec3iStruct
+    {
+        
+        Vec3iStruct();
+        Vec3iStruct(int x, int y, int z);
+        union
+        {
+            int m_pos[3];
+            struct {int m_x,m_y,m_z;};
+        };
+        
+        inline int & operator[](unsigned int index) { return m_pos[index]; }
+        inline const int & operator[](unsigned int index) const { return m_pos[index]; }
+        
+    }Vec3i;
+    
+    //Data Types
+    typedef struct Vec2fStruct
+    {
+        
+        Vec2fStruct();
+        Vec2fStruct(float x, float y);
+        
+        union
+        {
+            float m_pos[2];
+            struct { float m_x,m_y; };
+        };
+        
+        inline float & operator[](unsigned int index) { return m_pos[index]; }
+        inline const float & operator[](unsigned int index) const { return m_pos[index]; }
+        
+    }Vec2f;
+    
+    
+    //NOTE: Ray direction must be normalized
+    typedef struct RayStruct
+    {
+        
+        RayStruct();
+        RayStruct(Vec3 &origin, Vec3 &direction);
+        
+        Vec3 m_origin;
+        Vec3 m_direction;
+        
+    } Ray;
+    
+    typedef struct TriangleStruct
+    {
+        
+        TriangleStruct(Vec3 &v1, Vec3 &v2,Vec3 &v3);
+        
+        Vec3 m_v1, m_v2, m_v3;
+        
+    } Triangle;
+    
+    //Some useful operations
+    //Compute v1 - v2
+    Vec3 operator-(const Vec3 &v1, const Vec3 &v2);
+    
+    Vec3 operator+(const Vec3 &v1, const Vec3 &v2);
+    
+    //Cross Product
+    Vec3 operator%(const Vec3 &v1, const Vec3 &v2);
+    
+    //Dot Product
+    double operator*(const Vec3 &v1, const Vec3 &v2);
+    
+    
+    //Grid structure for Voxels
+    typedef struct VoxelGridStruct
+    {
+        //Square voxels only
+        VoxelGridStruct(Vec3 lowerLeft, unsigned int dimX, unsigned int dimY, unsigned int dimZ, double spacing);
+        virtual ~VoxelGridStruct();
+        
+        inline bool & isInside(unsigned int i, unsigned int j, unsigned int k)
+        {
+            return m_insideArray[k*(m_numX*m_numY)+j*m_numY + i];
+        }
+        
+        unsigned int m_numX, m_numY, m_numZ; // Number of voxels in each direction
+        unsigned int m_size; // Total number of voxels
+        bool *m_insideArray; // An array of all voxels
+        double m_spacing; // spacing between voxels
+        Vec3 m_corner; // corner of the whole grid with minimum x, minimum y, and minimum z
+        
+    } VoxelGrid;
+}
+
+
+
+#endif

include/Mesh.h

@@ -0,0 +1,55 @@
+#ifndef MESH_H
+#define MESH_H
+
+#include <map>
+#include <vector>
+#include <fstream>
+#include "../include/CompFab.h"
+
+class Mesh{
+public:
+    std::vector<CompFab::Vec3>v;
+    std::vector<CompFab::Vec3>n;
+    std::vector<CompFab::Vec2f>tex;
+    std::vector<CompFab::Vec3i>texId;
+    ///@brief triangles
+    std::vector<CompFab::Vec3i>t;
+    
+    Mesh();
+    Mesh(const std::vector<CompFab::Vec3>&_v,
+         const std::vector<CompFab::Vec3i>&_t);
+    Mesh(const CompFab::Vec3 * _v, const CompFab::Vec3i * _t);
+    
+    Mesh(const char * filename,bool normalize);
+    virtual ~Mesh();
+    void load_mesh(const char * filename, bool normalize=true);
+    void save(const char * filename);
+    void save(std::ostream &out, std::vector<CompFab::Vec3> *vert=0);
+    void load(std::istream &in);
+    void read_ply(std::istream & f);
+    void read_obj(std::istream &f);
+    
+    void save_obj(const char * filename);
+    void load_tex(const char * filename);
+    
+    void compute_norm();
+    void rescale();
+    void append(const Mesh & m);
+    Mesh & operator= (const Mesh& m);
+    virtual void update();
+};
+void makeCube(Mesh & m, const CompFab::Vec3 & mn,
+              const CompFab::Vec3 mx);
+///@brief cube [0,1]^3
+extern Mesh UNIT_CUBE;
+void BBox(const Mesh & m, CompFab::Vec3 & mn,
+          CompFab::Vec3 & mx);
+
+void BBox(const std::vector<CompFab::Vec3> & v, CompFab::Vec3 & mn,
+          CompFab::Vec3 & mx);
+
+bool ptInBox(const CompFab::Vec3 & mn,
+             const CompFab::Vec3 mx, const CompFab::Vec3 & x);
+void adjlist(const Mesh & m, std::vector<std::vector<int> > & adjMat);
+
+#endif

source/CompFab.cpp

@@ -0,0 +1,149 @@
+//
+//  CompFab.cpp
+//  voxelizer
+//
+//  Created by David Levin on 2/3/14.
+//
+//
+
+#include "../include/CompFab.h"
+
+using namespace CompFab;
+
+
+CompFab::Vec3Struct::Vec3Struct()
+{
+    m_x = m_y = m_z = 0.0;
+}
+
+CompFab::Vec3Struct::Vec3Struct(double x, double y, double z)
+{
+    m_x = x;
+    m_y = y;
+    m_z = z;
+}
+
+void CompFab::Vec3Struct::normalize() {
+    
+    double magnitude = sqrt(m_x*m_x+m_y*m_y+m_z*m_z);
+    
+    if(magnitude > EPSILON)
+    {
+        m_x /= magnitude;
+        m_y /= magnitude;
+        m_z /= magnitude;
+    }
+}
+
+//Data Types
+CompFab::Vec3iStruct::Vec3iStruct()
+{
+    m_x = m_y = m_z = 0;
+}
+
+CompFab::Vec3iStruct::Vec3iStruct(int x, int y, int z)
+{
+    m_x = x;
+    m_y = y;
+    m_z = z;
+}
+
+CompFab::Vec2fStruct::Vec2fStruct()
+{
+    m_x = m_y = 0.0;
+}
+
+CompFab::Vec2fStruct::Vec2fStruct(float x, float y)
+{
+    m_x = x;
+    m_y = y;
+}
+
+CompFab::RayStruct::RayStruct()
+{
+    m_origin[0] = m_origin[1] = m_origin[2] = 0.0;
+    m_direction[0] = 1.0;
+    m_direction[1] = m_direction[2] = 0.0;
+}
+
+CompFab::RayStruct::RayStruct(Vec3 &origin, Vec3 &direction)
+{
+    m_origin = origin;
+    m_direction = direction;
+}
+
+CompFab::TriangleStruct::TriangleStruct(Vec3 &v1, Vec3 &v2,Vec3 &v3)
+{
+    m_v1 = v1;
+    m_v2 = v2;
+    m_v3 = v3;
+}
+
+CompFab::Vec3 CompFab::operator-(const Vec3 &v1, const Vec3 &v2)
+{
+    Vec3 v3;
+    v3[0] = v1[0] - v2[0];
+    v3[1] = v1[1] - v2[1];
+    v3[2] = v1[2] - v2[2];
+    
+    return v3;
+}
+
+CompFab::Vec3 CompFab::operator+(const Vec3 &v1, const Vec3 &v2)
+{
+    Vec3 v3;
+    v3[0] = v1[0] + v2[0];
+    v3[1] = v1[1] + v2[1];
+    v3[2] = v1[2] + v2[2];
+    
+    return v3;
+}
+
+
+//Cross Product
+Vec3 CompFab::operator%(const Vec3 &v1, const Vec3 &v2)
+{
+    Vec3 v3;
+    v3[0] = v1[1]*v2[2] - v1[2]*v2[1];
+    v3[1] = v1[2]*v2[0] - v1[0]*v2[2];
+    v3[2] = v1[0]*v2[1] - v1[1]*v2[0];
+    
+    return v3;
+}
+
+//Dot Product
+double CompFab::operator*(const Vec3 &v1, const Vec3 &v2)
+{
+    return v1.m_x*v2.m_x + v1.m_y*v2.m_y+v1.m_z*v2.m_z;
+}
+
+
+//Grid structure for Voxels
+CompFab::VoxelGridStruct::VoxelGridStruct(Vec3 corner, unsigned int numX, unsigned int numY, unsigned int numZ, double spacing)
+{
+    m_corner = corner;
+    m_numX = numX;
+    m_numY = numY;
+    m_numZ = numZ;
+    m_size = numX*numY*numZ;
+    m_spacing = spacing;
+    
+    //Allocate Memory
+    m_insideArray = new bool[m_size];
+    
+    for(unsigned int ii=0; ii<m_size; ++ii)
+    {
+        m_insideArray[ii] = false;
+    }
+    
+}
+
+CompFab::VoxelGridStruct::~VoxelGridStruct()
+{
+    delete[] m_insideArray;
+}
+
+
+
+
+

source/Mesh.cpp

@@ -0,0 +1,466 @@
+#include "../include/Mesh.h"
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+#ifdef _WIN32
+#define NOMINMAX //Stop errors with std::max
+#include <windows.h>
+#endif
+#include <stdio.h>
+#include <cstdlib>
+#include <utility>
+#include <map>
+#include <sstream>
+#include <string.h>
+//#include "util.h"
+
+typedef double real_t;
+
+///@brief is a point inside a box
+bool ptInBox(const CompFab::Vec3 & mn,
+             const CompFab::Vec3 mx, const CompFab::Vec3 & x)
+{
+    for(int dim = 0 ;dim<3;dim++){
+        if(x[dim]<mn[dim] || x[dim] > mx[dim]){
+            return false;
+        }
+    }
+    return true;
+}
+
+void makeCube(Mesh & m, const CompFab::Vec3 & mn,
+              const CompFab::Vec3 mx)
+{
+    CompFab::Vec3 ss = mx -mn;
+    m=UNIT_CUBE;
+    for(unsigned int ii = 0;ii<m.v.size();ii++){
+        m.v[ii][0] = mn[0] + ss[0]*m.v[ii][0];
+        m.v[ii][1] = mn[1] + ss[1]*m.v[ii][1];
+        m.v[ii][2] = mn[2] + ss[2]*m.v[ii][2];
+        
+    }
+}
+
+void Mesh::append(const Mesh & m)
+{
+    unsigned int offset = v.size();
+    unsigned int ot = t.size();
+    v.insert(v.end(),m.v.begin(),m.v.end());
+    t.insert(t.end(),m.t.begin(), m.t.end());
+    for(unsigned int ii = ot;ii<t.size();ii++){
+        for(int jj = 0 ;jj<3;jj++){
+            t[ii][jj] += offset;
+        }
+    }
+}
+
+Mesh & Mesh::operator= (const Mesh& m)
+{
+    v = m.v;
+    t = m.t;
+    n = m.n;
+    return *this;
+}
+
+///@brief cube [0,1]^3
+CompFab::Vec3 CUBE_VERT[8]={
+    CompFab::Vec3 (0, 0, 0),
+    CompFab::Vec3 (1, 0, 0),
+    CompFab::Vec3 (1, 1, 0),
+    CompFab::Vec3 (0, 1, 0),
+    CompFab::Vec3 (0, 0, 1),
+    CompFab::Vec3 (1, 0, 1),
+    CompFab::Vec3 (1, 1, 1),
+    CompFab::Vec3 (0, 1, 1)
+};
+
+CompFab::Vec3i CUBE_TRIG[12]={CompFab::Vec3i(0,3,1),
+    CompFab::Vec3i(1, 3, 2),
+    CompFab::Vec3i(5, 4, 0),
+    CompFab::Vec3i(5, 0, 1),
+    CompFab::Vec3i(6, 5, 1),
+    CompFab:: Vec3i(1, 2, 6),
+    CompFab:: Vec3i(3, 6, 2),
+    CompFab:: Vec3i(3, 7, 6),
+    CompFab:: Vec3i(4, 3, 0),
+    CompFab:: Vec3i(4, 7, 3),
+    CompFab:: Vec3i(7, 4, 5),
+    CompFab:: Vec3i(7, 5, 6)};
+Mesh UNIT_CUBE(CUBE_VERT,CUBE_TRIG);
+
+Mesh::Mesh():v(0),t(0){}
+
+Mesh::Mesh(const std::vector<CompFab::Vec3>&_v,
+           const std::vector<CompFab::Vec3i>&_t):v(_v),t(_t)
+{
+    compute_norm();
+}
+
+Mesh::Mesh(const CompFab::Vec3 * _v,
+           const CompFab::Vec3i * _t)
+{
+    v.assign(_v,_v+8);
+    t.assign(_t,_t+12);
+    
+    compute_norm();
+}
+
+Mesh::~Mesh()
+{
+    
+}
+
+void Mesh::save(std::ostream & out, std::vector<CompFab::Vec3> * vert)
+{
+    std::string vTok("v");
+    std::string fTok("f");
+    std::string texTok("vt");
+    char bslash='/';
+    std::string tok;
+    if(vert==0){
+        vert = &v;
+    }
+    for(size_t ii=0;ii<vert->size();ii++){
+        out<<vTok<<" "<<(*vert)[ii][0]<<" "<<(*vert)[ii][1]<<" "<<(*vert)[ii][2]<<"\n";
+    }
+    if(tex.size()>0){
+        for(size_t ii=0;ii<tex.size();ii++){
+            out<<texTok<<" "<<tex[ii][0]<<" "<<tex[ii][1]<<"\n";
+        }
+        for(size_t ii=0;ii<t.size();ii++){
+            out<<fTok<<" "<<t[ii][0]+1<<bslash<<texId[ii][0]+1<<" "
+            <<t[ii][1]+1<<bslash<<texId[ii][1]+1<<" "
+            <<t[ii][2]+1<<bslash<<texId[ii][2]+1<<"\n";
+        }
+    }else{
+        for(size_t ii=0;ii<t.size();ii++){
+            out<<fTok<<" "<<t[ii][0]+1<<" "<<
+            t[ii][1]+1<<" "<<t[ii][2]+1<<"\n";
+        }
+    }
+    out<<"#end\n";
+}
+
+void Mesh::save(const char * filename)
+{
+    std::ofstream out;
+    out.open(filename);
+    save(out);
+    out.close();
+}
+
+
+void Mesh::load(std::istream &in)
+{
+    read_obj(in);
+}
+
+void Mesh::read_obj(std::istream & f)
+{
+    std::string line;
+    std::string vTok("v");
+    std::string fTok("f");
+    std::string texTok("vt");
+    char bslash='/',space=' ';
+    std::string tok;
+    while(1) {
+        std::getline(f,line);
+        if(f.eof()) {
+            break;
+        }
+        if(line == "#end"){
+            break;
+        }
+        if(line.size()<3) {
+            continue;
+        }
+        if(line.at(0)=='#') {
+            continue;
+        }
+        std::stringstream ss(line);
+        ss>>tok;
+        if(tok==vTok) {
+            CompFab::Vec3 vec;
+            ss>>vec[0]>>vec[1]>>vec[2];
+            v.push_back(vec);
+        } else if(tok==fTok) {
+            bool hasTexture = false;
+            if (line.find(bslash) != std::string::npos) {
+                std::replace(line.begin(), line.end(), bslash, space);
+                hasTexture = true;
+            }
+            std::stringstream facess(line);
+            facess>>tok;
+            std::vector<int> vidx;
+            std::vector<int> texIdx;
+            int x;
+            while(facess>>x){
+                vidx.push_back(x);
+                if(hasTexture){
+                    facess>>x;
+                    texIdx.push_back(x);
+                }
+            }
+            texIdx.resize(vidx.size());
+            for(int ii = 0;ii<(int)vidx.size()-2;ii++){
+                CompFab::Vec3i trig, textureId;
+                trig[0] = vidx[0]-1;
+                textureId[0] = texIdx[0]-1;
+                for (int jj = 1; jj < 3; jj++) {
+                    trig[jj] = vidx[ii+jj]-1;
+                    textureId[jj] = texIdx[ii+jj]-1;
+                }
+                t.push_back(trig);
+                texId.push_back(textureId);
+            }
+        } else if(tok==texTok) {
+            CompFab::Vec2f texcoord;
+            ss>>texcoord[0];
+            ss>>texcoord[1];
+            tex.push_back(texcoord);
+        }
+    }
+    std::cout<<"Num Triangles: "<< t.size()<<"\n";
+}
+
+void Mesh::read_ply(std::istream & f)
+{
+    std::string line;
+    std::string vertLine("element vertex");
+    std::string faceLine("element face");
+    std::string texLine("property float s");
+    std::string endHeaderLine("end_header");
+    while(true) {
+        std::getline(f,line);
+        if(std::string::npos!=line.find(vertLine)) {
+            break;
+        }
+    }
+    std::string token;
+    std::stringstream ss(line);
+    ss>>token>>token;
+    int nvert;
+    ss>>nvert;
+    bool hasTex=false;
+    while(true) {
+        std::getline(f,line);
+        if(std::string::npos!=line.find(faceLine)) {
+            break;
+        }
+        if(std::string::npos!=line.find(texLine)) {
+            hasTex=true;
+        }
+    }
+    std::stringstream ss1(line);
+    ss1>>token>>token;
+    int nface;
+    ss1>>nface;
+    while(true) {
+        std::getline(f,line);
+        if(std::string::npos!=line.find(endHeaderLine)) {
+            break;
+        }
+    }
+    
+    v.resize(nvert);
+    t.resize(nface);
+    if(hasTex) {
+        tex.resize(nvert);
+    }
+    for (int ii =0; ii<nvert; ii++) {
+        for (int jj=0; jj<3; jj++) {
+            f>>v[ii][jj];
+        }
+        if(hasTex) {
+            for (int jj=0; jj<2; jj++) {
+                f>>tex[ii][jj];
+            }
+            tex[ii][1]=1-tex[ii][1];;
+        }
+    }
+    for (int ii =0; ii<nface; ii++) {
+        int nidx;
+        f>>nidx;
+        for (int jj=0; jj<3; jj++) {
+            f>>t[ii][jj];
+        }
+    }
+}
+
+void Mesh::save_obj(const char * filename)
+{
+    std::ofstream out(filename);
+    if(!out.good()){
+        std::cout<<"cannot open output file"<<filename<<"\n";
+        return;
+    }
+    save(out);
+    out.close();
+}
+
+void Mesh::update()
+{}
+
+Mesh::Mesh(const char * filename,bool normalize)
+{
+    load_mesh(filename,normalize);
+}
+
+
+void Mesh::load_mesh(const char * filename, bool normalize)
+{
+    std::ifstream f ;
+    f.open(filename);
+    if(!f.good()) {
+        std::cout<<"Error: cannot open mesh "<<filename<<"\n";
+        return;
+    }
+    switch(filename[strlen(filename)-1]) {
+        case 'y':
+            read_ply(f);
+            break;
+        case 'j':
+            read_obj(f);
+            break;
+        default:
+            break;
+    }
+    if(normalize){
+        rescale();
+    }
+    compute_norm();
+    
+    f.close();
+}
+
+void Mesh::rescale()
+{
+    if(v.size()==0){
+        std::cout<<"empty mesh\n";
+        return;
+    }
+    CompFab::Vec3 mn=v[0],mx=v[0];
+    
+    //scale and translate to [0 , 1]
+    for (unsigned int dim = 0; dim<3; dim++) {
+        for( size_t ii=0; ii<v.size(); ii++) {
+            mn[dim]= std::min(v[ii][dim],mn[dim]);
+            mx[dim] = std::max(v[ii][dim],mx[dim]);
+        }
+        real_t translate = -mn[dim];
+        for(size_t ii=0; ii<v.size(); ii++) {
+            v[ii][dim]=(v[ii][dim]+translate);
+        }
+    }
+    
+    real_t scale = 1/(mx[0]-mn[0]);
+    for(unsigned int dim=1; dim<3; dim++) {
+        scale=std::min(1/(mx[dim]-mn[dim]),scale);
+    }
+    
+    for(size_t ii=0; ii<v.size(); ii++) {
+        for (unsigned int dim = 0; dim<3; dim++) {
+            v[ii][dim]=v[ii][dim]*scale;
+        }
+    }
+}
+
+void Mesh::compute_norm()
+{
+    CompFab::Vec3 ZERO;
+    
+    n.resize(v.size(), ZERO);
+    for(unsigned int ii=0; ii<t.size(); ii++) {
+        CompFab::Vec3 a = v[t[ii][1]] - v[t[ii][0]];
+        CompFab::Vec3 b = v[t[ii][2]] - v[t[ii][0]];
+        b=a%b;
+        b.normalize();
+        for(int jj=0; jj<3; jj++) {
+            n[t[ii][jj]]+=b;
+            if(t[ii][jj]>=(int)n.size() || t[ii][jj]<0){
+                std::cout<<ii<<" "<<jj<<" "<<t[ii][jj]<<" normal computation error\n";
+            }
+        }
+    }
+    for(unsigned int ii=0; ii<v.size(); ii++) {
+        n[ii].normalize();
+    }
+}
+
+void BBox(const Mesh & m,
+          CompFab::Vec3 & mn, CompFab::Vec3 & mx)
+{
+    BBox(m.v, mn, mx);
+}
+
+bool is_nbr(const CompFab::Vec3i & a, const CompFab::Vec3i&b, int vert)
+{
+    for (int ii=0; ii<3; ii++) {
+        
+        int va=a[ii];
+        if(va<=vert) {
+            continue;
+        }
+        
+        for (unsigned int jj=0; jj<3; jj++) {
+            int vb=b[jj];
+            if(vb<=vert) {
+                continue;
+            }
+            if(va==vb) {
+                return true;
+            }
+        }
+    }
+    return false;
+}
+
+
+void adjlist(const Mesh & m, std::vector<std::vector<int> > & adjMat)
+{
+    if(adjMat.size()==m.t.size()) {
+        return;
+    }
+    std::vector<std::vector<int> >trigList;
+    trigList.resize(m.v.size());
+    for (unsigned int ii=0; ii<m.t.size(); ii++) {
+        for (unsigned int jj=0; jj<3; jj++) {
+            int vidx=m.t[ii][jj];
+            trigList[vidx].push_back(ii);
+        }
+    }
+    adjMat.resize(m.t.size());
+    for (unsigned int ii=0; ii<m.v.size(); ii++) {
+        int n_nbr=trigList[ii].size();
+        for (int jj=0; jj<n_nbr; jj++) {
+            int tj=trigList[ii][jj];
+            for (int kk=(jj+1); kk<n_nbr; kk++) {
+                int tk=trigList[ii][kk];
+                if(is_nbr(m.t[tj],m.t[tk],ii)) {
+                    adjMat[tj].push_back(tk);
+                    adjMat[tk].push_back(tj);
+                }
+                
+            }
+        }
+    }
+}
+
+
+void BBox(const std::vector<CompFab::Vec3 >& v,
+          CompFab::Vec3 & mn, CompFab::Vec3 & mx)
+{
+    mn = v[0];
+    mx = v[0];
+    for(unsigned int ii = 1 ;ii<v.size();ii++){
+        for(int dim = 0 ; dim<3;dim++){
+            if(v[ii][dim]<mn[dim]){
+                mn[dim] = v[ii][dim];
+            }
+            if(v[ii][dim]>mx[dim]){
+                mx[dim] = v[ii][dim];
+            }
+        }
+    }
+}
+

source/main.cpp

@@ -0,0 +1,292 @@
+//Computational Fabrication Assignment #1
+// Created by David Levin 2014
+// Modified by Shinjiro Sueda 2016
+
+#include <iostream>
+#include <vector>
+#include "../include/CompFab.h"
+#include "../include/Mesh.h"
+
+//Triangle list (global)
+typedef std::vector<CompFab::Triangle> TriangleList;
+
+TriangleList g_triangleList;
+CompFab::VoxelGrid *g_voxelGrid;
+
+bool loadMesh(char *filename, unsigned int num)
+{
+    g_triangleList.clear();
+    
+    Mesh *tempMesh = new Mesh(filename, true);
+    
+    CompFab::Vec3 v1, v2, v3;
+    
+    //copy triangles to global list
+    for(unsigned int tri =0; tri<tempMesh->t.size(); ++tri)
+    {
+        v1 = tempMesh->v[tempMesh->t[tri][0]];
+        v2 = tempMesh->v[tempMesh->t[tri][1]];
+        v3 = tempMesh->v[tempMesh->t[tri][2]];
+        g_triangleList.push_back(CompFab::Triangle(v1,v2,v3));
+    }
+    
+    //Create Voxel Grid
+    CompFab::Vec3 bbMax, bbMin;
+    BBox(*tempMesh, bbMin, bbMax);
+    
+    //Build Voxel Grid
+    double bbX = bbMax[0] - bbMin[0];
+    double bbY = bbMax[1] - bbMin[1];
+    double bbZ = bbMax[2] - bbMin[2];
+    double spacing;
+    
+    if(bbX > bbY && bbX > bbZ)
+    {
+        spacing = bbX/(num-1);
+    } else if(bbY > bbX && bbY > bbZ) {
+        spacing = bbY/(num-1);
+    } else {
+        spacing = bbZ/(num-1);
+    }
+    
+    g_voxelGrid = new CompFab::VoxelGrid(bbMin, num, num, num, spacing);
+    
+    delete tempMesh;
+    
+    return true;
+    
+}
+
+void saveVoxelsToObj(const char * outfile)
+{
+    Mesh box;
+    Mesh mout;
+    int nx = g_voxelGrid->m_numX;
+    int ny = g_voxelGrid->m_numY;
+    int nz = g_voxelGrid->m_numZ;
+    double spacing = g_voxelGrid->m_spacing;
+    
+    CompFab::Vec3 hspacing(0.5*spacing, 0.5*spacing, 0.5*spacing);
+    
+    for (int ii = 0; ii < nx; ii++) {
+        for (int jj = 0; jj < ny; jj++) {
+            for (int kk = 0; kk < nz; kk++) {
+                if(!g_voxelGrid->isInside(ii,jj,kk)){
+                    continue;
+                }
+                CompFab::Vec3 coord(0.5f + ((double)ii)*spacing, 0.5f + ((double)jj)*spacing, 0.5f+((double)kk)*spacing);
+                CompFab::Vec3 box0 = coord - hspacing;
+                CompFab::Vec3 box1 = coord + hspacing;
+                makeCube(box, box0, box1);
+                mout.append(box);
+            }
+        }
+    }
+    
+    mout.save_obj(outfile);
+}
+
+// Code referenced is sourced in README.txt
+bool triangleIntersection(CompFab::Vec3 V1, CompFab::Vec3 V2, CompFab::Vec3 V3, CompFab::Vec3 O, CompFab::Vec3 D)
+{
+	CompFab::Vec3 e1, e2;  //Edge1, Edge2
+	CompFab::Vec3 P, Q, T;
+	double det, inv_det, u, v;
+	double t;
+
+	//Find vectors for two edges sharing V1
+	e1 = V2 - V1;
+	e2 = V3 - V1;
+	//Begin calculating determinant - also used to calculate u parameter
+	P = D % e2;
+	//if determinant is near zero, ray lies in plane of triangle or ray is parallel to plane of triangle
+	det = e1 * P;
+	//NOT CULLING
+	if (det > -EPSILON && det < EPSILON) return false;
+	inv_det = 1.f / det;
+
+	//calculate distance from V1 to ray origin
+	T = O - V1;
+
+	//Calculate u parameter and test bound
+	u = (T * P) * inv_det;
+	//The intersection lies outside of the triangle
+	if (u < 0.f || u > 1.f) return false;
+
+	//Prepare to test v parameter
+	Q = T % e1;
+
+	//Calculate V parameter and test bound
+	v = (D * Q) * inv_det;
+	//The intersection lies outside of the triangle
+	if (v < 0.f || u + v  > 1.f) return false;
+
+	t = (e2 * Q) * inv_det;
+
+	if (t > EPSILON) { //ray intersection
+		return true;
+	}
+
+	// No hit, no win
+	return false;
+}
+
+int castRays(CompFab::Vec3 t_origin)
+{
+
+	// 6 separate counter variables for the 6 raycasts. Additional coutner to count how many of them counted a voxel as being "inside"
+	int rc = 0, c1 = 0, c2 = 0, c3 = 0, c4 = 0, c5 = 0, c6 = 0;
+
+	// Calculating triangle intersections for each triangle with all 6 raycasts.
+	for (int i = 0; i < g_triangleList.size(); i++)
+	{
+		if (triangleIntersection(g_triangleList.at(i).m_v1, g_triangleList.at(i).m_v2, g_triangleList.at(i).m_v3, t_origin, CompFab::Vec3(1, 0, 0)))
+		{
+			++c1;
+		}
+		if (triangleIntersection(g_triangleList.at(i).m_v1, g_triangleList.at(i).m_v2, g_triangleList.at(i).m_v3, t_origin, CompFab::Vec3(-1, 0, 0)))
+		{
+			++c2;
+		}
+		if (triangleIntersection(g_triangleList.at(i).m_v1, g_triangleList.at(i).m_v2, g_triangleList.at(i).m_v3, t_origin, CompFab::Vec3(0, 1, 0)))
+		{
+			++c3;
+		}
+		if (triangleIntersection(g_triangleList.at(i).m_v1, g_triangleList.at(i).m_v2, g_triangleList.at(i).m_v3, t_origin, CompFab::Vec3(0, -1, 0)))
+		{
+			++c4;
+		}
+		if (triangleIntersection(g_triangleList.at(i).m_v1, g_triangleList.at(i).m_v2, g_triangleList.at(i).m_v3, t_origin, CompFab::Vec3(0, 0, 1)))
+		{
+			++c5;
+		}
+		if (triangleIntersection(g_triangleList.at(i).m_v1, g_triangleList.at(i).m_v2, g_triangleList.at(i).m_v3, t_origin, CompFab::Vec3(0, 0, -1)))
+		{
+			++c6;
+		}
+	}
+
+	// Determining which triangle intersections were valid.
+	if (c1 % 2 == 1) ++rc;
+	if (c2 % 2 == 1) ++rc;
+	if (c3 % 2 == 1) ++rc;
+	if (c4 % 2 == 1) ++rc;
+	if (c5 % 2 == 1) ++rc;
+	if (c6 % 2 == 1) ++rc;
+
+	return rc;
+}
+
+void findLW(Mesh &m, double &l, double &w)
+{
+    double minl, maxl, minw, maxw;
+    for(int i = 0; i < m.t.size(); i++)
+    {
+        if(m.v[i].m_x < minl)
+        {
+            minl = m.v[i].m_x;
+        }
+        if(m.v[i].m_x > maxl)
+        {
+            maxl = m.v[i].m_x;
+        }
+        if(m.v[i].m_y < minw)
+        {
+            minw = m.v[i].m_x;
+        }
+        if(m.v[i].m_y > maxw)
+        {
+            maxw = m.v[i].m_x;
+        }
+    }
+    l = maxl - minl;
+    w = maxw - minw;
+}
+
+int main(int argc, char **argv)
+{
+    unsigned int num = 16; //number of voxels (e.g. 16x16x16)
+    
+    if(argc < 4)
+    {
+        Mesh *test = new Mesh(argv[1], false);
+        Mesh *output = new Mesh(test->v, test->t);
+        double l = 0, w = 0;
+        double *length = &l, *width = &w;
+        findLW(*test, *length, *width);
+        for(int i = 0; i < test->v.size(); i++)
+        {
+            output->v.push_back(*new CompFab::Vec3(test->v[i].m_x + *length + 1, test->v[i].m_y, test->v[i].m_z));
+        }
+        for(int k = 0; k < test->t.size(); k++)
+        {
+            output->t.push_back(*new CompFab::Vec3i(test->t[k].m_x + test->v.size(), test->t[k].m_y + test->v.size(), test->t[k].m_z + test->v.size()));
+        }
+        for(int j = 0; j < output->v.size(); j++)
+        {
+            std::cout << output->v[j].m_x << " " << output->v[j].m_y << " " << output->v[j].m_z << std::endl;
+            std::cout << output->t[j].m_x << " " << output->t[j].m_y << " " << output->t[j].m_z << std::endl;
+            std::cout << std::endl;
+        }
+        output->save(argv[2]);
+    }
+    else
+    {
+        num = atoi(argv[3]);
+        
+        // The loadMesh() function loads the mesh and then creates:
+        // - g_triangleList: The list of triangles in the input OBJ mesh
+        // - g_voxelGrid: The VoxelGrid object, with all voxels marked as unoccupied
+        std::cout<<"Load Mesh : "<<argv[1]<<"\n";
+        loadMesh(argv[1], num);
+        
+        std::cout<<"Voxelizing into "<<num<<"x"<<num<<"x"<<num<<"\n";
+        
+        // Below, write a triple for-loop to go through all the voxels in X, Y, Z.
+        //   g_voxelGrid->m_numX is the number of voxels in the X direction.
+        //   g_voxelGrid->m_numY is the number of voxels in the Y direction.
+        //   g_voxelGrid->m_numZ is the number of voxels in the Z direction.
+        //   g_voxelGrid->m_spacing is the size of each voxel.
+        //
+        // Inside the triple for-loop, check if the voxel is inside or outside the
+        // mesh. Use the g_voxelGrid->isInside(...) method to set whether that voxel
+        // is inside or outside. E.g.,
+        //   g_voxelGrid->isInside(0,0,0) = false;
+
+        // <TRIPLE FOR-LOOP HERE>
+
+        // Temp varaible for the current corner
+        CompFab::Vec3 t_corner = g_voxelGrid->m_corner;
+
+        // Temp varaible for the current center
+        CompFab::Vec3 t_origin;
+
+        for (unsigned int i_x = 0; i_x < g_voxelGrid->m_numX; i_x++)
+        {
+            for (unsigned int i_y = 0; i_y < g_voxelGrid->m_numY; i_y++)
+            {
+                for (unsigned int i_z = 0; i_z < g_voxelGrid->m_numZ; i_z++)
+                {
+
+                    // Finding the corner of the current voxel.
+                    t_corner.m_x = g_voxelGrid->m_corner.m_x + i_x*g_voxelGrid->m_spacing;
+                    t_corner.m_y = g_voxelGrid->m_corner.m_y + i_y*g_voxelGrid->m_spacing;
+                    t_corner.m_z = g_voxelGrid->m_corner.m_z + i_z*g_voxelGrid->m_spacing;
+
+                    // Finding the origin of the current voxel.
+                    t_origin.m_x = t_corner.m_x + 0.5*g_voxelGrid->m_spacing;
+                    t_origin.m_y = t_corner.m_y + 0.5*g_voxelGrid->m_spacing;
+                    t_origin.m_z = t_corner.m_z + 0.5*g_voxelGrid->m_spacing;
+
+                    // Will say "inside" if 4 or more of the raycasts determine a voxel to be inside.
+                    g_voxelGrid->isInside(i_x, i_y, i_z) = (castRays(t_origin) >= 4);
+                }
+            }
+        }
+        
+        // Write out voxel data as obj
+        saveVoxelsToObj(argv[2]);
+    }
+    
+    delete g_voxelGrid;
+}