Simplify.h

/////////////////////////////////////////////
//
// Mesh Simplification Tutorial
//
// (C) by Sven Forstmann in 2014
//
// License : MIT
// http://opensource.org/licenses/MIT
//
//https://github.com/sp4cerat/Fast-Quadric-Mesh-Simplification
//
// 5/2016: Chris Rorden created minimal version for OSX/Linux/Windows compile

//#include <iostream>
//#include <stddef.h>
//#include <functional>
//#include <sys/stat.h>
//#include <stdbool.h>
#include <string.h>
//#include <ctype.h>
//#include <float.h>
#include <stdio.h>
#include <stdlib.h>
#include <map>
#include <vector>
#include <string>
#include <math.h>
#include <float.h> //FLT_EPSILON, DBL_EPSILON


struct vector3
{
double x, y, z;
};

struct vec3f
{
    double x, y, z;

    inline vec3f( void ) {}

    //inline vec3f operator =( vector3 a )
    // { vec3f b ; b.x = a.x; b.y = a.y; b.z = a.z; return b;}

    inline vec3f( vector3 a )
     { x = a.x; y = a.y; z = a.z; }

    inline vec3f( const vec3f & a )
    { x = a.x; y = a.y; z = a.z; }

    inline vec3f( const double X, const double Y, const double Z )
    { x = X; y = Y; z = Z; }

    inline vec3f operator + ( const vec3f& a ) const
    { return vec3f( x + a.x, y + a.y, z + a.z ); }

    inline vec3f operator += ( const vec3f& a ) const
    { return vec3f( x + a.x, y + a.y, z + a.z ); }

    inline vec3f operator * ( const double a ) const
    { return vec3f( x * a, y * a, z * a ); }

    inline vec3f operator * ( const vec3f a ) const
    { return vec3f( x * a.x, y * a.y, z * a.z ); }

    inline vec3f v3 () const
    { return vec3f( x , y, z ); }

    inline vec3f operator = ( const vector3 a )
    { x=a.x;y=a.y;z=a.z;return *this; }

    inline vec3f operator = ( const vec3f a )
    { x=a.x;y=a.y;z=a.z;return *this; }

    inline vec3f operator / ( const vec3f a ) const
    { return vec3f( x / a.x, y / a.y, z / a.z ); }

    inline vec3f operator - ( const vec3f& a ) const
    { return vec3f( x - a.x, y - a.y, z - a.z ); }

    inline vec3f operator / ( const double a ) const
    { return vec3f( x / a, y / a, z / a ); }

    inline double dot( const vec3f& a ) const
    { return a.x*x + a.y*y + a.z*z; }

    inline vec3f cross( const vec3f& a , const vec3f& b )
    {
        x = a.y * b.z - a.z * b.y;
        y = a.z * b.x - a.x * b.z;
        z = a.x * b.y - a.y * b.x;
        return *this;
    }

    inline double angle( const vec3f& v )
    {
        vec3f a = v , b = *this;
        double dot = v.x*x + v.y*y + v.z*z;
        double len = a.length() * b.length();
        if(len==0)len=0.00001f;
        double input = dot  / len;
        if (input<-1) input=-1;
        if (input>1) input=1;
        return (double) acos ( input );
    }

    inline double angle2( const vec3f& v , const vec3f& w )
    {
        vec3f a = v , b= *this;
        double dot = a.x*b.x + a.y*b.y + a.z*b.z;
        double len = a.length() * b.length();
        if(len==0)len=1;

        vec3f plane; plane.cross( b,w );

        if ( plane.x * a.x + plane.y * a.y + plane.z * a.z > 0 )
            return (double) -acos ( dot  / len );

        return (double) acos ( dot  / len );
    }

