#include "Cyl.h" #include "util/Solver.h" #include #include 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 _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), m_transform.transform_normal(Vector(0, 0, -1)))); } } } /* 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), m_transform.transform_normal(Vector(0, 0, 1)))); } } } /* * 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) { Vector normal; if ( FP_EQUAL(m_slope, 0.0) ) { normal = Vector(isect_point[0], isect_point[1], 0.0); } else { double x = isect_point[0]; double y = isect_point[1]; double dist = sqrt( isect_point[0] * isect_point[0] + isect_point[1] * isect_point[1] ); if (dist > 0.0) { double scale = 1.0 / dist; x *= scale; y *= scale; } normal = Vector(x, y, m_inv_slope); } normal.normalize(); res.add(Intersection(_this, m_transform.transform_point(isect_point), m_transform.transform_normal(normal))); } } } return res; }