parapiped.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 <cassert>
#include <cmath>
#include <iostream>
#include <vector>

#include <inttypes.h>

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

#include <jwsc/base/bits.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++/base/exception.h>
#include <jwsc++/graphics/polygon.h>
#include <jwsc++/graphics/parapiped.h>

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


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


static void create_adjacency(jwsc::Matrix_u32** adjacency_out)
{
    using namespace jwsc;

    create_matrix_u32(adjacency_out, 6, 4);
    Matrix_u32* adjacency = *adjacency_out;

    adjacency->elts[5][0] = 1;
    adjacency->elts[5][1] = 2;
    adjacency->elts[5][2] = 0;
    adjacency->elts[5][3] = 3;

    adjacency->elts[3][0] = 0;
    adjacency->elts[3][1] = 4;
    adjacency->elts[3][2] = 1;
    adjacency->elts[3][3] = 5;

    adjacency->elts[0][0] = 2;
    adjacency->elts[0][1] = 4;
    adjacency->elts[0][2] = 3;
    adjacency->elts[0][3] = 5;

    adjacency->elts[2][0] = 0;
    adjacency->elts[2][1] = 5;
    adjacency->elts[2][2] = 1;
    adjacency->elts[2][3] = 4;

    adjacency->elts[4][0] = 2;
    adjacency->elts[4][1] = 1;
    adjacency->elts[4][2] = 3;
    adjacency->elts[4][3] = 0;

    adjacency->elts[1][0] = 3;
    adjacency->elts[1][1] = 4;
    adjacency->elts[1][2] = 2;
    adjacency->elts[1][3] = 5;
}


Parapiped_f::Parapiped_f
(
    float x1, float y1, float z1,
    float x2, float y2, float z2,
    float x3, float y3, float z3,
    float x4, float y4, float z4
)
: faces(6, (Polygon_f*)0), points(4, (jwsc::Vector_f*)0)
{
    using namespace jwsc;

    center = 0;
    adjacency = 0;

    faces[5] = new Polygon_f(4);
    faces[5]->add_point(x1, y1, z1);
    faces[5]->add_point(x1 + (x3 - x2), y1 + (y3 - y2), z1 + (z3 - z2));
    faces[5]->add_point(x3, y3, z3);
    faces[5]->add_point(x2, y2, z2);

    faces[3] = new Polygon_f(4);
    faces[3]->add_point(x2, y2, z2);
    faces[3]->add_point(x2 + (x4 - x3), y2 + (y4 - y3), z2 + (z4 - z3));
    faces[3]->add_point(x1 + (x4 - x3), y1 + (y4 - y3), z1 + (z4 - z3));
    faces[3]->add_point(x1, y1, z1);

    faces[0] = new Polygon_f(4);
    faces[0]->add_point(x3, y3, z3);
    faces[0]->add_point(x4, y4, z4);
    faces[0]->add_point(x2 + (x4 - x3), y2 + (y4 - y3), z2 + (z4 - z3));
    faces[0]->add_point(x2, y2, z2);

    faces[2] = new Polygon_f(4);
    faces[2]->add_point(x4, y4, z4);
    faces[2]->add_point(x3, y3, z3);
    faces[2]->add_point(x3 + (x1 - x2), y3 + (y1 - y2), z3 + (z1 - z2));
    faces[2]->add_point(x4 + (x1 - x2), y4 + (y1 - y2), z4 + (z1 - z2));

    faces[4] = new Polygon_f(4);
    faces[4]->add_point(x4, y4, z4);
    faces[4]->add_point(x4 + (x1 - x2), y4 + (y1 - y2), z4 + (z1 - z2));
    faces[4]->add_point(x1 + (x4 - x3), y1 + (y4 - y3), z1 + (z4 - z3));
    faces[4]->add_point(x2 + (x4 - x3), y2 + (y4 - y3), z2 + (z4 - z3));

    faces[1] = new Polygon_f(4);
    faces[1]->add_point(x1, y1, z1);
    faces[1]->add_point(x1 + (x4 - x3), y1 + (y4 - y3), z1 + (z4 - z3));
    faces[1]->add_point(x4 + (x1 - x2), y4 + (y1 - y2), z4 + (z1 - z2));
    faces[1]->add_point(x3 + (x1 - x2), y3 + (y1 - y2), z3 + (z1 - z2));

    create_init_vector_f(&(points[0]), 4, 1);
    points[0]->elts[0] = x1;
    points[0]->elts[1] = y1;
    points[0]->elts[2] = z1;

    create_init_vector_f(&(points[1]), 4, 1);
    points[1]->elts[0] = x2;
    points[1]->elts[1] = y2;
    points[1]->elts[2] = z2;

    create_init_vector_f(&(points[2]), 4, 1);
    points[2]->elts[0] = x3;
    points[2]->elts[1] = y3;
    points[2]->elts[2] = z3;

    create_init_vector_f(&(points[3]), 4, 1);
    points[3]->elts[0] = x4;
    points[3]->elts[1] = y4;
    points[3]->elts[2] = z4;

    create_zero_vector_f(&center, 4);
    create_adjacency(&adjacency);

    for (size_t i = 0; i < 6; i++)
    {
        const Vector_f* centroid = faces[ i ]->get_centroid();

        // Ignore the homogeneous coordinate, since it must be 1.
        center->elts[0] += centroid->elts[0];
        center->elts[1] += centroid->elts[1];
        center->elts[2] += centroid->elts[2];
    }

    center->elts[0] /= 6.0f;
    center->elts[1] /= 6.0f;
    center->elts[2] /= 6.0f;
    center->elts[3]  = 1.0f;
}


