fart/shapes/Cyl.cc

#include "Cyl.h"
#include "util/Solver.h"
#include <math.h>
#include <iostream>
using namespace std;

#define FP_EQUAL(x,y) (fabs((x)-(y)) < 0.000001)

Cyl::Cyl(double bottom_radius, double top_radius, double height)
{
    m_bottom_radius = fabs(bottom_radius);
    m_bottom_radius_2 = bottom_radius * bottom_radius;
    m_top_radius = fabs(top_radius);
    m_top_radius_2 = top_radius * top_radius;
    m_height = fabs(height);
    if (m_height == 0.0)
        m_height = 1.0;
    m_slope = (m_top_radius - m_bottom_radius) / m_height; /* rise over run */
    if ( ! FP_EQUAL(m_slope, 0.0) )
        m_inv_slope = -1.0 / m_slope;
}

Shape::IntersectionList Cyl::intersect(refptr<Shape> _this, const Ray & ray)
{
    Ray ray_inv = m_inverse.transform_ray(ray);
    IntersectionList res;

    /* First intersect with the bottom plane, if it has positive area */
    if (m_bottom_radius > 0.0)
    {
        LinearSolver solver(-ray_inv.getDirection()[2],
                            -ray_inv.getOrigin()[2]);
        Solver::Result solutions = solver.solve();
        if (solutions.numResults > 0 && solutions.results[0] > 0.0)
        {
            Vector isect_point = ray_inv[solutions.results[0]];
            if (isect_point[0]*isect_point[0] + isect_point[1]*isect_point[1]
                    <= m_bottom_radius_2)
            {
                res.add(Intersection(_this,
                                     m_transform.transform_point(isect_point)));
            }
        }
    }

    /* Same for the top plane */
    if (m_top_radius > 0.0)
    {
        LinearSolver solver(ray_inv.getDirection()[2],
                            ray_inv.getOrigin()[2] - m_height);
        Solver::Result solutions = solver.solve();
        if (solutions.numResults > 0 && solutions.results[0] > 0.0)
        {
            Vector isect_point = ray_inv[solutions.results[0]];
            if (isect_point[0]*isect_point[0] + isect_point[1]*isect_point[1]
                    <= m_top_radius_2)
            {
                res.add(Intersection(_this,
                                     m_transform.transform_point(isect_point)));
            }
        }
    }

    /*
     * Now see if the ray hit the side of the cylinder/cone thingy
     * Ray equation:   R = R0 + tRd
     *                 x = R0x + tRdx
     *                 y = R0y + tRdy
     *                 z = R0z + tRdz
     * Side equation:  x^2 + y^2 = (m*z + b)^2
     * Combined:       (R0x+t*Rdx)^2 + (R0y+t*Rdy)^2 = (m*(R0z+t*Rdz) + b)^2
     * Expanded:       (R0x*R0x + (t*t)*Rdx*Rdx + 2*R0x*t*Rdx)
     *                 + (R0y*R0y + (t*t)*Rdy*Rdy + 2*R0y*t*Rdy)
     *                 = ((m*R0z+m*t*Rdz)^2 + b*b + 2*m*(R0z+t*Rdz)*b)
     *                 = ((m*R0z*m*R0z + m*m*t*t*Rdz*Rdz + 2*m*R0z*m*t*Rdz)
     *                  + b*b + 2*m*R0z*b + 2*m*t*Rdz*b)
     * Quadratic form: (t*t)(Rdx*Rdx + Rdy*Rdy - m*m*Rdz*Rdz)
     *                 + (t)(2*R0x*Rdx + 2*R0y*Rdy - 2*m*R0z*m*Rdz - 2*m*Rdz*b)
     *                 + (R0x*R0x + R0y*R0y - m*m*R0z*R0z - b*b - 2*m*R0z*b)
     *                 = 0
     */
    double Rdx = ray_inv.getDirection()[0];
    double Rdy = ray_inv.getDirection()[1];
    double Rdz = ray_inv.getDirection()[2];
    double R0x = ray_inv.getOrigin()[0];
    double R0y = ray_inv.getOrigin()[1];
    double R0z = ray_inv.getOrigin()[2];
    double m = m_slope;
    double m2 = m*m;
    double b = m_bottom_radius;
    QuadraticSolver solver(Rdx * Rdx + Rdy * Rdy - m2 * Rdz * Rdz,
                           2.0 * (R0x*Rdx + R0y*Rdy - m2*R0z*Rdz - m*Rdz*b),
                           R0x*R0x + R0y*R0y - m2*R0z*R0z - b*b - 2*m*R0z*b);

    Solver::Result solutions = solver.solve();
    for (int i = 0; i < solutions.numResults; i++)
    {
        if (solutions.results[i] >= 0.0)
        {
            Vector isect_point = ray_inv[solutions.results[i]];
            if (isect_point[2] > 0.0 && isect_point[2] < m_height)
            {
                res.add(Intersection(_this,
                                     m_transform.transform_point(isect_point)));
            }
        }
    }

    return res;
}

Vector Cyl::getNormalAt(const Vector & pt)
{
    Vector local_pt = m_inverse.transform_point(pt);

    Vector normal;

    if      ( FP_EQUAL(local_pt[2], 0.0) && m_bottom_radius > 0.0 )
    {
        normal = Vector(0, 0, -1);
    }
    else if ( FP_EQUAL(local_pt[2], m_height) && m_top_radius > 0.0 )
    {
        normal = Vector(0, 0, 1);
    }
    else
    {
        if ( FP_EQUAL(m_slope, 0.0) )
        {
            normal = Vector(local_pt[0], local_pt[1], 0.0);
        }
        else
        {
            double x = local_pt[0];
            double y = local_pt[1];
            double dist = sqrt(  local_pt[0] * local_pt[0]
                               + local_pt[1] * local_pt[1] );
            if (dist > 0.0)
            {
                double scale = 1.0 / dist;
                x *= scale;
                y *= scale;
            }
            normal = Vector(x, y, m_inv_slope);
        }
        normal.normalize();
    }

    return m_transform.transform_normal(normal);
}