alternaria_viewer.cpp

Go to the documentation of this file.

/*
 * This work is licensed under a Creative Commons 
 * Attribution-Noncommercial-Share Alike 3.0 United States License.
 * 
 *    http://creativecommons.org/licenses/by-nc-sa/3.0/us/
 * 
 * You are free:
 * 
 *    to Share - to copy, distribute, display, and perform the work
 *    to Remix - to make derivative works
 * 
 * Under the following conditions:
 * 
 *    Attribution. You must attribute the work in the manner specified by the
 *    author or licensor (but not in any way that suggests that they endorse you
 *    or your use of the work).
 * 
 *    Noncommercial. You may not use this work for commercial purposes.
 * 
 *    Share Alike. If you alter, transform, or build upon this work, you may
 *    distribute the resulting work only under the same or similar license to
 *    this one.
 * 
 * For any reuse or distribution, you must make clear to others the license
 * terms of this work. The best way to do this is by including this header.
 * 
 * Any of the above conditions can be waived if you get permission from the
 * copyright holder.
 * 
 * Apart from the remix rights granted under this license, nothing in this
 * license impairs or restricts the author's moral rights.
 */


#include <config.h>

#include <cstdlib>
#include <cstdio>
#include <cstring>
#include <cassert>
#include <iostream>
#include <sstream>
#include <vector>

#include <inttypes.h>

#if defined ALTERNARIA_HAVE_OPENGL_FRAMEWORK
#include <OpenGL/gl.h>
#include <OpenGL/glu.h>
#elif defined ALTERNARIA_HAVE_OPENGL
#include <GL/gl.h>
#include <GL/glu.h>
#endif

#if defined ALTERNARIA_HAVE_GLUT_FRAMEWORK
#include <GLUT/glut.h>
#elif defined ALTERNARIA_HAVE_GLUT
#include <GL/glut.h>
#endif

#include <jwsc/math/constants.h>
#include <jwsc/vector/vector.h>
#include <jwsc/vector/vector_math.h>
#include <jwsc/matrix/matrix.h>
#include <jwsc/matrix/matrix_math.h>
#include <jwsc/image/image.h>
#include <jwsc/image/image_io.h>
#include <jwsc/imgblock/imgblock.h>
#include <jwsc/imgblock/imgblock_io.h>
#include <jwsc/imgblock/surface.h>

#include <jwsc++/base/option.h>
#include <jwsc++/base/exception.h>
#include <jwsc++/graphics/camera.h>

#include "imaging_model.h"
#include "psf_model.h"
#include "structure.h"
#include "spore.h"
#include "hypha.h"
#include "alternaria_model.h"


/* Anaglyph glasses filter types. */
#define REDBLUE  1
#define REDGREEN 2
#define REDCYAN  3
#define BLUERED  4
#define GREENRED 5
#define CYANRED  6

/* Program defaults. */
#define  VIEW_WIDTH      800
#define  VIEW_HEIGHT     800
#define  APERTURE        60.0f
#define  CLIP_NEAR       10.0f
#define  CLIP_FAR        3000.0f
#define  GLASSES_TYPE    REDCYAN
#define  CAPTURE_FMT     "capture-%04u.tiff"
#define  MODEL_COMPARE_FNAME  0
#define  SURF_COMPARE_FNAME   0
#define  SURF_PATCH_SIZE 2.0
#define  OBJECT_SCALE    1.0

#define  CENTROID_MIN    -5000.0f
#define  CENTROID_MAX    5000.0f
#define  LENGTH_MU       40.0f
#define  LENGTH_SIGMA    4.0f
#define  LENGTH_MIN      1.f
#define  LENGTH_MAX      1000.0f
#define  WIDTH_MU        15.0f
#define  WIDTH_SIGMA     1.0f
#define  WIDTH_MIN       1.0f
#define  WIDTH_MAX       1000.0f
#define  DWIDTH_SIGMA    0.05f
#define  DWIDTH_MIN      -100.0f
#define  DWIDTH_MAX      100.0f
#define  THETA_MU        0
#define  THETA_SIGMA     0.2f
#define  THETA_MIN       0
#define  THETA_MAX       JWSC_PI
#define  PSI_MIN         -JWSC_PI
#define  PSI_MAX         JWSC_PI
#define  OPACITY_MU      0.9f
#define  OPACITY_SIGMA   0.2f
#define  OPACITY_MIN     0.0f
#define  OPACITY_MAX     1.0f
#define  DIST_MU         0.5f
#define  DIST_SIGMA      0.1f
#define  DIST_MIN        0
#define  DIST_MAX        1.0f



using std::cout;
using std::cerr;
using std::ostringstream;
using std::sscanf;
using std::sprintf;
using std::vector;
using namespace jwsc;
using jwscxx::base::Options;
using jwscxx::base::Exception;
using jwscxx::base::Arg_error;
using jwscxx::graphics::Stereo_camera_f;



Options* opts = 0;

uint32_t view_width = VIEW_WIDTH;

uint32_t view_height = VIEW_HEIGHT;

float aperture = APERTURE;

float clip_near = CLIP_NEAR;

float clip_far = CLIP_FAR;

int glasses_type = GLASSES_TYPE;

const char* capture_fmt = CAPTURE_FMT;

const char* model_compare_fname = MODEL_COMPARE_FNAME;

const char* surface_compare_fname = SURF_COMPARE_FNAME;

float object_scale = OBJECT_SCALE;

Stereo_camera_f* camera = 0;

