visualization.C

Go to the documentation of this file.

#ifdef HAVE_CONFIG_H
#include <tumble-conf.h>
#endif /* HAVE_CONFIG_H */

#include "visualization.h"

#include <iostream>
#include <cstdio>
#include <cstdlib>

#include <GL/gl.h>
#include <GL/glu.h>
#include <GL/freeglut_std.h>

#include "beziermesh.h"
#include "boundarymesh.h"
#include "hash_map.h"
#include "simulation.h"
#include "spline.h"
#include "timehelp.h"

using namespace std;
using namespace hashers;

Visualization::Visualization(Simulation *sim, int _screen_width, int _screen_height)
{
    if ( !sim ){
        cout<<"Visualization  --- received bad simulation!"<<endl;
        exit(-1);
    }
    initialize(sim->get_bezier_mesh(), sim->get_boundary_mesh(), _screen_width, _screen_height);
}

Visualization::Visualization(Simulation *sim)
{
    if ( !sim ){
        cout<<"Visualization  --- received bad simulation!"<<endl;
        exit(-1);
    }
    initialize(sim->get_bezier_mesh(), sim->get_boundary_mesh());
}

Visualization::Visualization(BezierMesh *_bezier_mesh, BoundaryMesh *_bdry_mesh, int _screen_width, int _screen_height)
{
    initialize(_bezier_mesh, _bdry_mesh, _screen_width, _screen_height);
}

Visualization::Visualization(BezierMesh *_bezier_mesh, BoundaryMesh *_bdry_mesh)
{
    initialize(_bezier_mesh, _bdry_mesh);
}

Visualization::Visualization(BezierMesh *_bezier_mesh, BoundaryMesh *_bdry_mesh, bool _project)
{
  initialize(_bezier_mesh, _bdry_mesh, -1, -1, true);
}

void Visualization::initialize(BezierMesh *_bezier_mesh, BoundaryMesh *_bdry_mesh, int _screen_width, int _screen_height, bool _project)
{
    if ( !_bezier_mesh || !_bdry_mesh ){
        cout<<"Visualization::initialize  --- received bad mesh!"<<endl;
        exit(-1);
    }
    bezier_mesh = _bezier_mesh;
    bdry_mesh = _bdry_mesh;

    if(_screen_width <= 0 ) _screen_width = 400;
    if(_screen_height <= 0 ) _screen_height = 400;

    if (!_project)
      {
        compute_world_size();
        view_top = world_top;
        view_bot = world_bot;
      }

    /* initialize color mapping */
    default_opengl_color_map( colormap );
    edge_draw_acc = 8;
    face_draw_acc = 12;

    init_zoom();
    resize(_screen_width, _screen_height);
}


Visualization::~Visualization()
{
    return;
}

Pic* Visualization::dump_pixels(Pic* &p)
{
    p = pic_alloc(screen_width,screen_height,3,NULL);
    if(!p) return p;
    for(int i= screen_height - 1; i >= 0; i--){
        glReadPixels(0,screen_height - i, screen_width, 1,
                     GL_RGB, GL_UNSIGNED_BYTE,
                     &p->pix[i * p->nx * p->bpp]);
    }
    return p;
}

int Visualization::map_color(int idx, const Color &color)
{
    if( !color.is_valid() ) return -1;
    colormap[idx] = color;
    return 0;
}


int Visualization::get_norms(vector<double> &min, vector<double> &max)
{
    int len = bezier_mesh->data_length();

    min.clear();
    max.clear();

    BezierMesh::Vertex_Hash_T::iterator vi;
    for(vi = bezier_mesh->get_vertices_begin();
        vi != bezier_mesh->get_vertices_end(); ++vi) {
        const LinearData& d = (*vi)->get_dp();
        for(int i = 0; i<len; i++){
            double data = d[i];
            if((signed)min.size() < i+1) min.push_back(data);
            else if(data < min[i]) min[i] = data;

            if((signed)max.size() < i+1) max.push_back(data);
            else if(data > max[i]) max[i] = data;
        }
    }
    return len;
}

int Visualization::write_jpeg(char *filename)
{
    int val;
    if(!filename) return FALSE;
    char file[MAX_FILENAME];
    if( add_filename_extension(filename, "jpg", file, MAX_FILENAME) ) return -1;

    Pic *p=NULL;
    val = jpeg_write(file, dump_pixels(p));
    pic_free(p);
    return val;
}

int Visualization::write_tiff(char *filename)
{
    int val;
    if(!filename) return FALSE;
    char file[MAX_FILENAME];
    if( add_filename_extension(filename, "tiff", file, MAX_FILENAME) ) return -1;
    Pic *p=NULL;
    val = tiff_write(file, dump_pixels(p));
    pic_free(p);
    return val;
}

int Visualization::write_ppm(char *filename)
{
    int val;
    if(!filename) return FALSE;
    char file[MAX_FILENAME];

    if( add_filename_extension(filename, "ppm", file, MAX_FILENAME) ) return -1;

    Pic *p=NULL;
    val = ppm_write(file, dump_pixels(p));
    pic_free(p);

    return val;
}

void Visualization::draw( const Color &bdry_color, const Color &edge_color){
    BezierMesh::Bezier_Edge_Hash_T::const_iterator ze;
    for(ze = bezier_mesh->get_edges_begin(); ze != bezier_mesh->get_edges_end(); ++ze) {
        draw_bezier_edge(*ze, edge_color);
    }

    BoundaryMesh::Boundary_Edge_Hash_T::const_iterator be;
    for(be = bdry_mesh->get_edges_begin(); be != bdry_mesh->get_edges_end(); ++be) {
        draw_boundary_edge(*be, bdry_color);
    }
}

void Visualization::draw_bezier( const Color &edge_color, const Color &tri_color)
{
    BezierMesh::Bezier_Face_Hash_T::const_iterator zt;
    for(zt = bezier_mesh->get_faces_begin(); zt != bezier_mesh->get_faces_end(); ++zt) {
        draw_bezier_triangle(*zt, tri_color);
    }

    BezierMesh::Bezier_Edge_Hash_T::const_iterator ze;
    for(ze = bezier_mesh->get_edges_begin(); ze != bezier_mesh->get_edges_end(); ++ze) {
        draw_bezier_edge(*ze, edge_color);
    }

}

