camera.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 <jwsc++/config.h>

#include <iostream>
#include <sstream>
#include <typeinfo>
#include <cstring>
#include <cmath>
#include <cassert>

#if defined JWSCXX_HAVE_OPENGL
#include <GL/gl.h>
#elif defined JWSCXX_HAVE_OPENGL_FRAMEWORK
#include <OpenGL/gl.h>
#endif

#include <jwsc/base/file_io.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/exception.h>
#include <jwsc++/base/readable.h>
#include <jwsc++/graphics/camera.h>

#ifdef JWSCXX_HAVE_DMALLOC
#include <dmalloc.h>
#endif


using namespace jwsc;
using namespace jwscxx::base;
using namespace jwscxx::graphics;


Camera_f::Camera_f
(
    float focal_pt_x,
    float focal_pt_y,
    float focal_pt_z,
    float ref_pt_x,
    float ref_pt_y,
    float ref_pt_z,
    float up_x,
    float up_y,
    float up_z,
    float aperture, 
    float clip_near, 
    float clip_far
)
throw (Arg_error)
{
    if (aperture <= 0)
    {
        throw Arg_error("Aperture <= 0");
    }
    this->aperture = aperture;

    if (clip_near < 0 || clip_far <= 0 || clip_near > clip_far)
    {
        throw Arg_error("Invalid clipping planes");
    }
    near = clip_near;
    far  = clip_far;

    prp     = 0;
    vrp     = 0;
    n       = 0;
    v       = 0;
    u       = 0;
    R_phi   = 0;
    R_theta = 0;
    R_psi   = 0;
    R       = 0;

    look_at(focal_pt_x, focal_pt_y, focal_pt_z, ref_pt_x, ref_pt_y, ref_pt_z,
            up_x, up_y, up_z);

    float tan_theta = std::tan(JWSC_DEG_TO_RAD * aperture / 2.0f);
    half_near_height = clip_near * tan_theta;
}


Camera_f::Camera_f
(
    const jwsc::Vector_f* focal_pt,
    const jwsc::Vector_f* ref_pt,
    const jwsc::Vector_f* up,
    float           aperture, 
    float           clip_near, 
    float           clip_far
)
throw (Arg_error)
{
    if (aperture <= 0)
    {
        throw Arg_error("Aperture <= 0");
    }
    this->aperture = aperture;

    if (clip_near < 0 || clip_far <= 0 || clip_near > clip_far)
    {
        throw Arg_error("Invalid clipping planes");
    }
    near = clip_near;
    far  = clip_far;

    prp     = 0;
    vrp     = 0;
    n       = 0;
    v       = 0;
    u       = 0;
    R_phi   = 0;
    R_theta = 0;
    R_psi   = 0;
    R       = 0;

    look_at(focal_pt, ref_pt, up);

    float tan_theta = std::tan(JWSC_DEG_TO_RAD * aperture / 2.0f);
    half_near_height = clip_near * tan_theta;
}


Camera_f::Camera_f(const Camera_f& c)
{
    prp     = 0;
    vrp     = 0;
    n       = 0;
    v       = 0;
    u       = 0;
    R_phi   = 0;
    R_theta = 0;
    R_psi   = 0;
    R       = 0;

    *this = c;
}


Camera_f::Camera_f(const char* fname) throw (jwscxx::base::Arg_error,
        jwscxx::base::IO_error)
{
    vrp     = 0;
    prp     = 0;
    u       = 0;
    v       = 0;
    n       = 0;
    R_phi   = 0;
    R_theta = 0;
    R_psi   = 0;
    R       = 0;

    jwscxx::base::Readable::read(fname);
}


Camera_f::Camera_f(std::istream& in) throw (jwscxx::base::Arg_error,
        jwscxx::base::IO_error)
{
    vrp     = 0;
    prp     = 0;
    u       = 0;
    v       = 0;
    n       = 0;
    R_phi   = 0;
    R_theta = 0;
    R_psi   = 0;
    R       = 0;

    read(in);
}


Camera_f::~Camera_f()
{
    free_vector_f(vrp);
    free_vector_f(prp);
    free_vector_f(u);
    free_vector_f(v);
    free_vector_f(n);
    free_matrix_f(R_phi);
    free_matrix_f(R_theta);
    free_matrix_f(R_psi);
    free_matrix_f(R);
}


Camera_f& Camera_f::operator=(const Camera_f &c)
{
    copy_vector_f(&vrp, c.vrp);
    copy_vector_f(&prp, c.prp);
    copy_vector_f(&u, c.u);
    copy_vector_f(&v, c.v);
    copy_vector_f(&n, c.n);

    aperture = c.aperture;
    focal_length = c.focal_length;
    half_focal_height = c.half_focal_height;
    half_near_height = c.half_near_height;
    near = c.near;
    far = c.far;

    if (c.R_phi   != 0)  copy_matrix_f(&R_phi, c.R_phi);
    if (c.R_theta != 0)  copy_matrix_f(&R_theta, c.R_theta);
    if (c.R_psi   != 0)  copy_matrix_f(&R_psi, c.R_psi);
    if (c.R       != 0)  copy_matrix_f(&R, c.R);

    return *this;
}


Camera_f* Camera_f::clone() const
{
    return new Camera_f(*this);
}


const jwsc::Vector_f* Camera_f::get_u() const 
{ 
    return (const jwsc::Vector_f*)u; 
}


const jwsc::Vector_f* Camera_f::get_v() const 
{ 
    return (const jwsc::Vector_f*)v; 
}


const jwsc::Vector_f* Camera_f::get_n() const 
{ 
    return (const jwsc::Vector_f*)n; 
}


const jwsc::Vector_f* Camera_f::get_prp() const 
{ 
    return (const jwsc::Vector_f*)prp; 
}


const jwsc::Vector_f* Camera_f::get_vrp() const 
{ 
    return (const jwsc::Vector_f*)vrp; 
}


float Camera_f::get_clipping_near() const 
{ 
    return near; 
}


float Camera_f::get_clipping_far() const 
{ 
    return far; 
}


float Camera_f::get_focal_length() const 
{ 
    return focal_length; 
}


float Camera_f::get_aperture() const 
{ 
    return aperture; 
}


float Camera_f::get_half_focal_height() const
{
    return half_focal_height;
}       


void Camera_f::set_aperture(float aperture) throw (Arg_error)
{
    if (aperture <= 0)
    {
        throw Arg_error("Aperture <= 0");
    }
    this->aperture = aperture;

    float tan_theta = std::tan(JWSC_DEG_TO_RAD * aperture / 2.0f);
    half_near_height = near * tan_theta;
    half_focal_height = focal_length * tan_theta;
}


void Camera_f::set_f_with_new_aperture_using_prp(float f) throw (Arg_error)
{
    if (f <= 0)
    {
        throw Arg_error("Focal length <= 0");
    }

    focal_length = f;

    multiply_vector_by_scalar_f(&prp, n, focal_length);
    add_vectors_f(&prp, prp, vrp);

    half_near_height = near * half_focal_height / focal_length;

    float theta = std::atan(half_focal_height / focal_length);
    aperture = JWSC_RAD_TO_DEG * 2.0f * theta;
}


void Camera_f::set_f_with_same_aperture_using_prp(float f) throw (Arg_error)
{
    if (f <= 0)
    {
        throw Arg_error("Focal length <= 0");
    }

    focal_length = f;

    multiply_vector_by_scalar_f(&prp, n, focal_length);
    add_vectors_f(&prp, prp, vrp);

    float tan_theta = std::tan(JWSC_DEG_TO_RAD * aperture / 2.0f);
    half_focal_height = focal_length * tan_theta;
}


void Camera_f::set_f_with_new_aperture_using_vrp(float f) throw (Arg_error)
{
    if (f <= 0)
    {
        throw Arg_error("Focal length <= 0");
    }

    focal_length = f;

    multiply_vector_by_scalar_f(&vrp, n, -focal_length);
    add_vectors_f(&vrp, vrp, prp);

    half_near_height = near * half_focal_height / focal_length;

    float theta = std::atan(half_focal_height / focal_length);
    aperture = JWSC_RAD_TO_DEG * 2.0f * theta;
}


void Camera_f::set_f_with_same_aperture_using_vrp(float f) throw (Arg_error)
{
    if (f <= 0)
    {
        throw Arg_error("Focal length <= 0");
    }

    focal_length = f;

    multiply_vector_by_scalar_f(&vrp, n, -focal_length);
    add_vectors_f(&vrp, vrp, prp);

    float tan_theta = std::tan(JWSC_DEG_TO_RAD * aperture / 2.0f);
    half_focal_height = focal_length * tan_theta;
}


void Camera_f::set_clipping(float clip_near, float clip_far) throw (Arg_error)
{
    if (clip_near < 0 || clip_far <= 0 || clip_near > clip_far)
    {
        throw Arg_error("Invalid clipping planes");
    }
    near = clip_near;
    far  = clip_far;

    float tan_theta   = std::tan(JWSC_DEG_TO_RAD * aperture / 2.0f);
    half_near_height  = near * tan_theta;
}


void Camera_f::rotate(float phi, float x, float y, float z) throw (Arg_error)
{
    float cos_phi, cos_theta, cos_psi;
    float sin_phi, sin_theta, sin_psi;
    float mag_xyz, mag_xy;

    mag_xyz = std::sqrt(x*x + y*y + z*z);

    if (mag_xyz < 1.0e-8)
    {
        throw Arg_error("Magnitude of rotation vector too small");
    }


    /* psi used to rotate around z-axis to put vector on the yz-plane. */
    mag_xy = std::sqrt(x*x + y*y);
    if (mag_xy > 1.0e-8)
    {
        cos_psi = y / mag_xy;
        sin_psi = x / mag_xy;
    }
    else
    {
        cos_psi = 0.0;
        sin_psi = 1.0;
    }

    /* theta used to rotate around x-axis to put vector onto the z-axis. */
    cos_theta = z / mag_xyz;
    sin_theta = mag_xy / mag_xyz;

    /* phi used to rotate around z-axis. This is actually the angle to rotate
     * around the vector. */
    cos_phi = std::cos(phi);
    sin_phi = std::sin(phi);

    create_zero_matrix_f(&R_phi, 3, 3);
    R_phi->elts[ 0 ][ 0 ] = cos_phi;
    R_phi->elts[ 0 ][ 1 ] = sin_phi;
    R_phi->elts[ 1 ][ 0 ] = -sin_phi;
    R_phi->elts[ 1 ][ 1 ] = cos_phi;
    R_phi->elts[ 2 ][ 2 ] = 1;

    create_zero_matrix_f(&R_theta, 3, 3);
    R_theta->elts[ 0 ][ 0 ] = 1;
    R_theta->elts[ 1 ][ 1 ] = cos_theta;
    R_theta->elts[ 1 ][ 2 ] = sin_theta;
    R_theta->elts[ 2 ][ 1 ] = -sin_theta;
    R_theta->elts[ 2 ][ 2 ] = cos_theta;

    create_zero_matrix_f(&R_psi, 3, 3);
    R_psi->elts[ 0 ][ 0 ] = cos_psi;
    R_psi->elts[ 0 ][ 1 ] = sin_psi;
    R_psi->elts[ 1 ][ 0 ] = -sin_psi;
    R_psi->elts[ 1 ][ 1 ] = cos_psi;
    R_psi->elts[ 2 ][ 2 ] = 1;

    /* R = R_psi * R_theta * R_phi * R_theta^T * R_psi^T */
    multiply_matrices_f(&R, R_psi, R_theta);
    multiply_matrices_f(&R, R, R_phi);
    transpose_matrix_f(&R_theta, R_theta);
    transpose_matrix_f(&R_psi, R_psi);
    multiply_matrices_f(&R, R, R_theta);
    multiply_matrices_f(&R, R, R_psi);

    multiply_matrix_and_vector_f(&vrp, R, vrp);
    multiply_matrix_and_vector_f(&prp, R, prp);
    multiply_matrix_and_vector_f(&u, R, u);
    multiply_matrix_and_vector_f(&v, R, v);
    multiply_matrix_and_vector_f(&n, R, n);
}