GLUquadric* quad_fill = 0;

GLUquadric* quad_wire = 0;

PSF_model* psf = 0;

Alternaria* alternaria = 0;

Alternaria* model_compare = 0;

list<Surface_patch_f*>* surface_compare = 0;

Imaging_scale* surface_scale = 0;

float surface_patch_size = SURF_PATCH_SIZE;

Imaging_model_density* imaging_d = 0;

PSF_model_density* psf_d = 0;

Apical_hypha_density* apical_hypha_d = 0;

Lateral_hypha_density* lateral_hypha_d = 0;

Spore_density* spore_d = 0;

Alternaria_density* alt_d = 0;

Structure* highlight = 0;

struct Alternaria_viewer_modes
{
    bool translate_camera;

    bool rotate_camera;

    bool capture_view;

    bool anaglyph;

    bool render_focal_plane;

    bool render_world_axis;

    bool focal_prp;

    bool focal_vrp;

    bool clip_near;
}
modes;

void init_modes()
{
    modes.translate_camera   = false;
    modes.rotate_camera      = true;
    modes.capture_view       = false;
    modes.anaglyph           = false;
    modes.render_focal_plane = false;
    modes.render_world_axis  = true;
    modes.focal_prp          = false;
    modes.focal_vrp          = true;
    modes.clip_near          = true;
}

void init_camera() throw (Arg_error)
{
    delete camera;

    // Generic
    camera = new Stereo_camera_f(0, 50, 300, 0, 50, 0, 0, 1, 0,
            aperture, clip_near, clip_far);

    // B
    //camera = new Stereo_camera_f(0,340,750, 0,340,450, 0, 1, 0,
    //        aperture, clip_near, clip_far);

    // A
    //camera = new Stereo_camera_f(-20,240,-490, -20,240,-190, 0, 1, 0,
}

void clean_up()
{
    if (surface_compare)
    {
        list<Surface_patch_f*>::iterator it;
        for (it = surface_compare->begin(); it != surface_compare->end(); it++)
        {
            delete *it;
        }
        delete surface_compare;
    }

    delete surface_scale;
    delete alternaria;
    delete model_compare;
    delete psf;
    delete imaging_d;
    delete psf_d;
    delete spore_d;
    delete apical_hypha_d;
    delete lateral_hypha_d;
    delete alt_d;
    delete camera;
    delete opts;
    gluDeleteQuadric(quad_fill);
    gluDeleteQuadric(quad_wire);
}

void reset()
{
    object_scale = OBJECT_SCALE;
    aperture     = APERTURE;
    highlight    = 0;
    init_camera();
    init_modes();
}

void render_world_axis()
{
    static GLfloat amb_dif[4] = {0, 0, 0, 1.0f};

    glMatrixMode(GL_MODELVIEW);

    glPushMatrix();

    float axis_length  = 100.0f;
    float arrow_radius = axis_length / 10.0f;
    float arrow_height = axis_length / 5.0f;

    /* X */
    amb_dif[0] = 0.9f;  amb_dif[1] = 1.0f;  amb_dif[2] = 0.3f;
    glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, amb_dif);
    glPushMatrix();
      glRotatef(90, 0, 1.0f, 0);
      gluCylinder(quad_fill, 1, 1, axis_length - 0.5f*arrow_height, 8, 1);
    glPopMatrix();
    glPushMatrix();
      glTranslatef(axis_length - arrow_height, 0, 0);
      glRotatef(90, 0, 1.0f, 0);
      gluCylinder(quad_fill, arrow_radius, 0, arrow_height, 16, 8);
    glPopMatrix();

    /* Y */
    amb_dif[0] = 0.9f;  amb_dif[1] = 1.0f;  amb_dif[2] = 0.3f;
    glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, amb_dif);
    glPushMatrix();
      glRotatef(-90, 1.0f, 0, 0);
      gluCylinder(quad_fill, 1, 1, axis_length - 0.5f*arrow_height, 8, 1);
    glPopMatrix();
    glPushMatrix();
      glTranslatef(0, axis_length - arrow_height, 0);
      glRotatef(-90, 1.0f, 0, 0);
      gluCylinder(quad_fill, arrow_radius, 0, arrow_height, 16, 8);
    glPopMatrix();

    /* Z */
    amb_dif[0] = 0.9f;  amb_dif[1] = 1.0f;  amb_dif[2] = 0.3f;
    glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, amb_dif);
    gluCylinder(quad_fill, 1, 1, axis_length - 0.5f*arrow_height, 8, 1);
    glPushMatrix();
      glTranslatef(0, 0, axis_length - arrow_height);
      gluCylinder(quad_fill, arrow_radius, 0, arrow_height, 16, 8);
    glPopMatrix();

    glPopMatrix();
}

void render_focal_plane()
{
    static GLfloat amb_dif[4] = {0.7f, 0.6f, 0.5f, 0.2f};
    GLboolean prev_blend;
    GLint     prev_blend_src;
    GLint     prev_blend_dst;

    if ((prev_blend = glIsEnabled(GL_BLEND)))
    {
        glGetIntegerv(GL_BLEND_SRC, &prev_blend_src);
        glGetIntegerv(GL_BLEND_DST, &prev_blend_dst);
    }
    else
    {
        glEnable(GL_BLEND);
    }
    glBlendFunc(GL_SRC_ALPHA, GL_SRC_ALPHA);

    glMatrixMode(GL_MODELVIEW);

    glPushMatrix();
    glLoadIdentity();

    glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, amb_dif);

    float x = 200.0f;
    float y = 200.0f;
    float z = -camera->get_focal_length();

    glBegin(GL_POLYGON);
    glNormal3f( 0,  0, 1);
    glVertex3f( x,  y, z);
    glVertex3f(-x,  y, z);
    glVertex3f(-x, -y, z);
    glVertex3f( x, -y, z);
    glEnd();

    glPopMatrix();

    if (prev_blend)
    {
        glBlendFunc(prev_blend_src, prev_blend_dst);
    }
    else
    {
        glDisable(GL_BLEND);
    }
}