    inline vec3f rot_x( double a )
    {
        double yy = cos ( a ) * y + sin ( a ) * z;
        double zz = cos ( a ) * z - sin ( a ) * y;
        y = yy; z = zz;
        return *this;
    }
    inline vec3f rot_y( double a )
    {
        double xx = cos ( -a ) * x + sin ( -a ) * z;
        double zz = cos ( -a ) * z - sin ( -a ) * x;
        x = xx; z = zz;
        return *this;
    }
    inline void clamp( double min, double max )
    {
        if (x<min) x=min;
        if (y<min) y=min;
        if (z<min) z=min;
        if (x>max) x=max;
        if (y>max) y=max;
        if (z>max) z=max;
    }
    inline vec3f rot_z( double a )
    {
        double yy = cos ( a ) * y + sin ( a ) * x;
        double xx = cos ( a ) * x - sin ( a ) * y;
        y = yy; x = xx;
        return *this;
    }
    inline vec3f invert()
    {
        x=-x;y=-y;z=-z;return *this;
    }
    inline vec3f frac()
    {
        return vec3f(
            x-double(int(x)),
            y-double(int(y)),
            z-double(int(z))
            );
    }

    inline vec3f integer()
    {
        return vec3f(
            double(int(x)),
            double(int(y)),
            double(int(z))
            );
    }

    inline double length() const
    {
        return (double)sqrt(x*x + y*y + z*z);
    }

    inline vec3f normalize()
    {
        double square = sqrt(x*x + y*y + z*z);
        x/=square;y/=square;z/=square;

        return *this;
    }
    static vec3f normalize( vec3f a );

    static void random_init();
    static double random_double();
    static vec3f random();

    static int random_number;

    double random_double_01(double a){
        double rnf=a*14.434252+a*364.2343+a*4213.45352+a*2341.43255+a*254341.43535+a*223454341.3523534245+23453.423412;
        int rni=((int)rnf)%100000;
        return double(rni)/(100000.0f-1.0f);
    }

    vec3f random01_fxyz(){
        x=(double)random_double_01(x);
        y=(double)random_double_01(y);
        z=(double)random_double_01(z);
        return *this;
    }

};

vec3f barycentric(const vec3f &p, const vec3f &a, const vec3f &b, const vec3f &c){
    vec3f v0 = b-a;
    vec3f v1 = c-a;
    vec3f v2 = p-a;
    double d00 = v0.dot(v0);
    double d01 = v0.dot(v1);
    double d11 = v1.dot(v1);
    double d20 = v2.dot(v0);
    double d21 = v2.dot(v1);
    double denom = d00*d11-d01*d01;
    double v = (d11 * d20 - d01 * d21) / denom;
    double w = (d00 * d21 - d01 * d20) / denom;
    double u = 1.0 - v - w;
    return vec3f(u,v,w);
}

vec3f interpolate(const vec3f &p, const vec3f &a, const vec3f &b, const vec3f &c, const vec3f attrs[3])
{
    vec3f bary = barycentric(p,a,b,c);
    vec3f out = vec3f(0,0,0);
    out = out + attrs[0] * bary.x;
    out = out + attrs[1] * bary.y;
    out = out + attrs[2] * bary.z;
    return out;
}

double min(double v1, double v2) {
    return fmin(v1,v2);
}


class SymetricMatrix {

    public:

    // Constructor

    SymetricMatrix(double c=0) { for (size_t i = 0; i < 10; ++i) { m[i] = c; } }

    SymetricMatrix(    double m11, double m12, double m13, double m14,
                        double m22, double m23, double m24,
                                    double m33, double m34,
                                                double m44) {
             m[0] = m11;  m[1] = m12;  m[2] = m13;  m[3] = m14;
                          m[4] = m22;  m[5] = m23;  m[6] = m24;
                                       m[7] = m33;  m[8] = m34;
                                                    m[9] = m44;
    }

    // Make plane

    SymetricMatrix(double a,double b,double c,double d)
    {
        m[0] = a*a;  m[1] = a*b;  m[2] = a*c;  m[3] = a*d;
                     m[4] = b*b;  m[5] = b*c;  m[6] = b*d;
                                  m[7 ] =c*c; m[8 ] = c*d;
                                               m[9 ] = d*d;
    }

    double operator[](int c) const { return m[c]; }

    // Determinant