void Camera_f::rotate(float phi, const jwsc::Vector_f* r) throw (Arg_error)
{
    if (r->num_elts < 3)
    {
        throw Arg_error("Vector to rotate around must be in R^3");
    }

    rotate(phi, r->elts[0], r->elts[1], r->elts[2]);
}


void Camera_f::translate(float x, float y, float z)
{
    prp->elts[0] += x;  vrp->elts[0] += x;
    prp->elts[1] += y;  vrp->elts[1] += y;
    prp->elts[2] += z;  vrp->elts[2] += z;
}


void Camera_f::translate(const jwsc::Vector_f* t) throw (Arg_error)
{
    if (t->num_elts < 3)
    {
        throw Arg_error("Vector to translate by must be in R^3");
    }

    add_vectors_f(&vrp, vrp, t);
    add_vectors_f(&prp, prp, t);
}


void Camera_f::look_at
(
    float focal_pt_x,
    float focal_pt_y,
    float focal_pt_z,
    float ref_pt_x,
    float ref_pt_y,
    float ref_pt_z,
    float up_x,
    float up_y,
    float up_z
)
throw (jwscxx::base::Arg_error)
{
    float n_x = focal_pt_x - ref_pt_x;
    float n_y = focal_pt_y - ref_pt_y;
    float n_z = focal_pt_z - ref_pt_z;

    float n_mag = sqrt(n_x*n_x + n_y*n_y + n_z*n_z);
    if (n_mag < 1.0e-8)
    {
        throw Arg_error("Focal length is zero");
    }

    float up_mag = sqrt(up_x*up_x + up_y*up_y + up_z*up_z);
    if (up_mag < 1.0e-8)
    {
        throw Arg_error("Up vector has no magnitude");
    }

    float dp = up_x*n_x + up_y*n_y + up_z*n_z;
    if ((n_mag*up_mag - fabs(dp)) < 1.0e-8)
    {
        throw Arg_error("Up vector parallel to view direction");
    }

    create_vector_f(&prp, 3);
    prp->elts[0] = focal_pt_x;
    prp->elts[1] = focal_pt_y;
    prp->elts[2] = focal_pt_z;

    create_vector_f(&vrp, 3);
    vrp->elts[0] = ref_pt_x;
    vrp->elts[1] = ref_pt_y;
    vrp->elts[2] = ref_pt_z;

    focal_length = n_mag;
    n_mag = 1.0f / n_mag;
    create_vector_f(&n, 3);
    n->elts[0] = n_x * n_mag;
    n->elts[1] = n_y * n_mag;
    n->elts[2] = n_z * n_mag;

    create_vector_f(&v, 3);
    v->elts[0] = up_x / up_mag;
    v->elts[1] = up_y / up_mag;
    v->elts[2] = up_z / up_mag;

    calc_3d_vector_cross_prod_f(&u, v, n);

    calc_3d_vector_cross_prod_f(&v, n, u);

    float tan_theta = std::tan(JWSC_DEG_TO_RAD * aperture / 2.0f);
    half_focal_height = focal_length * tan_theta;
}


void Camera_f::look_at
(
    const jwsc::Vector_f* focal_pt,
    const jwsc::Vector_f* ref_pt,
    const jwsc::Vector_f* up
)
throw (Arg_error)
{
    if (focal_pt->num_elts != 3 ||
        ref_pt->num_elts   != 3 ||
        up->num_elts       != 3)
    {
        throw Arg_error("Camera look vectors not in R^3");
    }

    look_at(focal_pt->elts[0], focal_pt->elts[1], focal_pt->elts[2],
            ref_pt->elts[0], ref_pt->elts[1], ref_pt->elts[2],
            up->elts[0], up->elts[1], up->elts[2]);
}


void Camera_f::set_gl_modelview() const
{
    GLfloat M[16] = {0};

    glMatrixMode(GL_MODELVIEW);

    glLoadIdentity();

    M[0] = u->elts[0];  M[4] = u->elts[1];  M[8]  = u->elts[2];
    M[1] = v->elts[0];  M[5] = v->elts[1];  M[9]  = v->elts[2];
    M[2] = n->elts[0];  M[6] = n->elts[1];  M[10] = n->elts[2];
    M[15] = 1;
    glMultMatrixf(M);

    glTranslatef(-prp->elts[0], -prp->elts[1], -prp->elts[2]);
}


void Camera_f::set_gl_projection() const
{
    glMatrixMode(GL_PROJECTION);

    glLoadIdentity();

    float aspect = get_gl_viewport_aspect();
    float left   = -aspect*half_near_height;
    float right  =  aspect*half_near_height;
    float top    =  half_near_height;
    float bottom = -half_near_height;

    glFrustum(left, right, bottom, top, near, far);
}


void Camera_f::get_image_pt_on_focal_plane
(
    jwsc::Vector_f** pt_out, 
    float            x, 
    float            y
)
{
    float half_view_height = 0.5f*get_gl_viewport_height();
    float focal_x = x * (half_focal_height / half_view_height);
    float focal_y = y * (half_focal_height / half_view_height);

    create_vector_f(pt_out, 3);
    Vector_f* pt = *pt_out;

    pt->elts[0] = focal_x;
    pt->elts[1] = focal_y;
    pt->elts[2] = -focal_length;

    create_matrix_f(&R, 3, 3);
    R->elts[0][0] = u->elts[0];
    R->elts[1][0] = u->elts[1];
    R->elts[2][0] = u->elts[2];

    R->elts[0][1] = v->elts[0];
    R->elts[1][1] = v->elts[1];
    R->elts[2][1] = v->elts[2];

    R->elts[0][2] = n->elts[0];
    R->elts[1][2] = n->elts[1];
    R->elts[2][2] = n->elts[2];

    multiply_matrix_and_vector_f(&pt, R, pt);

    add_vectors_f(&pt, pt, prp);
}


void Camera_f::get_image_pt_on_focal_plane
(
    jwsc::Vector_f**      pt_out, 
    const jwsc::Vector_f* pt_in
)
throw (Arg_error)
{
    if (pt_in->num_elts != 2)
    {
        throw Arg_error("Invalid image pt vector");
    }

    get_image_pt_on_focal_plane(pt_out, pt_in->elts[0], pt_in->elts[1]);
}


void Camera_f::get_image_pt_on_near_clipping_plane
(
    jwsc::Vector_f** pt_out, 
    float            x, 
    float            y
)
{
    float half_view_height = 0.5f*get_gl_viewport_height();
    float near_x = x * (half_near_height / half_view_height);
    float near_y = y * (half_near_height / half_view_height);

    create_vector_f(pt_out, 3);
    Vector_f* pt = *pt_out;

    pt->elts[0] = near_x;
    pt->elts[1] = near_y;
    pt->elts[2] = -near;

    create_matrix_f(&R, 3, 3);
    R->elts[0][0] = u->elts[0];
    R->elts[1][0] = u->elts[1];
    R->elts[2][0] = u->elts[2];

    R->elts[0][1] = v->elts[0];
    R->elts[1][1] = v->elts[1];
    R->elts[2][1] = v->elts[2];

    R->elts[0][2] = n->elts[0];
    R->elts[1][2] = n->elts[1];
    R->elts[2][2] = n->elts[2];

    multiply_matrix_and_vector_f(&pt, R, pt);

    add_vectors_f(&pt, pt, prp);
}


void Camera_f::get_image_pt_on_near_clipping_plane
(
    jwsc::Vector_f**      pt_out,
    const jwsc::Vector_f* pt_in
)
throw (Arg_error)
{
    if (pt_in->num_elts != 2)
    {
        throw Arg_error("invalid image pt vector");
    }

    get_image_pt_on_near_clipping_plane(pt_out, pt_in->elts[0], pt_in->elts[1]);
}


void Camera_f::get_image_pt_on_far_clipping_plane
(
    jwsc::Vector_f** pt_out, 
    float            x, 
    float            y
)
{
    float half_view_height = 0.5f*get_gl_viewport_height();
    float half_far_height = far * half_near_height / near;
    float far_x = x * (half_far_height / half_view_height);
    float far_y = y * (half_far_height / half_view_height);

    create_vector_f(pt_out, 3);
    Vector_f* pt = *pt_out;

    pt->elts[0] = far_x;
    pt->elts[1] = far_y;
    pt->elts[2] = -far;

    create_matrix_f(&R, 3, 3);
    R->elts[0][0] = u->elts[0];
    R->elts[1][0] = u->elts[1];
    R->elts[2][0] = u->elts[2];

    R->elts[0][1] = v->elts[0];
    R->elts[1][1] = v->elts[1];
    R->elts[2][1] = v->elts[2];

    R->elts[0][2] = n->elts[0];
    R->elts[1][2] = n->elts[1];
    R->elts[2][2] = n->elts[2];

    multiply_matrix_and_vector_f(&pt, R, pt);

    add_vectors_f(&pt, pt, prp);
}


void Camera_f::get_image_pt_on_far_clipping_plane
(
    jwsc::Vector_f**      pt_out,
    const jwsc::Vector_f* pt_in
)
throw (Arg_error)
{
    if (pt_in->num_elts != 2)
    {
        throw Arg_error("invalid image pt vector");
    }

    get_image_pt_on_far_clipping_plane(pt_out, pt_in->elts[0], pt_in->elts[1]);
}


float Camera_f::get_gl_viewport_aspect()
{
    GLfloat viewport[4] = {0};

    glGetFloatv(GL_VIEWPORT, viewport);
    return viewport[2] / viewport[3];
}


float Camera_f::get_gl_viewport_width()
{
    GLfloat viewport[4] = {0};

    glGetFloatv(GL_VIEWPORT, viewport);
    return viewport[2];
}


float Camera_f::get_gl_viewport_height()
{
    GLfloat viewport[4] = {0};

    glGetFloatv(GL_VIEWPORT, viewport);
    return viewport[3];
}


void Camera_f::capture_gl_view(jwsc::Image_f** img_out)
{
    uint32_t num_cols = (uint32_t)Camera_f::get_gl_viewport_width();
    uint32_t num_rows = (uint32_t)Camera_f::get_gl_viewport_height();

    float* img_buf = (float*)malloc(3*num_cols*num_rows*sizeof(float));
    assert(img_buf);

    glPixelStorei(GL_PACK_ALIGNMENT, 1);

    glReadPixels(0, 0, num_cols, num_rows, GL_RGB, GL_FLOAT, img_buf);

    create_image_f(img_out, num_rows, num_cols);
    Image_f* img = *img_out;

    for (uint32_t row = 0; row < num_rows; row++)
    {
        for (uint32_t col = 0; col < 3*num_cols; col += 3)
        {
            uint32_t r = num_rows - row - 1;

            img->pxls[ r ][ col/3 ].r = img_buf[ col+0 + row*3*num_cols ];
            img->pxls[ r ][ col/3 ].g = img_buf[ col+1 + row*3*num_cols ];
            img->pxls[ r ][ col/3 ].b = img_buf[ col+2 + row*3*num_cols ];
        }
    }

    free(img_buf);
}