void render_model_compare()
{
    if (!model_compare)
    {
        return;
    }

    model_compare->get_root()->recursively_draw_in_opengl(quad_wire, 
            object_scale);
}

void render_surface_compare()
{
    static GLfloat amb_dif[4] = {0, 0, 0, 1.0f};

    if (!surface_compare)
    {
        return;
    }

    glMatrixMode(GL_MODELVIEW);

    glPushMatrix();

    glDisable(GL_LIGHT0);
    glEnable(GL_LIGHT2);

    /* Render the surface patches. */
    list<Surface_patch_f*>::iterator it;
    for (it = surface_compare->begin(); it != surface_compare->end(); it++)
    {
        const Surface_patch_f* patch = *it;

        amb_dif[0] = patch->albedo*0.78f; 
        amb_dif[1] = patch->albedo*0.6f; 
        amb_dif[2] = patch->albedo*0.2f;

        glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, amb_dif);

        glNormal3f(patch->normal->elts[0], patch->normal->elts[1],
                patch->normal->elts[2]);

        glBegin(GL_POLYGON);
        for (int i = 0; i < 4; i++)
        {
            const Vector_f* point = patch->pts[ i ];
            glVertex3f(point->elts[0], point->elts[1], point->elts[2]);
        }
        glEnd();
    }

    glDisable(GL_LIGHT2);
    glEnable(GL_LIGHT0);

    glPopMatrix();
}

void render_scene()
{
    glMatrixMode(GL_MODELVIEW);

    glPushMatrix();

    float x = object_scale*alternaria->get_root()->get_centroid()->elts[0];
    float y = object_scale*alternaria->get_root()->get_centroid()->elts[1];
    float z = object_scale*alternaria->get_root()->get_centroid()->elts[2];
    glTranslatef(-x, -y, -z);

    if (highlight)
    {
        glDisable(GL_LIGHT0);
        glEnable(GL_LIGHT1);
        highlight->draw_in_opengl(quad_fill, object_scale);
        glDisable(GL_LIGHT1);
        glEnable(GL_LIGHT0);
    }

    alternaria->get_root()->recursively_draw_in_opengl(quad_fill, object_scale);

    render_model_compare();

    render_surface_compare();

    glPopMatrix();
}

void display_glut_anaglyph()
{
    /* Set the buffer for reading and writing. Using double buffering. */
    glDrawBuffer(GL_BACK);
    glReadBuffer(GL_BACK);

    /* Clear things */
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT | GL_ACCUM_BUFFER_BIT);

    /* Left camera filter. */
    glColorMask(GL_TRUE, GL_TRUE, GL_TRUE, GL_TRUE);
    switch (glasses_type) 
    {
        case REDBLUE:
        case REDGREEN:
        case REDCYAN:
            glColorMask(GL_TRUE, GL_FALSE, GL_FALSE, GL_TRUE);
            break;
        case BLUERED:
            glColorMask(GL_FALSE, GL_FALSE, GL_TRUE, GL_TRUE);
            break;
        case GREENRED:
            glColorMask(GL_FALSE, GL_TRUE, GL_FALSE, GL_TRUE);
            break;
        case CYANRED:
            glColorMask(GL_FALSE, GL_TRUE, GL_TRUE, GL_TRUE);
            break;
    } 

    /* Set the left camera projection. */
    camera->set_left_gl_projection();

    /* Set the camera model view. */
    camera->set_left_gl_modelview();

    /* Render the world axis. */
    if (modes.render_world_axis)
    {
        render_world_axis();
    }

    render_scene();

    /* Render the focal plane if in that mode. */
    if (modes.render_focal_plane)
    {
        render_focal_plane();
    }

    glFlush();

    glAccum(GL_LOAD, 0.5f);
    glClear(GL_DEPTH_BUFFER_BIT | GL_COLOR_BUFFER_BIT);

    /* Right camera filter. */
    switch (glasses_type) 
    {
        case REDBLUE:
            glColorMask(GL_FALSE, GL_FALSE, GL_TRUE, GL_TRUE);
            break;
        case REDGREEN:
            glColorMask(GL_FALSE, GL_TRUE, GL_FALSE, GL_TRUE);
            break;
        case REDCYAN:
            glColorMask(GL_FALSE, GL_TRUE, GL_TRUE, GL_TRUE);
            break;
        case BLUERED:
        case GREENRED:
        case CYANRED:
            glColorMask(GL_TRUE, GL_FALSE, GL_FALSE, GL_TRUE);
            break;
    } 

    /* Set the right camera projection. */
    camera->set_right_gl_projection();

    /* Set the camera model view. */
    camera->set_right_gl_modelview();

    /* Render the world axis. */
    if (modes.render_world_axis)
    {
        render_world_axis();
    }

    render_scene();

    /* Render the focal plane if in that mode. */
    if (modes.render_focal_plane)
    {
        render_focal_plane();
    }

    glFlush();

    glColorMask(GL_TRUE, GL_TRUE, GL_TRUE, GL_TRUE);

    glAccum(GL_ACCUM, 0.5f);
    glAccum(GL_RETURN, 2.0f);
}