    double det(    int a11, int a12, int a13,
                int a21, int a22, int a23,
                int a31, int a32, int a33)
    {
        double det =  m[a11]*m[a22]*m[a33] + m[a13]*m[a21]*m[a32] + m[a12]*m[a23]*m[a31]
                    - m[a13]*m[a22]*m[a31] - m[a11]*m[a23]*m[a32]- m[a12]*m[a21]*m[a33];
        return det;
    }

    const SymetricMatrix operator+(const SymetricMatrix& n) const
    {
        return SymetricMatrix( m[0]+n[0],   m[1]+n[1],   m[2]+n[2],   m[3]+n[3],
                                            m[4]+n[4],   m[5]+n[5],   m[6]+n[6],
                                                         m[ 7]+n[ 7], m[ 8]+n[8 ],
                                                                      m[ 9]+n[9 ]);
    }

    SymetricMatrix& operator+=(const SymetricMatrix& n)
    {
         m[0]+=n[0];   m[1]+=n[1];   m[2]+=n[2];   m[3]+=n[3];
         m[4]+=n[4];   m[5]+=n[5];   m[6]+=n[6];   m[7]+=n[7];
         m[8]+=n[8];   m[9]+=n[9];
        return *this;
    }

    double m[10];
};
///////////////////////////////////////////

namespace Simplify
{
    // Global Variables & Strctures
    enum Attributes {
        NONE,
        NORMAL = 2,
        TEXCOORD = 4,
        COLOR = 8
    };
    struct Triangle { int v[3];double err[4];int deleted,dirty,attr;vec3f n;vec3f uvs[3];int material; };
    struct Vertex { vec3f p;int tstart; size_t tcount;SymetricMatrix q;int border;};
    struct Ref { int tid,tvertex; };
    std::vector<Triangle> triangles;
    std::vector<Vertex> vertices;
    std::vector<Ref> refs;
    std::string mtllib;
    std::vector<std::string> materials;

    // Helper functions

    double vertex_error(SymetricMatrix q, double x, double y, double z);
    double calculate_error(int id_v1, int id_v2, vec3f &p_result);
    bool flipped(vec3f p,int i0,int i1,Vertex &v0,Vertex &v1,std::vector<int> &deleted);
    void update_uvs(int i0,const Vertex &v,const vec3f &p,std::vector<int> &deleted);
    void update_triangles(int i0,Vertex &v,std::vector<int> &deleted,int &deleted_triangles);
    void update_mesh(int iteration);
    void compact_mesh();
    //
    // Main simplification function
    //
    // target_count  : target nr. of triangles
    // agressiveness : sharpness to increase the threshold.
    //                 5..8 are good numbers
    //                 more iterations yield higher quality
    //

