spore.cpp

Go to the documentation of this file.

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


#include <config.h>

#include <cmath>
#include <cstring>
#include <cassert>
#include <iostream>
#include <sstream>
#include <list>
#include <vector>

#include <inttypes.h>

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

#include <jwsc/base/limits.h>
#include <jwsc/math/constants.h>
#include <jwsc/prob/pdf.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/matblock/matblock.h>

#include "density.h"
#include "printable.h"
#include "structure.h"
#include "dd_spore.h"
#include "hypha.h"
#include "spore.h"


#define  CENTROID_X_MIN  0
#define  CENTROID_X_MAX  0

#define  CENTROID_Y_MIN  0
#define  CENTROID_Y_MAX  0

#define  CENTROID_Z_MIN  0
#define  CENTROID_Z_MAX  0

#define  LENGTH_MU       40.0f
#define  LENGTH_SIGMA    4.0f
#define  LENGTH_MIN      10.f
#define  LENGTH_MAX      100.0f

#define  WIDTH_MU        15.0f
#define  WIDTH_SIGMA     1.0f
#define  WIDTH_MIN       2.0f
#define  WIDTH_MAX       30.0f

#define  THETA_MU        0
#define  THETA_SIGMA     0.2f
#define  THETA_MIN       0
#define  THETA_MAX       JWSC_PI_2

#define  PSI_MIN         -JWSC_PI
#define  PSI_MAX         JWSC_PI

#define  OPACITY_MU      0.9f
#define  OPACITY_SIGMA   0.2f
#define  OPACITY_MIN     0.1f
#define  OPACITY_MAX     1.0f


using std::istream;
using std::ostringstream;
using std::strncmp;
using std::strlen;
using std::list;
using std::vector;
using namespace jwsc;
using jwscxx::base::Arg_error;
using jwscxx::base::IO_error;


/* ----------------------------  Spore_density  ----------------------------- */

Spore_density::Spore_density(const DD_spore_set* dd_spores) throw (Arg_error)
{
    float mu, sigma;
    float min, max;

    min = CENTROID_X_MIN;
    max = CENTROID_X_MAX;
    set_centroid_x(min, max);

    min = CENTROID_Y_MIN;
    max = CENTROID_Y_MAX;
    set_centroid_y(min, max);

    min = CENTROID_Z_MIN;
    max = CENTROID_Z_MAX;
    set_centroid_z(min, max);

    mu    = LENGTH_MU;
    sigma = LENGTH_SIGMA;
    min   = LENGTH_MIN;
    max   = LENGTH_MAX;
    set_length(mu, sigma, min, max);

    mu    = WIDTH_MU;
    sigma = WIDTH_SIGMA;
    min   = WIDTH_MIN;
    max   = WIDTH_MAX;
    set_width(mu, sigma, min, max);

    mu    = THETA_MU;
    sigma = THETA_SIGMA;
    min   = THETA_MIN;
    max   = THETA_MAX;
    set_theta(mu, sigma, min, max);

    min = PSI_MIN;
    max = PSI_MAX;
    set_psi(min, max);

    mu    = OPACITY_MU;
    sigma = OPACITY_SIGMA;
    min   = OPACITY_MIN;
    max   = OPACITY_MAX;
    set_opacity(mu, sigma, min, max);

    this->dd_spores = dd_spores;
}


