alternaria_gen.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 <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/prob/pdf.h>
#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++/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  OBJECT_SCALE    1.0
#define  SEED            2537984
#define  IMG_IN_FNAME    0
#define  IMG_OUT_FNAME   "imaging.model"
#define  PSF_IN_FNAME    0
#define  PSF_OUT_FNAME   "psf.model"
#define  ALT_IN_FNAME    0
#define  ALT_OUT_FNAME   "alternaria.model"
#define  ALT_SCENE_FMT   0
#define  PSF_H_FMT       0
#define  NUM_SCENE_COLS  128
#define  NUM_SCENE_ROWS  128
#define  NUM_SCENE_IMGS  128
#define  CENTROID_MIN    -5000.0f
#define  CENTROID_MAX    5000.0f
#define  NUM_STRUCTS     300
#define  BASE_LEVEL      0


using std::cout;
using std::cerr;
using std::ostringstream;
using std::sscanf;
using std::vector;
using namespace jwsc;
using jwscxx::base::Options;
using jwscxx::base::Exception;
using jwscxx::base::Arg_error;
using jwscxx::base::IO_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;

float object_scale = OBJECT_SCALE;

uint32_t seed = SEED;

const char* img_in_fname = IMG_IN_FNAME;

const char* img_out_fname = IMG_OUT_FNAME;

const char* psf_in_fname = PSF_IN_FNAME;

const char* psf_out_fname = PSF_OUT_FNAME;

const char* alt_in_fname = ALT_IN_FNAME;

const char* alt_out_fname = ALT_OUT_FNAME;

const char* alt_scene_fmt = ALT_SCENE_FMT;

const char* psf_h_fmt = PSF_H_FMT;

uint32_t num_structs = NUM_STRUCTS;

uint32_t base_level = BASE_LEVEL;

Imaging_window* img_window = 0;

Imaging_scale* img_scale = 0;

Imaging_model* imaging = 0;

PSF_padding* padding = 0;

PSF_model* psf = 0;

Alternaria* alternaria = 0;

Imaging_model_density* imaging_d = 0;

PSF_model_density* psf_d = 0;

Apical_hypha_density* apical_hypha_d = 0;

Lateral_hypha_density* lateral_hypha_1_d = 0;

Lateral_hypha_density* lateral_hypha_n_d = 0;

Spore_density* spore_d = 0;

Alternaria_density* alt_d = 0;

bool interactive = 0;

Stereo_camera_f* camera = 0;

GLUquadric* quad_fill = 0;

GLUquadric* quad_wire = 0;

Structure* highlight = 0;

struct Alternaria_gen_modes
{
    bool translate_camera;

    bool rotate_camera;

    bool modify_structure;

    bool capture_view;

    bool anaglyph;

    bool render_focal_plane;

    bool render_world_axis;

    bool focal_prp;

    bool focal_vrp;
}
modes;

void init_modes()
{
    modes.translate_camera   = false;
    modes.rotate_camera      = true;
    modes.modify_structure   = false;
    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;
}

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

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