void display_glut_regular()
{
    /* Set the buffer for writing. Using double buffering. */
    glDrawBuffer(GL_BACK);
    glReadBuffer(GL_BACK);

    /* Clear things */
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    /* Set the camera projection. */
    camera->set_gl_projection();

    /* Set the camera model view. */
    camera->set_gl_modelview();

    /* Render the world axis. */
    if (modes.render_world_axis)
    {
        render_world_axis();
    }

    render_scene();

    /* Render the focal plane if in that mode. */
    if (modes.render_focal_plane)
    {
        render_focal_plane();
    }

    glFlush();
}

void display_glut()
{
    static uint32_t N = 1;

    if (modes.anaglyph)
    {
        display_glut_anaglyph();
    }
    else
    {
        display_glut_regular();
    }

    glutSwapBuffers();

    if (modes.capture_view)
    {
        Stereo_camera_f::capture_gl_view(capture_fmt, N++);
        modes.capture_view = 0;
    }
}

void reshape_glut(int w, int h)
{
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
    glViewport(0, 0, (GLsizei)w, (GLsizei)h); 
    view_width = w;
    view_height = h;
}

void keyboard_glut(unsigned char key, int a, int b)
{
    float translate_delta = 10;
    float rotate_delta    = 10*JWSC_DEG_TO_RAD;
    float focal_length;
    float theta;

    Vector_f* v_tmp = 0;

    switch (key) {
        case 't':
            modes.translate_camera = 1;
            modes.rotate_camera = 0;
            break;
        case 'r':
            modes.rotate_camera = 1;
            modes.translate_camera = 0;
            break;
        case 'l':
            if (modes.translate_camera)
            {
                const Vector_f* right = camera->get_u();
                camera->translate(
                        translate_delta * right->elts[0],
                        translate_delta * right->elts[1],
                        translate_delta * right->elts[2]);
            }
            else if (modes.rotate_camera)
            {
                camera->rotate(-rotate_delta, 0, 1, 0);
            }
            glutPostRedisplay();
            break;
        case 'h':
            if (modes.translate_camera)
            {
                const Vector_f* right = camera->get_u();
                camera->translate(
                        -translate_delta * right->elts[0],
                        -translate_delta * right->elts[1],
                        -translate_delta * right->elts[2]);
            }
            else if (modes.rotate_camera)
            {
                camera->rotate(rotate_delta, 0, 1, 0);
            }
            glutPostRedisplay();
            break;
        case 'L':
            if (modes.translate_camera)
            {
                const Vector_f* right = camera->get_u();
                camera->translate(
                        translate_delta * right->elts[0],
                        translate_delta * right->elts[1],
                        translate_delta * right->elts[2]);
            }
            else if (modes.rotate_camera)
            {
                copy_vector_f(&v_tmp, camera->get_prp());
                multiply_vector_by_scalar_f(&v_tmp, v_tmp, -1);
                camera->translate(v_tmp);
                camera->rotate(rotate_delta, camera->get_v());
                multiply_vector_by_scalar_f(&v_tmp, v_tmp, -1);
                camera->translate(v_tmp);
            }
            glutPostRedisplay();
            break;
        case 'H':
            if (modes.translate_camera)
            {
                const Vector_f* right = camera->get_u();
                camera->translate(
                        -translate_delta * right->elts[0],
                        -translate_delta * right->elts[1],
                        -translate_delta * right->elts[2]);
            }
            else if (modes.rotate_camera)
            {
                copy_vector_f(&v_tmp, camera->get_prp());
                multiply_vector_by_scalar_f(&v_tmp, v_tmp, -1);
                camera->translate(v_tmp);
                camera->rotate(-rotate_delta, camera->get_v());
                multiply_vector_by_scalar_f(&v_tmp, v_tmp, -1);
                camera->translate(v_tmp);
            }
            glutPostRedisplay();
            break;
        case 'k':
            if (modes.translate_camera)
            {
                const Vector_f* up = camera->get_v();
                camera->translate(
                        translate_delta * up->elts[0],
                        translate_delta * up->elts[1],
                        translate_delta * up->elts[2]);
            }
            else if (modes.rotate_camera)
            {
                camera->rotate(rotate_delta, camera->get_u());
            }
            glutPostRedisplay();
            break;
        case 'j':
            if (modes.translate_camera)
            {
                const Vector_f* up = camera->get_v();
                camera->translate(
                        -translate_delta * up->elts[0],
                        -translate_delta * up->elts[1],
                        -translate_delta * up->elts[2]);
            }
            else if (modes.rotate_camera)
            {
                camera->rotate(-rotate_delta, camera->get_u());
            }
            glutPostRedisplay();
            break;
        case 'K':
            if (modes.translate_camera)
            {
                const Vector_f* view = camera->get_n();
                camera->translate(
                        -translate_delta * view->elts[0],
                        -translate_delta * view->elts[1],
                        -translate_delta * view->elts[2]);
            }
            else if (modes.rotate_camera)
            {
                copy_vector_f(&v_tmp, camera->get_prp());
                multiply_vector_by_scalar_f(&v_tmp, v_tmp, -1);
                camera->translate(v_tmp);
                camera->rotate(-rotate_delta, camera->get_u());
                multiply_vector_by_scalar_f(&v_tmp, v_tmp, -1);
                camera->translate(v_tmp);
            }
            glutPostRedisplay();
            break;
        case 'J':
            if (modes.translate_camera)
            {
                const Vector_f* view = camera->get_n();
                camera->translate(
                        translate_delta * view->elts[0],
                        translate_delta * view->elts[1],
                        translate_delta * view->elts[2]);
            }
            else if (modes.rotate_camera)
            {
                copy_vector_f(&v_tmp, camera->get_prp());
                multiply_vector_by_scalar_f(&v_tmp, v_tmp, -1);
                camera->translate(v_tmp);
                camera->rotate(rotate_delta, camera->get_u());
                multiply_vector_by_scalar_f(&v_tmp, v_tmp, -1);
                camera->translate(v_tmp);
            }
            glutPostRedisplay();
            break;
        case 'F':
            if (modes.focal_prp)
            {
                focal_length = camera->get_focal_length() + 20;
                camera->set_f_with_same_aperture_using_prp(focal_length);
            }
            else if (modes.focal_vrp)
            {
                focal_length = camera->get_focal_length() + 20;
                camera->set_f_with_same_aperture_using_vrp(focal_length);
            }
            glutPostRedisplay();
            break;
        case 'f':
            if (modes.focal_prp)
            {
                focal_length = camera->get_focal_length();
                if (focal_length >= 22)
                {
                    focal_length -= 20;
                    camera->set_f_with_same_aperture_using_prp(focal_length);
                }
            }
            else if (modes.focal_vrp)
            {
                focal_length = camera->get_focal_length();
                if (focal_length >= 22)
                {
                    focal_length -= 20;
                    camera->set_f_with_same_aperture_using_vrp(focal_length);
                }
            }
            glutPostRedisplay();
            break;
        case 'S':
            if (object_scale >= 2.99f)
            {
                object_scale = 3.0f;
            }
            else
            {
                object_scale += 0.01f;
            }
            glutPostRedisplay();
            break;
        case 's':
            if (object_scale > 0.01f)
            {
                object_scale -= 0.01f;
            }
            glutPostRedisplay();
            break;
        case 'A':
            aperture += JWSC_PI_8;
            camera->set_aperture(aperture);
            glutPostRedisplay();
            break;
        case 'a':
            theta = aperture - JWSC_PI_8;
            if (theta >= 0)
            {
                aperture = theta;
                camera->set_aperture(theta);
            }
            glutPostRedisplay();
            break;
        case 'P':
            if (modes.clip_near)
            {
                clip_near += 1;
            }
            else
            {
                clip_far += 1;
            }
            camera->set_clipping(clip_near, clip_far);
            glutPostRedisplay();
            break;
        case 'p':
            if (clip_near > 0 && clip_far > 0)
            {
                if (modes.clip_near)
                {
                    clip_near -= 1;
                }
                else
                {
                    clip_far -= 1;
                }
                camera->set_clipping(clip_near, clip_far);
                glutPostRedisplay();
            }
            break;
        case 'c':
            modes.capture_view = 1;
            glutPostRedisplay();
            break;
        case 'o':
            reset();
            glutPostRedisplay();
            break;
        case 'q':
            std::cerr << "VRP: " 
                      << camera->get_vrp()->elts[0] << ", "
                      << camera->get_vrp()->elts[1] << ", "
                      << camera->get_vrp()->elts[2] << '\n';
            std::cerr << "PRP: " 
                      << camera->get_prp()->elts[0] << ", "
                      << camera->get_prp()->elts[1] << ", "
                      << camera->get_prp()->elts[2] << '\n';
            //free_vector_f(v_tmp);
            //clean_up();
            exit(EXIT_SUCCESS);
    }

    free_vector_f(v_tmp);
}