    void simplify_mesh(int target_count, double agressiveness=7, bool verbose=false)
    {
        // init
        for(Triangle& t: triangles) { t.deleted=0; }

        // main iteration loop
        int deleted_triangles=0;
        std::vector<int> deleted0,deleted1;
        int triangle_count=triangles.size();
        //int iteration = 0;
        //loop(iteration,0,100)
        for (int iteration = 0; iteration < 100; iteration ++)
        {
            if(triangle_count-deleted_triangles<=target_count)break;

            // update mesh once in a while
            if(iteration%5==0)
            {
                update_mesh(iteration);
            }

            // clear dirty flag
            for(Triangle& triangle: triangles) { triangle.dirty=0; }

            //
            // All triangles with edges below the threshold will be removed
            //
            // The following numbers works well for most models.
            // If it does not, try to adjust the 3 parameters
            //
            double threshold = 0.000000001*pow(double(iteration+3),agressiveness);

            // target number of triangles reached ? Then break
            if ((verbose) && (iteration%5==0)) {
                printf("iteration %d - triangles %d threshold %g\n",iteration,triangle_count-deleted_triangles, threshold);
            }

            // remove vertices & mark deleted triangles
            for(Triangle& t: triangles)
            {
                if(t.err[3]>threshold) continue;
                if(t.deleted) continue;
                if(t.dirty) continue;

                for(size_t j: {0, 1, 2})
                {
                    if(t.err[j] >= threshold)
                    {
                        continue;
                    }

                    int i0=t.v[ j     ]; Vertex &v0 = vertices[i0];
                    int i1=t.v[(j+1)%3]; Vertex &v1 = vertices[i1];
                    // Border check
                    if(v0.border != v1.border)  continue;

                    // Compute vertex to collapse to
                    vec3f p;
                    calculate_error(i0,i1,p);
                    deleted0.resize(v0.tcount); // normals temporarily
                    deleted1.resize(v1.tcount); // normals temporarily
                    // don't remove if flipped
                    if( flipped(p,i0,i1,v0,v1,deleted0) ) continue;

                    if( flipped(p,i1,i0,v1,v0,deleted1) ) continue;

                    if ( (t.attr & TEXCOORD) == TEXCOORD  )
                    {
                        update_uvs(i0,v0,p,deleted0);
                        update_uvs(i0,v1,p,deleted1);
                    }

                    // not flipped, so remove edge
                    v0.p=p;
                    v0.q=v1.q+v0.q;
                    int tstart=refs.size();

                    update_triangles(i0,v0,deleted0,deleted_triangles);
                    update_triangles(i0,v1,deleted1,deleted_triangles);

                    size_t tcount = refs.size() - tstart;

                    if(tcount<=v0.tcount)
                    {
                        // save ram
                        if(tcount)memcpy(&refs[v0.tstart],&refs[tstart],tcount*sizeof(Ref));
                    }
                    else
                        // append
                        v0.tstart=tstart;

                    v0.tcount=tcount;
                    break;
                }
                // done?
                if(triangle_count-deleted_triangles<=target_count)break;
            }
        }
        // clean up mesh
        compact_mesh();
    } //simplify_mesh()

    void simplify_mesh_lossless(bool verbose=false)
    {
        // init
        for(Triangle& t: triangles) { t.deleted=0; }

        // main iteration loop
        int deleted_triangles=0;
        std::vector<int> deleted0,deleted1;
        //int iteration = 0;
        //loop(iteration,0,100)
        for (int iteration = 0; iteration < 9999; iteration ++)
        {
            // update mesh constantly
            update_mesh(iteration);
            // clear dirty flag
            for(Triangle& t: triangles) { t.dirty=0; }
            //
            // All triangles with edges below the threshold will be removed
            //
            // The following numbers works well for most models.
            // If it does not, try to adjust the 3 parameters
            //
            double threshold = DBL_EPSILON; //1.0E-3 EPS;
            if (verbose) {
                printf("lossless iteration %d\n", iteration);
            }

            // remove vertices & mark deleted triangles
            for(Triangle& t: triangles)
            {
                if(t.err[3]>threshold) continue;
                if(t.deleted) continue;
                if(t.dirty) continue;

                for(size_t j: {0, 1, 2})
                {
                    if(t.err[j] >= threshold)
                    {
                        continue;
                    }
                    int i0=t.v[ j     ]; Vertex &v0 = vertices[i0];
                    int i1=t.v[(j+1)%3]; Vertex &v1 = vertices[i1];

                    // Border check
                    if(v0.border != v1.border)  continue;

                    // Compute vertex to collapse to
                    vec3f p;
                    calculate_error(i0,i1,p);

                    deleted0.resize(v0.tcount); // normals temporarily
                    deleted1.resize(v1.tcount); // normals temporarily

                    // don't remove if flipped
                    if( flipped(p,i0,i1,v0,v1,deleted0) ) continue;
                    if( flipped(p,i1,i0,v1,v0,deleted1) ) continue;

                    if ( (t.attr & TEXCOORD) == TEXCOORD )
                    {
                        update_uvs(i0,v0,p,deleted0);
                        update_uvs(i0,v1,p,deleted1);
                    }

                    // not flipped, so remove edge
                    v0.p=p;
                    v0.q=v1.q+v0.q;
                    size_t tstart = refs.size();

                    update_triangles(i0,v0,deleted0,deleted_triangles);
                    update_triangles(i0,v1,deleted1,deleted_triangles);

                    size_t tcount = refs.size() - tstart;

                    if(tcount <= v0.tcount)
                    {
                        // save ram
                        if(tcount)memcpy(&refs[v0.tstart],&refs[tstart],tcount*sizeof(Ref));
                    }
                    else
                        // append
                        v0.tstart=tstart;

                    v0.tcount=tcount;
                    break;
                }
            }
            if(deleted_triangles<=0)break;
            deleted_triangles=0;
        } //for each iteration
        // clean up mesh
        compact_mesh();
    } //simplify_mesh_lossless()