void Visualization::draw_solid( const Color &bdry_color, const Color &edge_color){
    BezierMesh::Bezier_Face_Hash_T::const_iterator zt;
    for(zt = bezier_mesh->get_faces_begin(); zt != bezier_mesh->get_faces_end(); ++zt) {
        draw_bezier_triangle(*zt);
    }

    BezierMesh::Bezier_Edge_Hash_T::const_iterator ze;
    for(ze = bezier_mesh->get_edges_begin(); ze != bezier_mesh->get_edges_end(); ++ze) {
        draw_bezier_edge(*ze, edge_color);
    }

    BoundaryMesh::Boundary_Edge_Hash_T::const_iterator be;
    for(be = bdry_mesh->get_edges_begin(); be != bdry_mesh->get_edges_end(); ++be) {
        draw_boundary_edge(*be, bdry_color);
    }
}

void Visualization::draw_data( const Color &bdry_color, const Color &edge_color, unsigned red_pos, unsigned green_pos, unsigned blue_pos){
    BezierMesh::Bezier_Face_Hash_T::const_iterator zt;
    for(zt = bezier_mesh->get_faces_begin(); zt != bezier_mesh->get_faces_end(); ++zt) {
        draw_bezier_triangle(*zt, red_pos, green_pos, blue_pos);
    }

    BezierMesh::Bezier_Edge_Hash_T::const_iterator ze;
    for(ze = bezier_mesh->get_edges_begin(); ze != bezier_mesh->get_edges_end(); ++ze) {
        draw_bezier_edge(*ze, edge_color);
    }

    BoundaryMesh::Boundary_Edge_Hash_T::const_iterator be;
    for(be = bdry_mesh->get_edges_begin(); be != bdry_mesh->get_edges_end(); ++be) {
        draw_boundary_edge(*be, bdry_color);
    }
}

void Visualization::draw_debug( const Color &bdry_vert_color, const Color &bdry_edge_color, const Color &mesh_color ){
    BezierMesh::Bezier_Face_Hash_T::const_iterator zt;
    for(zt = bezier_mesh->get_faces_begin(); zt != bezier_mesh->get_faces_end(); ++zt) {
        draw_bezier_triangle(*zt);
    }

    BezierMesh::Bezier_Edge_Hash_T::const_iterator ze;
    for(ze = bezier_mesh->get_edges_begin(); ze != bezier_mesh->get_edges_end(); ++ze) {
        draw_bezier_edge(*ze, mesh_color);
    }

    glPointSize(3.0);
    BezierMesh::Bezier_Vertex_Hash_T::const_iterator zv;
    for(zv = bezier_mesh->get_vertices_begin(); zv != bezier_mesh->get_vertices_end(); ++zv) {
        draw_bezier_vertex(*zv, mesh_color);
    }
    glPointSize(1.0);

    BoundaryMesh::Boundary_Edge_Hash_T::const_iterator be;
    for(be = bdry_mesh->get_edges_begin(); be != bdry_mesh->get_edges_end(); ++be) {
        draw_boundary_edge(*be, bdry_edge_color);
    }

    glPointSize(3.0);
    BoundaryMesh::Boundary_Vertex_Hash_T::const_iterator bv;
    for(bv = bdry_mesh->get_vertices_begin(); bv != bdry_mesh->get_vertices_end(); ++bv) {
        draw_boundary_vertex(*bv, bdry_vert_color);
    }
    glPointSize(1.0);
}

void Visualization::draw_debug( const Color &bdry_vert_color, const Color &bdry_edge_color, const Color &mesh_color, const Color &tri_color ){
    BezierMesh::Bezier_Face_Hash_T::const_iterator zt;
    for(zt = bezier_mesh->get_faces_begin(); zt != bezier_mesh->get_faces_end(); ++zt) {
        draw_bezier_triangle(*zt, tri_color);
    }

    BezierMesh::Bezier_Edge_Hash_T::const_iterator ze;
    for(ze = bezier_mesh->get_edges_begin(); ze != bezier_mesh->get_edges_end(); ++ze) {
        draw_bezier_edge(*ze, mesh_color);
    }

    glPointSize(3.0);
    BezierMesh::Bezier_Vertex_Hash_T::const_iterator zv;
    for(zv = bezier_mesh->get_vertices_begin(); zv != bezier_mesh->get_vertices_end(); ++zv) {
        draw_bezier_vertex(*zv, mesh_color);
    }
    glPointSize(1.0);

    BoundaryMesh::Boundary_Edge_Hash_T::const_iterator be;
    for(be = bdry_mesh->get_edges_begin(); be != bdry_mesh->get_edges_end(); ++be) {
        draw_boundary_edge(*be, bdry_edge_color);
    }

    glPointSize(3.0);
    BoundaryMesh::Boundary_Vertex_Hash_T::const_iterator bv;
    for(bv = bdry_mesh->get_vertices_begin(); bv != bdry_mesh->get_vertices_end(); ++bv) {
        draw_boundary_vertex(*bv, bdry_vert_color);
    }
    glPointSize(1.0);
}

void Visualization::draw_bdry_debug( const Color &vert_color, const Color &edge_color,
                      const Color &control_color, const Color &deboor_color)
{
    glPointSize(1.0);

    BoundaryMesh::Boundary_Edge_Hash_T::const_iterator be;
    for(be = bdry_mesh->get_edges_begin(); be != bdry_mesh->get_edges_end(); ++be) {
        draw_boundary_edge_debug(*be, edge_color, control_color, deboor_color);
    }

    glPointSize(4.0);
    BoundaryMesh::Boundary_Vertex_Hash_T::const_iterator bv;
    for(bv = bdry_mesh->get_vertices_begin(); bv != bdry_mesh->get_vertices_end(); ++bv) {
        draw_boundary_vertex(*bv, vert_color);
    }
    glPointSize(1.0);
}

/* methods for tumble-view.rb */
void Visualization::draw_skeleton( const Color &bdry_color, const Color &edge_color, bool draw_edges)
{
    if(draw_edges){
        BezierMesh::Bezier_Edge_Hash_T::const_iterator ze;
        for(ze = bezier_mesh->get_edges_begin(); ze != bezier_mesh->get_edges_end(); ++ze) {
            draw_bezier_edge(*ze, edge_color);
        }
    }

    BoundaryMesh::Boundary_Edge_Hash_T::const_iterator be;
    for(be = bdry_mesh->get_edges_begin(); be != bdry_mesh->get_edges_end(); ++be) {
        draw_boundary_edge(*be, bdry_color);
    }
}