void Camera_f::capture_gl_view(const char* fname) throw (jwscxx::base::IO_error)
{
    Image_f* img = 0;
    capture_gl_view(&img);

    Error* err;
    if ((err = write_image_f(img, fname)) != 0)
    {
        std::ostringstream ost(err->msg);
        free_error(err);
        throw IO_error(ost.str());
    }

    free_image_f(img);
}


void Camera_f::capture_gl_view(const char* fname_fmt, uint32_t N)
    throw (jwscxx::base::IO_error)
{
    char fname[256] = {0};

    snprintf(fname, 255, fname_fmt, N);

    capture_gl_view(fname);
}


void Camera_f::read(std::istream& in) throw (jwscxx::base::IO_error,
        jwscxx::base::Arg_error)
{
    using std::ostringstream;
    using std::istringstream;

    const char* type_name = typeid(*this).name();
    const char* field_value;

    // Type
    if (!(field_value = read_field_value(in, "Type")))
    {
        throw Arg_error("Missing Type field");
    }
    if (strncmp(field_value, type_name, strlen(type_name)) != 0)
    {
        ostringstream ost;
        ost << "Tried to read a '" << field_value << "' as a '"
            << type_name << "'";
        throw Arg_error(ost.str());
    }


    // VRP
    if (!(field_value = read_field_value(in, "VRP")))
    {
        throw Arg_error("Missing VRP fields");
    }
    istringstream ist(field_value);
    ist >> vrp->elts[0] >> vrp->elts[1] >> vrp->elts[2];
    if (ist.fail())
    {
        throw Arg_error("Invalid VRP");
    }
    ist.clear(std::ios_base::goodbit);


    // PRP
    if (!(field_value = read_field_value(in, "PRP")))
    {
        throw Arg_error("Missing PRP fields");
    }
    ist.str(field_value);
    ist >> prp->elts[0] >> prp->elts[1] >> prp->elts[2];
    if (ist.fail())
    {
        throw Arg_error("Invalid PRP");
    }
    ist.clear(std::ios_base::goodbit);


    // U
    if (!(field_value = read_field_value(in, "U")))
    {
        throw Arg_error("Missing U fields");
    }
    ist.str(field_value);
    ist >> u->elts[ 0 ] >> u->elts[1] >> u->elts[2];
    if (ist.fail())
    {
        throw Arg_error("Invalid U");
    }
    ist.clear(std::ios_base::goodbit);


    // V
    if (!(field_value = read_field_value(in, "V")))
    {
        throw Arg_error("Missing V fields");
    }
    ist.str(field_value);
    ist >> v->elts[ 0 ] >> v->elts[1] >> v->elts[2];
    if (ist.fail())
    {
        throw Arg_error("Invalid V");
    }
    ist.clear(std::ios_base::goodbit);


    // N
    if (!(field_value = read_field_value(in, "N")))
    {
        throw Arg_error("Missing N fields");
    }
    ist.str(field_value);
    ist >> n->elts[ 0 ] >> n->elts[1] >> n->elts[2];
    if (ist.fail())
    {
        throw Arg_error("Invalid N");
    }
    ist.clear(std::ios_base::goodbit);


    // F
    if (!(field_value = read_field_value(in, "F")))
    {
        throw Arg_error("Missing F field");
    }
    ist.str(field_value);
    ist >> focal_length;
    if (ist.fail())
    {
        throw Arg_error("Invalid focal length");
    }
    ist.clear(std::ios_base::goodbit);


    // Aperture
    if (!(field_value = read_field_value(in, "Aperture")))
    {
        throw Arg_error("Missing Aperture field");
    }
    ist.str(field_value);
    ist >> aperture;
    if (ist.fail())
    {
        throw Arg_error("Invalid aperture");
    }
    ist.clear(std::ios_base::goodbit);
    if (aperture <= 0)
    {
        throw Arg_error("Aperture <= 0");
    }


    // Near
    if (!(field_value = read_field_value(in, "Near")))
    {
        throw Arg_error("Missing Near field");
    }
    ist.str(field_value);
    ist >> near;
    if (ist.fail())
    {
        throw Arg_error("Invalid near clipping distance");
    }
    ist.clear(std::ios_base::goodbit);


    // Far
    if (!(field_value = read_field_value(in, "Far")))
    {
        throw Arg_error("Missing Far field");
    }
    ist.str(field_value);
    ist >> far;
    if (ist.fail())
    {
        throw Arg_error("Invalid far clipping distance");
    }
    ist.clear(std::ios_base::goodbit);
    if (near < 0 || far <= 0 || near > far)
    {
        throw Arg_error("Invalid clipping planes");
    }


    // Derived...
    float tan_theta = std::tan(JWSC_DEG_TO_RAD * aperture / 2.0f);
    half_near_height = near * tan_theta;
    half_focal_height = focal_length * tan_theta;
}


void Camera_f::write(std::ostream& out) const throw (jwscxx::base::IO_error)
{
    out << "     Type: " << typeid(*this).name() << '\n'
        << "      VRP: " << vrp->elts[0] << ' '
                         << vrp->elts[1] << ' '
                         << vrp->elts[2] << '\n'
        << "      PRP: " << prp->elts[0] << ' '
                         << prp->elts[1] << ' '
                         << prp->elts[2] << '\n'
        << "        U: " << u->elts[0] << ' '
                         << u->elts[1] << ' '
                         << u->elts[2] << '\n'
        << "        V: " << v->elts[0] << ' '
                         << v->elts[1] << ' '
                         << v->elts[2] << '\n'
        << "        N: " << n->elts[0] << ' '
                         << n->elts[1] << ' '
                         << n->elts[2] << '\n'
        << "        F: " << focal_length << '\n'
        << " Aperture: " << aperture << '\n'
        << "     Near: " << near << '\n'
        << "      Far: " << far << '\n';
}


Camera_d::Camera_d
(
    double focal_pt_x,
    double focal_pt_y,
    double focal_pt_z,
    double ref_pt_x,
    double ref_pt_y,
    double ref_pt_z,
    double up_x,
    double up_y,
    double up_z,
    double aperture, 
    double clip_near, 
    double clip_far
)
throw (Arg_error)
{
    if (aperture <= 0)
    {
        throw Arg_error("Aperture <= 0");
    }
    this->aperture = aperture;

    if (clip_near < 0 || clip_far <= 0 || clip_near > clip_far)
    {
        throw Arg_error("Invalid clipping planes");
    }
    near = clip_near;
    far  = clip_far;

    prp     = 0;
    vrp     = 0;
    n       = 0;
    v       = 0;
    u       = 0;
    R_phi   = 0;
    R_theta = 0;
    R_psi   = 0;
    R       = 0;

    look_at(focal_pt_x, focal_pt_y, focal_pt_z, ref_pt_x, ref_pt_y, ref_pt_z,
            up_x, up_y, up_z);

    double tan_theta = std::tan(JWSC_DEG_TO_RAD * aperture / 2.0);
    half_near_height = clip_near * tan_theta;
}


Camera_d::Camera_d
(
    const jwsc::Vector_d* focal_pt,
    const jwsc::Vector_d* ref_pt,
    const jwsc::Vector_d* up,
    double                aperture, 
    double                clip_near, 
    double                clip_far
)
throw (Arg_error)
{
    if (aperture <= 0)
    {
        throw Arg_error("Aperture <= 0");
    }
    this->aperture = aperture;

    if (clip_near < 0 || clip_far <= 0 || clip_near > clip_far)
    {
        throw Arg_error("Invalid clipping planes");
    }
    near = clip_near;
    far  = clip_far;

    prp     = 0;
    vrp     = 0;
    n       = 0;
    v       = 0;
    u       = 0;
    R_phi   = 0;
    R_theta = 0;
    R_psi   = 0;
    R       = 0;

    look_at(focal_pt, ref_pt, up);

    double tan_theta = std::tan(JWSC_DEG_TO_RAD * aperture / 2.0);
    half_near_height = clip_near * tan_theta;
}


Camera_d::Camera_d(const Camera_d& c)
{
    prp     = 0;
    vrp     = 0;
    n       = 0;
    v       = 0;
    u       = 0;
    R_phi   = 0;
    R_theta = 0;
    R_psi   = 0;
    R       = 0;

    *this = c;
}


Camera_d::Camera_d(const char* fname) throw (jwscxx::base::Arg_error,
        jwscxx::base::IO_error)
{
    vrp     = 0;
    prp     = 0;
    u       = 0;
    v       = 0;
    n       = 0;
    R_phi   = 0;
    R_theta = 0;
    R_psi   = 0;
    R       = 0;

    jwscxx::base::Readable::read(fname);
}


Camera_d::Camera_d(std::istream& in) throw (jwscxx::base::Arg_error,
        jwscxx::base::IO_error)
{
    vrp     = 0;
    prp     = 0;
    u       = 0;
    v       = 0;
    n       = 0;
    R_phi   = 0;
    R_theta = 0;
    R_psi   = 0;
    R       = 0;

    read(in);
}


Camera_d::~Camera_d()
{
    free_vector_d(vrp);
    free_vector_d(prp);
    free_vector_d(u);
    free_vector_d(v);
    free_vector_d(n);
    free_matrix_d(R_phi);
    free_matrix_d(R_theta);
    free_matrix_d(R_psi);
    free_matrix_d(R);
}


Camera_d& Camera_d::operator=(const Camera_d &c)
{
    copy_vector_d(&vrp, c.vrp);
    copy_vector_d(&prp, c.prp);
    copy_vector_d(&u, c.u);
    copy_vector_d(&v, c.v);
    copy_vector_d(&n, c.n);

    aperture = c.aperture;
    focal_length = c.focal_length;
    half_focal_height = c.half_focal_height;
    half_near_height = c.half_near_height;
    near = c.near;
    far = c.far;

    if (c.R_phi   != 0)  copy_matrix_d(&R_phi, c.R_phi);
    if (c.R_theta != 0)  copy_matrix_d(&R_theta, c.R_theta);
    if (c.R_psi   != 0)  copy_matrix_d(&R_psi, c.R_psi);
    if (c.R       != 0)  copy_matrix_d(&R, c.R);

    return *this;
}


Camera_d* Camera_d::clone() const
{
    return new Camera_d(*this);
}


const jwsc::Vector_d* Camera_d::get_u() const 
{ 
    return (const jwsc::Vector_d*)u; 
}


const jwsc::Vector_d* Camera_d::get_v() const 
{ 
    return (const jwsc::Vector_d*)v; 
}


const jwsc::Vector_d* Camera_d::get_n() const 
{ 
    return (const jwsc::Vector_d*)n; 
}


const jwsc::Vector_d* Camera_d::get_prp() const 
{ 
    return (const jwsc::Vector_d*)prp; 
}


const jwsc::Vector_d* Camera_d::get_vrp() const 
{ 
    return (const jwsc::Vector_d*)vrp; 
}


double Camera_d::get_clipping_near() const 
{ 
    return near; 
}


double Camera_d::get_clipping_far() const 
{ 
    return far; 
}


double Camera_d::get_focal_length() const 
{ 
    return focal_length; 
}