    // Check if a triangle flips when this edge is removed

    bool flipped(vec3f p,int i0,int i1,Vertex &v0,Vertex &v1,std::vector<int> &deleted)
    {
        (void)i0; // unused
        (void)v1; // unused

        for(size_t k = 0; k < v0.tcount; ++k)
        {
            Triangle &t=triangles[refs[v0.tstart+k].tid];
            if(t.deleted)continue;

            int s=refs[v0.tstart+k].tvertex;
            int id1=t.v[(s+1)%3];
            int id2=t.v[(s+2)%3];

            if(id1==i1 || id2==i1) // delete ?
            {

                deleted[k]=1;
                continue;
            }
            vec3f d1 = vertices[id1].p-p; d1.normalize();
            vec3f d2 = vertices[id2].p-p; d2.normalize();
            if(fabs(d1.dot(d2))>0.999) return true;
            vec3f n;
            n.cross(d1,d2);
            n.normalize();
            deleted[k]=0;
            if(n.dot(t.n)<0.2) return true;
        }
        return false;
    }

    // update_uvs

    void update_uvs(int i0,const Vertex &v,const vec3f &p,std::vector<int> &deleted)
    {
        (void)i0; // unused
        for(size_t k = 0; k < v.tcount; ++k)
        {
            Ref &r=refs[v.tstart+k];
            Triangle &t=triangles[r.tid];
            if(t.deleted)continue;
            if(deleted[k])continue;
            vec3f p1=vertices[t.v[0]].p;
            vec3f p2=vertices[t.v[1]].p;
            vec3f p3=vertices[t.v[2]].p;
            t.uvs[r.tvertex] = interpolate(p,p1,p2,p3,t.uvs);
        }
    }

    // Update triangle connections and edge error after a edge is collapsed

    void update_triangles(int i0,Vertex &v,std::vector<int> &deleted,int &deleted_triangles)
    {
        vec3f p;
        for(size_t k = 0; k < v.tcount; ++k)
        {
            Ref &r=refs[v.tstart+k];
            Triangle &t=triangles[r.tid];
            if(t.deleted)continue;
            if(deleted[k])
            {
                t.deleted=1;
                deleted_triangles++;
                continue;
            }
            t.v[r.tvertex]=i0;
            t.dirty=1;
            t.err[0]=calculate_error(t.v[0],t.v[1],p);
            t.err[1]=calculate_error(t.v[1],t.v[2],p);
            t.err[2]=calculate_error(t.v[2],t.v[0],p);
            t.err[3]=min(t.err[0],min(t.err[1],t.err[2]));
            refs.push_back(r);
        }
    }

    // compact triangles, compute edge error and build reference list