void mouse_glut(int button, int state, int x, int y)
{
    Vector_f* near_pt  = 0;
    Vector_f* far_pt   = 0;
    Vector_f* isect_pt = 0;

    if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN)
    {
        float u = (float)x - (float)view_width / 2.0f;
        float v = (float)view_height / 2.0f - (float)y;
        camera->get_image_pt_on_near_clipping_plane(&near_pt, u, v);
        camera->get_image_pt_on_far_clipping_plane(&far_pt, u, v);

        // Undo the transformation of the view in render_scene().
        multiply_vector_by_scalar_f(&near_pt, near_pt, 1.0f/object_scale);
        multiply_vector_by_scalar_f(&far_pt, far_pt, 1.0f/object_scale);

        float x = alternaria->get_root()->get_centroid()->elts[0];
        float y = alternaria->get_root()->get_centroid()->elts[1];
        float z = alternaria->get_root()->get_centroid()->elts[2];
        Vector_f* V = 0;
        create_vector_f(&V, 3);
        V->elts[0] = x; V->elts[1] = y; V->elts[2] = z;
        add_vectors_f(&near_pt, near_pt, V);
        add_vectors_f(&far_pt, far_pt, V);
        free_vector_f(V);

        highlight = 0;
        float t_min = 1.0f;

        vector<Structure*>* structs = alternaria->get_root()->get_descendants();

        for (size_t i = 0; i < structs->size(); i++)
        {
            float t = structs->at(i)->get_intersection(&isect_pt, 1, near_pt, 
                    far_pt);
            if (t >= 0 && t <= t_min)
            {
                t_min = t;
                highlight = structs->at(i);
                highlight->print();
            }
        }

        delete structs;

        glutPostRedisplay();
    }

    free_vector_f(near_pt);
    free_vector_f(far_pt);
    free_vector_f(isect_pt);
}