void Visualization::draw_solid( const Color &bdry_color, const Color &edge_color, bool draw_edges)
{
    BezierMesh::Bezier_Face_Hash_T::const_iterator zt;
    for(zt = bezier_mesh->get_faces_begin(); zt != bezier_mesh->get_faces_end(); ++zt) {
        draw_bezier_triangle(*zt);
    }

    if(draw_edges){
        BezierMesh::Bezier_Edge_Hash_T::const_iterator ze;
        for(ze = bezier_mesh->get_edges_begin(); ze != bezier_mesh->get_edges_end(); ++ze) {
            draw_bezier_edge(*ze, edge_color);
        }
    }

    BoundaryMesh::Boundary_Edge_Hash_T::const_iterator be;
    for(be = bdry_mesh->get_edges_begin(); be != bdry_mesh->get_edges_end(); ++be) {
        draw_boundary_edge(*be, bdry_color);
    }
}

void Visualization::draw_data_channel( const Color &bdry_color,
                                       const Color &edge_color,
                                       unsigned channel,
                                       double min,
                                       double max,
                                       bool draw_edges)
{
    BezierMesh::Bezier_Face_Hash_T::const_iterator zt;
    for(zt = bezier_mesh->get_faces_begin(); zt != bezier_mesh->get_faces_end(); ++zt) {
        draw_bezier_triangle(*zt,channel, min, max);
    }

    if(draw_edges){
        BezierMesh::Bezier_Edge_Hash_T::const_iterator ze;
        for(ze = bezier_mesh->get_edges_begin(); ze != bezier_mesh->get_edges_end(); ++ze) {
            draw_bezier_edge(*ze, edge_color);
        }
    }

    BoundaryMesh::Boundary_Edge_Hash_T::const_iterator be;
    glLineWidth(3.0);
    for(be = bdry_mesh->get_edges_begin(); be != bdry_mesh->get_edges_end(); ++be) {
        draw_boundary_edge(*be, bdry_color);
    }
    glLineWidth(1.0);
}

void Visualization::draw_point(BezierTriangle *t, double u, double v, const Color& tc, const Color &vc){
    draw_bezier_triangle(t, tc);
    Point2D pt = t->eval(u,v);
    draw_point(pt, vc, 5.0);
}

void Visualization::draw_point(BezierEdge *e, double u, const Color& ec, const Color& vc){
    draw_bezier_edge(e, ec);
    Point2D pt = e->eval(u);
    draw_point(pt, vc, 5.0);
}

void Visualization::draw_point(const Point2D &p, const Color &color, double size)
{
    glPointSize(size);
    glColor3fv(color.v);
    glBegin(GL_POINTS);
    glVertex2dv( p.coords );
    glEnd();
    glPointSize(1.0);
}

void Visualization::draw_circle(const Color &color, const Point2D &center, double radius)
{
    int i,acc=16;
    double angle=0.0, inc = 2.0*PI/(double)acc;
    Point2D p;
    glColor3fv(color.v);
    glBegin(GL_LINE_STRIP);
    for(i=0; i < acc+1; i++){
        p.assign(center.x()+cos(angle)*radius, center.y()+sin(angle)*radius);
        glVertex2dv(p.coords);
        angle+=inc;
    }
    glEnd();

}


void Visualization::draw_tuple(const BezierTuple &tup, const Color &vc, const Color &ec, const Color &fc){
    draw_bezier_triangle(tup.f, fc);
    draw_bezier_edge(tup.e, ec);
    glPointSize(5.0);
    draw_bezier_vertex(tup.v, vc);
    glPointSize(1.0);
}

void Visualization::draw_boundary_vertex(BoundaryVertex *bv, const Color &color) {
    if( !is_visible(bv) ) return;
    glColor3fv(color.v);
    glBegin(GL_POINTS);
    glVertex2dv( (bv->get_cp()).coords );
    glEnd();
}

void Visualization::draw_boundary_edge(const BoundaryEdge *be, const Color &c) {
    unsigned i,j;
    const QBSpline *spline = be->get_spline();
    glColor3fv(c.v);
    glBegin(GL_LINE_STRIP);
    for(i=0; i < ((spline->b.size()-3)/2+1); i++){
        for(j=0; j<=(unsigned)edge_draw_acc; j++){
            glVertex2dv(spline->evaluate_segment(i,double(j)/edge_draw_acc).coords);
        }
    }
    glEnd();
}

void Visualization::draw_boundary_edge_debug(const BoundaryEdge *be, const Color &edge,
                              const Color &control, const Color &deboor)
{
    unsigned i,j;
    const QBSpline *spline = be->get_spline();
    /* draw control */
    glColor3fv(control.v);
    glBegin(GL_LINE_STRIP);
    for(i=0; i < (unsigned)spline->get_num_deboor(); i++){
        glVertex2dv(spline->get_deboor(i)->coords);
    }
    glEnd();

    /* draw edge */
    glColor3fv(edge.v);
    glLineWidth(2.0);
    glBegin(GL_LINE_STRIP);
    for(i=0; i < ((spline->b.size()-3)/2+1); i++){
        for(j=0; j<=(unsigned)edge_draw_acc; j++){
            glVertex2dv(spline->evaluate_segment(i,double(j)/edge_draw_acc).coords);
        }
    }
    glEnd();
    glLineWidth(1.0);

    /* draw deboors */
    glColor3fv(deboor.v);
    glPointSize(4.0);
    glBegin(GL_POINTS);
    for(i=1; i< (unsigned)spline->get_num_deboor()-1; i++){
        glVertex2dv( spline->get_deboor(i)->coords );
    }
    glEnd();
    glPointSize(1.0);

}

void Visualization::draw_bezier_vertex(BezierVertex *zv, const Color &c){
    if(!is_visible(zv)) return;
    glColor3fv(c.v);
    glBegin(GL_POINTS);
        glVertex2dv( zv->get_cp().coords );
    glEnd();
}