    void update_mesh(int iteration)
    {
        if(iteration>0) // compact triangles
        {
            int dst=0;
            for (Triangle &t: triangles)
            {
                if(t.deleted)
                {
                    continue;
                }
                triangles[dst] = t;
                dst++;
            }
            triangles.resize(dst);
        }
        //
        // Init Quadrics by Plane & Edge Errors
        //
        // required at the beginning ( iteration == 0 )
        // recomputing during the simplification is not required,
        // but mostly improves the result for closed meshes
        //
        if( iteration == 0 )
        {
            for(Vertex& v: vertices)
            {
                v.q = SymetricMatrix(0.0);
            }

            for(Triangle& t: triangles)
            {
                vec3f n,p[3];
                for(size_t j: {0, 1, 2})
                {
                    p[j] = vertices[t.v[j]].p;
                }
                n.cross(p[1]-p[0],p[2]-p[0]);
                n.normalize();
                t.n=n;
                for(size_t j: {0, 1, 2})
                {
                    vertices[t.v[j]].q = vertices[t.v[j]].q+SymetricMatrix(n.x,n.y,n.z,-n.dot(p[0]));
                }
            }
            for(Triangle& t: triangles)
            {
                // Calc Edge Error
                vec3f p;
                for(size_t j: {0, 1, 2})
                {
                    t.err[j] = calculate_error(t.v[j],t.v[(j+1)%3],p);
                }
                t.err[3]=min(t.err[0],min(t.err[1],t.err[2]));
            }
        }

        // Init Reference ID list
        for(Vertex& v: vertices)
        {
            v.tstart=0;
            v.tcount=0;
        }
        for(Triangle& t: triangles)
        {
            for(size_t j: {0, 1, 2}) { vertices[t.v[j]].tcount++; }
        }
        int tstart=0;
        for(Vertex& v: vertices)
        {
            v.tstart=tstart;
            tstart+=v.tcount;
            v.tcount=0;
        }

        // Write References
        refs.resize(triangles.size()*3);
        for (size_t i = 0; i < triangles.size(); ++i)
        {
            Triangle &t = triangles[i];
            for(size_t j: {0, 1, 2})
            {
                Vertex &v=vertices[t.v[j]];
                refs[v.tstart+v.tcount].tid=i;
                refs[v.tstart+v.tcount].tvertex=j;
                v.tcount++;
            }
        }

        // Identify boundary : vertices[].border=0,1
        if( iteration == 0 )
        {
            std::vector<int> vcount,vids;

            for(Vertex& v: vertices)
            {
                v.border=0;
            }

            for(Vertex& v: vertices)
            {
                vcount.clear();
                vids.clear();
                for(size_t j = 0; j < v.tcount; ++j)
                {
                    int k=refs[v.tstart+j].tid;
                    Triangle &t=triangles[k];
                    for(size_t k: {0, 1, 2})
                    {
                        size_t ofs=0;
                        int id=t.v[k];
                        while(ofs<vcount.size())
                        {
                            if(vids[ofs]==id)break;
                            ofs++;
                        }
                        if(ofs==vcount.size())
                        {
                            vcount.push_back(1);
                            vids.push_back(id);
                        }
                        else
                            vcount[ofs]++;
                    }
                }
                for(size_t j = 0; j < vcount.size(); ++j)
                {
                    if(vcount[j] == 1)
                    {
                        vertices[vids[j]].border = 1;
                    }
               }
            }
        }
    }

    // Finally compact mesh before exiting

    void compact_mesh()
    {
        int dst=0;
        for(Vertex& v: vertices)
        {
            v.tcount = 0;
        }
        for(Triangle& t: triangles)
        {
            if(t.deleted)
            {
                continue;
            }
            triangles[dst++]=t;
            for(size_t j: {0, 1, 2}) { vertices[t.v[j]].tcount=1; }
        }
        triangles.resize(dst);
        dst=0;
        for(Vertex& v: vertices)
        {
            if(v.tcount == 0)
            {
                continue;
            }
            v.tstart=dst;
            vertices[dst].p=v.p;
            dst++;
        }
        for(Triangle& t: triangles)
        {
            for(size_t j: {0, 1, 2}) { t.v[j]=vertices[t.v[j]].tstart; }
        }
        vertices.resize(dst);
    }

    // Error between vertex and Quadric

    double vertex_error(SymetricMatrix q, double x, double y, double z)
    {
         return   q[0]*x*x + 2*q[1]*x*y + 2*q[2]*x*z + 2*q[3]*x + q[4]*y*y
              + 2*q[5]*y*z + 2*q[6]*y + q[7]*z*z + 2*q[8]*z + q[9];
    }

    // Error for one edge