double Camera_d::get_aperture() const 
{ 
    return aperture; 
}


double Camera_d::get_half_focal_height() const
{
    return half_focal_height;
}       


void Camera_d::set_aperture(double aperture) throw (Arg_error)
{
    if (aperture <= 0)
    {
        throw Arg_error("Aperture <= 0");
    }
    this->aperture = aperture;

    double tan_theta = std::tan(JWSC_DEG_TO_RAD * aperture / 2.0);
    half_near_height = near * tan_theta;
    half_focal_height = focal_length * tan_theta;
}


void Camera_d::set_f_with_new_aperture_using_prp(double f) throw (Arg_error)
{
    if (f <= 0)
    {
        throw Arg_error("Focal length <= 0");
    }

    focal_length = f;

    multiply_vector_by_scalar_d(&prp, n, focal_length);
    add_vectors_d(&prp, prp, vrp);

    half_near_height = near * half_focal_height / focal_length;

    double theta = std::atan(half_focal_height / focal_length);
    aperture = JWSC_RAD_TO_DEG * 2.0 * theta;
}


void Camera_d::set_f_with_same_aperture_using_prp(double f) throw (Arg_error)
{
    if (f <= 0)
    {
        throw Arg_error("Focal length <= 0");
    }

    focal_length = f;

    multiply_vector_by_scalar_d(&prp, n, focal_length);
    add_vectors_d(&prp, prp, vrp);

    double tan_theta = std::tan(JWSC_DEG_TO_RAD * aperture / 2.0);
    half_focal_height = focal_length * tan_theta;
}


void Camera_d::set_f_with_new_aperture_using_vrp(double f) throw (Arg_error)
{
    if (f <= 0)
    {
        throw Arg_error("Focal length <= 0");
    }

    focal_length = f;

    multiply_vector_by_scalar_d(&vrp, n, -focal_length);
    add_vectors_d(&vrp, vrp, prp);

    half_near_height = near * half_focal_height / focal_length;

    double theta = std::atan(half_focal_height / focal_length);
    aperture = JWSC_RAD_TO_DEG * 2.0 * theta;
}


void Camera_d::set_f_with_same_aperture_using_vrp(double f) throw (Arg_error)
{
    if (f <= 0)
    {
        throw Arg_error("Focal length <= 0");
    }

    focal_length = f;

    multiply_vector_by_scalar_d(&vrp, n, -focal_length);
    add_vectors_d(&vrp, vrp, prp);

    double tan_theta = std::tan(JWSC_DEG_TO_RAD * aperture / 2.0);
    half_focal_height = focal_length * tan_theta;
}


void Camera_d::set_clipping(double clip_near, double clip_far) throw (Arg_error)
{
    if (clip_near < 0 || clip_far <= 0 || clip_near > clip_far)
    {
        throw Arg_error("Invalid clipping planes");
    }
    near = clip_near;
    far  = clip_far;

    double tan_theta = std::tan(JWSC_DEG_TO_RAD * aperture / 2.0);
    half_near_height = near * tan_theta;
}


void Camera_d::rotate(double phi, double x, double y, double z) 
    throw (Arg_error)
{
    double cos_phi, cos_theta, cos_psi;
    double sin_phi, sin_theta, sin_psi;
    double mag_xyz, mag_xy;

    mag_xyz = std::sqrt(x*x + y*y + z*z);

    if (mag_xyz < 1.0e-8)
    {
        throw Arg_error("Magnitude of rotation vector too small");
    }


    /* psi used to rotate around z-axis to put vector on the yz-plane. */
    mag_xy = std::sqrt(x*x + y*y);
    if (mag_xy > 1.0e-8)
    {
        cos_psi = y / mag_xy;
        sin_psi = x / mag_xy;
    }
    else
    {
        cos_psi = 0.0;
        sin_psi = 1.0;
    }

    /* theta used to rotate around x-axis to put vector onto the z-axis. */
    cos_theta = z / mag_xyz;
    sin_theta = mag_xy / mag_xyz;

    /* phi used to rotate around z-axis. This is actually the angle to rotate
     * around the vector. */
    cos_phi = std::cos(phi);
    sin_phi = std::sin(phi);

    create_zero_matrix_d(&R_phi, 3, 3);
    R_phi->elts[ 0 ][ 0 ] = cos_phi;
    R_phi->elts[ 0 ][ 1 ] = sin_phi;
    R_phi->elts[ 1 ][ 0 ] = -sin_phi;
    R_phi->elts[ 1 ][ 1 ] = cos_phi;
    R_phi->elts[ 2 ][ 2 ] = 1;

    create_zero_matrix_d(&R_theta, 3, 3);
    R_theta->elts[ 0 ][ 0 ] = 1;
    R_theta->elts[ 1 ][ 1 ] = cos_theta;
    R_theta->elts[ 1 ][ 2 ] = sin_theta;
    R_theta->elts[ 2 ][ 1 ] = -sin_theta;
    R_theta->elts[ 2 ][ 2 ] = cos_theta;

    create_zero_matrix_d(&R_psi, 3, 3);
    R_psi->elts[ 0 ][ 0 ] = cos_psi;
    R_psi->elts[ 0 ][ 1 ] = sin_psi;
    R_psi->elts[ 1 ][ 0 ] = -sin_psi;
    R_psi->elts[ 1 ][ 1 ] = cos_psi;
    R_psi->elts[ 2 ][ 2 ] = 1;

    /* R = R_psi * R_theta * R_phi * R_theta^T * R_psi^T */
    multiply_matrices_d(&R, R_psi, R_theta);
    multiply_matrices_d(&R, R, R_phi);
    transpose_matrix_d(&R_theta, R_theta);
    transpose_matrix_d(&R_psi, R_psi);
    multiply_matrices_d(&R, R, R_theta);
    multiply_matrices_d(&R, R, R_psi);

    multiply_matrix_and_vector_d(&vrp, R, vrp);
    multiply_matrix_and_vector_d(&prp, R, prp);
    multiply_matrix_and_vector_d(&u, R, u);
    multiply_matrix_and_vector_d(&v, R, v);
    multiply_matrix_and_vector_d(&n, R, n);
}


void Camera_d::rotate(double phi, const jwsc::Vector_d* r) throw (Arg_error)
{
    if (r->num_elts < 3)
    {
        throw Arg_error("Vector to rotate around must be in R^3");
    }

    rotate(phi, r->elts[0], r->elts[1], r->elts[2]);
}


void Camera_d::translate(double x, double y, double z)
{
    prp->elts[0] += x;  vrp->elts[0] += x;
    prp->elts[1] += y;  vrp->elts[1] += y;
    prp->elts[2] += z;  vrp->elts[2] += z;
}


void Camera_d::translate(const jwsc::Vector_d* t) throw (Arg_error)
{
    if (t->num_elts < 3)
    {
        throw Arg_error("Vector to translate by must be in R^3");
    }

    add_vectors_d(&vrp, vrp, t);
    add_vectors_d(&prp, prp, t);
}


void Camera_d::look_at
(
    double focal_pt_x,
    double focal_pt_y,
    double focal_pt_z,
    double ref_pt_x,
    double ref_pt_y,
    double ref_pt_z,
    double up_x,
    double up_y,
    double up_z
)
throw (jwscxx::base::Arg_error)
{
    double n_x = focal_pt_x - ref_pt_x;
    double n_y = focal_pt_y - ref_pt_y;
    double n_z = focal_pt_z - ref_pt_z;

    double n_mag = sqrt(n_x*n_x + n_y*n_y + n_z*n_z);
    if (n_mag < 1.0e-8)
    {
        throw Arg_error("Focal length is zero");
    }

    double up_mag = sqrt(up_x*up_x + up_y*up_y + up_z*up_z);
    if (up_mag < 1.0e-8)
    {
        throw Arg_error("Up vector has no magnitude");
    }

    double dp = up_x*n_x + up_y*n_y + up_z*n_z;
    if ((n_mag*up_mag - fabs(dp)) < 1.0e-8)
    {
        throw Arg_error("Up vector parallel to view direction");
    }

    create_vector_d(&prp, 3);
    prp->elts[0] = focal_pt_x;
    prp->elts[1] = focal_pt_y;
    prp->elts[2] = focal_pt_z;

    create_vector_d(&vrp, 3);
    vrp->elts[0] = ref_pt_x;
    vrp->elts[1] = ref_pt_y;
    vrp->elts[2] = ref_pt_z;

    focal_length = n_mag;
    n_mag = 1.0 / n_mag;
    create_vector_d(&n, 3);
    n->elts[0] = n_x * n_mag;
    n->elts[1] = n_y * n_mag;
    n->elts[2] = n_z * n_mag;

    create_vector_d(&v, 3);
    v->elts[0] = up_x / up_mag;
    v->elts[1] = up_y / up_mag;
    v->elts[2] = up_z / up_mag;

    calc_3d_vector_cross_prod_d(&u, v, n);

    calc_3d_vector_cross_prod_d(&v, n, u);

    double tan_theta = std::tan(JWSC_DEG_TO_RAD * aperture / 2.0);
    half_focal_height = focal_length * tan_theta;
}


void Camera_d::look_at
(
    const jwsc::Vector_d* focal_pt,
    const jwsc::Vector_d* ref_pt,
    const jwsc::Vector_d* up
)
throw (Arg_error)
{
    if (focal_pt->num_elts != 3 ||
        ref_pt->num_elts   != 3 ||
        up->num_elts       != 3)
    {
        throw Arg_error("Camera look vectors not in R^3");
    }

    look_at(focal_pt->elts[0], focal_pt->elts[1], focal_pt->elts[2],
            ref_pt->elts[0], ref_pt->elts[1], ref_pt->elts[2],
            up->elts[0], up->elts[1], up->elts[2]);
}


void Camera_d::set_gl_modelview() const
{
    GLdouble M[16] = {0};

    glMatrixMode(GL_MODELVIEW);

    glLoadIdentity();

    M[0] = u->elts[0];  M[4] = u->elts[1];  M[8]  = u->elts[2];
    M[1] = v->elts[0];  M[5] = v->elts[1];  M[9]  = v->elts[2];
    M[2] = n->elts[0];  M[6] = n->elts[1];  M[10] = n->elts[2];
    M[15] = 1;
    glMultMatrixd(M);

    glTranslated(-prp->elts[0], -prp->elts[1], -prp->elts[2]);
}


void Camera_d::set_gl_projection() const
{
    glMatrixMode(GL_PROJECTION);

    glLoadIdentity();

    double aspect = get_gl_viewport_aspect();
    double left   = -aspect*half_near_height;
    double right  =  aspect*half_near_height;
    double top    =  half_near_height;
    double bottom = -half_near_height;

    glFrustum(left, right, bottom, top, near, far);
}


void Camera_d::get_image_pt_on_focal_plane
(
    jwsc::Vector_d** pt_out, 
    double           x, 
    double           y
)
{
    double half_view_height = 0.5*get_gl_viewport_height();
    double focal_x = x * (half_focal_height / half_view_height);
    double focal_y = y * (half_focal_height / half_view_height);

    create_vector_d(pt_out, 3);
    Vector_d* pt = *pt_out;

    pt->elts[0] = focal_x;
    pt->elts[1] = focal_y;
    pt->elts[2] = -focal_length;

    create_matrix_d(&R, 3, 3);
    R->elts[0][0] = u->elts[0];
    R->elts[1][0] = u->elts[1];
    R->elts[2][0] = u->elts[2];

    R->elts[0][1] = v->elts[0];
    R->elts[1][1] = v->elts[1];
    R->elts[2][1] = v->elts[2];

    R->elts[0][2] = n->elts[0];
    R->elts[1][2] = n->elts[1];
    R->elts[2][2] = n->elts[2];

    multiply_matrix_and_vector_d(&pt, R, pt);

    add_vectors_d(&pt, pt, prp);
}


