Ellipsoid tesselation

267

Author

June 09, 2009 08:41 AM

I've simply adapted the alogrithm of sphere tesselation, when starting with an octahedron. Which means replacing the vertices, (a,0,0),(0,a,0),(0,0,a) with (a,0,0),(0,b,0),(0,0,c) and their opposites accordingly. I subdivide a triangle whose vertices are already on the ellipse, by four triangles that have the same area( they each contain two points that are the middle of two edges of the triangle being subdivided ). I compute the points lying on the ellipse whith the following formulas: r = sqrt( x * x + y * y + z * z ); new_x = a * x / r; new_y = b * y / r; new_z = c * z / r; x,y,z being the coordinates of the middle of an edge of a triangle from the previous iteration. All the triangles have their edges correctly ordered( their normals point outside of the shape generated ) when a==b and b==c( sphere ) but some( or many ) of them don't have when an ellipsoid mesh is generated... This isn't apparently a floating point precision problem. Here's the code: main.cpp


#include <fstream>
#include <GL/gl.h>
#include <GL/glu.h>
#include <iostream>
#include <SDL/SDL.h>
#include <vector>
#include "../include/Octahedron.h"

void draw( std::vector< Vector3 > const& vertices,
           std::vector< Triple > const& triples,
           std::vector< Real > const& pos )
{
    glClear( GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT );
    glMatrixMode( GL_MODELVIEW );
    glLoadIdentity( );
    gluLookAt((double)pos[0],(double)pos[1],(double)pos[2],0,0,0,0,0,1);

    size_t n_triples = triples.size();

    glBegin( GL_TRIANGLES );

    for( size_t i = 0; i < n_triples; ++i )
    {
        glColor3ub(255,000,000);
        glVertex3d( (double)vertices[triples.a].x,
                    (double)vertices[triples.a].y,
                    (double)vertices[triples.a].z );

        glColor3ub(000,255,000);
        glVertex3d( (double)vertices[triples.b].x,
                    (double)vertices[triples.b].y,
                    (double)vertices[triples.b].z );

        glColor3ub(000,000,255);
        glVertex3d( (double)vertices[triples.c].x,
                    (double)vertices[triples.c].y,
                    (double)vertices[triples.c].z );
    }

    glEnd();

    glFlush();
    SDL_GL_SwapBuffers();
}

Real compute_area( Triple const& triple, std::vector< Vector3 > const& vertices )
{
    return 0.5*length((vertices[triple.b]-vertices[triple.a]).cross(vertices[triple.c]-vertices[triple.a]));
}

bool const correctly_oriented( Triple const& triple, std::vector< Vector3 > const& vertices )
{
    Vector3 center = vertices[triple.a];
    center += vertices[triple.b];
    center += vertices[triple.c];
    center /= ( Real )3.0;

    Vector3 d1 = vertices[triple.b] - vertices[triple.a];
    Vector3 d2 = vertices[triple.c] - vertices[triple.a];

    return ( dot( d1.cross(d2), center ) > (Real)0.0 );
}

void tesselate( Real a,
                Real b,
                Real c,
                std::vector< Triple >& triples,
                std::vector< Vector3 >& vertices )
{
    Vector3 vec;
    Real inv_r;
    int count = -1;
    size_t n_triples = triples.size();
    std::vector< Triple > temp( n_triples * 4 );

    for( size_t i = 0; i < n_triples; ++i )
    {
        //assert( compute_area( triples, vertices ) > (Real)1.0e-6 );

        vec = (vertices[triples.a]+vertices[triples.b])*(Real)0.5;
        inv_r = (Real)1.0/length(vec);
        vec.x = vec.x*a*inv_r;
        vec.y = vec.y*b*inv_r;
        vec.z = vec.z*c*inv_r;
        vertices.push_back(vec);//d -3

        vec = (vertices[triples.b]+vertices[triples.c])*(Real)0.5;
        inv_r =(Real)1.0/length(vec);
        vec.x = vec.x*a*inv_r;
        vec.y = vec.y*b*inv_r;
        vec.z = vec.z*c*inv_r;
        vertices.push_back(vec);//e -2

        vec = (vertices[triples.c]+vertices[triples.a])*(Real)0.5;
        inv_r =(Real)1.0/length(vec);
        vec.x = vec.x*a*inv_r;
        vec.y = vec.y*b*inv_r;
        vec.z = vec.z*c*inv_r;
        vertices.push_back(vec);//f -1

        size_t n_vertices = vertices.size();
        assert( n_vertices > 2 );

        size_t d = n_vertices-3;
        size_t e = n_vertices-2;
        size_t f = n_vertices-1;

        temp[++count] = Triple(triples.a,d,f);
        //assert( correctly_oriented( temp[count], vertices ) );
        //std::cout << "temp[" << count << "] = " << temp[count] << "\n";

        temp[++count] = Triple(d,triples.b,e);
        //assert( correctly_oriented( temp[count], vertices ) );
        //std::cout << "temp[" << count << "] = " << temp[count] << "\n";

        temp[++count] = Triple(d,e,f);
        //assert( correctly_oriented( temp[count], vertices ) );
        //std::cout << "temp[" << count << "] = " << temp[count] << "\n";

        temp[++count] = Triple(f,e,triples.c);
        //assert( correctly_oriented( temp[count], vertices ) );
        //std::cout << "temp[" << count << "] = " << temp[count] << "\n";
    }

    triples.swap( temp );
}

