Sign in to follow this  
johnstanp

Ellipsoid tesselation

Recommended Posts

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[i].a].x,
                    (double)vertices[triples[i].a].y,
                    (double)vertices[triples[i].a].z );

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

        glColor3ub(000,000,255);
        glVertex3d( (double)vertices[triples[i].c].x,
                    (double)vertices[triples[i].c].y,
                    (double)vertices[triples[i].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[i], vertices ) > (Real)1.0e-6 );

        vec = (vertices[triples[i].a]+vertices[triples[i].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[i].b]+vertices[triples[i].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[i].c]+vertices[triples[i].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[i].a,d,f);
        //assert( correctly_oriented( temp[count], vertices ) );
        //std::cout << "temp[" << count << "] = " << temp[count] << "\n";

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

        return out;
    }






param.txt
500
400
300
5
500
500
500





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

Share this post


Link to post
Share on other sites
Since I can't delete this topic...Well, I simply tesselate a unit sphere then scale it to obtain an ellipsoid mesh. It works fine.

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this