void Camera_d::get_image_pt_on_focal_plane
(
    jwsc::Vector_d**      pt_out, 
    const jwsc::Vector_d* pt_in
)
throw (Arg_error)
{
    if (pt_in->num_elts != 2)
    {
        throw Arg_error("Invalid image pt vector");
    }

    get_image_pt_on_focal_plane(pt_out, pt_in->elts[0], pt_in->elts[1]);
}


void Camera_d::get_image_pt_on_near_clipping_plane
(
    jwsc::Vector_d** pt_out, 
    double           x, 
    double           y
)
{
    double half_view_height = 0.5*get_gl_viewport_height();
    double near_x = x * (half_near_height / half_view_height);
    double near_y = y * (half_near_height / half_view_height);

    create_vector_d(pt_out, 3);
    Vector_d* pt = *pt_out;

    pt->elts[0] = near_x;
    pt->elts[1] = near_y;
    pt->elts[2] = -near;

    create_matrix_d(&R, 3, 3);
    R->elts[0][0] = u->elts[0];
    R->elts[1][0] = u->elts[1];
    R->elts[2][0] = u->elts[2];

    R->elts[0][1] = v->elts[0];
    R->elts[1][1] = v->elts[1];
    R->elts[2][1] = v->elts[2];

    R->elts[0][2] = n->elts[0];
    R->elts[1][2] = n->elts[1];
    R->elts[2][2] = n->elts[2];

    multiply_matrix_and_vector_d(&pt, R, pt);

    add_vectors_d(&pt, pt, prp);
}


void Camera_d::get_image_pt_on_near_clipping_plane
(
    jwsc::Vector_d**      pt_out,
    const jwsc::Vector_d* pt_in
)
throw (Arg_error)
{
    if (pt_in->num_elts != 2)
    {
        throw Arg_error("invalid image pt vector");
    }

    get_image_pt_on_near_clipping_plane(pt_out, pt_in->elts[0], pt_in->elts[1]);
}


void Camera_d::get_image_pt_on_far_clipping_plane
(
    jwsc::Vector_d** pt_out, 
    double           x, 
    double           y
)
{
    double half_view_height = 0.5*get_gl_viewport_height();
    double half_far_height = far * half_near_height / near;
    double far_x = x * (half_far_height / half_view_height);
    double far_y = y * (half_far_height / half_view_height);

    create_vector_d(pt_out, 3);
    Vector_d* pt = *pt_out;

    pt->elts[0] = far_x;
    pt->elts[1] = far_y;
    pt->elts[2] = -far;

    create_matrix_d(&R, 3, 3);
    R->elts[0][0] = u->elts[0];
    R->elts[1][0] = u->elts[1];
    R->elts[2][0] = u->elts[2];

    R->elts[0][1] = v->elts[0];
    R->elts[1][1] = v->elts[1];
    R->elts[2][1] = v->elts[2];

    R->elts[0][2] = n->elts[0];
    R->elts[1][2] = n->elts[1];
    R->elts[2][2] = n->elts[2];

    multiply_matrix_and_vector_d(&pt, R, pt);

    add_vectors_d(&pt, pt, prp);
}


void Camera_d::get_image_pt_on_far_clipping_plane
(
    jwsc::Vector_d**      pt_out,
    const jwsc::Vector_d* pt_in
)
throw (Arg_error)
{
    if (pt_in->num_elts != 2)
    {
        throw Arg_error("invalid image pt vector");
    }

    get_image_pt_on_far_clipping_plane(pt_out, pt_in->elts[0], pt_in->elts[1]);
}


double Camera_d::get_gl_viewport_aspect()
{
    GLdouble viewport[4] = {0};

    glGetDoublev(GL_VIEWPORT, viewport);
    return viewport[2] / viewport[3];
}


double Camera_d::get_gl_viewport_width()
{
    GLdouble viewport[4] = {0};

    glGetDoublev(GL_VIEWPORT, viewport);
    return viewport[2];
}


double Camera_d::get_gl_viewport_height()
{
    GLdouble viewport[4] = {0};

    glGetDoublev(GL_VIEWPORT, viewport);
    return viewport[3];
}


void Camera_d::capture_gl_view(jwsc::Image_f** img_out)
{
    uint32_t num_cols = (uint32_t)Camera_d::get_gl_viewport_width();
    uint32_t num_rows = (uint32_t)Camera_d::get_gl_viewport_height();

    float* img_buf = (float*)malloc(3*num_cols*num_rows*sizeof(float));
    assert(img_buf);

    glPixelStorei(GL_PACK_ALIGNMENT, 1);

    glReadPixels(0, 0, num_cols, num_rows, GL_RGB, GL_FLOAT, img_buf);

    create_image_f(img_out, num_rows, num_cols);
    Image_f* img = *img_out;

    for (uint32_t row = 0; row < num_rows; row++)
    {
        for (uint32_t col = 0; col < 3*num_cols; col += 3)
        {
            uint32_t r = num_rows - row - 1;

            img->pxls[ r ][ col/3 ].r = img_buf[ col+0 + row*3*num_cols ];
            img->pxls[ r ][ col/3 ].g = img_buf[ col+1 + row*3*num_cols ];
            img->pxls[ r ][ col/3 ].b = img_buf[ col+2 + row*3*num_cols ];
        }
    }

    free(img_buf);
}


void Camera_d::capture_gl_view(const char* fname) throw (jwscxx::base::IO_error)
{
    Image_f* img = 0;
    capture_gl_view(&img);

    Error* err;
    if ((err = write_image_f(img, fname)) != 0)
    {
        std::ostringstream ost(err->msg);
        free_error(err);
        throw IO_error(ost.str());
    }

    free_image_f(img);
}


void Camera_d::capture_gl_view(const char* fname_fmt, uint32_t N)
    throw (jwscxx::base::IO_error)
{
    char fname[256] = {0};

    snprintf(fname, 255, fname_fmt, N);

    capture_gl_view(fname);
}


void Camera_d::read(std::istream& in) throw (jwscxx::base::IO_error,
        jwscxx::base::Arg_error)
{
    using std::ostringstream;
    using std::istringstream;

    const char* type_name = typeid(*this).name();
    const char* field_value;

    // Type
    if (!(field_value = read_field_value(in, "Type")))
    {
        throw Arg_error("Missing Type field");
    }
    if (strncmp(field_value, type_name, strlen(type_name)) != 0)
    {
        ostringstream ost;
        ost << "Tried to read a '" << field_value << "' as a '"
            << type_name << "'";
        throw Arg_error(ost.str());
    }


    // VRP
    if (!(field_value = read_field_value(in, "VRP")))
    {
        throw Arg_error("Missing VRP fields");
    }
    istringstream ist(field_value);
    ist >> vrp->elts[0] >> vrp->elts[1] >> vrp->elts[2];
    if (ist.fail())
    {
        throw Arg_error("Invalid VRP");
    }
    ist.clear(std::ios_base::goodbit);


    // PRP
    if (!(field_value = read_field_value(in, "PRP")))
    {
        throw Arg_error("Missing PRP fields");
    }
    ist.str(field_value);
    ist >> prp->elts[0] >> prp->elts[1] >> prp->elts[2];
    if (ist.fail())
    {
        throw Arg_error("Invalid PRP");
    }
    ist.clear(std::ios_base::goodbit);


    // U
    if (!(field_value = read_field_value(in, "U")))
    {
        throw Arg_error("Missing U fields");
    }
    ist.str(field_value);
    ist >> u->elts[ 0 ] >> u->elts[1] >> u->elts[2];
    if (ist.fail())
    {
        throw Arg_error("Invalid U");
    }
    ist.clear(std::ios_base::goodbit);


    // V
    if (!(field_value = read_field_value(in, "V")))
    {
        throw Arg_error("Missing V fields");
    }
    ist.str(field_value);
    ist >> v->elts[ 0 ] >> v->elts[1] >> v->elts[2];
    if (ist.fail())
    {
        throw Arg_error("Invalid V");
    }
    ist.clear(std::ios_base::goodbit);


    // N
    if (!(field_value = read_field_value(in, "N")))
    {
        throw Arg_error("Missing N fields");
    }
    ist.str(field_value);
    ist >> n->elts[ 0 ] >> n->elts[1] >> n->elts[2];
    if (ist.fail())
    {
        throw Arg_error("Invalid N");
    }
    ist.clear(std::ios_base::goodbit);


    // F
    if (!(field_value = read_field_value(in, "F")))
    {
        throw Arg_error("Missing F field");
    }
    ist.str(field_value);
    ist >> focal_length;
    if (ist.fail())
    {
        throw Arg_error("Invalid focal length");
    }
    ist.clear(std::ios_base::goodbit);


    // Aperture
    if (!(field_value = read_field_value(in, "Aperture")))
    {
        throw Arg_error("Missing Aperture field");
    }
    ist.str(field_value);
    ist >> aperture;
    if (ist.fail())
    {
        throw Arg_error("Invalid aperture");
    }
    ist.clear(std::ios_base::goodbit);
    if (aperture <= 0)
    {
        throw Arg_error("Aperture <= 0");
    }


    // Near
    if (!(field_value = read_field_value(in, "Near")))
    {
        throw Arg_error("Missing Near field");
    }
    ist.str(field_value);
    ist >> near;
    if (ist.fail())
    {
        throw Arg_error("Invalid near clipping distance");
    }
    ist.clear(std::ios_base::goodbit);


    // Far
    if (!(field_value = read_field_value(in, "Far")))
    {
        throw Arg_error("Missing Far field");
    }
    ist.str(field_value);
    ist >> far;
    if (ist.fail())
    {
        throw Arg_error("Invalid far clipping distance");
    }
    ist.clear(std::ios_base::goodbit);
    if (near < 0 || far <= 0 || near > far)
    {
        throw Arg_error("Invalid clipping planes");
    }


    // Derived...
    double tan_theta = std::tan(JWSC_DEG_TO_RAD * aperture / 2.0);
    half_near_height = near * tan_theta;
    half_focal_height = focal_length * tan_theta;
}


void Camera_d::write(std::ostream& out) const throw (jwscxx::base::IO_error)
{
    out << "     Type: " << typeid(*this).name() << '\n'
        << "      VRP: " << vrp->elts[0] << ' '
                         << vrp->elts[1] << ' '
                         << vrp->elts[2] << '\n'
        << "      PRP: " << prp->elts[0] << ' '
                         << prp->elts[1] << ' '
                         << prp->elts[2] << '\n'
        << "        U: " << u->elts[0] << ' '
                         << u->elts[1] << ' '
                         << u->elts[2] << '\n'
        << "        V: " << v->elts[0] << ' '
                         << v->elts[1] << ' '
                         << v->elts[2] << '\n'
        << "        N: " << n->elts[0] << ' '
                         << n->elts[1] << ' '
                         << n->elts[2] << '\n'
        << "        F: " << focal_length << '\n'
        << " Aperture: " << aperture << '\n'
        << "     Near: " << near << '\n'
        << "      Far: " << far << '\n';
}





