hypha.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 <cassert>
#include <cmath>
#include <cstring>
#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 <jwsc++/base/exception.h>

#include "density.h"
#include "printable.h"
#include "structure.h"
#include "spore.h"
#include "hypha.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  APICAL_LENGTH_MU       15.0f
#define  APICAL_LENGTH_SIGMA    1.0f
#define  APICAL_LENGTH_MIN      10.f
#define  APICAL_LENGTH_MAX      100.0f

#define  LATERAL_LENGTH_MU      15.0f
#define  LATERAL_LENGTH_SIGMA   1.0f
#define  LATERAL_LENGTH_MIN     10.f
#define  LATERAL_LENGTH_MAX     100.0f

#define  APICAL_WIDTH_MU        6.0f
#define  APICAL_WIDTH_SIGMA     0.3f
#define  APICAL_WIDTH_MIN       1.0f
#define  APICAL_WIDTH_MAX       20.0f

#define  APICAL_DWIDTH_SIGMA    0.05f
#define  APICAL_DWIDTH_MIN      -2.0f
#define  APICAL_DWIDTH_MAX      2.0f

#define  LATERAL_WIDTH_MU       6.0f
#define  LATERAL_WIDTH_SIGMA    0.3f
#define  LATERAL_WIDTH_MIN      1.0f
#define  LATERAL_WIDTH_MAX      20.0f

#define  APICAL_THETA_MU        0
#define  APICAL_THETA_SIGMA     0.2f
#define  APICAL_THETA_MIN       0
#define  APICAL_THETA_MAX       JWSC_PI

#define  LATERAL_THETA_MU       JWSC_PI_4
#define  LATERAL_THETA_SIGMA    0.2f
#define  LATERAL_THETA_MIN      0
#define  LATERAL_THETA_MAX      JWSC_PI

#define  PSI_MIN                -JWSC_PI
#define  PSI_MAX                JWSC_PI

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

#define  LAT_DIST_MU            0.5f
#define  LAT_DIST_SIGMA         0.1f
#define  LAT_DIST_MIN           0
#define  LAT_DIST_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;


/* -----------------------  Apical_hypha_density  --------------------------- */

Apical_hypha_density::Apical_hypha_density() 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    = APICAL_LENGTH_MU;
    sigma = APICAL_LENGTH_SIGMA;
    min   = APICAL_LENGTH_MIN;
    max   = APICAL_LENGTH_MAX;
    set_length(mu, sigma, min, max);

    mu    = APICAL_WIDTH_MU;
    sigma = APICAL_WIDTH_SIGMA;
    min   = APICAL_WIDTH_MIN;
    max   = APICAL_WIDTH_MAX;
    set_width(mu, sigma, min, max);

    sigma = APICAL_DWIDTH_SIGMA;
    min   = APICAL_DWIDTH_MIN;
    max   = APICAL_DWIDTH_MAX;
    set_dwidth(sigma, min, max);

    mu    = APICAL_THETA_MU;
    sigma = APICAL_THETA_SIGMA;
    min   = APICAL_THETA_MIN;
    max   = APICAL_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);
}


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

    Apical_hypha_density* h  = new Apical_hypha_density();

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

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

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

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

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

    sigma = dwidth.get_sigma();
    min   = dwidth.get_min();
    max   = dwidth.get_max();
    h->set_dwidth(sigma, min, max);

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

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

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

    return h;
}


void Apical_hypha_density::set_dwidth
(
    float sigma, 
    float min, 
    float max
) 
throw (Arg_error)
{
    dwidth.set(0, sigma, min, max);
}


Apical_hypha* Apical_hypha_density::sample
(
    float base_theta, 
    float base_psi,
    size_t level
) 
const
{
    return new Apical_hypha(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);
}


Apical_hypha* Apical_hypha_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 Apical_hypha(centroid_x, centroid_y, centroid_z, sample_length(),
            sample_width(), base_theta, base_psi, sample_theta(), sample_psi(), 
            sample_opacity(), level, this);
}