void compute_ellipsoid_mesh( Real a,
                             Real b,
                             Real c,
                             size_t n_iterations,
                             std::vector< Triple >& triples,
                             std::vector< Vector3 >& vertices )
{
    for( size_t i = 0; i < n_iterations; ++i )
    {
        tesselate( a, b, c, triples, vertices );
    }
}

int main()
{
    std::ifstream in( "/home/stanp/codeblocks/tesselation/data/param.txt" );
    assert( in );
    Real a;
    Real b;
    Real c;
    size_t n_iterations;
    in >> a;
    in >> b;
    in >> c;
    in >> n_iterations;

    Octahedron oh( a, b, c );
    std::vector< Triple > mesh_triples( oh.triples );
    std::vector< Vector3 > mesh_vertices( oh.vertices );

    compute_ellipsoid_mesh( a,
                            b,
                            c,
                            n_iterations,
                            mesh_triples,
                            mesh_vertices );

    SDL_Init(SDL_INIT_VIDEO);
    SDL_WM_SetCaption("SDL GL Application", NULL);
    SDL_SetVideoMode(640, 480, 32, SDL_OPENGL);
    glMatrixMode( GL_PROJECTION );
    glLoadIdentity( );
    gluPerspective(70,(double)640/480,1,1000);
    glEnable(GL_DEPTH_TEST);

    SDL_Event event;
    std::vector< Real > pos(3);
    in >> pos[0];
    in >> pos[1];
    in >> pos[2];

    for (;;)
    {
        while (SDL_PollEvent(&event))
        {
            switch(event.type)
            {
                case SDL_QUIT:
                    exit(0);
                    break;
            }
        }

        draw(mesh_vertices,mesh_triples,pos);
        SDL_Delay( 40 );
    }

    return 0;
}

Octahedron.h


#ifndef OCTAHEDRON_H
#define OCTAHEDRON_H

#include "prerequisites.h"
#include "Triple.h"
#include "Vector3.h"
#include <cassert>
#include <iostream>
#include <vector>

struct Octahedron
{
    Real a, b, c;
    std::vector< Vector3 > vertices;
    std::vector< Triple > triples;

