83 lines
2.6 KiB
C++
83 lines
2.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;
|
|
/*
|
|
* 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_m;
|
|
double m2 = m*m;
|
|
QuadraticSolver solver(Rdx * Rdx + Rdy * Rdy - m2 * Rdz * Rdz,
|
|
2.0 * (R0x*Rdx + R0y*Rdy - m2*R0z*Rdz - m*Rdz*m_b),
|
|
R0x*R0x + R0y*R0y
|
|
- m2*R0z*R0z - m_b*m_b - 2*m*R0z*m_b
|
|
);
|
|
|
|
Solver::Results results = solver.solve();
|
|
|
|
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
|
|
{
|
|
}
|
|
|
|
return m_transform.transform_normal(normal);
|
|
}
|