void main_menu_glut(int id)
{
    switch (id)
    {
        case 0:
            modes.capture_view = 1;
            glutPostRedisplay();
            break;
        case 1:
            reset();
            glutPostRedisplay();
            break;
        case 2:
            exit(EXIT_SUCCESS);
    }
}

void option_menu_glut(int id)
{
    switch (id)
    {
        case 0:
            modes.anaglyph = !modes.anaglyph;
            glutPostRedisplay();
            break;
        case 1:
            modes.render_focal_plane = !modes.render_focal_plane;
            glutPostRedisplay();
            break;
        case 2:
            modes.render_world_axis = !modes.render_world_axis;
            glutPostRedisplay();
            break;
    }
}

void navigation_menu_glut(int id)
{
    switch (id)
    {
        case 0:
            modes.rotate_camera    = true;
            modes.translate_camera = false;
            break;
        case 1:
            modes.rotate_camera    = false;
            modes.translate_camera = true;
            break;
        case 2:
            modes.rotate_camera    = false;
            modes.translate_camera = false;
            break;
    }
}

void focal_length_menu_glut(int id)
{
    switch (id)
    {
        case 0:
            modes.focal_prp = true;
            modes.focal_vrp = false;
            break;
        case 1:
            modes.focal_prp = false;
            modes.focal_vrp = true;
            break;
    }
}

void clipping_menu_glut(int id)
{
    switch (id)
    {
        case 0:
            modes.clip_near = true;
            break;
        case 1:
            modes.clip_near = false;
            break;
    }
}

void init_glut(int argc, char** argv)
{
    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA | GLUT_DEPTH | GLUT_ACCUM
                                    | GLUT_ALPHA);
    glutInitWindowSize(view_width, view_height);
    glutCreateWindow("Alternaria Structure Tester");
    glutDisplayFunc(display_glut); 
    glutReshapeFunc(reshape_glut);
    glutKeyboardFunc(keyboard_glut);
    glutMouseFunc(mouse_glut);

    int navigation = glutCreateMenu(navigation_menu_glut);
    glutAddMenuEntry("Rotate Camera", 0);
    glutAddMenuEntry("Translate Camera", 1);
    glutAddMenuEntry("Modify Structure", 2);

    int focal_length = glutCreateMenu(focal_length_menu_glut);
    glutAddMenuEntry("PRP", 0);
    glutAddMenuEntry("VRP", 1);

    int clipping = glutCreateMenu(clipping_menu_glut);
    glutAddMenuEntry("Near", 0);
    glutAddMenuEntry("Far", 1);

    int options = glutCreateMenu(option_menu_glut);
    glutAddMenuEntry("Anaglyph", 0);
    glutAddMenuEntry("Focal Plane", 1);
    glutAddMenuEntry("Axis", 2);

    glutCreateMenu(main_menu_glut);
    glutAddSubMenu("Navigation", navigation);
    glutAddSubMenu("Focal Length", focal_length);
    glutAddSubMenu("Clipping", clipping);
    glutAddSubMenu("Options", options);
    glutAddMenuEntry("Capture", 1);
    glutAddMenuEntry("Reset", 2);
    glutAddMenuEntry("Quit", 3);

    glutAttachMenu(GLUT_RIGHT_BUTTON);
}

void init_gl() 
{
    glEnable(GL_DEPTH_TEST);
    glEnable(GL_NORMALIZE);
    glEnable(GL_LIGHTING);

    static GLfloat base_amb[] = {0.0f, 0.0f, 0.0f, 0.0f};
    glLightModelfv(GL_LIGHT_MODEL_AMBIENT, base_amb);

    static GLfloat ambient0  [] = {0.2f, 0.2f, 0.2f, 1.0f};
    static GLfloat diffuse0  [] = {0.8f, 0.8f, 0.8f, 1.0f};
    static GLfloat specular0 [] = {0.9f, 0.9f, 0.9f, 1.0f};
    static GLfloat direction0[] = {500.0f, 500.0f, 2000.0f, 1.0f};

    glEnable(GL_LIGHT0);
    glLightfv(GL_LIGHT0, GL_AMBIENT, ambient0);
    glLightfv(GL_LIGHT0, GL_DIFFUSE, diffuse0);
    glLightfv(GL_LIGHT0, GL_SPECULAR, specular0);
    glLightfv(GL_LIGHT0, GL_POSITION, direction0);

    static GLfloat ambient1  [] = {0.1f, 0.1f, 0.1f, 1.0f};
    static GLfloat diffuse1  [] = {0.4f, 0.4f, 0.4f, 1.0f};
    static GLfloat specular1 [] = {0.5f, 0.5f, 0.5f, 1.0f};
    static GLfloat direction1[] = {500.0f, 500.0f, 2000.0f, 1.0f};

    glDisable(GL_LIGHT1);
    glLightfv(GL_LIGHT1, GL_AMBIENT, ambient1);
    glLightfv(GL_LIGHT1, GL_DIFFUSE, diffuse1);
    glLightfv(GL_LIGHT1, GL_SPECULAR, specular1);
    glLightfv(GL_LIGHT1, GL_POSITION, direction1);

    static GLfloat ambient2  [] = {0.2f, 0.2f, 0.2f, 1.0f};
    static GLfloat diffuse2  [] = {0.9f, 0.9f, 0.9f, 1.0f};
    static GLfloat specular2 [] = {0.0f, 0.0f, 0.0f, 1.0f};
    static GLfloat direction2[] = {500.0f, 500.0f, 2000.0f, 1.0f};

    glDisable(GL_LIGHT2);
    glLightfv(GL_LIGHT2, GL_AMBIENT, ambient2);
    glLightfv(GL_LIGHT2, GL_DIFFUSE, diffuse2);
    glLightfv(GL_LIGHT2, GL_SPECULAR, specular2);
    glLightfv(GL_LIGHT2, GL_POSITION, direction2);


    if (!quad_fill)
    {
        quad_fill = gluNewQuadric();
        gluQuadricDrawStyle(quad_fill, GLU_FILL);
    }

    if (!quad_wire)
    {
        quad_wire = gluNewQuadric();
        gluQuadricDrawStyle(quad_wire, GLU_LINE);
    }

    glLineWidth(1.0);
    glPointSize(1.0);
    glClearColor(0.0, 0.0, 0.0, 0.0);
    glClearAccum(0.0, 0.0, 0.0, 0.0);

    glViewport(0, 0, view_width, view_height);
}