Apical_hypha* Apical_hypha_density::sample
(
    Structure* parent,
    float      base_theta,
    float      base_psi
) 
const throw (Arg_error)
{
    if (dynamic_cast<Apical_hypha*>(parent) || 
        dynamic_cast<Lateral_hypha*>(parent))
    {
        return new Apical_hypha(parent, sample_length(), 
                parent->get_width() + sample_dwidth(), base_theta, base_psi,
                sample_theta(), sample_psi(), sample_opacity(), this);
    }

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


float Apical_hypha_density::sample_dwidth() const
{
    return (float)sample_gaussian_pdf_d(dwidth.get_mu(), dwidth.get_sigma(), 
            dwidth.get_min(), dwidth.get_max());
}

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






/* ----------------------  Lateral_hypha_density  --------------------------- */

Lateral_hypha_density::Lateral_hypha_density() 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    = LATERAL_LENGTH_MU;
    sigma = LATERAL_LENGTH_SIGMA;
    min   = LATERAL_LENGTH_MIN;
    max   = LATERAL_LENGTH_MAX;
    set_length(mu, sigma, min, max);

    mu    = LATERAL_WIDTH_MU;
    sigma = LATERAL_WIDTH_SIGMA;
    min   = LATERAL_WIDTH_MIN;
    max   = LATERAL_WIDTH_MAX;
    set_width(mu, sigma, min, max);

    mu    = LATERAL_THETA_MU;
    sigma = LATERAL_THETA_SIGMA;
    min   = LATERAL_THETA_MIN;
    max   = LATERAL_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);

    mu    = LAT_DIST_MU;
    sigma = LAT_DIST_SIGMA;
    min   = LAT_DIST_MIN;
    max   = LAT_DIST_MAX;
    set_lat_dist(mu, sigma, min, max);
}


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

    Lateral_hypha_density* h  = new Lateral_hypha_density();

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

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

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

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

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

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

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

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

    mu    = lat_dist.get_mu();
    sigma = lat_dist.get_sigma();
    min   = lat_dist.get_min();
    max   = lat_dist.get_max();
    h->set_lat_dist(mu, sigma, min, max);

    return h;
}


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


Lateral_hypha* Lateral_hypha_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 Lateral_hypha(centroid_x, centroid_y, centroid_z, 
            sample_lat_dist(), sample_length(), sample_width(), base_theta,
            base_psi, sample_theta(), sample_psi(), sample_opacity(), level,
            this);
}


Lateral_hypha* Lateral_hypha_density::sample
(
    Structure* parent,
    float      base_theta,
    float      base_psi
) 
const throw (Arg_error)
{
    return new Lateral_hypha(parent, sample_lat_dist(), sample_length(), 
            sample_width(), base_theta, base_psi, sample_theta(), sample_psi(), 
            sample_opacity(), this);
}

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







/* ----------------------------  Apical_hypha  ------------------------------ */

Apical_hypha::Apical_hypha
(
    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 Apical_hypha_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();
}


Apical_hypha::Apical_hypha
(
    Structure* parent,
    float      length,
    float      width,
    float      base_theta,
    float      base_psi,
    float      theta,
    float      psi,
    float      opacity,
    const class Apical_hypha_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();
}


Apical_hypha::Apical_hypha(const Apical_hypha& h) : Apical_structure(h)
{
    density = h.density;
}


Apical_hypha::Apical_hypha
(
    const class Spore&          s, 
    float                       width,
    const Apical_hypha_density* density
)
throw (Arg_error) : Apical_structure(s)
{
    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();
}


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


Apical_hypha& Apical_hypha::operator= (const Apical_hypha& h)
{
    Apical_structure::operator=(h);

    density = h.density;

    return *this;
}


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


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

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

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