Spore_density* Spore_density::clone() const
{
    float mu, sigma;
    float min, max;

    Spore_density* s  = new Spore_density();

    min = centroid_x.get_min();
    max = centroid_x.get_max();
    s->set_centroid_x(min, max);

    min = centroid_y.get_min();
    max = centroid_y.get_max();
    s->set_centroid_y(min, max);

    min = centroid_z.get_min();
    max = centroid_z.get_max();
    s->set_centroid_z(min, max);

    mu    = length.get_mu();
    sigma = length.get_sigma();
    min   = length.get_min();
    max   = length.get_max();
    s->set_length(mu, sigma, min, max);

    mu    = width.get_mu();
    sigma = width.get_sigma();
    min   = width.get_min();
    max   = width.get_max();
    s->set_width(mu, sigma, min, max);

    mu    = theta.get_mu();
    sigma = theta.get_sigma();
    min   = theta.get_min();
    max   = theta.get_max();
    s->set_theta(mu, sigma, min, max);

    min = psi.get_min();
    max = psi.get_max();
    s->set_psi(min, max);

    mu    = opacity.get_mu();
    sigma = opacity.get_sigma();
    min   = opacity.get_min();
    max   = opacity.get_max();
    s->set_opacity(mu, sigma, min, max);

    return s;
}


Spore* Spore_density::sample(float base_theta, float base_psi, size_t level) 
    const
{
    return new Spore(sample_centroid_x(), sample_centroid_y(), 
            sample_centroid_z(), sample_length(), sample_width(), base_theta, 
            base_psi, sample_theta(), sample_psi(), sample_opacity(), level,
            this);
}


Spore* Spore_density::sample
(
    float centroid_x, 
    float centroid_y, 
    float centroid_z,
    float base_theta,
    float base_psi,
    size_t level
) 
const throw (Arg_error)
{
    return new Spore(centroid_x, centroid_y, centroid_z, sample_length(), 
            sample_width(), base_theta, base_psi, sample_theta(), sample_psi(), 
            sample_opacity(), level, this);
}


Spore* Spore_density::sample
(
    Structure* parent,
    float      base_theta,
    float      base_psi
) 
const throw (Arg_error)
{
    if (dd_spores)
    {
        return dd_spores->propose_similar_spore(parent, base_theta, base_psi, 
                this, 0);
    }

    return new Spore(parent, sample_length(), sample_width(), base_theta,
            base_psi, sample_theta(), sample_psi(), sample_opacity(), this);
}

/* -------------------------------------------------------------------------- */






/* --------------------------------  Spore  --------------------------------- */

Spore::Spore
(
    float centroid_x, 
    float centroid_y,
    float centroid_z,
    float length,
    float width,
    float base_theta,
    float base_psi,
    float theta,
    float psi,
    float opacity,
    size_t level,
    const class Spore_density* density
) 
throw (Arg_error)
    : Apical_structure(centroid_x, centroid_y, centroid_z, length, width, 
            base_theta, base_psi, theta, psi, opacity, level, density)
{
    this->density = density;

    update_ancestry_log_prob();
}


Spore::Spore
(
    Structure* parent,
    float      length,
    float      width,
    float      base_theta,
    float      base_psi,
    float      theta,
    float      psi,
    float      opacity,
    const class Spore_density* density
)
throw (Arg_error, Apical_error)
    : Apical_structure(parent, length, width, base_theta, base_psi, theta, psi, 
            opacity, density)
{
    this->density = density;

    update_ancestry_log_prob();
}


Spore::Spore
(
    Structure* parent,
    float      centroid_x, 
    float      centroid_y,
    float      centroid_z,
    float      length,
    float      width,
    float      base_theta,
    float      base_psi,
    float      opacity,
    const class Spore_density* density
)
throw (Arg_error, Apical_error)
    : Apical_structure(parent, centroid_x, centroid_y, centroid_z, length, 
            width, base_theta, base_psi, opacity, density)
{
    this->density = density;

    update_ancestry_log_prob();
}


Spore::Spore(const Spore& s) : Apical_structure(s)
{
    density = s.density;
}


Spore::Spore
(
    const class Apical_hypha& h, 
    float                     width, 
    const Spore_density*      density
) 
throw (Arg_error) : Apical_structure(h)
{
    Structure::density        = density;
    Apical_structure::density = density;
    this->density             = density;

    this->width               = width;

    try
    {
        check_width();
        check_length();
        check_theta();
        check_psi();
        check_opacity();
    }
    catch (Arg_error e)
    {
        apical = 0;
        laterals.clear();
        throw e;
    }

    update_log_prob();
}


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