void clean_up()
{
    delete img_window;
    delete img_scale;
    delete padding;
    delete alternaria;
    delete psf;
    delete imaging;
    delete imaging_d;
    delete psf_d;
    delete spore_d;
    delete apical_hypha_d;
    delete lateral_hypha_1_d;
    delete lateral_hypha_n_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_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)
    {
        glEnable(GL_LIGHT1);
        highlight->draw_in_opengl(quad_fill, object_scale);
        glDisable(GL_LIGHT1);
    }

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

    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;
    Alternaria* alt;

    switch (key) {
        case 't':
            modes.translate_camera = true;
            modes.rotate_camera = false;
            modes.modify_structure = false;
            break;
        case 'r':
            modes.rotate_camera = true;
            modes.translate_camera = false;
            modes.modify_structure = false;
            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);
            }
            else if (modes.modify_structure && highlight)
            {
                float width = highlight->get_width();
                try 
                {
                    highlight->set_width(width + 0.5);
                }
                catch (Arg_error e)
                {
                    e.print_details();
                }
            }
            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);
            }
            else if (modes.modify_structure && highlight)
            {
                float width = highlight->get_width();
                try 
                {
                    highlight->set_width(width - 0.5);
                }
                catch (Arg_error e)
                {
                    e.print_details();
                }
            }
            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);
            }
            else if (modes.modify_structure && highlight)
            {
                float psi = highlight->get_psi() + 0.1f*JWSC_PI_2;
                if (psi > spore_d->get_psi().get_max() 
                        && (psi - JWSC_2_PI) >= spore_d->get_psi().get_min())
                {
                    psi -= JWSC_2_PI;
                }
                try 
                {
                    highlight->set_psi(psi);
                }
                catch (Arg_error e)
                {
                    e.print_details();
                }
            }
            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);
            }
            else if (modes.modify_structure && highlight)
            {
                float psi = highlight->get_psi() - 0.1f*JWSC_PI_2;
                if (psi < spore_d->get_psi().get_min() 
                        && (psi + JWSC_2_PI) <= spore_d->get_psi().get_max())
                {
                    psi += JWSC_2_PI;
                }
                try 
                {
                    highlight->set_psi(psi);
                }
                catch (Arg_error e)
                {
                    e.print_details();
                }
            }
            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());
            }
            else if (modes.modify_structure && highlight)
            {
                float length = highlight->get_length();
                try 
                {
                    highlight->set_length(length + 1);
                }
                catch (Arg_error e)
                {
                    e.print_details();
                }
            }
            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());
            }
            else if (modes.modify_structure && highlight)
            {
                float length = highlight->get_length();
                try 
                {
                    highlight->set_length(length - 1);
                }
                catch (Arg_error e)
                {
                    e.print_details();
                }
            }
            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);
            }
            else if (modes.modify_structure && highlight)
            {
                float theta = highlight->get_theta();
                try 
                {
                    highlight->set_theta(theta + 0.1f*JWSC_PI_2);
                }
                catch (Arg_error e)
                {
                    e.print_details();
                }
            }
            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);
            }
            else if (modes.modify_structure && highlight)
            {
                float theta = highlight->get_theta();
                try 
                {
                    highlight->set_theta(theta - 0.1f*JWSC_PI_2);
                }
                catch (Arg_error e)
                {
                    e.print_details();
                }
            }
            glutPostRedisplay();
            break;
        case 'v':
            if (modes.modify_structure && highlight 
                    && dynamic_cast<Apical_structure*>(highlight))
            {
                Apical_structure* apical = (Apical_structure*)highlight;
                const Vector_f* centroid = apical->get_centroid();
                float x = centroid->elts[0] + 
                          sample_gaussian_pdf_d(0, 1, -10, 10);
                float y = centroid->elts[1] +
                          sample_gaussian_pdf_d(0, 1, -10, 10);
                float z = centroid->elts[2] +
                          sample_gaussian_pdf_d(0, 1, -10, 10);
                try 
                {
                    apical->set_position(x, y, z);
                }
                catch (Arg_error e)
                {
                    e.print_details();
                }
                glutPostRedisplay();
            }
            break;
        case 'V':
            break;
        case 'b':
            if (modes.modify_structure && highlight 
                    && dynamic_cast<Lateral_structure*>(highlight))
            {
                Lateral_structure* lat = (Lateral_structure*)highlight;
                float lat_dist = lat->get_lat_dist();
                try 
                {
                    lat->set_lat_dist(lat_dist - 0.1);
                }
                catch (Arg_error e)
                {
                    Structure* parent = 0;
                    float max = lat->get_density()->get_lat_dist().get_max();

                    if (dynamic_cast<Apical_structure*>(lat->get_parent()) &&
                            (parent = lat->get_parent()->get_parent()))
                    {
                        try
                        {
                            lat->set_parent(parent, max - 0.1);
                        }
                        catch (Arg_error e)
                        {
                            e.print_details();
                        }
                    }
                    else
                    {
                        e.print_details();
                    }
                }
                glutPostRedisplay();
            }
            break;
        case 'B':
            if (modes.modify_structure && highlight 
                    && dynamic_cast<Lateral_structure*>(highlight))
            {
                Lateral_structure* lat = (Lateral_structure*)highlight;
                float lat_dist = lat->get_lat_dist();
                try 
                {
                    lat->set_lat_dist(lat_dist + 0.1);
                }
                catch (Arg_error e)
                {
                    Apical_structure* apical = 0;
                    float min = lat->get_density()->get_lat_dist().get_min();

                    if (lat->get_parent() && 
                            (apical = lat->get_parent()->get_apical()))
                    {
                        try
                        {
                            lat->set_parent(apical, min + 0.1);
                        }
                        catch (Arg_error e)
                        {
                            e.print_details();
                        }
                    }
                    else
                    {
                        e.print_details();
                    }
                }
                glutPostRedisplay();
            }
            break;
        case 'm':
            if (modes.modify_structure && highlight 
                    && dynamic_cast<Apical_hypha*>(highlight))
            {
                Apical_hypha* rvals_1 = 0;
                Apical_hypha* rvals_2 = 0;

                try
                {
                    float width_1 = highlight->get_width() + 
                        apical_hypha_d->sample_dwidth();
                    float width_2 = highlight->get_width() + 
                        apical_hypha_d->sample_dwidth();
                    rvals_1 = apical_hypha_d->sample();
                    rvals_1->set_width(width_1);
                    rvals_2 = apical_hypha_d->sample();
                    rvals_2->set_width(width_2);
                    alternaria->split_apical_hypha_into_apical(
                            (Apical_hypha*)highlight, rvals_1, rvals_2);
                    highlight = 0;
                    delete rvals_1;
                    delete rvals_2;
                }
                catch (Arg_error e)
                {
                    delete rvals_1;
                    delete rvals_2;
                    e.print_details();
                }
            }
            else if (modes.modify_structure && highlight 
                    && dynamic_cast<Spore*>(highlight))
            {
                Spore* rvals_1 = spore_d->sample();
                Spore* rvals_2 = spore_d->sample();

                float width = 0.65f*highlight->get_width();
                float sigma = spore_d->get_width().get_sigma();
                float min   = -6.0f*sigma;
                float max   = 6.0f*sigma;

                try
                {
                    width += sample_gaussian_pdf_d(0, sigma, min, max);
                    rvals_1->set_width(width);
                    width += sample_gaussian_pdf_d(0, sigma, min, max);
                    rvals_1->set_width(width);
                    alternaria->split_spore_into_spore(
                            (Spore*)highlight, rvals_1, rvals_2);
                    highlight = 0;
                    delete rvals_1;
                    delete rvals_2;
                }
                catch (Arg_error e)
                {
                    delete rvals_1;
                    delete rvals_2;
                    e.print_details();
                }
            }
            glutPostRedisplay();
            break;
        case 'M':
            if (modes.modify_structure && highlight 
                    && dynamic_cast<Apical_hypha*>(highlight))
            {
                Apical_hypha* rvals = 0;

                try
                {
                    float width = highlight->get_width() + 
                        apical_hypha_d->sample_dwidth();
                    rvals = apical_hypha_d->sample();
                    rvals->set_width(width);
                    alternaria->merge_apical_hypha_with_apical(
                            (Apical_hypha*)highlight, rvals);
                    highlight = 0;
                    delete rvals;
                }
                catch (Arg_error e)
                {
                    delete rvals;
                    e.print_details();
                }
            }
            else if (modes.modify_structure && highlight 
                    && highlight->get_apical()
                    && dynamic_cast<Spore*>(highlight))
            {
                Spore* rvals = spore_d->sample();

                float width = 0.65f*highlight->get_width() + 
                              0.65f*highlight->get_apical()->get_width();
                float sigma = spore_d->get_width().get_sigma();
                float min   = -6.0f*sigma;
                float max   = 6.0f*sigma;

                try
                {
                    width += sample_gaussian_pdf_d(0, sigma, min, max);
                    rvals->set_width(width);
                    alternaria->merge_spore_with_spore((Spore*)highlight,rvals);
                    highlight = 0;
                    delete rvals;
                }
                catch (Arg_error e)
                {
                    delete rvals;
                    e.print_details();
                }
            }
            glutPostRedisplay();
            break;
        case 'n':
            if (modes.modify_structure && highlight 
                    && dynamic_cast<Apical_hypha*>(highlight))
            {
                Apical_hypha* rvals_1 = 0;
                Lateral_hypha* rvals_2 = 0;

                try
                {
                    float width = highlight->get_width() + 
                        apical_hypha_d->sample_dwidth();
                    rvals_1 = apical_hypha_d->sample();
                    rvals_1->set_width(width);
                    if (highlight->get_level() == 0)
                    {
                        rvals_2 = lateral_hypha_1_d->sample();
                    }
                    else
                    {
                        rvals_2 = lateral_hypha_n_d->sample();
                    }
                    alternaria->split_apical_hypha_into_lateral(
                            (Apical_hypha*)highlight, rvals_1, rvals_2);
                    delete rvals_1;
                    delete rvals_2;
                }
                catch (Arg_error e)
                {
                    delete rvals_1;
                    delete rvals_2;
                    e.print_details();
                }
                glutPostRedisplay();
            }
            break;
        case 'N':
            if (modes.modify_structure && highlight && highlight->get_parent()
                    && dynamic_cast<Apical_hypha*>(highlight->get_parent())
                    && dynamic_cast<Lateral_hypha*>(highlight))
            {
                Apical_hypha* rvals = 0;

                try
                {
                    float width = highlight->get_width() + 
                        apical_hypha_d->sample_dwidth();
                    rvals = apical_hypha_d->sample();
                    rvals->set_width(width);
                    alternaria->merge_apical_hypha_with_lateral(
                            (Apical_hypha*)highlight->get_parent(), 
                            (Lateral_hypha*)highlight, rvals);
                    highlight = 0;
                    delete rvals;
                }
                catch (Arg_error e)
                {
                    delete rvals;
                    e.print_details();
                }
                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 'c':
            modes.capture_view = 1;
            glutPostRedisplay();
            break;
        case 'C':
            alt = alternaria->clone();
            delete alternaria;
            alternaria = alt;
            highlight = 0;
            glutPostRedisplay();
            break;
        case 'o':
            reset();
            glutPostRedisplay();
            break;
        case 'w':
            imaging->print(img_out_fname);
            psf->print(psf_out_fname);
            alternaria->print(alt_out_fname);
            break;
        case 'W':
            if (psf_h_fmt)
            {
                psf->write_h(psf_h_fmt);
            }
            if (alt_scene_fmt)
            {
                Alternaria_model* model = new Alternaria_model(*alternaria, 
                        psf, imaging);
                model->write_scene(alt_scene_fmt, psf);
                delete model;
            }
            break;
        case 'q':
            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);
            }
        }

        delete structs;

        if (highlight)
        {
            highlight->print();
        }

        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:
            clean_up();
            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;
            modes.modify_structure = false;
            break;
        case 1:
            modes.rotate_camera    = false;
            modes.translate_camera = true;
            modes.modify_structure = false;
            break;
        case 2:
            modes.rotate_camera    = false;
            modes.translate_camera = false;
            modes.modify_structure = true;
            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 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 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("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);

    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-gen OPTIONS\n";
    cerr << "where OPTIONS is one or more of the following:\n";
    opts->print(cerr, 33);
}

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

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

