2 This file is part of the n3dmath2 library.
4 Copyright (c) 2004, 2005 John Tsiombikas <nuclear@siggraph.org>
6 The n3dmath2 library is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
11 The n3dmath2 library is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
16 You should have received a copy of the GNU General Public License
17 along with the n3dmath2 library; if not, write to the Free Software
18 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21 #include "n3dmath2_qdr.hpp"
23 Surface::Surface(const Vector3 &pos) {
27 Surface::~Surface() {}
29 void Surface::set_position(const Vector3 &pos) {
33 Vector3 Surface::get_position() const {
37 void Surface::set_rotation(const Quaternion &rot) {
41 Quaternion Surface::get_rotation() const {
45 //////////////// sphere implementation ///////////////
47 Sphere::Sphere(const Vector3 &pos, scalar_t rad) : Surface(pos) {
53 void Sphere::set_radius(scalar_t rad) {
57 scalar_t Sphere::get_radius() const {
61 Vector2 Sphere::inv_map(const Vector3 &pt) const {
62 Vector3 normal = (pt - pos) / radius;
63 Vector3 pole = Vector3(0, 1, 0).transformed(rot);
64 Vector3 equator = Vector3(0, 0, 1).transformed(rot);
67 scalar_t phi = acos(dot_product(normal, pole));
70 if(imap.y < xsmall_number || 1.0 - imap.y < xsmall_number) {
75 scalar_t theta = acos(dot_product(equator, normal) / sin(phi)) / two_pi;
77 imap.x = (dot_product(cross_product(pole, equator), normal) < 0.0) ? theta : 1.0 - theta;
82 bool Sphere::check_intersection(const Ray &ray) const {
83 return find_intersection(ray, 0);
86 bool Sphere::find_intersection(const Ray &ray, SurfPoint *isect) const {
87 // find terms of the quadratic equation
88 scalar_t a = SQ(ray.dir.x) + SQ(ray.dir.y) + SQ(ray.dir.z);
89 scalar_t b = 2.0 * ray.dir.x * (ray.origin.x - pos.x) +
90 2.0 * ray.dir.y * (ray.origin.y - pos.y) +
91 2.0 * ray.dir.z * (ray.origin.z - pos.z);
92 scalar_t c = SQ(pos.x) + SQ(pos.y) + SQ(pos.z) +
93 SQ(ray.origin.x) + SQ(ray.origin.y) + SQ(ray.origin.z) +
94 2.0 * (-pos.x * ray.origin.x - pos.y * ray.origin.y - pos.z * ray.origin.z) - SQ(radius);
96 // find the discriminant
97 scalar_t d = SQ(b) - 4.0 * a * c;
98 if(d < 0.0) return false;
101 scalar_t sqrt_d = sqrt(d);
102 scalar_t t1 = (-b + sqrt_d) / (2.0 * a);
103 scalar_t t2 = (-b - sqrt_d) / (2.0 * a);
105 if(t1 < error_margin && t2 < error_margin) return false;
108 if(t1 < error_margin) t1 = t2;
109 if(t2 < error_margin) t2 = t1;
110 isect->t = t1 < t2 ? t1 : t2;
111 isect->pos = ray.origin + ray.dir * isect->t;
113 isect->normal = (isect->pos - pos) / radius;
114 isect->pre_ior = ray.ior;
115 //isect->post_ior = mat.ior;
122 Plane::Plane(const Vector3 &pos, const Vector3 &normal) : Surface(pos) {
123 this->normal = normal;
126 Plane::Plane(const Vector3 &p1, const Vector3 &p2, const Vector3 &p3) : Surface(p1)
128 normal = cross_product(p2 - p1, p3 - p1).normalized();
133 void Plane::set_normal(const Vector3 &normal) {
134 this->normal = normal;
137 Vector3 Plane::get_normal() const {
141 Vector2 Plane::inv_map(const Vector3 &pt) const {
142 static int dbg; dbg++;
143 if(dbg == 1) std::cerr << "inverse mapping for planes not implemented yet!\n";
144 return Vector2(0, 0);
147 bool Plane::check_intersection(const Ray &ray) const {
148 return find_intersection(ray, 0);
151 bool Plane::find_intersection(const Ray &ray, SurfPoint *isect) const {
152 scalar_t normal_dot_dir = dot_product(normal, ray.dir);
153 if(fabs(normal_dot_dir) < error_margin) return false;
155 // TODO: this is only correct if pos is the projection of the origin on the plane
156 scalar_t d = pos.length();
158 scalar_t t = -(dot_product(normal, ray.origin) + d) / normal_dot_dir;
160 if(t < error_margin) return false;
163 isect->pos = ray.origin + ray.dir * t;
164 isect->normal = normal;
166 isect->pre_ior = ray.ior;
167 //isect->post_ior = mat.ior;
173 Box::Box(const Vector3 &min_vec, const Vector3 &max_vec) {
174 verts[0] = Vector3(min_vec.x, max_vec.y, min_vec.z);
175 verts[1] = Vector3(max_vec.x, max_vec.y, min_vec.z);
176 verts[2] = Vector3(max_vec.x, min_vec.y, min_vec.z);
177 verts[3] = Vector3(min_vec.x, min_vec.y, min_vec.z);
179 verts[4] = Vector3(min_vec.x, max_vec.y, max_vec.z);
180 verts[5] = Vector3(max_vec.x, max_vec.y, max_vec.z);
181 verts[6] = Vector3(max_vec.x, min_vec.y, max_vec.z);
182 verts[7] = Vector3(min_vec.x, min_vec.y, max_vec.z);
185 Box::Box(const Vector3 &v0, const Vector3 &v1, const Vector3 &v2, const Vector3 &v3,
186 const Vector3 &v4, const Vector3 &v5, const Vector3 &v6, const Vector3 &v7) {
197 Box::Box(const Vector3 *array) {
198 memcpy(verts, array, sizeof verts);
201 Vector2 Box::inv_map(const Vector3 &pt) const {
202 return Vector2(); // TODO: implement
205 bool Box::check_intersection(const Ray &ray) const {
206 return false; // TODO: implement
209 bool Box::find_intersection(const Ray &ray, SurfPoint *isect) const {
210 return false; // TODO: implement
215 * PointOverPlane (MG)
216 * returns true if the point is in the positive side of the
219 bool point_over_plane(const Plane &plane, const Vector3 &point)
221 if (dot_product(plane.get_position() - point, plane.get_normal()) < 0)
230 static inline bool check_correct_winding(const Vector3 &v0, const Vector3 &v1, const Vector3 &v2, const Vector3 &normal) {
231 Vector3 tri_normal = cross_product(v1 - v0, v2 - v0);
232 return dot_product(tri_normal, normal) > 0.0;
236 bool check_tri_ray_intersection(const Vector3 &v0, const Vector3 &v1, const Vector3 &v2, const Ray &ray) {
237 // find the cosine of the angle between the normal and the line
238 Vector3 normal = cross_product(v1 - v0, v2 - v0);
239 scalar_t dot = dot_product(normal, ray.dir);
241 if(fabs(dot) < small_number) {
245 // further testing to verify intersection in the area of the triangle
246 scalar_t t = -(normal.x*(ray.origin.x - v0.x) + normal.y*(ray.origin.y - v0.y) + normal.z*(ray.origin.z - v0.z)) / (normal.x * ray.dir.x + normal.y * ray.dir.y + normal.z * ray.dir.z);
247 if(t > small_number) {
248 Vector3 intersect(ray.origin.x + ray.dir.x*t, ray.origin.y + ray.dir.y*t, ray.origin.z + ray.dir.z*t);
250 if(check_correct_winding(v0, v1, intersect, normal)
251 && check_correct_winding(v1, v2, intersect, normal)
252 && check_correct_winding(v2, v0, intersect, normal)) {