141 lines
4.6 KiB
C++
141 lines
4.6 KiB
C++
|
|
#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 */
|
|
}
|
|
|
|
Shape::IntersectList Cyl::intersect(const Ray & ray)
|
|
{
|
|
Ray ray_inv = m_inverse.transform_ray(ray);
|
|
IntersectList 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)
|
|
{
|
|
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.push_back(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)
|
|
{
|
|
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.push_back(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.push_back(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 (m_slope == 0.0)
|
|
{
|
|
normal = Vector(local_pt[0], local_pt[1], 0.0);
|
|
normal.normalize();
|
|
}
|
|
else
|
|
{
|
|
double norm_slope = -1.0 / m_slope;
|
|
double z = sqrt(local_pt[0]*local_pt[0] + local_pt[1]*local_pt[1])
|
|
* norm_slope;
|
|
normal = Vector(local_pt[0], local_pt[1], z);
|
|
normal.normalize();
|
|
}
|
|
}
|
|
|
|
return m_transform.transform_normal(normal);
|
|
}
|