void process_options_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'options' requires an argument");
    }
    try
    {
        opts->process(arg);
    }
    catch (IO_error e)
    {
        throw Arg_error(e.get_msg());
    }
}

void process_seed_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'seed' requires an argument");
    }
    if (sscanf(arg, "%u", &seed) != 1)
    {
        throw Arg_error("Option 'seed' must be a non-negative integer");
    }
}

void process_apical_hypha_length_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'apical-hypha-length' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'apical-hypha-length' has format mu,sigma,min,max");
    }
    try
    {
        apical_hypha_d->set_length(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'apical-hypha-length': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_lateral_hypha_length_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'lateral-hypha-length' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'lateral-hypha-length' has format mu,sigma,min,max");
    }
    try
    {
        lateral_hypha_1_d->set_length(mu, sigma, min, max);
        lateral_hypha_n_d->set_length(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'lateral-hypha-length': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_apical_hypha_width_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'apical-hypha-width' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'apical-hypha-width' has format mu,sigma,min,max");
    }
    try
    {
        apical_hypha_d->set_width(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'apical-hypha-width': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_lateral_hypha_width_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'lateral-hypha-width' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'lateral-hypha-width' has format mu,sigma,min,max");
    }
    try
    {
        lateral_hypha_1_d->set_width(mu, sigma, min, max);
        lateral_hypha_n_d->set_width(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'lateral-hypha-width': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_apical_hypha_theta_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'apical-hypha-theta' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'apical-hypha-theta' has format mu,sigma,min,max");
    }
    try
    {
        apical_hypha_d->set_theta(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'apical-hypha-theta': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_lateral_hypha_1_theta_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'lateral-hypha-1-theta' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'lateral-hypha-1-theta' has format mu,sigma,min,max");
    }
    try
    {
        lateral_hypha_1_d->set_theta(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'lateral-hypha-1-theta': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_lateral_hypha_n_theta_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'lateral-hypha-n-theta' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'lateral-hypha-n-theta' has format mu,sigma,min,max");
    }
    try
    {
        lateral_hypha_n_d->set_theta(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'lateral-hypha-n-theta': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_hypha_psi_opt(const char* arg) throw (Arg_error)
{
    float min, max;

    if (!arg)
    {
        throw Arg_error("Option 'hypha-psi' requires an argument");
    }
    if (sscanf(arg, "%f,%f", &min, &max) != 2)
    {
        throw Arg_error("Option 'hypha-psi' has format min,max");
    }
    try
    {
        apical_hypha_d->set_psi(min, max);
        lateral_hypha_1_d->set_psi(min, max);
        lateral_hypha_n_d->set_psi(min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'hypha-psi': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_hypha_opacity_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'hypha-opacity' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'hypha-opacity' has format mu,sigma,min,max");
    }
    try
    {
        apical_hypha_d->set_opacity(mu, sigma, min, max);
        lateral_hypha_1_d->set_opacity(mu, sigma, min, max);
        lateral_hypha_n_d->set_opacity(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'hypha-opacity': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_lateral_hypha_dist_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'lateral-hypha-dist' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'lateral-hypha-dist' has format mu,sigma,min,max");
    }
    try
    {
        lateral_hypha_1_d->set_lat_dist(mu, sigma, min, max);
        lateral_hypha_n_d->set_lat_dist(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'lateral-hypha-dist': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_spore_length_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'spore-length' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'spore-length' has format mu,sigma,min,max");
    }
    try
    {
        spore_d->set_length(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'spore-length': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_spore_width_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'spore-width' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'spore-width' has format mu,sigma,min,max");
    }
    try
    {
        spore_d->set_width(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'spore-width': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_spore_theta_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'spore-theta' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'spore-theta' has format mu,sigma,min,max");
    }
    try
    {
        spore_d->set_theta(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'spore-theta': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_spore_psi_opt(const char* arg) throw (Arg_error)
{
    float min, max;

    if (!arg)
    {
        throw Arg_error("Option 'spore-psi' requires an argument");
    }
    if (sscanf(arg, "%f,%f", &min, &max) != 2)
    {
        throw Arg_error("Option 'spore-psi' has format min,max");
    }
    try
    {
        spore_d->set_psi(min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'spore-psi': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_spore_opacity_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'spore-opacity' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'spore-opacity' has format mu,sigma,min,max");
    }
    try
    {
        spore_d->set_opacity(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'spore-opacity': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_psf_alpha_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'psf-alpha' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'psf-alpha' has format mu,sigma,min,max");
    }
    try
    {
        psf_d->set_alpha(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'psf-alpha': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_psf_beta_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'psf-beta' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'psf-beta' has format mu,sigma,min,max");
    }
    try
    {
        psf_d->set_beta(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'psf-beta': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_psf_gamma_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'psf-gamma' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'psf-gamma' has format mu,sigma,min,max");
    }
    try
    {
        psf_d->set_gamma(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'psf-gamma': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_img_bg_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'img-bg' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'img-bg' has format mu,sigma,min,max");
    }
    try
    {
        imaging_d->set_bg(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'img-bg': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_levels_p_opt(const char* arg) throw (Arg_error)
{
    float p;

    if (!arg)
    {
        throw Arg_error("Option 'levels-p' requires an argument");
    }
    if (sscanf(arg, "%f", &p) != 1)
    {
        throw Arg_error("Option 'levels-p' must be in a number in (0,1)");
    }
    try
    {
        alt_d->set_num_levels_p(p);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'levels-p': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_structs_p_opt(const char* arg) throw (Arg_error)
{
    float spore_p, apical_hypha_p, lateral_hypha_1_p, lateral_hypha_n_p;

    if (!arg)
    {
        throw Arg_error("Option 'structs-p' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &spore_p, &apical_hypha_p, &lateral_hypha_1_p, &lateral_hypha_n_p) != 4)
    {
        throw Arg_error("Option 'structs-p' has format spore,apical-hypha,lateral-hypha-1,lateral-hypha-n");
    }
    try
    {
        alt_d->set_structure_p(spore_p, apical_hypha_p, lateral_hypha_1_p,
                lateral_hypha_n_p);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'structs-p': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_img_in_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'img-in' requires an argument");
    }
    img_in_fname = arg;
}

void process_img_out_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'img-out' requires an argument");
    }
    img_out_fname = arg;
}

void process_psf_in_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'psf-in' requires an argument");
    }
    psf_in_fname = arg;
}

void process_psf_out_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'psf-out' requires an argument");
    }
    psf_out_fname = arg;
}

void process_alt_in_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'alt-in' requires an argument");
    }
    alt_in_fname = arg;
}

void process_alt_out_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'alt-out' requires an argument");
    }
    alt_out_fname = arg;
}

void process_scene_fmt_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'scene-fmt' requires an argument");
    }
    alt_scene_fmt = arg;
}

void process_scene_size_opt(const char* arg) throw (Arg_error)
{
    uint32_t ncols, nrows, nimgs;

    if (!arg)
    {
        throw Arg_error("Option 'scene-size' requires an argument");
    }
    if (sscanf(arg, "%u,%u,%u", &ncols, &nrows, &nimgs) != 3)
    {
        throw Arg_error("Option 'scene-size' has format ncols,nrows,nimgs");
    }
    try
    {
        img_window->set(0, 0, 0, ncols, nrows, nimgs);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'scene-size': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

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

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

void process_psf_fmt_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'psf-fmt' requires an argument");
    }
    psf_h_fmt = arg;
}

void process_num_structs_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'num-structs' requires an argument");
    }
    if (sscanf(arg, "%u", &num_structs) != 1)
    {
        throw Arg_error("Option 'num-structs' must be an integer");
    }
}

void process_interactive_opt() throw (Arg_error)
{
    interactive = true;
}

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_base_level_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'base-level' requires an argument");
    }
    if (sscanf(arg, "%u", &base_level) != 1)
    {
        throw Arg_error("Option 'base-level' must be >= 0");
    }
}

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 = "options";
    desc   = "Process a file of program options.";
    opts->add(new Option_with_arg(0, l_name, desc, process_options_opt));

    l_name = "interactive";
    desc   = "Show the generated Alternaria in an interactive GLUT/GL window.";
    opts->add(new Option_no_arg('i', l_name, desc, process_interactive_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 = "seed";
    desc   = "Random seed to use in srand(). Must be a non-negative integer.";
    opts->add(new Option_with_arg(0, l_name, desc, process_seed_opt));

    l_name = "img-in";
    desc   = "Imaging input file to read from.";
    opts->add(new Option_with_arg(0, l_name, desc, process_img_in_opt));

    l_name = "img-out";
    desc   = "Imaging output file to write to.";
    opts->add(new Option_with_arg(0, l_name, desc, process_img_out_opt));

    l_name = "psf-in";
    desc   = "PSF input file to read from.";
    opts->add(new Option_with_arg(0, l_name, desc, process_psf_in_opt));

    l_name = "psf-out";
    desc   = "PSF output file to write to.";
    opts->add(new Option_with_arg(0, l_name, desc, process_psf_out_opt));

    l_name = "alt-in";
    desc   = "Alternaria input file to read from.";
    opts->add(new Option_with_arg(0, l_name, desc, process_alt_in_opt));

    l_name = "alt-out";
    desc   = "Alternaria output file to write to.";
    opts->add(new Option_with_arg(0, l_name, desc, process_alt_out_opt));

    l_name = "scene-fmt";
    desc   = "Printf-formatted file name for the optically sectioned scene. Must contain one integer conversion";
    opts->add(new Option_with_arg(0, l_name, desc, process_scene_fmt_opt));

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

    l_name = "scene-size";
    desc   = "Width, height, and number of images to write the optically sectioned scene to. Has format ncols,nrows,nimgs.";
    opts->add(new Option_with_arg(0, l_name, desc, process_scene_size_opt));

    l_name = "psf-fmt";
    desc   = "Printf-formatted file name for the optically sectioned PSF. Must contain one integer conversion.";
    opts->add(new Option_with_arg(0, l_name, desc, process_psf_fmt_opt));

    l_name = "num-structs";
    desc   = "Number of structure elements to generate in the Alternaria. Must be an integer.";
    opts->add(new Option_with_arg(0, l_name, desc, process_num_structs_opt));

    l_name = "apical-hypha-length";
    desc   = "Length of an apical hypha. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, 
                process_apical_hypha_length_opt));

    l_name = "lateral-hypha-length";
    desc   = "Length of a lateral hypha. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, 
                process_lateral_hypha_length_opt));

    l_name = "apical-hypha-width";
    desc   = "Width of an apical hypha. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, 
                process_apical_hypha_width_opt));

    l_name = "lateral-hypha-width";
    desc   = "Width of a lateral hypha. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, 
                process_lateral_hypha_width_opt));

    l_name = "apical-hypha-theta";
    desc   = "Theta Euler angle of an apical hypha. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, 
                process_apical_hypha_theta_opt));

    l_name = "lateral-hypha-1-theta";
    desc   = "Theta Euler angle of a lateral hypha at level 1. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, 
                process_lateral_hypha_1_theta_opt));

    l_name = "lateral-hypha-n-theta";
    desc   = "Theta Euler angle of a lateral hypha at level > 1. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, 
                process_lateral_hypha_n_theta_opt));

    l_name = "hypha-psi";
    desc   = "Psi Euler angle of a hypha. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, process_hypha_psi_opt));

    l_name = "hypha-opacity";
    desc   = "Opacity of a hypha. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, process_hypha_opacity_opt));

    l_name = "lateral-hypha-dist";
    desc   = "Distance of a lateral hypha from its parent. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, 
                process_lateral_hypha_dist_opt));

    l_name = "spore-width";
    desc   = "Width of a spore. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, process_spore_width_opt));

    l_name = "spore-length";
    desc   = "Length of a spore. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, process_spore_length_opt));

    l_name = "spore-theta";
    desc   = "Theta Euler angle of a spore. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, process_spore_theta_opt));

    l_name = "spore-psi";
    desc   = "Psi Euler angle of a spore. Has argument format min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, process_spore_psi_opt));

    l_name = "spore-opacity";
    desc   = "Opacity of a spore. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, process_spore_opacity_opt));

    l_name = "psf-alpha";
    desc   = "PSF alpha parameter. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, process_psf_alpha_opt));

    l_name = "psf-beta";
    desc   = "PSF beta parameter. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, process_psf_beta_opt));

    l_name = "psf-gamma";
    desc   = "PSF gamma parameter. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, process_psf_gamma_opt));

    l_name = "img-bg";
    desc   = "Image background parameter. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, process_img_bg_opt));

    l_name = "levels-p";
    desc   = "Geometric distribution parameter over growth levels. Must be in (0,1)";
    opts->add(new Option_with_arg(0, l_name, desc, process_levels_p_opt));

    l_name = "structs-p";
    desc   = "Distribution over structure type. Has argument format spore,apical-hypha,lateral-hypha. All must sum to 1.";
    opts->add(new Option_with_arg(0, l_name, desc, process_structs_p_opt));

    l_name = "base-level";
    desc   = "Base level for the root Structure. Must be >= 0.";
    opts->add(new Option_with_arg(0, l_name, desc, process_base_level_opt));
}