void Visualization::draw_bezier_edge(BezierEdge *ze, const Color& c) {
    if( !is_visible(ze) ) return;
    glColor3fv(c.v);
    glBegin(GL_LINE_STRIP);
    for(int i=0; i <= edge_draw_acc; i++){
        glVertex2dv( ze->eval(double(i)/edge_draw_acc).coords );
    }
    glEnd();
}

void Visualization::draw_bezier_triangle(BezierTriangle *zt)
{
    draw_bezier_triangle(zt, colormap[ zt->get_color() ]);
}

void Visualization::draw_bezier_triangle(BezierTriangle *zt, const Color &c)
{
    if( !is_visible(zt) ) return;
    double d = 1.0 / face_draw_acc;
    glColor3fv(c.v);
    glBegin(GL_TRIANGLES);
    int i,j;
    Point2D p1,p2,p3;
    for(i=0; i <= face_draw_acc; i++ ) {
        for(j=0; j+i < face_draw_acc; j++ ) {
            p1 = zt->eval( i*d,j*d );
            p2 = zt->eval( i*d, (j+1)*d );
            if( i < face_draw_acc ) {
                p3 = zt->eval( (i+1)*d, j*d );
                glVertex2dv( p1.coords );
                glVertex2dv( p2.coords );
                glVertex2dv( p3.coords );
            }
            if(i>0){
                p3 = zt->eval( (i-1)*d, (j+1)*d );
                glVertex2dv( p1.coords );
                glVertex2dv( p2.coords );
                glVertex2dv( p3.coords );
            }
        }
    }
    glEnd();
}

void Visualization::draw_bezier_triangle(BezierTriangle *zt, unsigned r, unsigned g, unsigned b) {
    if( !is_visible(zt) ) return;
    double d = 1.0 / face_draw_acc;
    glBegin(GL_TRIANGLES);
    int i,j;
    Point2D p1,p2,p3;
    LinearData d1,d2,d3;
    for(i=0; i <= face_draw_acc; i++ ) {
        for(j=0; j+i < face_draw_acc; j++ ) {
            p1 = zt->eval( i*d,j*d );
            p2 = zt->eval( i*d, (j+1)*d );
            d1 = zt->dataeval(i*d, j*d);
            d2 = zt->dataeval(i*d,(j+1)*d);

            if( i < face_draw_acc ) {
                p3 = zt->eval( (i+1)*d, j*d );
                d3 = zt->dataeval ((i+1)*d, j*d );
                glColor3f(d1[r],d1[g],d1[b]);
                glVertex2dv( p1.coords );
                glColor3f(d2[r],d2[g],d2[b]);
                glVertex2dv( p2.coords );
                glColor3f(d3[r],d3[g],d3[b]);
                glVertex2dv( p3.coords );
            }
            if(i>0){
                p3 = zt->eval( (i-1)*d, (j+1)*d );
                d3 = zt->dataeval( (i-1)*d, (j+1)*d );
                glColor3f(d1[r],d1[g],d1[b]);
                glVertex2dv( p1.coords );
                glColor3f(d2[r],d2[g],d2[b]);
                glVertex2dv( p2.coords );
                glColor3f(d3[r],d3[g],d3[b]);
                glVertex2dv( p3.coords );
            }
        }
    }
    glEnd();
}

void Visualization::draw_bezier_triangle(BezierTriangle *zt, unsigned channel, double min, double max)
{
    if( !is_visible(zt) ) return;
    double d = 1.0 / face_draw_acc;
    glBegin(GL_TRIANGLES);
    int i,j;

    Point2D p1,p2,p3;
    LinearData d1,d2,d3;
    for(i=0; i <= face_draw_acc; i++ ) {
        for(j=0; j+i < face_draw_acc; j++ ) {
            p1 = zt->eval( i*d,j*d );
            p2 = zt->eval( i*d, (j+1)*d );
            d1 = zt->dataeval(i*d, j*d);
            d2 = zt->dataeval(i*d,(j+1)*d);

            if( i < face_draw_acc ) {
                p3 = zt->eval( (i+1)*d, j*d );
                d3 = zt->dataeval ((i+1)*d, j*d );

                glColor3fv( spectrum_color(d1[channel], min, max).v );
                glVertex2dv( p1.coords );

                glColor3fv( spectrum_color(d2[channel], min, max).v );
                glVertex2dv( p2.coords );

                glColor3fv( spectrum_color(d3[channel], min, max).v );
                glVertex2dv( p3.coords );
            }
            if(i>0){
                p3 = zt->eval( (i-1)*d, (j+1)*d );
                d3 = zt->dataeval( (i-1)*d, (j+1)*d );

                glColor3fv( spectrum_color(d1[channel], min, max).v );
                glVertex2dv( p1.coords );

                glColor3fv( spectrum_color(d2[channel], min, max).v );
                glVertex2dv( p2.coords );

                glColor3fv( spectrum_color(d3[channel], min, max).v );
                glVertex2dv( p3.coords );
            }
        }
    }
    glEnd();
}


void Visualization::draw_control_net(BezierTriangle *t, const Color &color)
{
    if( !is_visible(t) ) return;
    glColor3fv(color.v);

    /* draw perimiter */
    glBegin(GL_LINE_STRIP);
    glVertex2dv( t->get_cp(0).coords );
    glVertex2dv( t->get_cp(3).coords );
    glVertex2dv( t->get_cp(5).coords );
    glVertex2dv( t->get_cp(4).coords );
    glVertex2dv( t->get_cp(2).coords );
    glVertex2dv( t->get_cp(1).coords );
    glVertex2dv( t->get_cp(0).coords );
    glEnd();

    /*craw inside triangle */
    glBegin(GL_LINE_STRIP);
    glVertex2dv( t->get_cp(1).coords );
    glVertex2dv( t->get_cp(3).coords );
    glVertex2dv( t->get_cp(4).coords );
    glVertex2dv( t->get_cp(1).coords );
    glEnd();
}

void Visualization::draw_polygon(Point2D ps[], int len, const Color &color)
{
    glColor3fv(color.v);
    glBegin(GL_LINE_STRIP);
    for(int i=0; i< len; ++i) glVertex2dv( ps[i].coords );
    glVertex2dv( ps[0].coords );
    glEnd();
}

void Visualization::draw_flip_debug(BezierEdge *e, Point2D ps[], const Point2D &primary, const Point2D & secondary)
{
    draw_smooth_debug(e,ps,primary,secondary);
}