    Octahedron( Real a, Real b, Real c )
        :a(a), b(b), c(c)
    {
        assert( a > 0 );
        assert( b > 0 );
        assert( c > 0 );

        vertices.resize(6);
        vertices[0] = Vector3(a,0,0);
        vertices[1] = Vector3(0,b,0);
        vertices[2] = Vector3(0,0,c);
        vertices[3] = Vector3(-a,0,0);
        vertices[4] = Vector3(0,-b,0);
        vertices[5] = Vector3(0,0,-c);

        triples.resize(8);
        triples[0] = Triple(0,1,2);
        triples[1] = Triple(1,3,2);
        triples[2] = Triple(2,3,4);
        triples[3] = Triple(0,2,4);
        triples[4] = Triple(0,5,1);
        triples[5] = Triple(1,5,3);
        triples[6] = Triple(3,5,4);
        triples[7] = Triple(0,4,5);

//        std::vector< Vector3 > centers(8);
//
//        for( size_t i = 0; i < 8; ++i )
//        {
//            centers = vertices[triples.a];
//            //std::cout << "centers[" << i << "] =\n" << centers << "\n";
//            centers += vertices[triples.b];
//            //std::cout << "centers[" << i << "] =\n" << centers << "\n";
//            centers += vertices[triples.c];
//            //std::cout << "centers[" << i << "] =\n" << centers << "\n";
//
//            centers /= ( Real )3.0;
//            //std::cout << "centers[" << i << "] =\n" << centers << "\n";
//        }
//
//        Vector3 d1;
//        Vector3 d2;
//        Vector3 v;
//        size_t temp;
//
//        for( size_t i = 0; i < 8; ++i )
//        {
//            d1 = vertices[triples.b] - vertices[triples.a];
//            d2 = vertices[triples.c] - vertices[triples.a];
//            v = d1.cross(d2);
//
//            //std::cout << "v =\n" << v << "\n";
//            //std::cout << "dot( v, centers[" << i << "] ) =\n" << dot( v, centers ) << "\n";
//
//            if( dot( v, centers ) < (Real)0.0 )
//            {
//                temp = triples.b;
//                triples.b = triples.c;
//                triples.c = temp;
//
//                //std::cout << "i = " << i << "\n";
//            }
//        }
    }
};


#endif // OCTAHEDRON_H

Triple.h


#ifndef TRIPLE_H
#define TRIPLE_H

#include <iostream>

struct Triple
{
    size_t a, b, c;

    Triple()
        :a(0), b(1), c(2)
    {
    }

    Triple( size_t a, size_t b, size_t c )
        :a(a), b(b), c(c)
    {
    }

    bool operator==( Triple const& triple )const
    {
        return ( this->a == triple.a &&
                 this->b == triple.b &&
                 this->c == triple.c );
    }

    bool operator!=( Triple const& triple )const
    {
        return !( this->operator==( triple ) );
    }
};

std::ostream& operator<<( std::ostream& out, Triple const& triple )
{
    out << triple.a << ", " << triple.b << ", " << triple.c;

    return out;
}

#endif // TRIPLE_H

Vector3.h


#ifndef VECTOR3_H
#define VECTOR3_H

#include "prerequisites.h"
#include <iostream>

    class Vector3
    {
        public:
            static Vector3 const ZERO;

            Real x;
            Real y;
            Real z;

            Vector3();
            Vector3( Real, Real, Real );
            explicit Vector3( Real* );

            Real operator [] ( size_t )const;
            Real& operator [] ( size_t );
            Vector3& operator += ( Vector3 const& );
            Vector3 const operator + ( Vector3 const& )const;
            Vector3& operator -= ( Vector3 const& );
            Vector3 const operator - ( Vector3 const& )const;
            Vector3& operator *= ( Real );
            Vector3 const operator * ( Real )const;
            Vector3& operator /= ( Real );
            Vector3 const operator / ( Real )const;
            bool const operator == ( Vector3 const& )const;
            bool const operator != ( Vector3 const& )const;

            Vector3 const cross( Vector3 const& )const;
            bool const is_equal( Vector3 const&, Real = 1.0e-6f )const;
    };

    Vector3 const abs( Vector3 const& );
    Real dot( Vector3 const&, Vector3 const& );
    Real length( Vector3 const& );
    Real squared_length( Vector3 const& );
    void normalize( Vector3& );
    Vector3 const operator* ( Real, Vector3 const& );
    std::ostream& operator<< ( std::ostream& out, Vector3 const& );

#endif // VECTOR3_H

prerequisites.h


#ifndef PREREQUISITES_H
#define PREREQUISITES_H

typedef long double Real;

#endif

Vector3.cpp