Stereo_camera_f::Stereo_camera_f
(
    float focal_pt_x,
    float focal_pt_y,
    float focal_pt_z,
    float ref_pt_x,
    float ref_pt_y,
    float ref_pt_z,
    float up_x,
    float up_y,
    float up_z,
    float aperture, 
    float clip_near, 
    float clip_far
)
throw (Arg_error)
    : Camera_f(focal_pt_x, focal_pt_y, focal_pt_z, ref_pt_x, ref_pt_y, ref_pt_z,
            up_x, up_y, up_z, aperture, clip_near, clip_far)
{
    left_prp  = 0;
    left_vrp  = 0;
    right_prp = 0;
    right_vrp = 0;

    look_at(focal_pt_x, focal_pt_y, focal_pt_z, ref_pt_x, ref_pt_y, ref_pt_z,
            up_x, up_y, up_z);
}


Stereo_camera_f::Stereo_camera_f
(
    const jwsc::Vector_f* focal_pt,
    const jwsc::Vector_f* ref_pt,
    const jwsc::Vector_f* up,
    float                 aperture, 
    float                 clip_near, 
    float                 clip_far
)
throw (Arg_error)
    : Camera_f(focal_pt, ref_pt, up, aperture, clip_near, clip_far)
{
    left_prp  = 0;
    left_vrp  = 0;
    right_prp = 0;
    right_vrp = 0;

    Camera_f::look_at(focal_pt, ref_pt, up);
}


Stereo_camera_f::Stereo_camera_f(const Stereo_camera_f& c)
    : Camera_f(c)
{
    left_prp  = 0;
    right_prp = 0;
    left_vrp  = 0;
    right_vrp = 0;

    *this = c;
}


Stereo_camera_f::Stereo_camera_f(const char* fname) 
    throw (jwscxx::base::Arg_error, jwscxx::base::IO_error)
    : Camera_f(0,0,0, 0,0,-1, 0,1,0, 60.0f, 0, 1)
{
    left_prp  = 0;
    right_prp = 0;
    left_vrp  = 0;
    right_vrp = 0;

    jwscxx::base::Readable::read(fname);
}


Stereo_camera_f::Stereo_camera_f(std::istream& in) 
    throw (jwscxx::base::Arg_error, jwscxx::base::IO_error)
    : Camera_f(0,0,0, 0,0,-1, 0,1,0, 60.0f, 0, 1)
{
    left_prp  = 0;
    right_prp = 0;
    left_vrp  = 0;
    right_vrp = 0;

    read(in);
}


Stereo_camera_f::~Stereo_camera_f()
{
    free_vector_f(left_vrp);
    free_vector_f(left_prp);
    free_vector_f(right_vrp);
    free_vector_f(right_prp);
}


Stereo_camera_f& Stereo_camera_f::operator=(const Stereo_camera_f &c)
{
    Camera_f::operator=(c);

    copy_vector_f(&left_vrp, c.left_vrp);
    copy_vector_f(&right_vrp, c.right_vrp);
    copy_vector_f(&left_prp, c.left_prp);
    copy_vector_f(&right_prp, c.right_prp);

    separation = c.separation;
    ndfl = c.ndfl;

    return *this;
}


Stereo_camera_f* Stereo_camera_f::clone() const
{
    return new Stereo_camera_f(*this);
}


const jwsc::Vector_f* Stereo_camera_f::get_left_prp() const 
{ 
    return (const jwsc::Vector_f*)left_prp; 
}


const jwsc::Vector_f* Stereo_camera_f::get_right_prp() const 
{ 
    return (const jwsc::Vector_f*)right_prp; 
}


const jwsc::Vector_f* Stereo_camera_f::get_left_vrp() const 
{ 
    return (const jwsc::Vector_f*)left_vrp; 
}


const jwsc::Vector_f* Stereo_camera_f::get_right_vrp() const 
{ 
    return (const jwsc::Vector_f*)right_vrp; 
}


void Stereo_camera_f::set_aperture(float aperture) throw (Arg_error)
{
    Camera_f::set_aperture(aperture);

    ndfl = near / focal_length;
}


void Stereo_camera_f::set_f_with_new_aperture_using_prp(float f) 
    throw (Arg_error)
{
    Camera_f::set_f_with_new_aperture_using_prp(f);

    separation = focal_length / 30.0f;
    ndfl = near / focal_length;

    assert(copy_vector_f(&left_prp, prp) == 0);
    assert(copy_vector_f(&left_vrp, vrp) == 0);

    assert(copy_vector_f(&right_prp, prp) == 0);
    assert(copy_vector_f(&right_vrp, vrp) == 0);

    float t_x = 0.5f*separation*u->elts[0];
    float t_y = 0.5f*separation*u->elts[1];
    float t_z = 0.5f*separation*u->elts[2];

    left_prp->elts[0] -= t_x;  left_vrp->elts[0] -= t_x;
    left_prp->elts[1] -= t_y;  left_vrp->elts[1] -= t_y;
    left_prp->elts[2] -= t_z;  left_vrp->elts[2] -= t_z;

    right_prp->elts[0] += t_x;  right_vrp->elts[0] += t_x;
    right_prp->elts[1] += t_y;  right_vrp->elts[1] += t_y;
    right_prp->elts[2] += t_z;  right_vrp->elts[2] += t_z;
}


void Stereo_camera_f::set_f_with_same_aperture_using_prp(float f) 
    throw (Arg_error)
{
    Camera_f::set_f_with_same_aperture_using_prp(f);

    separation = focal_length / 30.0f;
    ndfl = near / focal_length;

    assert(copy_vector_f(&left_prp, prp) == 0);
    assert(copy_vector_f(&left_vrp, vrp) == 0);

    assert(copy_vector_f(&right_prp, prp) == 0);
    assert(copy_vector_f(&right_vrp, vrp) == 0);

    float t_x = 0.5f*separation*u->elts[0];
    float t_y = 0.5f*separation*u->elts[1];
    float t_z = 0.5f*separation*u->elts[2];

    left_prp->elts[0] -= t_x;  left_vrp->elts[0] -= t_x;
    left_prp->elts[1] -= t_y;  left_vrp->elts[1] -= t_y;
    left_prp->elts[2] -= t_z;  left_vrp->elts[2] -= t_z;

    right_prp->elts[0] += t_x;  right_vrp->elts[0] += t_x;
    right_prp->elts[1] += t_y;  right_vrp->elts[1] += t_y;
    right_prp->elts[2] += t_z;  right_vrp->elts[2] += t_z;
}


void Stereo_camera_f::set_f_with_new_aperture_using_vrp(float f) 
    throw (Arg_error)
{
    Camera_f::set_f_with_new_aperture_using_vrp(f);

    separation = focal_length / 30.0f;
    ndfl = near / focal_length;

    assert(copy_vector_f(&left_prp, prp) == 0);
    assert(copy_vector_f(&left_vrp, vrp) == 0);

    assert(copy_vector_f(&right_prp, prp) == 0);
    assert(copy_vector_f(&right_vrp, vrp) == 0);

    float t_x = 0.5f*separation*u->elts[0];
    float t_y = 0.5f*separation*u->elts[1];
    float t_z = 0.5f*separation*u->elts[2];

    left_prp->elts[0] -= t_x;  left_vrp->elts[0] -= t_x;
    left_prp->elts[1] -= t_y;  left_vrp->elts[1] -= t_y;
    left_prp->elts[2] -= t_z;  left_vrp->elts[2] -= t_z;

    right_prp->elts[0] += t_x;  right_vrp->elts[0] += t_x;
    right_prp->elts[1] += t_y;  right_vrp->elts[1] += t_y;
    right_prp->elts[2] += t_z;  right_vrp->elts[2] += t_z;
}


void Stereo_camera_f::set_f_with_same_aperture_using_vrp(float f) 
    throw (Arg_error)
{
    Camera_f::set_f_with_same_aperture_using_vrp(f);

    separation = focal_length / 30.0f;
    ndfl = near / focal_length;

    assert(copy_vector_f(&left_prp, prp) == 0);
    assert(copy_vector_f(&left_vrp, vrp) == 0);

    assert(copy_vector_f(&right_prp, prp) == 0);
    assert(copy_vector_f(&right_vrp, vrp) == 0);

    float t_x = 0.5f*separation*u->elts[0];
    float t_y = 0.5f*separation*u->elts[1];
    float t_z = 0.5f*separation*u->elts[2];

    left_prp->elts[0] -= t_x;  left_vrp->elts[0] -= t_x;
    left_prp->elts[1] -= t_y;  left_vrp->elts[1] -= t_y;
    left_prp->elts[2] -= t_z;  left_vrp->elts[2] -= t_z;

    right_prp->elts[0] += t_x;  right_vrp->elts[0] += t_x;
    right_prp->elts[1] += t_y;  right_vrp->elts[1] += t_y;
    right_prp->elts[2] += t_z;  right_vrp->elts[2] += t_z;
}


void Stereo_camera_f::set_clipping(float clip_near, float clip_far) 
    throw (Arg_error)
{
    Camera_f::set_clipping(clip_near, clip_far);

    ndfl = near / focal_length;
}


void Stereo_camera_f::rotate(float phi, float x, float y, float z) 
    throw (Arg_error)
{
    Camera_f::rotate(phi, x, y, z);

    multiply_matrix_and_vector_f(&left_vrp, R, left_vrp);
    multiply_matrix_and_vector_f(&left_prp, R, left_prp);
    multiply_matrix_and_vector_f(&right_vrp, R, right_vrp);
    multiply_matrix_and_vector_f(&right_prp, R, right_prp);
}


void Stereo_camera_f::rotate(float angle, const jwsc::Vector_f* r) 
    throw (jwscxx::base::Arg_error) 
{ 
    Camera_f::rotate(angle, r); 
}


void Stereo_camera_f::translate(float x, float y, float z)
{
    Camera_f::translate(x, y, z);

    left_prp->elts[0] += x;  left_vrp->elts[0] += x;
    left_prp->elts[1] += y;  left_vrp->elts[1] += y;
    left_prp->elts[2] += z;  left_vrp->elts[2] += z;

    right_prp->elts[0] += x;  right_vrp->elts[0] += x;
    right_prp->elts[1] += y;  right_vrp->elts[1] += y;
    right_prp->elts[2] += z;  right_vrp->elts[2] += z;
}


void Stereo_camera_f::translate(const jwsc::Vector_f* t) throw (Arg_error)
{
    Camera_f::translate(t);

    add_vectors_f(&left_vrp, left_vrp, t);
    add_vectors_f(&left_prp, left_prp, t);

    add_vectors_f(&right_vrp, right_vrp, t);
    add_vectors_f(&right_prp, right_prp, t);
}


void Stereo_camera_f::look_at
(
    float focal_pt_x,
    float focal_pt_y,
    float focal_pt_z,
    float ref_pt_x,
    float ref_pt_y,
    float ref_pt_z,
    float up_x,
    float up_y,
    float up_z
)
throw (jwscxx::base::Arg_error)
{
    Camera_f::look_at(focal_pt_x, focal_pt_y, focal_pt_z,
            ref_pt_x, ref_pt_y, ref_pt_z, up_x, up_y, up_z);

    separation = focal_length / 30.0f;
    ndfl = near / focal_length;

    copy_vector_f(&left_prp, prp);
    copy_vector_f(&left_vrp, vrp);

    copy_vector_f(&right_prp, prp);
    copy_vector_f(&right_vrp, vrp);

    float t_x = 0.5f*separation*u->elts[0];
    float t_y = 0.5f*separation*u->elts[1];
    float t_z = 0.5f*separation*u->elts[2];

    left_prp->elts[0] -= t_x;  left_vrp->elts[0] -= t_x;
    left_prp->elts[1] -= t_y;  left_vrp->elts[1] -= t_y;
    left_prp->elts[2] -= t_z;  left_vrp->elts[2] -= t_z;

    right_prp->elts[0] += t_x;  right_vrp->elts[0] += t_x;
    right_prp->elts[1] += t_y;  right_vrp->elts[1] += t_y;
    right_prp->elts[2] += t_z;  right_vrp->elts[2] += t_z;
}