Apical_hypha* Apical_hypha::merge_with_apical
(
    class Spore**        spore_out,
    const class Spore*   rvals
) 
throw (Arg_error)
{
    assert(*spore_out == 0);

    if (!apical)
    {
        throw Arg_error("No apical to merge with");
    }

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

    Apical_hypha hypha(*rvals, width, density);

    Apical_hypha* merged = (Apical_hypha*)Apical_structure::merge_with_apical(
            &hypha);

    *spore_out = new Spore(*this, rvals->get_width(), rvals->get_density());

    try
    {
        this->replace(*spore_out);
    }
    catch (Arg_error e)
    {
        delete *spore_out;
        delete merged;
        throw e;
    }

    return merged;
}


Lateral_hypha* Apical_hypha::split_into_lateral
(
    const Apical_hypha*        rvals_1,
    const class Lateral_hypha* rvals_2
) 
throw (Arg_error)
{
    return (Lateral_hypha*)Apical_structure::split_into_lateral(rvals_1, 
            rvals_2);
}


Lateral_hypha* Apical_hypha::merge_with_lateral
(
    class Lateral_hypha* lateral,
    const Apical_hypha*  rvals
) 
throw (Arg_error)
{
    return (Lateral_hypha*)Apical_structure::merge_with_lateral(lateral,rvals);
}


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


bool Apical_hypha::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, num_rows, num_cols;
    uint32_t mat, row, col;

    float height, radius;
    float min_scale;
    float phi;
    float phi_step;
    float h;
    float h_step;
    float r;
    float r_step;
    float x_0, y_0, z_0;

    Vector_f* point = 0;

    num_mats = M_blk->num_mats;
    num_rows = M_blk->num_rows;
    num_cols = M_blk->num_cols;

    height = length;
    radius = 0.5*width;

    create_vector_f(&point, 3);

    if (x_scale < y_scale && x_scale < z_scale)
        min_scale = x_scale;
    else if (y_scale < z_scale)
        min_scale = y_scale;
    else
        min_scale = z_scale;

    h_step = 0.5*min_scale;
    r_step = (fill) ? 0.5*min_scale : radius;
    phi_step = 0.5*min_scale / radius;

    bool intersect = false;

    // Iterate over the points on the hypha parametrically and paint each
    // point in the matrix block.
    for (r = radius; r > 0; r -= r_step)
    {
        for (phi = 0; phi < 2*JWSC_PI; phi += phi_step)
        {
            x_0 = r*std::cos(phi);
            y_0 = r*std::sin(phi);

            for (h = 0; h <= height; h += h_step)
            {
                z_0 = h;

                point->elts[0] = x_0;
                point->elts[1] = y_0;
                point->elts[2] = z_0;
                multiply_matrix_and_vector_f(&point, R, point);
                add_vectors_f(&point, point, begin_pt);

                col = static_cast<uint32_t>(std::floor(
                            (point->elts[0] - x) / x_scale + 0.5));
                row = static_cast<uint32_t>(std::floor(
                            (point->elts[1] - y) / y_scale + 0.5));
                mat = static_cast<uint32_t>(std::floor(
                            (point->elts[2] - z) / z_scale + 0.5));

                if (mat >= 0 && mat < num_mats &&
                    row >= 0 && row < num_rows && 
                    col >= 0 && col < num_cols)
                {
                    uint32_t current_id = (uint32_t)floor(
                            M_blk->elts[ mat ][ row ][ col ]);

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

                    M_blk->elts[ mat ][ row ][ col ] = (float)id + opacity;
                }
            }
        }
    }

    free_vector_f(point);

    return intersect;
}


#if defined ALTERNARIA_HAVE_OPENGL_FRAMEWORK || defined ALTERNARIA_HAVE_OPENGL
void Apical_hypha::draw_in_opengl(GLUquadric* quad, float scale) const
{
    static GLfloat amb_dif[4] = {0.5f, 0.7f, 0.6f, 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*begin_pt->elts[0], scale*begin_pt->elts[1], 
            scale*begin_pt->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.0f;
    glMultMatrixf(r);

    s[0]  = scale*0.5f*width;
    s[5]  = scale*0.5f*width;
    s[10] = scale*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 < 1)
        stacks = 1;
    gluCylinder(quad, 1, 1, 1, slices, stacks);

    glPopMatrix();
}
#endif