Spore& Spore::operator= (const Spore& s)
{
    Apical_structure::operator=(s);

    density = s.density;

    return *this;
}


Spore* Spore::split_into_apical
(
    const Spore* rvals_1,
    const Spore* rvals_2
) 
throw (Arg_error)
{
    return (Spore*)Apical_structure::split_into_apical(rvals_1, rvals_2);
}


Spore* Spore::merge_with_apical
(
    const Spore* rvals
) 
throw (Arg_error)
{
    if (!apical)
    {
        throw Arg_error("No apical to merge with");
    }

    if (!dynamic_cast<Spore*>(apical))
    {
        throw Arg_error("Apical to merge with not a spore");
    }

    return (Spore*)Apical_structure::merge_with_apical(rvals);
}


class Apical_hypha* Spore::split_into_apical
(
    const Apical_hypha* rvals_1,
    const Apical_hypha* rvals_2
) 
throw (Arg_error)
{
    Apical_hypha* split = new Apical_hypha(*this, rvals_1->get_width(),
            rvals_1->get_density());

    try
    {
        this->replace(split);
    }
    catch (Arg_error e)
    {
        delete split;
        throw e;
    }

    try
    {
        split->split_into_apical(rvals_1, rvals_2);
    }
    catch (Arg_error e)
    {
        split->replace(this);
        delete split;
        throw e;
    }

    return split;
}


void Spore::replace(class Apical_hypha* hypha) throw (Arg_error)
{
    Apical_structure::replace(hypha);
}


bool Spore::draw_in_matrix_f
(
     jwsc::Matrix_f* M, 
     float     x, 
     float     y,
     float     z,
     float     x_scale,
     float     y_scale,
     float     z_scale,
     bool      fill
) 
const
{
    uint32_t num_rows;
    uint32_t num_cols;
    uint32_t row, col;
    int64_t  row2, col2;

    float center;
    float x_pad;
    float y_pad;

    Vector_f* C  = 0;
    Matrix_f* M2 = 0;

    /* Create a sub-matrix with just the spore in it. */
    center = (width > length) ? width : length;

    x_pad = 3*x_scale;
    y_pad = 3*y_scale;

    create_vector_f(&C, 3);
    C->elts[0] = center + x_pad;
    C->elts[1] = center + y_pad;
    C->elts[2] = centroid->elts[2];

    // Copy this spore.
    Spore spore(*this);

    // Position the copy in the center.
    subtract_vectors_f(&C, C, spore.centroid);
    add_vectors_f(&(spore.begin_pt), spore.begin_pt, C);
    add_vectors_f(&(spore.centroid), spore.centroid, C);
    add_vectors_f(&(spore.end_pt), spore.end_pt, C);

    free_vector_f(C);

    // Clear the laterals and the apical so deleting the copy doesn't mess
    // things up.
    spore.apical = 0;
    spore.laterals.clear();

    num_rows = static_cast<uint32_t>(std::floor((2 * (center+y_pad)) / 
            y_scale + 0.5));
    num_cols = static_cast<uint32_t>(std::floor((2 * (center+x_pad)) / 
            x_scale + 0.5));

    create_zero_matrix_f(&M2, num_rows, num_cols);

    spore.draw_in_zeroed_matrix_f(M2, 0, 0, z, x_scale, y_scale, z_scale, fill);

    bool intersect = false;

    /* Copy the sub-matrix back into the target matrix. */
    row2 = static_cast<int64_t>(floor(
                (centroid->elts[1] - (center+y_pad) - y) / y_scale + 0.5));

    for (row = 0; row < num_rows; row++)
    {
        col2 = static_cast<int64_t>(floor(
                    (centroid->elts[0] - (center+x_pad) - x) / x_scale + 0.5));

        for (col = 0; col < num_cols; col++)
        {
            if ((col2 >= 0) && (col2 < M->num_cols) && 
                (row2 >= 0) && (row2 < M->num_rows))
            {
                if (M2->elts[ row ][ col ] > 0)
                {
                    uint32_t current_id = (uint32_t)floor(
                            M->elts[ row2 ][ col2 ]);

                    if (current_id > 0 && current_id != id && 
                            !has_child_id(current_id))
                    {
                        intersect = true;
                    }

                    M->elts[ row2 ][ col2 ] = M2->elts[ row ][ col ];
                }
            }
            col2++;
        }
        row2++; 
    }

    free_matrix_f(M2);

    return intersect;
}