#include "../include/Vector3.h"
#include <cmath>
#include <cassert>
#include <iomanip>


    Vector3 const Vector3::ZERO( ( Real )0.0, ( Real )0.0, ( Real )0.0 );

    Vector3::Vector3()
    {
    }

    Vector3::Vector3( Real x, Real y, Real z )
        :x( x ), y( y ), z( z )
    {
    }

    Vector3::Vector3( Real* coords )
    {
        assert( coords );

        this->x = coords[0];
        this->y = coords[1];
        this->z = coords[2];
    }

    Real Vector3::operator [] ( size_t i )const
    {
        assert( i < 3 );

        return *( &this->x + i );
    }

    Real& Vector3::operator [] ( size_t i )
    {
        assert( i < 3 );

        return *( &this->x + i );
    }

    Vector3& Vector3::operator += ( Vector3 const& rhs )
    {
        this->x += rhs.x;
        this->y += rhs.y;
        this->z += rhs.z;

        return *this;
    }

    Vector3 const Vector3::operator + ( Vector3 const& rhs )const
    {
        Vector3 res( *this );

        return res += rhs;
    }

    Vector3& Vector3::operator -= ( Vector3 const& rhs )
    {
        this->x -= rhs.x;
        this->y -= rhs.y;
        this->z -= rhs.z;

        return *this;
    }

    Vector3 const Vector3::operator - ( Vector3 const& rhs )const
    {
        Vector3 res( *this );

        return res -= rhs;
    }

    Vector3& Vector3::operator *= ( Real rhs )
    {
        this->x *= rhs;
        this->y *= rhs;
        this->z *= rhs;

        return *this;
    }

    Vector3 const Vector3::operator * ( Real rhs )const
    {
        Vector3 res( *this );

        return res *= rhs;
    }

    Vector3& Vector3::operator /= ( Real rhs )
    {
        Real inv_rhs = ( Real )1.0 / rhs;
        this->x *= inv_rhs;
        this->y *= inv_rhs;
        this->z *= inv_rhs;

        return *this;
    }

    Vector3 const Vector3::operator / ( Real rhs )const
    {
        Vector3 res( *this );

        return res /= rhs;
    }

    bool const Vector3::operator == ( Vector3 const& rhs )const
    {
        return this->is_equal( rhs, ( Real )1.0e-6 );
    }

    bool const Vector3::operator != ( Vector3 const& rhs )const
    {
        return !( this->operator==( rhs ) );
    }

    Vector3 const Vector3::cross( Vector3 const& rhs )const
    {
        return Vector3( this->y * rhs.z - this->z * rhs.y,
                        this->z * rhs.x - this->x * rhs.z,
                        this->x * rhs.y - this->y * rhs.x );
    }

    Vector3 const abs( Vector3 const& vec )
    {
        return Vector3( fabs( vec.x ), fabs( vec.y ), fabs( vec.z ) );
    }

    Real dot( Vector3 const& lhs, Vector3 const& rhs )
    {
        return ( lhs.x * rhs.x + lhs.y * rhs.y + lhs.z * rhs.z );
    }

    bool const Vector3::is_equal( Vector3 const& rhs, Real tol )const
    {
        return ( std::fabs( this->x - rhs.x ) < tol &&
                 std::fabs( this->y - rhs.y ) < tol &&
                 std::fabs( this->z - rhs.z ) < tol );
    }

    Real length( Vector3 const& vec )
    {
        return std::sqrt( vec.x * vec.x + vec.y * vec.y + vec.z * vec.z );
    }

    Real squared_length( Vector3 const& vec )
    {
        return ( vec.x * vec.x + vec.y * vec.y + vec.z * vec.z );
    }

    void normalize( Vector3& vec )
    {
        Real inv_len = ( Real )1.0 / length( vec );
        vec.x *= inv_len;
        vec.y *= inv_len;
        vec.z *= inv_len;
    }

    Vector3 const operator* ( Real a, Vector3 const& vec )
    {
        Vector3 res( vec );

        return res.operator*=( a );
    }

    std::ostream& operator<< ( std::ostream& out, Vector3 const& vec )
    {
        for( size_t i = 0; i < 3; ++i )
        {
            out << std::setiosflags( std::ios::right | std::ios::fixed | std::ios::scientific );
            out << std::setw( 15 ) << std::setprecision( 6 )  << vec;
        }

        return out;
    }

param.txt

[Edited by - johnstanp on June 10, 2009 2:07:06 PM]

Ellipsoid tesselation

This topic is closed to new replies.

Popular Topics

Recommended Tutorials

Ellipsoid tesselation

This topic is closed to new replies.

Popular Topics

Recommended Tutorials

Reticulating splines