void Stereo_camera_f::set_left_gl_modelview()
{
    GLfloat M[16] = {0};

    glMatrixMode(GL_MODELVIEW);

    glLoadIdentity();

    M[0] = u->elts[0];  M[4] = u->elts[1];  M[8]  = u->elts[2];
    M[1] = v->elts[0];  M[5] = v->elts[1];  M[9]  = v->elts[2];
    M[2] = n->elts[0];  M[6] = n->elts[1];  M[10] = n->elts[2];
    M[15] = 1;
    glMultMatrixf(M);

    glTranslatef(-left_prp->elts[0], -left_prp->elts[1], -left_prp->elts[2]);
}


void Stereo_camera_f::set_right_gl_modelview()
{
    GLfloat M[16] = {0};

    glMatrixMode(GL_MODELVIEW);

    glLoadIdentity();

    M[0] = u->elts[0];  M[4] = u->elts[1];  M[8]  = u->elts[2];
    M[1] = v->elts[0];  M[5] = v->elts[1];  M[9]  = v->elts[2];
    M[2] = n->elts[0];  M[6] = n->elts[1];  M[10] = n->elts[2];
    M[15] = 1;
    glMultMatrixf(M);

    glTranslatef(-right_prp->elts[0], -right_prp->elts[1], -right_prp->elts[2]);
}


void Stereo_camera_f::set_left_gl_projection()
{
    glMatrixMode(GL_PROJECTION);

    glLoadIdentity();

    float aspect = get_gl_viewport_aspect();
    float left   = -aspect*half_near_height + 0.5f*separation*ndfl;
    float right  =  aspect*half_near_height + 0.5f*separation*ndfl;
    float top    =  half_near_height;
    float bottom = -half_near_height;

    glFrustum(left, right, bottom, top, near, far);
}


void Stereo_camera_f::set_right_gl_projection()
{
    glMatrixMode(GL_PROJECTION);

    glLoadIdentity();

    float aspect = get_gl_viewport_aspect();
    float left   = -aspect*half_near_height - 0.5f*separation*ndfl;
    float right  =  aspect*half_near_height - 0.5f*separation*ndfl;
    float top    =  half_near_height;
    float bottom = -half_near_height;

    glFrustum(left, right, bottom, top, near, far);
}


void Stereo_camera_f::read(std::istream& in) throw (jwscxx::base::IO_error,
        jwscxx::base::Arg_error)
{
    using std::ostringstream;
    using std::istringstream;

    const char* field_value;

    Camera_f::read(in);


    // Left VRP
    if (!(field_value = read_field_value(in, "Left VRP")))
    {
        throw Arg_error("Missing Left VRP fields");
    }
    istringstream ist(field_value);
    ist >> left_vrp->elts[0] >> left_vrp->elts[1] >> left_vrp->elts[2];
    if (ist.fail())
    {
        throw Arg_error("Invalid Left VRP");
    }
    ist.clear(std::ios_base::goodbit);


    // Right VRP
    if (!(field_value = read_field_value(in, "Right VRP")))
    {
        throw Arg_error("Missing Right VRP fields");
    }
    ist.str(field_value);
    ist >> right_vrp->elts[0] >> right_vrp->elts[1] >> right_vrp->elts[2];
    if (ist.fail())
    {
        throw Arg_error("Invalid Right VRP");
    }
    ist.clear(std::ios_base::goodbit);


    // Left PRP
    if (!(field_value = read_field_value(in, "Left PRP")))
    {
        throw Arg_error("Missing Left PRP fields");
    }
    ist.str(field_value);
    ist >> left_prp->elts[0] >> left_prp->elts[1] >> left_prp->elts[2];
    if (ist.fail())
    {
        throw Arg_error("Invalid Left PRP");
    }
    ist.clear(std::ios_base::goodbit);


    // Right PRP
    if (!(field_value = read_field_value(in, "Right PRP")))
    {
        throw Arg_error("Missing Right PRP fields");
    }
    ist.str(field_value);
    ist >> right_prp->elts[0] >> right_prp->elts[1] >> right_prp->elts[2];
    if (ist.fail())
    {
        throw Arg_error("Invalid Right PRP");
    }
    ist.clear(std::ios_base::goodbit);


    // Derived...
    separation = focal_length / 30.0f;
    ndfl = near / focal_length;
}


void Stereo_camera_f::write(std::ostream& out) const 
    throw (jwscxx::base::IO_error)
{
    Camera_f::write(out);

    out << "   Left VRP: " << left_vrp->elts[0] << ' '
                           << left_vrp->elts[1] << ' '
                           << left_vrp->elts[2] << '\n'
        << "  Right VRP: " << right_vrp->elts[0] << ' '
                           << right_vrp->elts[1] << ' '
                           << right_vrp->elts[2] << '\n'
        << "   Left PRP: " << left_prp->elts[0] << ' '
                           << left_prp->elts[1] << ' '
                           << left_prp->elts[2] << '\n'
        << "  Right PRP: " << right_prp->elts[0] << ' '
                           << right_prp->elts[1] << ' '
                           << right_prp->elts[2] << '\n';
}



Stereo_camera_d::Stereo_camera_d
(
    double focal_pt_x,
    double focal_pt_y,
    double focal_pt_z,
    double ref_pt_x,
    double ref_pt_y,
    double ref_pt_z,
    double up_x,
    double up_y,
    double up_z,
    double aperture, 
    double clip_near, 
    double clip_far
)
throw (Arg_error)
    : Camera_d(focal_pt_x, focal_pt_y, focal_pt_z, ref_pt_x, ref_pt_y, ref_pt_z,
            up_x, up_y, up_z, aperture, clip_near, clip_far)
{
    left_prp  = 0;
    left_vrp  = 0;
    right_prp = 0;
    right_vrp = 0;

    look_at(focal_pt_x, focal_pt_y, focal_pt_z, ref_pt_x, ref_pt_y, ref_pt_z,
            up_x, up_y, up_z);
}


Stereo_camera_d::Stereo_camera_d
(
    const jwsc::Vector_d* focal_pt,
    const jwsc::Vector_d* ref_pt,
    const jwsc::Vector_d* up,
    double                aperture, 
    double                clip_near, 
    double                clip_far
)
throw (Arg_error)
    : Camera_d(focal_pt, ref_pt, up, aperture, clip_near, clip_far)
{
    left_prp  = 0;
    left_vrp  = 0;
    right_prp = 0;
    right_vrp = 0;

    Camera_d::look_at(focal_pt, ref_pt, up);
}


Stereo_camera_d::Stereo_camera_d(const Stereo_camera_d& c)
    : Camera_d(c)
{
    left_prp  = 0;
    right_prp = 0;
    left_vrp  = 0;
    right_vrp = 0;

    *this = c;
}


Stereo_camera_d::Stereo_camera_d(const char* fname) 
    throw (jwscxx::base::Arg_error, jwscxx::base::IO_error)
    : Camera_d(0,0,0, 0,0,-1, 0,1,0, 60.0, 0, 1)
{
    left_prp  = 0;
    right_prp = 0;
    left_vrp  = 0;
    right_vrp = 0;

    jwscxx::base::Readable::read(fname);
}


Stereo_camera_d::Stereo_camera_d(std::istream& in) 
    throw (jwscxx::base::Arg_error, jwscxx::base::IO_error)
    : Camera_d(0,0,0, 0,0,-1, 0,1,0, 60.0, 0, 1)
{
    left_prp  = 0;
    right_prp = 0;
    left_vrp  = 0;
    right_vrp = 0;

    read(in);
}


Stereo_camera_d::~Stereo_camera_d()
{
    free_vector_d(left_vrp);
    free_vector_d(left_prp);
    free_vector_d(right_vrp);
    free_vector_d(right_prp);
}


Stereo_camera_d& Stereo_camera_d::operator=(const Stereo_camera_d &c)
{
    Camera_d::operator=(c);

    copy_vector_d(&left_vrp, c.left_vrp);
    copy_vector_d(&right_vrp, c.right_vrp);
    copy_vector_d(&left_prp, c.left_prp);
    copy_vector_d(&right_prp, c.right_prp);

    separation = c.separation;
    ndfl = c.ndfl;

    return *this;
}


Stereo_camera_d* Stereo_camera_d::clone() const
{
    return new Stereo_camera_d(*this);
}


const jwsc::Vector_d* Stereo_camera_d::get_left_prp() const 
{ 
    return (const jwsc::Vector_d*)left_prp; 
}


const jwsc::Vector_d* Stereo_camera_d::get_right_prp() const 
{ 
    return (const jwsc::Vector_d*)right_prp; 
}


const jwsc::Vector_d* Stereo_camera_d::get_left_vrp() const 
{ 
    return (const jwsc::Vector_d*)left_vrp; 
}


const jwsc::Vector_d* Stereo_camera_d::get_right_vrp() const 
{ 
    return (const jwsc::Vector_d*)right_vrp; 
}


void Stereo_camera_d::set_aperture(double aperture) throw (Arg_error)
{
    Camera_d::set_aperture(aperture);

    ndfl = near / focal_length;
}


void Stereo_camera_d::set_f_with_new_aperture_using_prp(double f) 
    throw (Arg_error)
{
    Camera_d::set_f_with_new_aperture_using_prp(f);

    separation = focal_length / 30.0;
    ndfl = near / focal_length;

    assert(copy_vector_d(&left_prp, prp) == 0);
    assert(copy_vector_d(&left_vrp, vrp) == 0);

    assert(copy_vector_d(&right_prp, prp) == 0);
    assert(copy_vector_d(&right_vrp, vrp) == 0);

    double t_x = 0.5*separation*u->elts[0];
    double t_y = 0.5*separation*u->elts[1];
    double t_z = 0.5*separation*u->elts[2];

    left_prp->elts[0] -= t_x;  left_vrp->elts[0] -= t_x;
    left_prp->elts[1] -= t_y;  left_vrp->elts[1] -= t_y;
    left_prp->elts[2] -= t_z;  left_vrp->elts[2] -= t_z;

    right_prp->elts[0] += t_x;  right_vrp->elts[0] += t_x;
    right_prp->elts[1] += t_y;  right_vrp->elts[1] += t_y;
    right_prp->elts[2] += t_z;  right_vrp->elts[2] += t_z;
}


