1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195
|
#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <vector>
#include <cmath>
#include <string>
using namespace std;
const double PI = 3.141592653589793;
const bool ANTICLOCKWISE = true; // Traversal of triangles, seen from OUTSIDE
//================================================
struct Pt // Coordinates of a vector or point
{
double x, y, z;
};
struct Triangle // Holds the INDICES and vector area of triangles
{
int v[3];
Pt area;
Pt normal;
};
Pt operator +( Pt p, Pt q ) { return { p.x + q.x, p.y + q.y, p.z + q.z }; }
Pt operator -( Pt p, Pt q ) { return { p.x - q.x, p.y - q.y, p.z - q.z }; }
Pt operator *( double r, Pt p ) { return { r * p.x, r * p.y, r * p.z }; }
Pt operator *( Pt p, double r ) { return r * p; }
Pt operator /( Pt p, double r ) { return { p.x / r, p.y / r, p.z / r }; }
double dot( Pt p, Pt q ) { return p.x * q.x + p.y * q.y + p.z * q.z; }
double abs( Pt p ) { return sqrt( dot( p, p ) ); }
Pt cross( Pt p, Pt q ) { return { p.y * q.z - p.z * q.y, p.z * q.x - p.x * q.z, p.x * q.y - p.y * q.x }; }
ostream &operator << ( ostream &strm, Pt p ) { return strm << p.x << " " << p.y << " " << p.z; }
istream &operator >> ( istream &strm, Pt &p ) { return strm >> p.x >> p.y >> p.z; }
//================================================
void readFile( string filename, vector<Pt> &vertices, vector<Triangle> &triangles )
{
vertices.clear(); // Comment out if adding to an existing collection
triangles.clear(); //
string line;
ifstream in( filename );
while ( getline( in, line ) ) // Read whole line
{
if ( line.find_first_of( "vVfF" ) == string::npos ) continue; // skip pointless lines
istringstream ss( line ); // Put line into a stream for input
string s;
ss >> s; // Get identifier
if ( s == "v" || s == "V" ) // VERTICES
{
Pt p;
ss >> p;
vertices.push_back( p );
}
else if ( s == "f" || s == "F" ) // FACES
{
string si, sj, sk;
ss >> si >> sj >> sk; // FIRST, read into individual strings
Triangle t;
t.v[0] = stoi( si ) - 1; // Get the FIRST integer from each string
t.v[1] = stoi( sj ) - 1; // NOTE: subtract 1 because arrays start from 0
t.v[2] = stoi( sk ) - 1;
triangles.push_back( t );
}
}
in.close();
}
//======================================================================
void calcFaces( const vector<Pt> &vertices, vector<Triangle> &triangles )
{
for ( int t = 0; t < triangles.size(); t++ )
{
int v0 = triangles[t].v[0]; // Global vertex index
int v1 = triangles[t].v[1];
int v2 = triangles[t].v[2];
Pt side1 = vertices[v1] - vertices[v0]; // Side vectors meeting at this vertex
Pt side2 = vertices[v2] - vertices[v0];
triangles[t].area = 0.5 * cross( side1, side2 ); // Triangle vector area is half the cross product of the sides
if ( !ANTICLOCKWISE ) triangles[t].area = -1.0 * triangles[t].area;
double A = abs( triangles[t].area );
triangles[t].normal = triangles[t].area / A; // Unit normal vector
}
}
//================================================
vector<Pt> calcNormals( const vector<Pt> &vertices, const vector<Triangle> &triangles )
{
vector<Pt> vnormal( vertices.size(), { 0.0, 0.0, 0.0 } );
for ( int t = 0; t < triangles.size(); t++ ) // Loop through triangle list
{
Pt norm = triangles[t].area;
norm = norm / abs( norm ); // Unit normal vector for this triangle
for ( int i = 0; i < 3; i++ ) // Update weighted sum at each vertex
{
int v0 = triangles[t].v[i]; // Global vertex index
int v1 = triangles[t].v[(i+1)%3];
int v2 = triangles[t].v[(i+2)%3];
Pt side1 = vertices[v1] - vertices[v0]; // Side vectors meeting at this vertex
Pt side2 = vertices[v2] - vertices[v0];
// double weight = acos( dot( side1, side2 ) / ( abs( side1 ) * abs( side2 ) ) ); // Angle enclosed at this vertex
double weight = 1.0; // Use this if you want unweighted averages
vnormal[v0] = vnormal[v0] + weight * norm; // Contribution to normal vector
}
}
for ( int v = 0; v < vertices.size(); v++ )
{
vnormal[v] = vnormal[v] / abs( vnormal[v] ); // Normalise to unit length
}
return vnormal;
}
//================================================
void output( ostream &out, const vector<Pt> &vertices, const vector<Pt> &vnormal )
{
int N = vertices.size();
// Formatting (adapt to your needs)
out.setf( ios::fixed ); // Fixed or scientific
out.precision( 6 ); // Decimal digits after point
#define SP << " " << setw( 12 ) << // Spacer and field width
for ( int i = 0; i < N; i++ )
{
out SP i
SP vertices[i].x SP vertices[i].y SP vertices[i].z
SP vnormal [i].x SP vnormal [i].y SP vnormal [i].z << '\n';
}
}
//================================================
void summary( const vector<Pt> &vertices, const vector<Triangle> &triangles )
{
cout << "\nVertices: " << vertices.size() << '\n';
for ( int i = 0; i < vertices.size(); i++ )
{
cout << i << ": " << vertices[i] << '\n';
}
cout << "\nTriangles: " << triangles.size() << '\n';
for ( int t = 0; t < triangles.size(); t++ )
{
cout << t << ": ";
for ( int i = 0; i < 3; i++ )
{
cout << triangles[t].v[i] << " ";
}
cout << triangles[t].normal << '\n';
}
}
//================================================
int main()
{
string filename = "input.obj";
vector<Pt> vertices;
vector<Triangle> triangles;
readFile( filename, vertices, triangles ); // Read vertices and triangles
calcFaces( vertices, triangles ); // Calculate area vectors for triangles
summary( vertices, triangles ); // Check
vector<Pt> vnormal = calcNormals( vertices, triangles ); // Calculate normals at vertices
cout << "\nVertices and vertex normals:\n";
output( cout, vertices, vnormal );
}
//================================================
|