    double calculate_error(int id_v1, int id_v2, vec3f &p_result)
    {
        // compute interpolated vertex

        SymetricMatrix q = vertices[id_v1].q + vertices[id_v2].q;
        bool   border = vertices[id_v1].border & vertices[id_v2].border;
        double error=0;
        double det = q.det(0, 1, 2, 1, 4, 5, 2, 5, 7);
        if ( det != 0 && !border )
        {

            // q_delta is invertible
            p_result.x = -1/det*(q.det(1, 2, 3, 4, 5, 6, 5, 7 , 8));    // vx = A41/det(q_delta)
            p_result.y =  1/det*(q.det(0, 2, 3, 1, 5, 6, 2, 7 , 8));    // vy = A42/det(q_delta)
            p_result.z = -1/det*(q.det(0, 1, 3, 1, 4, 6, 2, 5,  8));    // vz = A43/det(q_delta)

            error = vertex_error(q, p_result.x, p_result.y, p_result.z);
        }
        else
        {
            // det = 0 -> try to find best result
            vec3f p1=vertices[id_v1].p;
            vec3f p2=vertices[id_v2].p;
            vec3f p3=(p1+p2)/2;
            double error1 = vertex_error(q, p1.x,p1.y,p1.z);
            double error2 = vertex_error(q, p2.x,p2.y,p2.z);
            double error3 = vertex_error(q, p3.x,p3.y,p3.z);
            error = min(error1, min(error2, error3));
            if (error1 == error) p_result=p1;
            if (error2 == error) p_result=p2;
            if (error3 == error) p_result=p3;
        }
        return error;
    }

    char *trimwhitespace(char *str)
    {
        char *end;

        // Trim leading space
        while(isspace((unsigned char)*str)) str++;

        if(*str == 0)  // All spaces?
        return str;

        // Trim trailing space
        end = str + strlen(str) - 1;
        while(end > str && isspace((unsigned char)*end)) end--;

        // Write new null terminator
        *(end+1) = 0;

        return str;
    }

    //Option : Load OBJ
    void load_obj(const char* filename, bool process_uv=false){
        vertices.clear();
        triangles.clear();
        //printf ( "Loading Objects %s ... \n",filename);
        FILE* fn;
        if(filename==NULL)        return ;
        if((char)filename[0]==0)    return ;
        if ((fn = fopen(filename, "rb")) == NULL)
        {
            printf ( "File %s not found!\n" ,filename );
            return;
        }
        char line[1000];
        memset ( line,0,1000 );
        int vertex_cnt = 0;
        int material = -1;
        std::map<std::string, int> material_map;
        std::vector<vec3f> uvs;
        std::vector<std::vector<int> > uvMap;

        while(fgets( line, 1000, fn ) != NULL)
        {
            Vertex v;
            vec3f uv;

            if (strncmp(line, "mtllib", 6) == 0)
            {
                mtllib = trimwhitespace(&line[7]);
            }
            if (strncmp(line, "usemtl", 6) == 0)
            {
                std::string usemtl = trimwhitespace(&line[7]);
                if (material_map.find(usemtl) == material_map.end())
                {
                    material_map[usemtl] = materials.size();
                    materials.push_back(usemtl);
                }
                material = material_map[usemtl];
            }

            if ( line[0] == 'v' && line[1] == 't' )
            {
                if ( line[2] == ' ' )
                {
                    if(sscanf(line,"vt %lf %lf",
                        &uv.x,&uv.y)==2)
                    {
                        uv.z = 0;
                        uvs.push_back(uv);
                    } else
                    if(sscanf(line,"vt %lf %lf %lf",
                        &uv.x,&uv.y,&uv.z)==3)
                    {
                        uvs.push_back(uv);
                    }
                }
            }
            else if ( line[0] == 'v' )
            {
                if ( line[1] == ' ' )
                {
                    if(sscanf(line,"v %lf %lf %lf",
                        &v.p.x,    &v.p.y,    &v.p.z)==3)
                    {
                        vertices.push_back(v);
                    }
                }
            }
            int integers[9];
            if ( line[0] == 'f' )
            {
                Triangle t;
                bool tri_ok = false;
                bool has_uv = false;

                if(sscanf(line,"f %d %d %d",
                    &integers[0],&integers[1],&integers[2])==3)
                {
                    tri_ok = true;
                }else
                if(sscanf(line,"f %d// %d// %d//",
                    &integers[0],&integers[1],&integers[2])==3)
                {
                    tri_ok = true;
                }else
                if(sscanf(line,"f %d//%d %d//%d %d//%d",
                    &integers[0],&integers[3],
                    &integers[1],&integers[4],
                    &integers[2],&integers[5])==6)
                {
                    tri_ok = true;
                }else
                if(sscanf(line,"f %d/%d/%d %d/%d/%d %d/%d/%d",
                    &integers[0],&integers[6],&integers[3],
                    &integers[1],&integers[7],&integers[4],
                    &integers[2],&integers[8],&integers[5])==9)
                {
                    tri_ok = true;
                    has_uv = true;
                }
                else
                {
                    printf("unrecognized sequence\n");
                    printf("%s\n",line);
                    while(1);
                }
                if ( tri_ok )
                {
                    t.v[0] = integers[0]-1-vertex_cnt;
                    t.v[1] = integers[1]-1-vertex_cnt;
                    t.v[2] = integers[2]-1-vertex_cnt;
                    t.attr = 0;

                    if ( process_uv && has_uv )
                    {
                        std::vector<int> indices;
                        indices.push_back(integers[6]-1-vertex_cnt);
                        indices.push_back(integers[7]-1-vertex_cnt);
                        indices.push_back(integers[8]-1-vertex_cnt);
                        uvMap.push_back(indices);
                        t.attr |= TEXCOORD;
                    }

                    t.material = material;
                    //geo.triangles.push_back ( tri );
                    triangles.push_back(t);
                    //state_before = state;
                    //state ='f';
                }
            }
        }

        if ( process_uv && uvs.size() )
        {
            for(size_t i = 0; i < triangles.size(); ++i)
            {
                Triangle& t = triangles[i];
                for(size_t j: {0, 1, 2})
                {
                    t.uvs[j] = uvs[uvMap[i][j]];
                }
            }
        }

        fclose(fn);

        //printf("load_obj: vertices = %lu, triangles = %lu, uvs = %lu\n", vertices.size(), triangles.size(), uvs.size() );
    } // load_obj()