Parapiped_f::Parapiped_f
(
    const jwsc::Vector_f* p1, 
    const jwsc::Vector_f* p2, 
    const jwsc::Vector_f* p3, 
    const jwsc::Vector_f* p4
)
throw (jwscxx::base::Arg_error)
: faces(6, (Polygon_f*)0), points(4, (jwsc::Vector_f*)0)
{
    using namespace jwsc;

    center = 0;
    adjacency = 0;

    if (p1->num_elts != 3 || p2->num_elts != 3 ||
        p3->num_elts != 3 || p4->num_elts != 3)
    {
        throw Arg_error("Vectors for parapiped not 3D");
    }

    float x1 = p1->elts[0]; float y1 = p1->elts[1]; float z1 = p1->elts[2];
    float x2 = p2->elts[0]; float y2 = p2->elts[1]; float z2 = p2->elts[2];
    float x3 = p3->elts[0]; float y3 = p3->elts[1]; float z3 = p3->elts[2];
    float x4 = p4->elts[0]; float y4 = p4->elts[1]; float z4 = p4->elts[2];

    faces[5] = new Polygon_f(4);
    faces[5]->add_point(x1, y1, z1);
    faces[5]->add_point(x1 + (x3 - x2), y1 + (y3 - y2), z1 + (z3 - z2));
    faces[5]->add_point(x3, y3, z3);
    faces[5]->add_point(x2, y2, z2);

    faces[3] = new Polygon_f(4);
    faces[3]->add_point(x2, y2, z2);
    faces[3]->add_point(x2 + (x4 - x3), y2 + (y4 - y3), z2 + (z4 - z3));
    faces[3]->add_point(x1 + (x4 - x3), y1 + (y4 - y3), z1 + (z4 - z3));
    faces[3]->add_point(x1, y1, z1);

    faces[0] = new Polygon_f(4);
    faces[0]->add_point(x3, y3, z3);
    faces[0]->add_point(x4, y4, z4);
    faces[0]->add_point(x2 + (x4 - x3), y2 + (y4 - y3), z2 + (z4 - z3));
    faces[0]->add_point(x2, y2, z2);

    faces[2] = new Polygon_f(4);
    faces[2]->add_point(x4, y4, z4);
    faces[2]->add_point(x3, y3, z3);
    faces[2]->add_point(x3 + (x1 - x2), y3 + (y1 - y2), z3 + (z1 - z2));
    faces[2]->add_point(x4 + (x1 - x2), y4 + (y1 - y2), z4 + (z1 - z2));

    faces[4] = new Polygon_f(4);
    faces[4]->add_point(x4, y4, z4);
    faces[4]->add_point(x4 + (x1 - x2), y4 + (y1 - y2), z4 + (z1 - z2));
    faces[4]->add_point(x1 + (x4 - x3), y1 + (y4 - y3), z1 + (z4 - z3));
    faces[4]->add_point(x2 + (x4 - x3), y2 + (y4 - y3), z2 + (z4 - z3));

    faces[1] = new Polygon_f(4);
    faces[1]->add_point(x1, y1, z1);
    faces[1]->add_point(x1 + (x4 - x3), y1 + (y4 - y3), z1 + (z4 - z3));
    faces[1]->add_point(x4 + (x1 - x2), y4 + (y1 - y2), z4 + (z1 - z2));
    faces[1]->add_point(x3 + (x1 - x2), y3 + (y1 - y2), z3 + (z1 - z2));

    create_init_vector_f(&(points[0]), 4, 1);
    points[0]->elts[0] = x1;
    points[0]->elts[1] = y1;
    points[0]->elts[2] = z1;

    create_init_vector_f(&(points[1]), 4, 1);
    points[1]->elts[0] = x2;
    points[1]->elts[1] = y2;
    points[1]->elts[2] = z2;

    create_init_vector_f(&(points[2]), 4, 1);
    points[2]->elts[0] = x3;
    points[2]->elts[1] = y3;
    points[2]->elts[2] = z3;

    create_init_vector_f(&(points[3]), 4, 1);
    points[3]->elts[0] = x4;
    points[3]->elts[1] = y4;
    points[3]->elts[2] = z4;

    create_zero_vector_f(&center, 4);
    create_adjacency(&adjacency);

    for (size_t i = 0; i < 6; i++)
    {
        const Vector_f* centroid = faces[ i ]->get_centroid();

        // Ignore the homogeneous coordinate, since it must be 1.
        center->elts[0] += centroid->elts[0];
        center->elts[1] += centroid->elts[1];
        center->elts[2] += centroid->elts[2];
    }

    center->elts[0] /= 6.0f;
    center->elts[1] /= 6.0f;
    center->elts[2] /= 6.0f;
    center->elts[3]  = 1.0f;
}