float Apical_hypha::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 hypha 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 r = 0.5f*width;

    float a = i*i + j*j;
    float b = 2*x1*i + 2*y1*j;
    float c = x1*x1 + y1*y1 - r*r;

    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;

    float x = x1 + i*t;
    float y = y1 + j*t;
    float z = z1 + k*t;

    if (fabs(z) > 0.5f*length)
    {
        return -1;
    }

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

    isect->elts[0] = x;
    isect->elts[1] = y;
    isect->elts[2] = z;

    // 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 Apical_hypha::update_log_prob()
{
    Apical_structure::update_log_prob();

    if (parent && (dynamic_cast<Apical_hypha*>(parent) || 
                   dynamic_cast<Lateral_hypha*>(parent)))
    {
        float    mu, sigma;
        float    min, max;
        double   delta;
        double   p_1, p_2;

        /* Take out the width log prob. */
        mu        = density->get_width().get_mu();
        sigma     = density->get_width().get_sigma();
        min       = density->get_width().get_min();
        max       = density->get_width().get_max();
        delta     = (max - min) / 1000;
        p_1       = integrate_gaussian_pdf_d(mu, sigma, width, width+delta); 
        p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
        log_prob -= log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
        log_prob += log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

        float w = width - parent->get_width();

        /* Put in the differential width log prob. */
        mu        = density->get_dwidth().get_mu();
        sigma     = density->get_dwidth().get_sigma();
        min       = density->get_dwidth().get_min();
        max       = density->get_dwidth().get_max();
        delta     = (max - min) / 1000;
        p_1       = integrate_gaussian_pdf_d(mu, sigma, w, w+delta); 
        p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
        log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
        log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

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

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







/* -----------------------  Printed_apical_hypha  --------------------------- */

Printed_apical_hypha::Printed_apical_hypha(const Apical_hypha_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;
    level = 0;
    parent = 0;
    apical = 0;
    this->density = density;
}


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


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

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

    if (strncmp(member_value, Apical_hypha::type_str(), 
                strlen(Apical_hypha::type_str())) != 0)
    {
        ostringstream ost;
        ost << "Tried to read a '" << member_value << "' as a '"
            << Apical_hypha::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);
}


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

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

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







/* ---------------------------  Lateral_hypha  ------------------------------ */

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

    update_ancestry_log_prob();
}


Lateral_hypha::Lateral_hypha
(
    Structure* parent,
    float      lat_dist,
    float      length,
    float      width,
    float      base_theta,
    float      base_psi,
    float      theta,
    float      psi,
    float      opacity,
    const class Lateral_hypha_density* density
)
throw (Arg_error, Apical_error)
    : Lateral_structure(parent, lat_dist, length, width, base_theta, base_psi,
            theta, psi, opacity, density)
{
    this->density = density;

    update_ancestry_log_prob();
}


Lateral_hypha::Lateral_hypha(const Lateral_hypha& h) : Lateral_structure(h)
{
    density = h.density;
}


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


Lateral_hypha& Lateral_hypha::operator= (const Lateral_hypha& h)
{
    Lateral_structure::operator=(h);

    density = h.density;

    return *this;
}