bool Spore::draw_in_matblock_f
(
     jwsc::Matblock_f* M_blk, 
     float       x, 
     float       y,
     float       z,
     float       x_scale,
     float       y_scale,
     float       z_scale,
     bool        fill
)
const
{
    uint32_t num_mats;
    uint32_t mat;

    double mat_z;

    Matrix_f* M = 0;

    num_mats  = M_blk->num_mats;
    mat_z     = z;

    bool intersect = false;

    for (mat = 0; mat < num_mats; mat++)
    {
        mat_z = z + mat*z_scale;

        copy_matrix_from_matblock_f(&M, M_blk, MATBLOCK_MAT_MATRIX, mat);
        if (draw_in_matrix_f(M, x, y, mat_z, x_scale, y_scale, z_scale, fill))
        {
            intersect = true;
        }
        copy_matrix_into_matblock_f(&M_blk, M_blk, M, MATBLOCK_MAT_MATRIX, mat);
    }

    free_matrix_f(M);

    return intersect;
}


void Spore::draw_in_zeroed_matrix_f
(
     jwsc::Matrix_f* M, 
     float     x, 
     float     y,
     float     z,
     float     x_scale,
     float     y_scale,
     float     z_scale,
     bool      fill
) 
const
{
    uint32_t num_rows;
    uint32_t num_cols;
    uint32_t row, col;
    uint32_t left_edge;
    uint32_t middle;
    uint32_t elements_to_fill;
    uint32_t begin_fill_col;
    uint32_t end_fill_col;

    float minor, major;
    float M_x, M_y, M_z;
    float d, d_prev;

    Matrix_f* R_inv       = 0;
    Vector_f* point_world = 0;
    Vector_f* point_spore = 0;

    num_rows = M->num_rows;
    num_cols = M->num_cols;

    minor = 0.5*width;
    major = 0.5*length;

    /* Inverse rotation matrix. */
    transpose_matrix_f(&R_inv, R);

    create_vector_f(&point_world, 3);
    create_vector_f(&point_spore, 3);

    for (row = 0; row < num_rows; row++)
    {
        /* Get coordinates of matrix element in world units. Note that we
         * are using the center of matrix elements for position. */
        point_world->elts[ 0 ] = x + x_scale*((-1) + 0.5);
        point_world->elts[ 1 ] = y + y_scale*(row + 0.5);
        point_world->elts[ 2 ] = z + z_scale*0.5;

        /* Translate matrix element and ... */
        point_world->elts[ 0 ] -= centroid->elts[0];
        point_world->elts[ 1 ] -= centroid->elts[1];
        point_world->elts[ 2 ] -= centroid->elts[2];

        /* Rotate matrix element to spore space. */
        multiply_matrix_and_vector_f(&point_spore, R_inv, point_world);

        M_x = point_spore->elts[ 0 ];
        M_y = point_spore->elts[ 1 ];
        M_z = point_spore->elts[ 2 ];

        d_prev = (M_x*M_x)/(minor*minor) + 
                 (M_y*M_y)/(minor*minor) + 
                 (M_z*M_z)/(major*major);

        for (col = 0; col < num_cols; col++)
        {
            /* Get coordinates of matrix element in world units. Note that we
             * are using the center of matrix elements for position. */
            point_world->elts[ 0 ] = x + x_scale*(col + 0.5);
            point_world->elts[ 1 ] = y + y_scale*(row + 0.5);
            point_world->elts[ 2 ] = z + z_scale*0.5;

            /* Translate matrix element and ... */
            point_world->elts[ 0 ] -= centroid->elts[0];
            point_world->elts[ 1 ] -= centroid->elts[1];
            point_world->elts[ 2 ] -= centroid->elts[2];

            /* Rotate matrix element to spore space. */
            multiply_matrix_and_vector_f(&point_spore, R_inv, point_world);

            M_x = point_spore->elts[ 0 ];
            M_y = point_spore->elts[ 1 ];
            M_z = point_spore->elts[ 2 ];

            d = (M_x*M_x)/(minor*minor) + 
                (M_y*M_y)/(minor*minor) + 
                (M_z*M_z)/(major*major);

            if (((d - 1.0) <= 0 && (d_prev - 1.0) >= 0) ||
                ((d - 1.0) >= 0 && (d_prev - 1.0) <= 0))
            {
                M->elts[ row ][ col ] = id + opacity;
            }

            d_prev = d;
        }
    }
    for (col = 0; col < num_cols; col++)
    {
        /* Get coordinates of matrix element in world units. Note that we
         * are using the center of matrix elements for position. */
        point_world->elts[ 0 ] = x + x_scale*(col + 0.5);
        point_world->elts[ 1 ] = y + y_scale*((-1) + 0.5);
        point_world->elts[ 2 ] = z + z_scale*0.5;

        /* Translate matrix element and ... */
        point_world->elts[ 0 ] -= centroid->elts[0];
        point_world->elts[ 1 ] -= centroid->elts[1];
        point_world->elts[ 2 ] -= centroid->elts[2];

        /* Rotate matrix element to spore space. */
        multiply_matrix_and_vector_f(&point_spore, R_inv, point_world);

        M_x = point_spore->elts[ 0 ];
        M_y = point_spore->elts[ 1 ];
        M_z = point_spore->elts[ 2 ];

        d_prev = (M_x*M_x)/(minor*minor) + 
                 (M_y*M_y)/(minor*minor) + 
                 (M_z*M_z)/(major*major);

        for (row = 0; row < num_rows; row++)
        {
            /* Get coordinates of matrix element in world units. Note that we
             * are using the center of matrix elements for position. */
            point_world->elts[ 0 ] = x + x_scale*(col + 0.5);
            point_world->elts[ 1 ] = y + y_scale*(row + 0.5);
            point_world->elts[ 2 ] = z + z_scale*0.5;

            /* Translate matrix element and ... */
            point_world->elts[ 0 ] -= centroid->elts[0];
            point_world->elts[ 1 ] -= centroid->elts[1];
            point_world->elts[ 2 ] -= centroid->elts[2];

            /* Rotate matrix element to spore space. */
            multiply_matrix_and_vector_f(&point_spore, R_inv, point_world);

            M_x = point_spore->elts[ 0 ];
            M_y = point_spore->elts[ 1 ];
            M_z = point_spore->elts[ 2 ];

            d = (M_x*M_x)/(minor*minor) + 
                (M_y*M_y)/(minor*minor) + 
                (M_z*M_z)/(major*major);

            if (((d - 1.0) <= 0 && (d_prev - 1.0) >= 0) ||
                ((d - 1.0) >= 0 && (d_prev - 1.0) <= 0))
            {
                M->elts[ row ][ col ] = id + opacity;
            }

            d_prev = d;
        }
    }

    /* Fill in the spore. */
    if (fill)
    {
        for (row = 0; row < num_rows; row++)
        {
            left_edge        = 0;
            middle           = 0;
            elements_to_fill = 0;
            begin_fill_col   = 0;
            end_fill_col     = 0;

            for (col = 0; col < num_cols; col++)
            {
                if (left_edge && (M->elts[ row ][ col ] == 0))
                {
                    left_edge = 0;
                    middle = 1;
                    begin_fill_col = col;
                }
                else if (middle && (M->elts[ row ][ col ] > 0))
                {
                    elements_to_fill = 1;
                    end_fill_col = col;
                }
                else if (M->elts[ row ][ col ] > 0)
                {
                    left_edge = 1;
                }
            }

            if (elements_to_fill)
            {
                for (col = begin_fill_col; col < end_fill_col; col++)
                {
                    M->elts[ row ][ col ] = id + opacity;
                }
            }
        }
    }

    free_vector_f(point_world);
    free_vector_f(point_spore);
    free_matrix_f(R_inv);
}