void print_usage(void)
{
    cerr << "usage: alternaria_viewer OPTIONS [model]\n";
    cerr << "where OPTIONS is one or more of the following:\n";
    opts->print(cerr, 30);
}

void process_help_opt(void)
{   
    print_usage();
    exit(EXIT_SUCCESS);
}

void process_version_opt(void)
{
    cout << ALTERNARIA_PACKAGE_STRING << "\n";
    exit(EXIT_SUCCESS);
}

void process_view_size_opt(const char* arg) throw (Arg_error)
{
    int32_t w, h;

    if (!arg)
    {
        throw Arg_error("Option 'view-size' requires an argument");
    }
    if (sscanf(arg, "%d,%d", &w, &h) != 2)
    {
        throw Arg_error("Option 'view-size' has format width,height");
    }
    if (w <= 0 || h <= 0)
    {
        throw Arg_error("Option 'view-size' must have width,height > 0");
    }
    view_width = (uint32_t)w;
    view_height = (uint32_t)h;
}

void process_glasses_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'glasses' requires an argument");
    }
    if (strcmp("red-blue", arg) == 0)
        glasses_type = REDBLUE;
    else if (strcmp("red-green", arg) == 0)
        glasses_type = REDGREEN;
    else if (strcmp("red-cyan", arg) == 0)
        glasses_type = REDCYAN;
    else if (strcmp("blue-red", arg) == 0)
        glasses_type = BLUERED;
    else if (strcmp("green-red", arg) == 0)
        glasses_type = GREENRED;
    else if (strcmp("cyan-red", arg) == 0)
        glasses_type = CYANRED;
    else
    {
        throw Arg_error("Invalid 'glasses' type");
    }
}

void process_capture_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'capture' requires an argument");
    }
    capture_fmt = arg;
}

void process_surf_compare_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'surf-compare' requires an argument");
    }
    surface_compare_fname = arg;
}

void process_surf_scale_opt(const char* arg) throw (Arg_error)
{
    float x, y, z;

    if (!arg)
    {
        throw Arg_error("Option 'surf-scale' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f", &x, &y, &z) != 3)
    {
        throw Arg_error("Option 'surf-scale' has format x,y,z");
    }
    try
    {
        surface_scale->set(x, y, z);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'surf-scale': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_patch_size_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'patch-size' requires an argument");
    }
    if (sscanf(arg, "%f", &surface_patch_size) != 1 || surface_patch_size <= 0)
    {
        throw Arg_error("Option 'patch-size' must be > 0");
    }
}

void process_model_compare_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'model-compare' requires an argument");
    }
    model_compare_fname = arg;
}

void init_options()
{
    using jwscxx::base::Option_no_arg;
    using jwscxx::base::Option_with_arg;

    const char* l_name;
    const char* desc;

    opts = new Options();

    l_name = "help";
    desc   = "Program usage";
    opts->add(new Option_no_arg('h', l_name, desc, process_help_opt));

    l_name = "version";
    desc   = "Program version.";
    opts->add(new Option_no_arg('v', l_name, desc, process_version_opt));

    l_name = "view-size";
    desc   = "Interactive window View size with format width,height. Default is 400,400.";
    opts->add(new Option_with_arg(0, l_name, desc, process_view_size_opt));

    l_name = "glasses";
    desc   = "Type of colored glasses for anaglyph mode in the interactive window. Possible left-right filter types include include red-blue, red-green, red-cyan, blue-red, green-red, cyan-red. Default is cyan-red.";
    opts->add(new Option_with_arg(0, l_name, desc, process_glasses_opt));

    l_name = "capture";
    desc   = "Printf-formatted file write captured view images to.";
    opts->add(new Option_with_arg(0, l_name, desc, process_capture_opt));

    l_name = "model-compare";
    desc   = "Alternaria model to compare with. Rendered as a wireframe.";
    opts->add(new Option_with_arg(0, l_name, desc, process_model_compare_opt));

    l_name = "surf-compare";
    desc   = "Alternaria surface to compare with.";
    opts->add(new Option_with_arg(0, l_name, desc, process_surf_compare_opt));

    l_name = "surf-scale";
    desc   = "Number of world units per image pixel. Has format x,y,z.";
    opts->add(new Option_with_arg(0, l_name, desc, process_surf_scale_opt));

    l_name = "patch-size";
    desc   = "Surface patch size. Must be > 0.";
    opts->add(new Option_with_arg(0, l_name, desc, process_patch_size_opt));
}