bool Lateral_hypha::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, num_rows, num_cols;
    uint32_t mat, row, col;

    float height, radius;
    float min_scale;
    float phi;
    float phi_step;
    float h;
    float h_step;
    float r;
    float r_step;
    float x_0, y_0, z_0;

    Vector_f* point = 0;

    num_mats = M_blk->num_mats;
    num_rows = M_blk->num_rows;
    num_cols = M_blk->num_cols;

    height = length;
    radius = 0.5*width;

    create_vector_f(&point, 3);

    if (x_scale < y_scale && x_scale < z_scale)
        min_scale = x_scale;
    else if (y_scale < z_scale)
        min_scale = y_scale;
    else
        min_scale = z_scale;

    h_step = 0.5*min_scale;
    r_step = (fill) ? 0.5*min_scale : radius;
    phi_step = 0.5*min_scale / radius;

    bool intersect = false;

    // Iterate over the points on the hypha parametrically and paint each
    // point in the matrix block.
    for (r = radius; r > 0; r -= r_step)
    {
        for (phi = 0; phi < 2*JWSC_PI; phi += phi_step)
        {
            x_0 = r*std::cos(phi);
            y_0 = r*std::sin(phi);

            for (h = 0; h <= height; h += h_step)
            {
                z_0 = h;

                point->elts[0] = x_0;
                point->elts[1] = y_0;
                point->elts[2] = z_0;
                multiply_matrix_and_vector_f(&point, R, point);
                add_vectors_f(&point, point, begin_pt);

                col = static_cast<uint32_t>(std::floor(
                            (point->elts[0] - x) / x_scale + 0.5));
                row = static_cast<uint32_t>(std::floor(
                            (point->elts[1] - y) / y_scale + 0.5));
                mat = static_cast<uint32_t>(std::floor(
                            (point->elts[2] - z) / z_scale + 0.5));


                if (mat >= 0 && mat < num_mats &&
                    row >= 0 && row < num_rows && 
                    col >= 0 && col < num_cols)
                {
                    uint32_t current_id = (uint32_t)floor(
                            M_blk->elts[ mat ][ row ][ col ]);

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

                    M_blk->elts[ mat ][ row ][ col ] = id + opacity;
                }
            }
        }
    }

    free_vector_f(point);

    return intersect;
}


#if defined ALTERNARIA_HAVE_OPENGL_FRAMEWORK || defined ALTERNARIA_HAVE_OPENGL
void Lateral_hypha::draw_in_opengl(GLUquadric* quad, float scale) const
{
    static GLfloat amb_dif[4] = {0.7f, 0.8f, 0.9f, 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*begin_pt->elts[0], scale*begin_pt->elts[1], 
            scale*begin_pt->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.0f;
    glMultMatrixf(r);

    s[0]  = scale*0.5f*width;
    s[5]  = scale*0.5f*width;
    s[10] = scale*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 < 1)
        stacks = 1;
    gluCylinder(quad, 1, 1, 1, slices, stacks);

    glPopMatrix();
}
#endif


float Lateral_hypha::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 hypha 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 r = 0.5f*width;

    float a = i*i + j*j;
    float b = 2*x1*i + 2*y1*j;
    float c = x1*x1 + y1*y1 - r*r;

    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;

    float x = x1 + i*t;
    float y = y1 + j*t;
    float z = z1 + k*t;

    if (fabs(z) > 0.5f*length)
    {
        return -1;
    }

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

    isect->elts[0] = x;
    isect->elts[1] = y;
    isect->elts[2] = z;

    // 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;
}

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








/* -----------------------  Printed_lateral_hypha  -------------------------- */

Printed_lateral_hypha::Printed_lateral_hypha
(
    const Lateral_hypha_density* level_1_density,
    const Lateral_hypha_density* level_n_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;
    level = 0;
    parent = 0;
    apical = 0;
    this->level_1_density = level_1_density;
    this->level_n_density = level_n_density;
}


Printed_lateral_hypha::Printed_lateral_hypha
(
    std::istream& in,
    const Lateral_hypha_density* level_1_density,
    const Lateral_hypha_density* level_n_density
)
throw (IO_error, Arg_error)
{
    this->level_1_density = level_1_density;
    this->level_n_density = level_n_density;
    read(in);
}


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

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

    if (strncmp(member_value, Lateral_hypha::type_str(), 
                strlen(Lateral_hypha::type_str())) != 0)
    {
        ostringstream ost;
        ost << "Tried to read a '" << member_value << "' as a '"
            << Lateral_hypha::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);
    read_lat_dist(in);
}


Lateral_hypha* Printed_lateral_hypha::convert(Structure* parent) const 
    throw (Arg_error)
{
    const Lateral_hypha_density* density = (level == 1) ? level_1_density 
                                                        : level_n_density;

    if (parent)
    {
        return new Lateral_hypha(parent, lat_dist, length, width, base_theta,
                base_psi, theta, psi, opacity, density);
    }

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

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