Parapiped_f::Parapiped_f(const Parapiped_f& p)
    : faces(6, (Polygon_f*)0), points(4, (jwsc::Vector_f*)0)
{
    center = 0;

    for (size_t i = 0; i < 6; i++)
    {
        faces[ i ] = p.faces[ i ]->clone();
    }

    for (size_t i = 0; i < 4; i++)
    {
        copy_vector_f(&(points[ i ]), p.points[ i ]);
    }

    copy_vector_f(&center, p.center);
    copy_matrix_u32(&adjacency, p.adjacency);
}


Parapiped_f::Parapiped_f(const char* fname) 
    throw (jwscxx::base::Arg_error, jwscxx::base::IO_error)
    : faces(6, (Polygon_f*)0), points(4, (jwsc::Vector_f*)0)
{
    center = 0;
    adjacency = 0;

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


Parapiped_f::Parapiped_f(std::istream& in)
    throw (jwscxx::base::Arg_error, jwscxx::base::IO_error)
    : faces(6, (Polygon_f*)0), points(4, (jwsc::Vector_f*)0)
{
    center = 0;
    adjacency = 0;

    read(in);
}


Parapiped_f::~Parapiped_f()
{
    for (size_t i = 0; i < 6; i++)
    {
        delete faces[ i ];
    }

    for (size_t i = 0; i < 4; i++)
    {
        free_vector_f(points[ i ]);
    }

    free_vector_f(center);
    free_matrix_u32(adjacency);
}


Parapiped_f& Parapiped_f::operator= (const Parapiped_f& p)
{
    // Don't try to assign to self.
    if (points[0] == p.points[0])
        return *this;

    for (size_t i = 0; i < 6; i++)
    {
        *(faces[ i ]) = *(p.faces[ i ]);
    }

    for (size_t i = 0; i < 4; i++)
    {
        copy_vector_f(&(points[ i ]), p.points[ i ]);
    }

    copy_vector_f(&center, p.center);
    copy_matrix_u32(&adjacency, p.adjacency);

    return *this;
}


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


size_t Parapiped_f::num_faces() const
{
    return 6;
}


const jwsc::Vector_f* Parapiped_f::get_point(size_t i) const
    throw (jwscxx::base::Arg_error)
{
    if (i >= 4)
    { 
        throw jwscxx::base::Arg_error("Invalid point index"); 
    }

    return static_cast<const jwsc::Vector_f*>(points[ i ]);
}


const jwsc::Vector_f* Parapiped_f::get_center() const 
{ 
    return center; 
}


const jwscxx::graphics::Polygon_f* Parapiped_f::get_face(size_t i) const 
    throw (jwscxx::base::Arg_error)
{ 
    if (i >= 6)
    { 
        throw jwscxx::base::Arg_error("Invalid face index"); 
    }

    return static_cast<const jwscxx::graphics::Polygon_f*>(faces[ i ]);
}


const std::vector<const jwscxx::graphics::Polygon_f*>& Parapiped_f::get_faces() 
    const
{
    return (const std::vector<const jwscxx::graphics::Polygon_f*>&)(faces);
}


void Parapiped_f::set_face_id(size_t i, uint32_t id)
    throw (jwscxx::base::Arg_error)
{
    if (i >= 6)
    {
        throw jwscxx::base::Arg_error("Invalid face index");
    }

    faces[ i ]->set_id(id);
}


void Parapiped_f::read(std::istream& in) throw (jwscxx::base::IO_error,
        jwscxx::base::Arg_error)
{
    using namespace jwsc;

    for (size_t i = 0; i < 6; i++)
    {
        if (faces[ i ] == 0)
        {
            faces[ i ] = new Polygon_f(4);
        }
        faces[ i ]->read(in);
    }

    assert(sizeof(float) == sizeof(uint32_t));

    float x, y, z, w;

    for (size_t i = 0; i < 4; i++)
    {
        in.read((char*)&x, sizeof(float));
        in.read((char*)&y, sizeof(float));
        in.read((char*)&z, sizeof(float));
        in.read((char*)&w, sizeof(float));
        if (in.fail() || in.eof())
        {
            throw IO_error("Could not read point");
        }

#ifndef JWSC_BIGENDIAN
        bswap_u32((uint32_t*)&(x));
        bswap_u32((uint32_t*)&(y));
        bswap_u32((uint32_t*)&(z));
        bswap_u32((uint32_t*)&(w));
#endif

        create_vector_f(&(points[ i ]), 4);
        points[ i ]->elts[0] = x;
        points[ i ]->elts[1] = y;
        points[ i ]->elts[2] = z;
        points[ i ]->elts[3] = w;
    }

    in.read((char*)&x, sizeof(float));
    in.read((char*)&y, sizeof(float));
    in.read((char*)&z, sizeof(float));
    in.read((char*)&w, sizeof(float));
    if (in.fail() || in.eof())
    {
        throw IO_error("Could not read center");
    }

#ifndef JWSC_BIGENDIAN
    bswap_u32((uint32_t*)&(x));
    bswap_u32((uint32_t*)&(y));
    bswap_u32((uint32_t*)&(z));
    bswap_u32((uint32_t*)&(w));
#endif

    create_vector_f(&center, 4);
    center->elts[0] = x;
    center->elts[1] = y;
    center->elts[2] = z;
    center->elts[3] = w;

    create_adjacency(&adjacency);
}


void Parapiped_f::write(std::ostream& out) const throw (jwscxx::base::IO_error)
{
    using namespace jwsc;

    for (size_t i = 0; i < 6; i++)
    {
        faces[ i ]->write(out);
    }

    assert(sizeof(float) == sizeof(uint32_t));

    for (size_t i = 0; i < 4; i++)
    {
        float x = points[ i ]->elts[0];
        float y = points[ i ]->elts[1];
        float z = points[ i ]->elts[2];
        float w = points[ i ]->elts[3];

#ifndef JWSC_BIGENDIAN
        bswap_u32((uint32_t*)&(x));
        bswap_u32((uint32_t*)&(y));
        bswap_u32((uint32_t*)&(z));
        bswap_u32((uint32_t*)&(w));
#endif

        out.write((char*)&x, sizeof(float));
        out.write((char*)&y, sizeof(float));
        out.write((char*)&z, sizeof(float));
        out.write((char*)&w, sizeof(float));
        if (out.fail() || out.eof())
        {
            throw IO_error("Could not write point");
        }
    }

    float x = center->elts[0];
    float y = center->elts[1];
    float z = center->elts[2];
    float w = center->elts[3];

#ifndef JWSC_BIGENDIAN
    bswap_u32((uint32_t*)&(x));
    bswap_u32((uint32_t*)&(y));
    bswap_u32((uint32_t*)&(z));
    bswap_u32((uint32_t*)&(w));
#endif

    out.write((char*)&x, sizeof(float));
    out.write((char*)&y, sizeof(float));
    out.write((char*)&z, sizeof(float));
    out.write((char*)&w, sizeof(float));
    if (out.fail() || out.eof())
    {
        throw IO_error("Could not write center");
    }
}


void Parapiped_f::transform(const jwsc::Matrix_f* M) 
    throw (jwscxx::base::Arg_error)
{
    for (size_t i = 0; i < 6; i++)
    {
        faces[ i ]->transform(M);
    }

    for (size_t i = 0; i < 4; i++)
    {
        multiply_matrix_and_vector_f(&(points[ i ]), M, points[ i ]);
    }

    // Transform the center.
    multiply_matrix_and_vector_f(&center, M, center);
}


void Parapiped_f::wire_render() const
{
    for (size_t i = 0; i < 6; i++)
    {
        faces[ i ]->wire_render();
    }
}


void Parapiped_f::wire_occlude_render() const
{
    for (size_t i = 0; i < 6; i++)
    {
        faces[ i ]->wire_occlude_render();
    }
}


void Parapiped_f::solid_render() const
{
    for (size_t i = 0; i < 6; i++)
    {
        faces[ i ]->solid_render();
    }
}


void Parapiped_f::project()
{
    for (size_t i = 0; i < 6; i++)
    {
        faces[ i ]->project();
    }
}


uint32_t Parapiped_f::adjacent_face(uint32_t i, uint32_t p) const 
    throw (jwscxx::base::Arg_error)
{
    if (i >= 6)
    { 
        throw jwscxx::base::Arg_error("Invalid face index"); 
    }
    if (p > 3)
    {
        throw jwscxx::base::Arg_error("Invalid face point index");
    }

    return adjacency->elts[ i ][ p ];
}


Parapiped_d::Parapiped_d
(
    double x1, double y1, double z1,
    double x2, double y2, double z2,
    double x3, double y3, double z3,
    double x4, double y4, double z4
)
: faces(6, (Polygon_d*)0), points(4, (jwsc::Vector_d*)0)
{
    using namespace jwsc;

    center = 0;
    adjacency = 0;

    faces[5] = new Polygon_d(4);
    faces[5]->add_point(x1, y1, z1);
    faces[5]->add_point(x1 + (x3 - x2), y1 + (y3 - y2), z1 + (z3 - z2));
    faces[5]->add_point(x3, y3, z3);
    faces[5]->add_point(x2, y2, z2);

    faces[3] = new Polygon_d(4);
    faces[3]->add_point(x2, y2, z2);
    faces[3]->add_point(x2 + (x4 - x3), y2 + (y4 - y3), z2 + (z4 - z3));
    faces[3]->add_point(x1 + (x4 - x3), y1 + (y4 - y3), z1 + (z4 - z3));
    faces[3]->add_point(x1, y1, z1);

    faces[0] = new Polygon_d(4);
    faces[0]->add_point(x3, y3, z3);
    faces[0]->add_point(x4, y4, z4);
    faces[0]->add_point(x2 + (x4 - x3), y2 + (y4 - y3), z2 + (z4 - z3));
    faces[0]->add_point(x2, y2, z2);

    faces[2] = new Polygon_d(4);
    faces[2]->add_point(x4, y4, z4);
    faces[2]->add_point(x3, y3, z3);
    faces[2]->add_point(x3 + (x1 - x2), y3 + (y1 - y2), z3 + (z1 - z2));
    faces[2]->add_point(x4 + (x1 - x2), y4 + (y1 - y2), z4 + (z1 - z2));

    faces[4] = new Polygon_d(4);
    faces[4]->add_point(x4, y4, z4);
    faces[4]->add_point(x4 + (x1 - x2), y4 + (y1 - y2), z4 + (z1 - z2));
    faces[4]->add_point(x1 + (x4 - x3), y1 + (y4 - y3), z1 + (z4 - z3));
    faces[4]->add_point(x2 + (x4 - x3), y2 + (y4 - y3), z2 + (z4 - z3));

    faces[1] = new Polygon_d(4);
    faces[1]->add_point(x1, y1, z1);
    faces[1]->add_point(x1 + (x4 - x3), y1 + (y4 - y3), z1 + (z4 - z3));
    faces[1]->add_point(x4 + (x1 - x2), y4 + (y1 - y2), z4 + (z1 - z2));
    faces[1]->add_point(x3 + (x1 - x2), y3 + (y1 - y2), z3 + (z1 - z2));

    create_init_vector_d(&(points[0]), 4, 1);
    points[0]->elts[0] = x1;
    points[0]->elts[1] = y1;
    points[0]->elts[2] = z1;

    create_init_vector_d(&(points[1]), 4, 1);
    points[1]->elts[0] = x2;
    points[1]->elts[1] = y2;
    points[1]->elts[2] = z2;

    create_init_vector_d(&(points[2]), 4, 1);
    points[2]->elts[0] = x3;
    points[2]->elts[1] = y3;
    points[2]->elts[2] = z3;

    create_init_vector_d(&(points[3]), 4, 1);
    points[3]->elts[0] = x4;
    points[3]->elts[1] = y4;
    points[3]->elts[2] = z4;

    create_zero_vector_d(&center, 4);
    create_adjacency(&adjacency);

    for (size_t i = 0; i < 6; i++)
    {
        const Vector_d* centroid = faces[ i ]->get_centroid();

        // Ignore the homogeneous coordinate, since it must be 1.
        center->elts[0] += centroid->elts[0];
        center->elts[1] += centroid->elts[1];
        center->elts[2] += centroid->elts[2];
    }

    center->elts[0] /= 6.0;
    center->elts[1] /= 6.0;
    center->elts[2] /= 6.0;
    center->elts[3]  = 1.0;
}


Parapiped_d::Parapiped_d
(
    const jwsc::Vector_d* p1, 
    const jwsc::Vector_d* p2, 
    const jwsc::Vector_d* p3, 
    const jwsc::Vector_d* p4
)
throw (jwscxx::base::Arg_error)
: faces(6, (Polygon_d*)0), points(4, (jwsc::Vector_d*)0)
{
    using namespace jwsc;

    center = 0;
    adjacency = 0;

    if (p1->num_elts != 3 || p2->num_elts != 3 ||
        p3->num_elts != 3 || p4->num_elts != 3)
    {
        throw Arg_error("Vectors for parapiped not 3D");
    }

    double x1 = p1->elts[0]; double y1 = p1->elts[1]; double z1 = p1->elts[2];
    double x2 = p2->elts[0]; double y2 = p2->elts[1]; double z2 = p2->elts[2];
    double x3 = p3->elts[0]; double y3 = p3->elts[1]; double z3 = p3->elts[2];
    double x4 = p4->elts[0]; double y4 = p4->elts[1]; double z4 = p4->elts[2];

    faces[5] = new Polygon_d(4);
    faces[5]->add_point(x1, y1, z1);
    faces[5]->add_point(x1 + (x3 - x2), y1 + (y3 - y2), z1 + (z3 - z2));
    faces[5]->add_point(x3, y3, z3);
    faces[5]->add_point(x2, y2, z2);

    faces[3] = new Polygon_d(4);
    faces[3]->add_point(x2, y2, z2);
    faces[3]->add_point(x2 + (x4 - x3), y2 + (y4 - y3), z2 + (z4 - z3));
    faces[3]->add_point(x1 + (x4 - x3), y1 + (y4 - y3), z1 + (z4 - z3));
    faces[3]->add_point(x1, y1, z1);

    faces[0] = new Polygon_d(4);
    faces[0]->add_point(x3, y3, z3);
    faces[0]->add_point(x4, y4, z4);
    faces[0]->add_point(x2 + (x4 - x3), y2 + (y4 - y3), z2 + (z4 - z3));
    faces[0]->add_point(x2, y2, z2);

    faces[2] = new Polygon_d(4);
    faces[2]->add_point(x4, y4, z4);
    faces[2]->add_point(x3, y3, z3);
    faces[2]->add_point(x3 + (x1 - x2), y3 + (y1 - y2), z3 + (z1 - z2));
    faces[2]->add_point(x4 + (x1 - x2), y4 + (y1 - y2), z4 + (z1 - z2));

    faces[4] = new Polygon_d(4);
    faces[4]->add_point(x4, y4, z4);
    faces[4]->add_point(x4 + (x1 - x2), y4 + (y1 - y2), z4 + (z1 - z2));
    faces[4]->add_point(x1 + (x4 - x3), y1 + (y4 - y3), z1 + (z4 - z3));
    faces[4]->add_point(x2 + (x4 - x3), y2 + (y4 - y3), z2 + (z4 - z3));

    faces[1] = new Polygon_d(4);
    faces[1]->add_point(x1, y1, z1);
    faces[1]->add_point(x1 + (x4 - x3), y1 + (y4 - y3), z1 + (z4 - z3));
    faces[1]->add_point(x4 + (x1 - x2), y4 + (y1 - y2), z4 + (z1 - z2));
    faces[1]->add_point(x3 + (x1 - x2), y3 + (y1 - y2), z3 + (z1 - z2));

    create_init_vector_d(&(points[0]), 4, 1);
    points[0]->elts[0] = x1;
    points[0]->elts[1] = y1;
    points[0]->elts[2] = z1;

    create_init_vector_d(&(points[1]), 4, 1);
    points[1]->elts[0] = x2;
    points[1]->elts[1] = y2;
    points[1]->elts[2] = z2;

    create_init_vector_d(&(points[2]), 4, 1);
    points[2]->elts[0] = x3;
    points[2]->elts[1] = y3;
    points[2]->elts[2] = z3;

    create_init_vector_d(&(points[3]), 4, 1);
    points[3]->elts[0] = x4;
    points[3]->elts[1] = y4;
    points[3]->elts[2] = z4;

    create_zero_vector_d(&center, 4);
    create_adjacency(&adjacency);

    for (size_t i = 0; i < 6; i++)
    {
        const Vector_d* centroid = faces[ i ]->get_centroid();

        // Ignore the homogeneous coordinate, since it must be 1.
        center->elts[0] += centroid->elts[0];
        center->elts[1] += centroid->elts[1];
        center->elts[2] += centroid->elts[2];
    }

    center->elts[0] /= 6.0;
    center->elts[1] /= 6.0;
    center->elts[2] /= 6.0;
    center->elts[3]  = 1.0;
}


Parapiped_d::Parapiped_d(const Parapiped_d& p)
    : faces(6, (Polygon_d*)0), points(4, (jwsc::Vector_d*)0)
{
    center = 0;
    adjacency = 0;

    for (size_t i = 0; i < 6; i++)
    {
        faces[ i ] = p.faces[ i ]->clone();
    }

    for (size_t i = 0; i < 4; i++)
    {
        copy_vector_d(&(points[ i ]), p.points[ i ]);
    }

    copy_vector_d(&center, p.center);
    copy_matrix_u32(&adjacency, p.adjacency);
}


Parapiped_d::Parapiped_d(const char* fname) 
    throw (jwscxx::base::Arg_error, jwscxx::base::IO_error)
    : faces(6, (Polygon_d*)0), points(4, (jwsc::Vector_d*)0)
{
    center = 0;
    adjacency = 0;

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


Parapiped_d::Parapiped_d(std::istream& in)
    throw (jwscxx::base::Arg_error, jwscxx::base::IO_error)
    : faces(6, (Polygon_d*)0), points(4, (jwsc::Vector_d*)0)
{
    center = 0;
    adjacency = 0;

    read(in);
}


Parapiped_d::~Parapiped_d()
{
    for (size_t i = 0; i < 6; i++)
    {
        delete faces[ i ];
    }

    for (size_t i = 0; i < 4; i++)
    {
        free_vector_d(points[ i ]);
    }

    free_vector_d(center);
    free_matrix_u32(adjacency);
}


Parapiped_d& Parapiped_d::operator= (const Parapiped_d& p)
{
    // Don't try to assign to self.
    if (points[0] == p.points[0])
        return *this;

    for (size_t i = 0; i < 6; i++)
    {
        *(faces[ i ]) = *(p.faces[ i ]);
    }

    for (size_t i = 0; i < 4; i++)
    {
        copy_vector_d(&(points[ i ]), p.points[ i ]);
    }

    copy_vector_d(&center, p.center);
    copy_matrix_u32(&adjacency, p.adjacency);

    return *this;
}


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


size_t Parapiped_d::num_faces() const
{
    return 6;
}


const jwsc::Vector_d* Parapiped_d::get_point(size_t i) const
    throw (jwscxx::base::Arg_error)
{
    if (i >= 4)
    { 
        throw jwscxx::base::Arg_error("Invalid point index"); 
    }

    return static_cast<const jwsc::Vector_d*>(points[ i ]);
}


const jwsc::Vector_d* Parapiped_d::get_center() const 
{ 
    return center; 
}


const jwscxx::graphics::Polygon_d* Parapiped_d::get_face(size_t i) const 
    throw (jwscxx::base::Arg_error)
{ 
    if (i >= 6)
    { 
        throw jwscxx::base::Arg_error("Invalid face index"); 
    }

    return static_cast<const jwscxx::graphics::Polygon_d*>(faces[ i ]);
}


const std::vector<const jwscxx::graphics::Polygon_d*>& Parapiped_d::get_faces() 
    const
{
    return (const std::vector<const jwscxx::graphics::Polygon_d*>&)(faces);
}


void Parapiped_d::set_face_id(size_t i, uint32_t id)
    throw (jwscxx::base::Arg_error)
{
    if (i >= 6)
    {
        throw jwscxx::base::Arg_error("Invalid face index");
    }

    faces[ i ]->set_id(id);
}


void Parapiped_d::read(std::istream& in) throw (jwscxx::base::IO_error,
        jwscxx::base::Arg_error)
{
    using namespace jwsc;

    for (size_t i = 0; i < 6; i++)
    {
        if (faces[ i ] == 0)
        {
            faces[ i ] = new Polygon_d(4);
        }
        faces[ i ]->read(in);
    }

    assert(sizeof(double) == sizeof(uint64_t));

    double x, y, z, w;

    for (size_t i = 0; i < 4; i++)
    {
        in.read((char*)&x, sizeof(double));
        in.read((char*)&y, sizeof(double));
        in.read((char*)&z, sizeof(double));
        in.read((char*)&w, sizeof(double));
        if (in.fail() || in.eof())
        {
            throw IO_error("Could not read point");
        }

#ifndef JWSC_BIGENDIAN
        bswap_u64((uint64_t*)&(x));
        bswap_u64((uint64_t*)&(y));
        bswap_u64((uint64_t*)&(z));
        bswap_u64((uint64_t*)&(w));
#endif

        create_vector_d(&(points[ i ]), 4);
        points[ i ]->elts[0] = x;
        points[ i ]->elts[1] = y;
        points[ i ]->elts[2] = z;
        points[ i ]->elts[3] = w;
    }

    in.read((char*)&x, sizeof(double));
    in.read((char*)&y, sizeof(double));
    in.read((char*)&z, sizeof(double));
    in.read((char*)&w, sizeof(double));
    if (in.fail() || in.eof())
    {
        throw IO_error("Could not read center");
    }

#ifndef JWSC_BIGENDIAN
    bswap_u64((uint64_t*)&(x));
    bswap_u64((uint64_t*)&(y));
    bswap_u64((uint64_t*)&(z));
    bswap_u64((uint64_t*)&(w));
#endif

    create_vector_d(&center, 4);
    center->elts[0] = x;
    center->elts[1] = y;
    center->elts[2] = z;
    center->elts[3] = w;

    create_adjacency(&adjacency);
}


void Parapiped_d::write(std::ostream& out) const throw (jwscxx::base::IO_error)
{
    using namespace jwsc;

    for (size_t i = 0; i < 6; i++)
    {
        faces[ i ]->write(out);
    }

    assert(sizeof(double) == sizeof(uint64_t));

    for (size_t i = 0; i < 4; i++)
    {
        double x = points[ i ]->elts[0];
        double y = points[ i ]->elts[1];
        double z = points[ i ]->elts[2];
        double w = points[ i ]->elts[3];

#ifndef JWSC_BIGENDIAN
        bswap_u64((uint64_t*)&(x));
        bswap_u64((uint64_t*)&(y));
        bswap_u64((uint64_t*)&(z));
        bswap_u64((uint64_t*)&(w));
#endif

        out.write((char*)&x, sizeof(double));
        out.write((char*)&y, sizeof(double));
        out.write((char*)&z, sizeof(double));
        out.write((char*)&w, sizeof(double));
        if (out.fail() || out.eof())
        {
            throw IO_error("Could not write point");
        }
    }

    double x = center->elts[0];
    double y = center->elts[1];
    double z = center->elts[2];
    double w = center->elts[3];

#ifndef JWSC_BIGENDIAN
    bswap_u64((uint64_t*)&(x));
    bswap_u64((uint64_t*)&(y));
    bswap_u64((uint64_t*)&(z));
    bswap_u64((uint64_t*)&(w));
#endif

    out.write((char*)&x, sizeof(double));
    out.write((char*)&y, sizeof(double));
    out.write((char*)&z, sizeof(double));
    out.write((char*)&w, sizeof(double));
    if (out.fail() || out.eof())
    {
        throw IO_error("Could not write center");
    }
}


void Parapiped_d::transform(const jwsc::Matrix_d* M) 
    throw (jwscxx::base::Arg_error)
{
    for (size_t i = 0; i < 6; i++)
    {
        faces[ i ]->transform(M);
    }

    for (size_t i = 0; i < 4; i++)
    {
        multiply_matrix_and_vector_d(&(points[ i ]), M, points[ i ]);
    }

    // Transform the center.
    multiply_matrix_and_vector_d(&center, M, center);
}


void Parapiped_d::wire_render() const
{
    for (size_t i = 0; i < 6; i++)
    {
        faces[ i ]->wire_render();
    }
}


void Parapiped_d::wire_occlude_render() const
{
    for (size_t i = 0; i < 6; i++)
    {
        faces[ i ]->wire_occlude_render();
    }
}


void Parapiped_d::solid_render() const
{
    for (size_t i = 0; i < 6; i++)
    {
        faces[ i ]->solid_render();
    }
}


void Parapiped_d::project()
{
    for (size_t i = 0; i < 6; i++)
    {
        faces[ i ]->project();
    }
}


uint32_t Parapiped_d::adjacent_face(uint32_t i, uint32_t p) const 
    throw (jwscxx::base::Arg_error)
{
    if (i >= 6)
    { 
        throw jwscxx::base::Arg_error("Invalid face index"); 
    }
    if (p > 3)
    {
        throw jwscxx::base::Arg_error("Invalid face point index");
    }

    return adjacency->elts[ i ][ p ];
}