void init_densities()
{
    spore_d->set_centroid_x(CENTROID_MIN, CENTROID_MAX);
    spore_d->set_centroid_y(CENTROID_MIN, CENTROID_MAX);
    spore_d->set_centroid_z(CENTROID_MIN, CENTROID_MAX);
    spore_d->set_length(LENGTH_MU, LENGTH_SIGMA, LENGTH_MIN, LENGTH_MAX);
    spore_d->set_width(WIDTH_MU, WIDTH_SIGMA, WIDTH_MIN, WIDTH_MAX);
    spore_d->set_theta(THETA_MU, THETA_SIGMA, THETA_MIN, THETA_MAX);
    spore_d->set_psi(PSI_MIN, PSI_MAX);
    spore_d->set_opacity(OPACITY_MU, OPACITY_SIGMA, OPACITY_MIN, OPACITY_MAX);

    apical_hypha_d->set_centroid_x(CENTROID_MIN, CENTROID_MAX);
    apical_hypha_d->set_centroid_y(CENTROID_MIN, CENTROID_MAX);
    apical_hypha_d->set_centroid_z(CENTROID_MIN, CENTROID_MAX);
    apical_hypha_d->set_length(LENGTH_MU, LENGTH_SIGMA, LENGTH_MIN, LENGTH_MAX);
    apical_hypha_d->set_width(WIDTH_MU, WIDTH_SIGMA, WIDTH_MIN, WIDTH_MAX);
    apical_hypha_d->set_theta(THETA_MU, THETA_SIGMA, THETA_MIN, THETA_MAX);
    apical_hypha_d->set_psi(PSI_MIN, PSI_MAX);
    apical_hypha_d->set_opacity(OPACITY_MU, OPACITY_SIGMA, OPACITY_MIN, OPACITY_MAX);
    apical_hypha_d->set_dwidth(DWIDTH_SIGMA, DWIDTH_MIN, DWIDTH_MAX);

    lateral_hypha_d->set_centroid_x(CENTROID_MIN, CENTROID_MAX);
    lateral_hypha_d->set_centroid_y(CENTROID_MIN, CENTROID_MAX);
    lateral_hypha_d->set_centroid_z(CENTROID_MIN, CENTROID_MAX);
    lateral_hypha_d->set_length(LENGTH_MU, LENGTH_SIGMA, LENGTH_MIN, LENGTH_MAX);
    lateral_hypha_d->set_width(WIDTH_MU, WIDTH_SIGMA, WIDTH_MIN, WIDTH_MAX);
    lateral_hypha_d->set_theta(THETA_MU, THETA_SIGMA, THETA_MIN, THETA_MAX);
    lateral_hypha_d->set_psi(PSI_MIN, PSI_MAX);
    lateral_hypha_d->set_opacity(OPACITY_MU, OPACITY_SIGMA, OPACITY_MIN, OPACITY_MAX);
    lateral_hypha_d->set_lat_dist(DIST_MU, DIST_SIGMA, DIST_MIN, DIST_MAX);
}


void read_model_compare() throw (Exception)
{
    if (!model_compare_fname)
    {
        return;
    }

    assert(model_compare == 0);

    model_compare = new Alternaria(model_compare_fname, alt_d, spore_d, 
            apical_hypha_d, lateral_hypha_d, lateral_hypha_d);
}


void read_surface_compare() throw (Exception)
{
    if (!surface_compare_fname)
    {
        return;
    }

    assert(surface_compare == 0);

    Surface_point_f* points = 0;
    uint32_t num_pts;
    Error* err;
    if ((err = read_surface_points_f(&points, &num_pts, 
                    surface_compare_fname)))
    {
        ostringstream ost;
        ost << surface_compare_fname << ": " << err->msg;
        free(err);
        throw Arg_error(ost.str());
    }

    if (num_pts == 0)
    {
        ostringstream ost;
        ost << surface_compare_fname << ": doesn't contain any surface points";
        throw Arg_error(ost.str());
    }


    surface_compare = new list<Surface_patch_f*>();

    for (uint32_t i = 0; i < num_pts; i++)
    {
        Surface_patch_f* patch = 0;
        create_surface_patch_f(&patch, &(points[ i ]), 
                surface_patch_size, surface_scale->get_x(), 
                surface_scale->get_y(), surface_scale->get_z());
        surface_compare->push_back(patch);
    }

    free(points);
}




int main(int argc, char** argv)
{
    init_options();

    try
    {
        surface_scale   = new Imaging_scale();
        spore_d         = new Spore_density();
        apical_hypha_d  = new Apical_hypha_density();
        lateral_hypha_d = new Lateral_hypha_density();
        alt_d           = new Alternaria_density();

        init_densities();

        int nargs = opts->process(argc, argv); 

        if (nargs < 2)
        {
            throw Arg_error("No input model");
        }

        alternaria = new Alternaria(argv[ argc - nargs + 1], alt_d, spore_d,
                apical_hypha_d, lateral_hypha_d, lateral_hypha_d);

        read_model_compare();
        read_surface_compare();

        init_glut(argc, argv);
        init_gl();
        init_camera();
        init_modes();

        glutMainLoop();
    }
    catch (Exception e)
    {
        cerr << "alternaria_viewer: ";
        e.print();
    }

    return EXIT_FAILURE;
}