void* operator new (size_t s) throw (std::bad_alloc)
{
    return malloc(s);
}


void operator delete (void* ptr) throw ()
{
    if (ptr)
        free(ptr);
}




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

    try
    {
        img_window        = new Imaging_window();
        img_scale         = new Imaging_scale();
        padding           = new PSF_padding();
        imaging_d         = new Imaging_model_density();
        psf_d             = new PSF_model_density();
        spore_d           = new Spore_density();
        apical_hypha_d    = new Apical_hypha_density();
        lateral_hypha_1_d = new Lateral_hypha_density();
        lateral_hypha_n_d = new Lateral_hypha_density();
        alt_d             = new Alternaria_density();

        img_window->set(0, 0, 0, NUM_SCENE_COLS, NUM_SCENE_ROWS, 
                NUM_SCENE_IMGS);

        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);

        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);

        lateral_hypha_1_d->set_centroid_x(CENTROID_MIN, CENTROID_MAX);
        lateral_hypha_1_d->set_centroid_y(CENTROID_MIN, CENTROID_MAX);
        lateral_hypha_1_d->set_centroid_z(CENTROID_MIN, CENTROID_MAX);

        lateral_hypha_n_d->set_centroid_x(CENTROID_MIN, CENTROID_MAX);
        lateral_hypha_n_d->set_centroid_y(CENTROID_MIN, CENTROID_MAX);
        lateral_hypha_n_d->set_centroid_z(CENTROID_MIN, CENTROID_MAX);

        opts->process(argc, argv); 

        std::srand(seed);

        // Sometimes the first few draws off the generator are always ~0.
        std::rand();
        std::rand();
        std::rand();

        if (img_in_fname)
        {
            imaging = new Imaging_model(img_in_fname, img_window, img_scale,
                    imaging_d);
        }
        else
        {
            imaging = imaging_d->sample(img_window, img_scale);
        }

        if (psf_in_fname)
        {
            psf = new PSF_model(psf_in_fname, padding, psf_d, imaging);
        }
        else
        {
            psf = psf_d->sample(padding, imaging);
        }

        if (alt_in_fname)
        {
            alternaria = new Alternaria(alt_in_fname, alt_d, spore_d,
                    apical_hypha_d, lateral_hypha_1_d, lateral_hypha_n_d);
        }
        else
        {
            float x = 0.5f*img_scale->get_x()*img_window->get_num_cols();
            float y = 0.5f*apical_hypha_d->sample_length();
            float z = 0.5f*img_scale->get_z()*img_window->get_num_imgs();
            alternaria = alt_d->sample(num_structs, x, y, z, JWSC_PI_2, 0, 
                    base_level, spore_d, apical_hypha_d, lateral_hypha_1_d, 
                    lateral_hypha_n_d);
        }

        if (interactive)
        {
            init_glut(argc, argv);
            init_gl();
            init_camera();
            init_modes();
            glutMainLoop();
        }
        else
        {
            if (img_out_fname)
            {
                imaging->print(img_out_fname);
            }
            if (psf_out_fname)
            {
                psf->print(psf_out_fname);
            }
            if (alt_out_fname)
            {
                alternaria->print(alt_out_fname);
            }
            if (psf_h_fmt)
            {
                psf->write_h(psf_h_fmt);
            }
            if (alt_scene_fmt)
            {
                Alternaria_model* model = new Alternaria_model(*alternaria, 
                        psf, imaging);
                model->write_scene(alt_scene_fmt, psf);
                delete model;
            }
        }
    }
    catch (Exception e)
    {
        cerr << "alternaria_gen: ";
        e.print();
    }

    clean_up();

    return EXIT_SUCCESS;
}