void Visualization::draw_smooth_debug(BezierEdge *e, Point2D ps[], const Point2D &primary, const Point2D & secondary)
{
    Point2D top, bot;
    BezierTriangle* t[2];
    assert(e->num_faces() == 2);
    t[0] = *e->begin_faces();
    t[1] = *++e->begin_faces();

    get_bbox(t,2,top,bot);
    zoomin(top,bot);
    draw_bezier_triangle(t[0], DIMYELLOW);
    draw_bezier_triangle(t[1], DIMCYAN);
    draw( BLUE, GREEN );
    glLineWidth(4.0);
    draw_control_net(t[0], YELLOW);
    draw_control_net(t[1], CYAN);
    glLineWidth(2.0);
    draw_polygon(ps, 6, RED);
    glLineWidth(1.0);
    draw_point(primary, MAGENTA, 7.0);
    draw_point(secondary, DARKRED, 7.0);

}

void Visualization::resize(int _screen_width, int _screen_height)
{
    screen_width = _screen_width;
    screen_height = _screen_height;
    screen_ratio = (double) screen_width / (double) screen_height;
    project();
}

void Visualization::compute_world_size( )
{
    world_top.assign(-1E6, -1E6);
    world_bot.assign(1E6,1E6);
    bdry_mesh->bbox( world_top, world_bot );
    //cout<<"World top: "<<world_top<<"World bot: "<<world_bot<<endl;
    if( (world_top.x() <= world_bot.x() ) || (world_top.y() <= world_bot.y() ) ){
        cout<<"Visualization  --- received mesh with bad width or height!"<<endl;
        cout<<"World top: "<<world_top<<" World bot: "<<world_bot<<endl;
        exit(-1);
    }

    world_width = world_top.x() - world_bot.x();
    world_height = world_top.y() - world_bot.y();

    double x_border = world_width/100.0;
    double y_border = world_height/100.0;
    world_top.coords[0] += x_border;
    world_top.coords[1] += y_border;
    world_bot.coords[0] -= x_border;
    world_bot.coords[1] -= y_border;

    world_width = world_top.x() - world_bot.x();
    world_height = world_top.y() - world_bot.y();
    world_ratio = world_width/world_height;
    world_center.coords[0] = world_bot.x() + world_width * 0.5;
    world_center.coords[1] = world_bot.y() + world_height * 0.5;
}

void Visualization::project()
{
    glViewport(0, 0, screen_width, screen_height);
    glMatrixMode( GL_PROJECTION );
    glLoadIdentity();

    compute_world_size();

    double new_x, new_y;
    Point2D image_top, image_bot, image_center;
    double image_ratio, image_width, image_height;

    /* If zoomed then adjust the zoom size to maintain same aspect ratio as the world */
    /* Then set the image varibles to be the new aspect corrected zoom, if zoomed,
     * or the normal world size if not zoomed
     */
    if(zoomed){
        new_x = zoom_height *  world_ratio * 0.5;
        image_top.assign(zoom_top.x()-new_x, zoom_top.y());
        image_bot.assign(zoom_bot.x()+new_x, zoom_bot.y());

        image_width  = image_top.x() - image_bot.x();
        image_height = image_top.y() - image_bot.y();
        image_ratio  = image_width / image_height;
        image_center.coords[0] = image_bot.x() + image_width * 0.5;
        image_center.coords[1] = image_bot.y() + image_height * 0.5;
        //assert( fabs(image_ratio - world_ratio) < 1E-4 );
    }
    else{
        image_ratio = world_ratio;
        image_width = world_width;
        image_height = world_height;
        image_top = world_top;
        image_bot = world_bot;
        image_center = world_center;
    }

    /* Now given the image to display, correct to maintain aspect ratio */
    if(image_ratio > screen_ratio){
        new_y = image_width * (1.0 / screen_ratio) * 0.5;
        view_top.assign(image_top.x(), image_center.y() + new_y);
        view_bot.assign(image_bot.x(), image_center.y() - new_y);
    }
    else{
        new_x = image_height * screen_ratio * 0.5;
        view_top.assign(image_center.x() + new_x, image_top.y());
        view_bot.assign(image_center.x() - new_x, image_bot.y());
    }

    /* With all aspect ratios maintained, and global view variables set,
     * now just project the image */
    //cout<<"View top: "<<view_top<<" View bot: "<<view_bot<<endl;
    gluOrtho2D(view_bot.x(), view_top.x(), view_bot.y(), view_top.y());
    pixel_area = fabs((view_top.x() - view_bot.x()) * (view_top.y() - view_bot.y())) / (screen_width * screen_height);
}

void Visualization::init_zoom() {
    zoomed = false;
}

void Visualization::zoomin(Point2D zt, Point2D zb) {
    zoomed=true;
    zoom_top = zt;
    zoom_bot = zb;
    zoom_height = zoom_top.y() - zoom_bot.y();
    zoom_width  = zoom_top.x() - zoom_bot.x();
    zoom_center.coords[0] = zoom_bot.x() + zoom_width * 0.5;
    zoom_center.coords[1] = zoom_bot.y() + zoom_height * 0.5;
    zoom_ratio = zoom_width/zoom_height;
    project();
}

void Visualization::zoomout(){
    zoomed=false;
    project();
}

Point2D Visualization::convert_from_screen_coords(int x, int y){
    Point2D p;
    p.coords[0] = double(x) * ((view_top.x() - view_bot.x())/(double)screen_width) + view_bot.x();
    p.coords[1] = double(screen_height - y) * ((view_top.y() - view_bot.y())/(double)screen_height) + view_bot.y();
    return p;
}

void Visualization::get_bbox(BezierTriangle* t[], int num, Point2D &top, Point2D &bot)
{
    int i,j;
    Point2D p;
    for(i=0; i<num; i++){
        for(j=0; j<6; j++){
            p = t[i]->get_cp(j);
            if(i==0 && j==0) {
                bot=p;
                top=p;
            }
            if(p.x() < bot.x()) bot.coords[0]=p.x();
            if(p.x() > top.x()) top.coords[0]=p.x();
            if(p.y() < bot.y()) bot.coords[1]=p.y();
            if(p.y() > top.y()) top.coords[1]=p.y();
        }
    }
}