#if defined ALTERNARIA_HAVE_OPENGL_FRAMEWORK || defined ALTERNARIA_HAVE_OPENGL
void Spore::draw_in_opengl(GLUquadric* quad, float scale) const
{
    static GLfloat amb_dif[4] = {0.7f, 0.6f, 0.2f, 1.0f};
    static GLfloat shininess[1] = {100.0f};
    static GLfloat r[16] = {0};
    static GLfloat s[16] = {0};


    glMatrixMode(GL_MODELVIEW);

    glPushMatrix();

    glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, amb_dif);
    glMaterialfv(GL_FRONT, GL_SPECULAR, amb_dif);
    glMaterialfv(GL_FRONT, GL_SHININESS, shininess);

    glTranslatef(scale*centroid->elts[0], scale*centroid->elts[1], 
            scale*centroid->elts[2]);

    const float* elts = *(R->elts);
    r[0] = elts[0];  r[4] = elts[1];  r[8]  = elts[2];
    r[1] = elts[3];  r[5] = elts[4];  r[9]  = elts[5];
    r[2] = elts[6];  r[6] = elts[7];  r[10] = elts[8];
    r[15] = 1;
    glMultMatrixf(r);

    s[0]  = scale*0.5f*width;
    s[5]  = scale*0.5f*width;
    s[10] = scale*0.5f*length;
    s[15] = 1.0f;
    glMultMatrixf(s);

    GLint slices = static_cast<GLint>(2.0f*width*scale);
    if (slices < 6)
        slices = 6;
    GLint stacks = static_cast<GLint>(length*scale);
    if (stacks < 3)
        stacks = 3;
    gluSphere(quad, 1, slices, stacks);

    glPopMatrix();
}
#endif