void Stereo_camera_d::set_f_with_same_aperture_using_prp(double f) 
    throw (Arg_error)
{
    Camera_d::set_f_with_same_aperture_using_prp(f);

    separation = focal_length / 30.0;
    ndfl = near / focal_length;

    assert(copy_vector_d(&left_prp, prp) == 0);
    assert(copy_vector_d(&left_vrp, vrp) == 0);

    assert(copy_vector_d(&right_prp, prp) == 0);
    assert(copy_vector_d(&right_vrp, vrp) == 0);

    double t_x = 0.5*separation*u->elts[0];
    double t_y = 0.5*separation*u->elts[1];
    double t_z = 0.5*separation*u->elts[2];

    left_prp->elts[0] -= t_x;  left_vrp->elts[0] -= t_x;
    left_prp->elts[1] -= t_y;  left_vrp->elts[1] -= t_y;
    left_prp->elts[2] -= t_z;  left_vrp->elts[2] -= t_z;

    right_prp->elts[0] += t_x;  right_vrp->elts[0] += t_x;
    right_prp->elts[1] += t_y;  right_vrp->elts[1] += t_y;
    right_prp->elts[2] += t_z;  right_vrp->elts[2] += t_z;
}


void Stereo_camera_d::set_f_with_new_aperture_using_vrp(double f) 
    throw (Arg_error)
{
    Camera_d::set_f_with_new_aperture_using_vrp(f);

    separation = focal_length / 30.0;
    ndfl = near / focal_length;

    assert(copy_vector_d(&left_prp, prp) == 0);
    assert(copy_vector_d(&left_vrp, vrp) == 0);

    assert(copy_vector_d(&right_prp, prp) == 0);
    assert(copy_vector_d(&right_vrp, vrp) == 0);

    double t_x = 0.5*separation*u->elts[0];
    double t_y = 0.5*separation*u->elts[1];
    double t_z = 0.5*separation*u->elts[2];

    left_prp->elts[0] -= t_x;  left_vrp->elts[0] -= t_x;
    left_prp->elts[1] -= t_y;  left_vrp->elts[1] -= t_y;
    left_prp->elts[2] -= t_z;  left_vrp->elts[2] -= t_z;

    right_prp->elts[0] += t_x;  right_vrp->elts[0] += t_x;
    right_prp->elts[1] += t_y;  right_vrp->elts[1] += t_y;
    right_prp->elts[2] += t_z;  right_vrp->elts[2] += t_z;
}


void Stereo_camera_d::set_f_with_same_aperture_using_vrp(double f) 
    throw (Arg_error)
{
    Camera_d::set_f_with_same_aperture_using_vrp(f);

    separation = focal_length / 30.0;
    ndfl = near / focal_length;

    assert(copy_vector_d(&left_prp, prp) == 0);
    assert(copy_vector_d(&left_vrp, vrp) == 0);

    assert(copy_vector_d(&right_prp, prp) == 0);
    assert(copy_vector_d(&right_vrp, vrp) == 0);

    double t_x = 0.5*separation*u->elts[0];
    double t_y = 0.5*separation*u->elts[1];
    double t_z = 0.5*separation*u->elts[2];

    left_prp->elts[0] -= t_x;  left_vrp->elts[0] -= t_x;
    left_prp->elts[1] -= t_y;  left_vrp->elts[1] -= t_y;
    left_prp->elts[2] -= t_z;  left_vrp->elts[2] -= t_z;

    right_prp->elts[0] += t_x;  right_vrp->elts[0] += t_x;
    right_prp->elts[1] += t_y;  right_vrp->elts[1] += t_y;
    right_prp->elts[2] += t_z;  right_vrp->elts[2] += t_z;
}


void Stereo_camera_d::set_clipping(double clip_near, double clip_far) 
    throw (Arg_error)
{
    Camera_d::set_clipping(clip_near, clip_far);

    ndfl = near / focal_length;
}


void Stereo_camera_d::rotate(double phi, double x, double y, double z) 
    throw (Arg_error)
{
    Camera_d::rotate(phi, x, y, z);

    multiply_matrix_and_vector_d(&left_vrp, R, left_vrp);
    multiply_matrix_and_vector_d(&left_prp, R, left_prp);
    multiply_matrix_and_vector_d(&right_vrp, R, right_vrp);
    multiply_matrix_and_vector_d(&right_prp, R, right_prp);
}


void Stereo_camera_d::rotate(double angle, const jwsc::Vector_d* r) 
    throw (jwscxx::base::Arg_error) 
{ 
    Camera_d::rotate(angle, r); 
}


void Stereo_camera_d::translate(double x, double y, double z)
{
    Camera_d::translate(x, y, z);

    left_prp->elts[0] += x;  left_vrp->elts[0] += x;
    left_prp->elts[1] += y;  left_vrp->elts[1] += y;
    left_prp->elts[2] += z;  left_vrp->elts[2] += z;

    right_prp->elts[0] += x;  right_vrp->elts[0] += x;
    right_prp->elts[1] += y;  right_vrp->elts[1] += y;
    right_prp->elts[2] += z;  right_vrp->elts[2] += z;
}


void Stereo_camera_d::translate(const jwsc::Vector_d* t) throw (Arg_error)
{
    Camera_d::translate(t);

    add_vectors_d(&left_vrp, left_vrp, t);
    add_vectors_d(&left_prp, left_prp, t);

    add_vectors_d(&right_vrp, right_vrp, t);
    add_vectors_d(&right_prp, right_prp, t);
}


void Stereo_camera_d::look_at
(
    double focal_pt_x,
    double focal_pt_y,
    double focal_pt_z,
    double ref_pt_x,
    double ref_pt_y,
    double ref_pt_z,
    double up_x,
    double up_y,
    double up_z
)
throw (jwscxx::base::Arg_error)
{
    Camera_d::look_at(focal_pt_x, focal_pt_y, focal_pt_z,
            ref_pt_x, ref_pt_y, ref_pt_z, up_x, up_y, up_z);

    separation = focal_length / 30.0;
    ndfl = near / focal_length;

    copy_vector_d(&left_prp, prp);
    copy_vector_d(&left_vrp, vrp);

    copy_vector_d(&right_prp, prp);
    copy_vector_d(&right_vrp, vrp);

    double t_x = 0.5*separation*u->elts[0];
    double t_y = 0.5*separation*u->elts[1];
    double t_z = 0.5*separation*u->elts[2];

    left_prp->elts[0] -= t_x;  left_vrp->elts[0] -= t_x;
    left_prp->elts[1] -= t_y;  left_vrp->elts[1] -= t_y;
    left_prp->elts[2] -= t_z;  left_vrp->elts[2] -= t_z;

    right_prp->elts[0] += t_x;  right_vrp->elts[0] += t_x;
    right_prp->elts[1] += t_y;  right_vrp->elts[1] += t_y;
    right_prp->elts[2] += t_z;  right_vrp->elts[2] += t_z;
}


void Stereo_camera_d::set_left_gl_modelview()
{
    GLdouble M[16] = {0};

    glMatrixMode(GL_MODELVIEW);

    glLoadIdentity();

    M[0] = u->elts[0];  M[4] = u->elts[1];  M[8]  = u->elts[2];
    M[1] = v->elts[0];  M[5] = v->elts[1];  M[9]  = v->elts[2];
    M[2] = n->elts[0];  M[6] = n->elts[1];  M[10] = n->elts[2];
    M[15] = 1;
    glMultMatrixd(M);

    glTranslated(-left_prp->elts[0], -left_prp->elts[1], -left_prp->elts[2]);
}


void Stereo_camera_d::set_right_gl_modelview()
{
    GLdouble M[16] = {0};

    glMatrixMode(GL_MODELVIEW);

    glLoadIdentity();

    M[0] = u->elts[0];  M[4] = u->elts[1];  M[8]  = u->elts[2];
    M[1] = v->elts[0];  M[5] = v->elts[1];  M[9]  = v->elts[2];
    M[2] = n->elts[0];  M[6] = n->elts[1];  M[10] = n->elts[2];
    M[15] = 1;
    glMultMatrixd(M);

    glTranslated(-right_prp->elts[0], -right_prp->elts[1], -right_prp->elts[2]);
}


void Stereo_camera_d::set_left_gl_projection()
{
    glMatrixMode(GL_PROJECTION);

    glLoadIdentity();

    double aspect = get_gl_viewport_aspect();
    double left   = -aspect*half_near_height + 0.5*separation*ndfl;
    double right  =  aspect*half_near_height + 0.5*separation*ndfl;
    double top    =  half_near_height;
    double bottom = -half_near_height;

    glFrustum(left, right, bottom, top, near, far);
}


void Stereo_camera_d::set_right_gl_projection()
{
    glMatrixMode(GL_PROJECTION);

    glLoadIdentity();

    double aspect = get_gl_viewport_aspect();
    double left   = -aspect*half_near_height - 0.5*separation*ndfl;
    double right  =  aspect*half_near_height - 0.5*separation*ndfl;
    double top    =  half_near_height;
    double bottom = -half_near_height;

    glFrustum(left, right, bottom, top, near, far);
}


void Stereo_camera_d::read(std::istream& in) throw (jwscxx::base::IO_error,
        jwscxx::base::Arg_error)
{
    using std::ostringstream;
    using std::istringstream;

    const char* field_value;

    Camera_d::read(in);


    // Left VRP
    if (!(field_value = read_field_value(in, "Left VRP")))
    {
        throw Arg_error("Missing Left VRP fields");
    }
    istringstream ist(field_value);
    ist >> left_vrp->elts[0] >> left_vrp->elts[1] >> left_vrp->elts[2];
    if (ist.fail())
    {
        throw Arg_error("Invalid Left VRP");
    }
    ist.clear(std::ios_base::goodbit);


    // Right VRP
    if (!(field_value = read_field_value(in, "Right VRP")))
    {
        throw Arg_error("Missing Right VRP fields");
    }
    ist.str(field_value);
    ist >> right_vrp->elts[0] >> right_vrp->elts[1] >> right_vrp->elts[2];
    if (ist.fail())
    {
        throw Arg_error("Invalid Right VRP");
    }
    ist.clear(std::ios_base::goodbit);


    // Left PRP
    if (!(field_value = read_field_value(in, "Left PRP")))
    {
        throw Arg_error("Missing Left PRP fields");
    }
    ist.str(field_value);
    ist >> left_prp->elts[0] >> left_prp->elts[1] >> left_prp->elts[2];
    if (ist.fail())
    {
        throw Arg_error("Invalid Left PRP");
    }
    ist.clear(std::ios_base::goodbit);


    // Right PRP
    if (!(field_value = read_field_value(in, "Right PRP")))
    {
        throw Arg_error("Missing Right PRP fields");
    }
    ist.str(field_value);
    ist >> right_prp->elts[0] >> right_prp->elts[1] >> right_prp->elts[2];
    if (ist.fail())
    {
        throw Arg_error("Invalid Right PRP");
    }
    ist.clear(std::ios_base::goodbit);


    // Derived...
    separation = focal_length / 30.0;
    ndfl = near / focal_length;
}


void Stereo_camera_d::write(std::ostream& out) const 
    throw (jwscxx::base::IO_error)
{
    Camera_d::write(out);

    out << "   Left VRP: " << left_vrp->elts[0] << ' '
                           << left_vrp->elts[1] << ' '
                           << left_vrp->elts[2] << '\n'
        << "  Right VRP: " << right_vrp->elts[0] << ' '
                           << right_vrp->elts[1] << ' '
                           << right_vrp->elts[2] << '\n'
        << "   Left PRP: " << left_prp->elts[0] << ' '
                           << left_prp->elts[1] << ' '
                           << left_prp->elts[2] << '\n'
        << "  Right PRP: " << right_prp->elts[0] << ' '
                           << right_prp->elts[1] << ' '
                           << right_prp->elts[2] << '\n';
}