bool Visualization::is_visible(const Point2D &p)
{
    return (p.x() >= view_bot.x()) && (p.x() <= view_top.x()) &&
           (p.y() >= view_bot.y()) && (p.y() <= view_top.y());
}

bool Visualization::is_visible(const BoundaryVertex *v)
{
    return is_visible(v->get_cp());
}

bool Visualization::is_visible(const BezierVertex *v)
{
    return is_visible(v->get_cp());
}

bool Visualization::is_visible(const BezierEdge *e)
{
    return is_visible(e->get_cp(0))
      || is_visible(e->get_cp(1)) || is_visible(e->get_cp(2));
}

bool Visualization::is_visible(const BezierTriangle *t)
{
    Point2D cps[6];
    int i;
    for(i=0; i<6; i++) cps[i] = t->get_cp(i);

    bool on_screen = is_visible(cps[0]) || is_visible(cps[1]) || is_visible(cps[2]) ||
                     is_visible(cps[3]) || is_visible(cps[4]) || is_visible(cps[5]);
    if(!on_screen) return false;

#ifndef DRAW_SMALL_TRIS
    /* otherwise check if we are less than a pixel */
    double minx=cps[0].x(), maxx=cps[0].x(), miny=cps[0].y(), maxy=cps[0].y();

    for(i=1; i<6; i++){
        if( cps[i].x() < minx ) minx = cps[i].x();
        if( cps[i].x() > maxx ) maxx = cps[i].x();
        if( cps[i].y() < miny ) miny = cps[i].y();
        if( cps[i].y() > maxy ) maxy = cps[i].y();
    }
    double tri_area = (maxx-minx) * (maxy-miny);
    if( tri_area < 3*pixel_area) return false;
#endif /* DRAW_SMALL_TRIS */

    return true;
}


void Visualization::start_draw()
{
  glutSwapBuffers();
  glClear(GL_COLOR_BUFFER_BIT);
  glutSwapBuffers();
  glFlush();
}

void Visualization::finish_draw()
{
  glutSwapBuffers();
  glFlush();
}

void Visualization::visual_switch(int snum, BezierTuple &tup)
{
    glClear(GL_COLOR_BUFFER_BIT);

    draw(GREEN, BLUE);
    tup = bezier_mesh->Switch(snum,tup);
    draw_tuple(tup,RED,WHITE,GRAY2);
    glFlush();
}

BezierTuple Visualization::visual_locate(Point2D pt)
{
    BezierTuple tup;
    double u,v;
    int cd;

    cd = bezier_mesh->blind_locate_point(pt, tup, u, v);
    draw_point(pt, GREEN, 5.0); /* draw the searched-for point in green */
    switch(cd) {
      case 0:
        draw_point(pt, RED, 5.0);
        break;
      case 1:
        draw_point(tup.e, u, YELLOW, RED);
        break;
      case 2:
        draw_point(tup.f, u, v, YELLOW, RED);
        break;
      default:
        cerr << "pointlocation error: returned " << cd << endl;
        tup = bezier_mesh->get_tuple();
    }
    return tup;
}

#if 0
double Simulation::get_jacobian(BezierTriangle *t, Point2D p){
    double u,v;
    if( ! t->newton_find(p,u,v)){
        cout<<"Simulation::get_jacobian - newton find failure"<<endl;
        return 0.0;
    }
    return t->jacobian(u,v);
}

void Simulation::check_for_inversions(){
    int num_inversions=0;
    glClear( GL_COLOR_BUFFER_BIT );
    glLineWidth(1.0);
    draw_debug(RED,GREEN,BLUE);
    BezierMesh::Bezier_Face_Hash_T::iterator i;
    for(i = bezier_mesh->get_faces_begin(); i != bezier_mesh->get_faces_end(); ++i){
        if( (*i)->is_inverted()) {
            (*i)->draw(RED);
            num_inversions++;
        //(*i)->print();
        }
    }
    glFlush();
    cout<<"Found "<<num_inversions<<" Inversions"<<endl;
}

void Simulation::draw_troubled_triangles(BezierTriangle *t1, BezierTriangle *t2, Point2D pt, Point2D *poly){
    Point2D top, bot, p;
    int i,j;
    BezierTriangle *t;
    BezierEdge *e;

    bot=top=*t1->cp[0];
    for(i=0; i<2; i++){
        if(i==0) t= t1;
        else t=t2;

        for(j=0; j<6; j++){
            p=*t->cp[j];
            if(p.x() < bot.x()) bot.coords[0]=p.x();
            if(p.x() > top.x()) top.coords[0]=p.x();
            if(p.y() < bot.y()) bot.coords[1]=p.y();
            if(p.y() > top.y()) top.coords[1]=p.y();
        }
    }
    double world_width = top.x() - bot.x();
    double world_height = top.y() - bot.y();
    double x_border = world_width/100.0;
    double y_border = world_height/100.0;
    top.coords[0] += x_border;
    top.coords[1] += y_border;
    bot.coords[0] -= x_border;
    bot.coords[1] -= y_border;
    world_width = top.x() - bot.x();
    world_height = top.y() - bot.y();
    double screen_ratio = (double) screen_width / (double) screen_height;
    double world_ratio=world_width/world_height;
    double y_center = bot.y()+(top.y() - bot.y()) /2.0;
    double x_center = bot.x()+(top.x() - bot.x()) /2.0;

    if(world_ratio > screen_ratio){
        double new_y=world_width*(1/screen_ratio)*0.5;
        top.assign(top.x(), y_center+new_y);
        bot.assign(bot.x(), y_center-new_y);
    }
    else{
        double new_x = world_height * screen_ratio * 0.5;
        top.assign(x_center+new_x, top.y());
        bot.assign(x_center-new_x, bot.y());
    }
    glClear( GL_COLOR_BUFFER_BIT );
    glLineWidth(1.0);
    draw_debug(RED,BLUE,GREEN);

    glLineWidth(3.0);
    for(i=0; i<2; i++){
        if(i==0) t= t1;
        else t=t2;

        for(j=0; j<3; j++){
            e=(BezierEdge*) t->edges[j];
            e->draw(RED);
        }
    }
    glColor3f(MAGENTA.R,MAGENTA.G,MAGENTA.B);
    glPointSize(6.0);
    glBegin(GL_POINTS);
    glVertex2dv( pt.coords );
    glEnd();
    glFlush();
    cout<<"Enconutered mesh geometry error!!!  Bad triangles drawn in red"<<endl;
    cout<<"press <enter> to zoom in on problem"<<endl;

    getchar();



    glViewport(0,0,screen_width,screen_height);
    glMatrixMode( GL_PROJECTION );
    glLoadIdentity();

    gluOrtho2D(bot.x(), top.x(), bot.y(), top.y());
    glClear( GL_COLOR_BUFFER_BIT );
    glLineWidth(1.0);
    draw_debug(RED,BLUE,GREEN);
    glLineWidth(3.0);
    for(i=0; i<2; i++){
        if(i==0) t= t1;
        else t=t2;

        for(j=0; j<3; j++){
            e=(BezierEdge*) t->edges[j];
            e->draw(RED);
        }
    }
    glColor3f(MAGENTA.R,MAGENTA.G,MAGENTA.B);
    glPointSize(6.0);
    glBegin(GL_POINTS);
    glVertex2dv( pt.coords );
    glEnd();

    glBegin(GL_LINE_STRIP);
    glColor3f(YELLOW.R,YELLOW.G,YELLOW.B);
    glLineWidth(2.0);
    for(i=0; i<6; i++) {
        glVertex2dv( poly[i].coords);
    }
    glVertex2dv( poly[0].coords);
    glEnd();
    glFlush();
    cout<<"Zoomed view of problem"<<endl;
    cout<<"press <enter> to continue meshing"<<endl;
    getchar();

    /* redraw old screen */
    project_to_screen(screen_width, screen_height);
    glClear( GL_COLOR_BUFFER_BIT );
    glLineWidth(1.0);
    draw_debug(RED,BLUE,GREEN);
    glFlush();
}