float Spore::get_intersection
(
    jwsc::Vector_f**      isect_out, 
    uint32_t              n,
    const jwsc::Vector_f* p1_in, 
    const jwsc::Vector_f* p2_in
)
throw (Arg_error)
{
    if (n != 1 && n != 2)
    {
        throw Arg_error("Intersection number is not 1 or 2");
    }

    if (p1_in->num_elts != 3 || p2_in->num_elts != 3)
    {
        throw Arg_error("Line end-points are not in R^3");
    }

    // Transform a copy of the vectors to bring them into spore coords.
    Vector_f* p1 = 0;
    Vector_f* p2 = 0;
    Matrix_f* M  = 0;
    subtract_vectors_f(&p1, p1_in, centroid);
    subtract_vectors_f(&p2, p2_in, centroid);
    transpose_matrix_f(&M, R);
    multiply_matrix_and_vector_f(&p1, M, p1);
    multiply_matrix_and_vector_f(&p2, M, p2);
    free_matrix_f(M);

    float x1 = p1->elts[0];
    float y1 = p1->elts[1];
    float z1 = p1->elts[2];

    free_vector_f(p1);

    float x2 = p2->elts[0];
    float y2 = p2->elts[1];
    float z2 = p2->elts[2];

    free_vector_f(p2);

    float i = x2 - x1;
    float j = y2 - y1;
    float k = z2 - z1;

    float A = (2.0f / width)*(2.0f / width);
    float B = A;
    float C = (2.0f / length)*(2.0f / length);

    float a = A*i*i + B*j*j + C*k*k;
    float b = 2*A*x1*i + 2*B*y1*j + 2*C*z1*k;
    float c = A*x1*x1 + B*y1*y1 + C*z1*z1 - 1;

    float discriminant = b*b - 4*a*c;
    float t;

    if (discriminant < 0.0f)
    {
        // no real roots
        return -1;
    }
    else if (discriminant > 1.0e-6f)
    {
        // two real roots
        float t1 = (-b + sqrt(discriminant)) / (2*a);
        float t2 = (-b - sqrt(discriminant)) / (2*a);

        if ((n == 1 && t1 < t2) || (n == 2 && t1 > t2))
        {
            t = t1;
        }
        else
        {
            t = t2;
        }
    }
    else
    {
        // one real root
        t = (-b + sqrt(discriminant)) / (2*a);
    }

    // Handle precision problems.
    if (t < 0)
        t = 0;
    else if (t > 1.0f)
        t = 1.0f;

    create_vector_f(isect_out, 3);
    Vector_f* isect = *isect_out;

    isect->elts[0] = x1 + i*t;
    isect->elts[1] = y1 + j*t;
    isect->elts[2] = z1 + k*t;

    // Transform the intersection point to put it back in world coords.
    multiply_matrix_and_vector_f(&isect, R, isect);
    add_vectors_f(&isect, isect, centroid);

    return t;
}