    // Optional : Store as OBJ

    void write_obj(const char* filename)
    {
        FILE *file=fopen(filename, "w");
        int cur_material = -1;
        bool has_uv = (triangles.size() && (triangles[0].attr & TEXCOORD) == TEXCOORD);

        if (!file)
        {
            printf("write_obj: can't write data file \"%s\".\n", filename);
            exit(0);
        }
        if (!mtllib.empty())
        {
            fprintf(file, "mtllib %s\n", mtllib.c_str());
        }
        for(Vertex& v: vertices)
        {
            //fprintf(file, "v %lf %lf %lf\n", v.p.x, v.p.y, v.p.z);
            fprintf(file, "v %g %g %g\n", v.p.x, v.p.y, v.p.z); //more compact: remove trailing zeros
        }
        if (has_uv)
        {
            for(Triangle& t: triangles)
            {
                if(t.deleted)
                {
                    continue;
                }
                fprintf(file, "vt %g %g\n", t.uvs[0].x, t.uvs[0].y);
                fprintf(file, "vt %g %g\n", t.uvs[1].x, t.uvs[1].y);
                fprintf(file, "vt %g %g\n", t.uvs[2].x, t.uvs[2].y);
            }
        }
        int uv = 1;
        for(Triangle& t: triangles)
        {
            if(t.deleted)
            {
                continue;
            }
            if (t.material != cur_material)
            {
                cur_material = t.material;
                fprintf(file, "usemtl %s\n", materials[t.material].c_str());
            }
            if (has_uv)
            {
                fprintf(file, "f %d/%d %d/%d %d/%d\n", t.v[0]+1, uv, t.v[1]+1, uv+1, t.v[2]+1, uv+2);
                uv += 3;
            }
            else
            {
                fprintf(file, "f %d %d %d\n", t.v[0]+1, t.v[1]+1, t.v[2]+1);
            }
            //fprintf(file, "f %d// %d// %d//\n", t.v[0]+1, t.v[1]+1, t.v[2]+1); //more compact: remove trailing zeros
        }
        fclose(file);
    }
};
///////////////////////////////////////////