void Simulation::draw_troubled_interp_triangles(BezierTriangle *t1, GhostTriangle **oldts, Point2D pt){
    Point2D top, bot, p;
    int i,j;
    GhostTriangle *t;
    BezierEdge *e;

    bot=top=*t1->cp[0];
    for(i=0; i<2; i++){
        if(i==0) t = oldts[0];
        else     t = oldts[1];

        for(j=0; j<6; j++){
            p = t->cp[j];
            if(p.x() < bot.x()) bot.coords[0]=p.x();
            if(p.x() > top.x()) top.coords[0]=p.x();
            if(p.y() < bot.y()) bot.coords[1]=p.y();
            if(p.y() > top.y()) top.coords[1]=p.y();
        }
    }
    double world_width = top.x() - bot.x();
    double world_height = top.y() - bot.y();
    double x_border = world_width/100.0;
    double y_border = world_height/100.0;
    top.coords[0] += x_border;
    top.coords[1] += y_border;
    bot.coords[0] -= x_border;
    bot.coords[1] -= y_border;
    world_width = top.x() - bot.x();
    world_height = top.y() - bot.y();
    double screen_ratio = (double) screen_width / (double) screen_height;
    double world_ratio=world_width/world_height;
    double y_center = bot.y()+(top.y() - bot.y()) /2.0;
    double x_center = bot.x()+(top.x() - bot.x()) /2.0;

    if(world_ratio > screen_ratio){
        double new_y=world_width*(1/screen_ratio)*0.5;
        top.assign(top.x(), y_center+new_y);
        bot.assign(bot.x(), y_center-new_y);
    }
    else{
        double new_x = world_height * screen_ratio * 0.5;
        top.assign(x_center+new_x, top.y());
        bot.assign(x_center-new_x, bot.y());
    }
    glClear( GL_COLOR_BUFFER_BIT );
    glLineWidth(1.0);
    draw_debug(RED,BLUE,GREEN);

    /* draw old triangles */
    glLineWidth(4.0);
    glColor3f(RED.R, RED.G, RED.B);
    double u;
    for(j=0; j<2; j++){
        t = oldts[j];
        /*draw edge 1 */
        glBegin(GL_LINE_STRIP);
        for(int i=0; i <= edge_draw_acc; i++){
            u=double(i)/edge_draw_acc;
            glVertex2dv( t->eval(u,0).coords );
        }
        glEnd();

        /*draw edge 2 */
        glBegin(GL_LINE_STRIP);
        for(int i=0; i <= edge_draw_acc; i++){
            u=double(i)/edge_draw_acc;
            glVertex2dv( t->eval(0,u).coords );
        }
        glEnd();

        /*draw edge 3 */
        glBegin(GL_LINE_STRIP);
        for(int i=0; i <= edge_draw_acc; i++){
            u=double(i)/edge_draw_acc;
            glVertex2dv( t->eval(u,1-u).coords );
        }
        glEnd();
    }
    /*draw new triangle */
    glLineWidth(2.0);
    for(j=0; j<3; j++){
        e=(BezierEdge*) t1->edges[j];
        e->draw(CYAN);
    }
    glColor3f(MAGENTA.R,MAGENTA.G,MAGENTA.B);
    glPointSize(6.0);
    glBegin(GL_POINTS);
    glVertex2dv( pt.coords );
    glEnd();
    glFlush();
    cout<<"Enconutered mesh geometry error!!!  Bad triangles drawn in red"<<endl;
    cout<<"press <enter> to zoom in on problem"<<endl;

    getchar();



    glViewport(0,0,screen_width,screen_height);
    glMatrixMode( GL_PROJECTION );
    glLoadIdentity();

    gluOrtho2D(bot.x(), top.x(), bot.y(), top.y());
    glClear( GL_COLOR_BUFFER_BIT );
    glLineWidth(0.5);
    draw_debug(RED,BLUE,GREEN);

    /* draw old triangles */
    glLineWidth(5.0);
    glColor3f(RED.R, RED.G, RED.B);
    for(j=0; j<2; j++){
        t = oldts[j];
        /*draw edge 1 */
        glBegin(GL_LINE_STRIP);
        for(int i=0; i <= edge_draw_acc; i++){
            u=double(i)/edge_draw_acc;
            glVertex2dv( t->eval(u,0).coords );
        }
        glEnd();

        /*draw edge 2 */
        glBegin(GL_LINE_STRIP);
        for(int i=0; i <= edge_draw_acc; i++){
            u=double(i)/edge_draw_acc;
            glVertex2dv( t->eval(0,u).coords );
        }
        glEnd();

        /*draw edge 3 */
        glBegin(GL_LINE_STRIP);
        for(int i=0; i <= edge_draw_acc; i++){
            u=double(i)/edge_draw_acc;
            glVertex2dv( t->eval(u,1-u).coords );
        }
        glEnd();
    }
    /*draw new triangle */
    glLineWidth(1.0);
    for(j=0; j<3; j++){
        e=(BezierEdge*) t1->edges[j];
        e->draw(CYAN);
    }
    glColor3f(MAGENTA.R,MAGENTA.G,MAGENTA.B);
    glPointSize(6.0);
    glBegin(GL_POINTS);
    glVertex2dv( pt.coords );
    glEnd();
    glFlush();
    cout<<"Zoomed view of problem"<<endl;
    cout<<"This error occured while trying to reinteroplate after a flip"<<endl;
    cout<<"The red triangles are the original triangles (old triangles)"<<endl;
    cout<<"The cyan triangle is on of the new triangles created"<<endl;
    cout<<"The magenta point is the point we were sure was in the blue triangle, but were not able to lacte in either red triangle"<<endl;
    cout<<"press <enter> to continue meshing"<<endl;
    getchar();

    /* redraw old screen */
    project_to_screen(screen_width, screen_height);
    glClear( GL_COLOR_BUFFER_BIT );
    glLineWidth(1.0);
    draw(BLUE,GREEN);
    glFlush();
}