void Spore::update_log_prob()
{
    Apical_structure::update_log_prob();

    const DD_spore_set* dd_spores = density->get_dd_spores();

    if (dd_spores)
    {
        log_prob += dd_spores->get_similar_spore(this)->get_log_prob();

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
        assert(!std::isinf(log_prob) && !std::isnan(log_prob));
#endif
    }
}

/* -------------------------------------------------------------------------- */







/* ---------------------------  Printed_spore  ------------------------------ */

Printed_spore::Printed_spore(const Spore_density* density)
{
    id = 0;
    centroid_x = 0;
    centroid_y = 0;
    centroid_z = 0;
    length = 0;
    width = 0;
    base_theta = 0;
    base_psi = 0;
    theta = 0;
    psi = 0;
    opacity = 0;
    parent = 0;
    apical = 0;
    this->density = density;
}


Printed_spore::Printed_spore(std::istream& in, const Spore_density* density) 
    throw (IO_error, Arg_error)
{
    this->density = density;
    read(in);
}


void Printed_spore::read(std::istream& in) throw (IO_error, Arg_error)
{
    const char* member_value;

    if (!(member_value = Spore::get_member_value("Structure", in)))
    {
        throw Arg_error("Missing structure type");
    }

    if (strncmp(member_value, Spore::type_str(), strlen(Spore::type_str())) != 0)
    {
        ostringstream ost;
        ost << "Tried to read a '" << member_value << "' as a '"
            << Spore::type_str() << "'";
        throw Arg_error(ost.str());
    }

    read_id(in);
    read_level(in);
    read_centroid(in);
    read_length(in);
    read_width(in);
    read_base_theta(in);
    read_base_psi(in);
    read_theta(in);
    read_psi(in);
    read_opacity(in);
    read_parent(in);
    read_apical(in);
    read_laterals(in);
    read_log_prob(in);
}


Spore* Printed_spore::convert(Structure* parent) const throw (Arg_error)
{
    if (parent)
    {
        return new Spore(parent, length, width, base_theta, base_psi, theta, 
                psi, opacity, density);
    }

    return new Spore(centroid_x, centroid_y, centroid_z, length, width,
            base_theta, base_psi, theta, psi, opacity, level, density);
}

/* -------------------------------------------------------------------------- */