/*
void Simulation::visual_link(BezierTuple tup){
    list<BezierVertex*> link;
    bezier_mesh->get_link(tup,link);
    draw_link(link,MAGENTA,RED,GREEN);
}
*/

/*
void Simulation::draw_link(list<BezierVertex*> &link, Color startc, Color midc, Color endc){
    list<BezierVertex*>::iterator i;
    if(link.size()<2){
        cout<<"Error link too small"<<endl;
        return;
}
    glPointSize(10.0);
    i=link.begin();
    (*i)->draw(startc);
    glPointSize(4.0);
    ++i;
    for(; i != --link.end(); ++i){
        (*i)->draw(midc);
}
    glPointSize(6.0);
    (*i)->draw(endc);
    glPointSize(1.0);
}
*/
#endif


void Visualization::draw_newton(){
    double red,green;
    double u,v;
    double alpha,beta;
    Point2D pt;
    int newton_iterations=20;
    int j=0;
    unsigned num_calls=0, num_failures=0;
    timeval start, end;
    int niter;
    int edge_num;
    double elapsed;
    BezierMesh::Bezier_Face_Hash_T::iterator i;

    glBegin(GL_POINTS);
    glPointSize(1.0);
    cout<<"$$$$$$$$Newton Convergence Map$$$$$$$$"<<endl;
    cout<<"Processing";
    cout.flush();
    gettimeofday(&start,NULL);
    for(i=bezier_mesh->get_faces_begin(); i != bezier_mesh->get_faces_end(); ++i){
        for(u=0; u<1.0; u+=0.025){
            for(v=0; v+u<1.0; v+=0.025){
                pt = (*i)->eval(u,v);
                if( (*i)->newton_find(pt, alpha, beta, niter, edge_num)){
                    if(niter==newton_iterations){
                        red=1.0;
                        green=1.0;
                    } else{
                        green=1.0-((double)niter/newton_iterations);
                        red=0.0;
                    }
                } else {
                    red=1.0;
                    green=0.0;
                    num_failures++;
                }
                glColor3f(red,green,0.0);
                glVertex2dv(pt.coords);
                num_calls++;
            }
        }
        if(j++ % 25==24) {cout<<"."; cout.flush(); }
    }
    gettimeofday(&end,NULL);
    glEnd();
    cout<<endl;
    elapsed = 1000.0 * (end.tv_sec  - start.tv_sec) + (end.tv_usec - start.tv_usec)/1000.0;
    cout<<"Newton calls: "<<num_calls<<" Num failures: "<<num_failures<<endl
            <<"  total_time(sec): "<<elapsed/1000.0
            <<"   per_call_time(ms): "<<elapsed/num_calls<<endl;
}

void Visualization::draw_inverted_classes()
{
    vector<BezierTriangle*> normal, inverted, unknown;
    bezier_mesh->classify_inverted_triangles(normal, inverted, unknown);
    cout<<"Inverted Classification:"<<endl;
    cout<<"normal: "<<normal.size()<<"  inverted: "<<inverted.size()<<"  unknown: "<<unknown.size()
        <<"  errors: "<<bezier_mesh->get_num_faces() - normal.size() - inverted.size() - unknown.size()<<endl<<endl;
    vector<BezierTriangle*>::iterator i;
    for(i=normal.begin();   i!=normal.end();   ++i) draw_bezier_triangle(*i, GREEN);
    for(i=inverted.begin(); i!=inverted.end(); ++i) draw_bezier_triangle(*i, RED);
    for(i=unknown.begin();  i!=unknown.end();  ++i) draw_bezier_triangle(*i, YELLOW);
}

int Visualization::draw_smoothable_edges()
{
    int unsmoothable=0;
    for(BezierMesh::Edge_Hash_T::iterator i = bezier_mesh->get_edges_begin();
        i != bezier_mesh->get_edges_end(); ++i){
        BezierEdge *e = *i;
        if (e->is_boundary()){
            draw_bezier_edge(e, BLUE);
        } else {
            if(bezier_mesh->can_smooth(e)) draw_bezier_edge(e, GREEN);
            else{
                draw_bezier_edge(e, RED);
                unsmoothable++;
            }
        }
    }
    return unsmoothable;
}

int Visualization::draw_flippable_edges()
{
    int unflippable=0;
    for(BezierMesh::Edge_Hash_T::iterator i = bezier_mesh->get_edges_begin();
        i != bezier_mesh->get_edges_end(); ++i){
        BezierEdge *e = *i;
        if (e->is_boundary()){
            draw_bezier_edge(e, BLUE);
        } else if(bezier_mesh->should_flip(e)) {
            if(bezier_mesh->can_flip(e)) draw_bezier_edge(e, CYAN);
            else{
                draw_bezier_edge(e, RED);
                unflippable++;
            }
        } else {
            draw_bezier_edge(e, GREEN);
        }
    }
    return unflippable;
}

---