git.mutantstargoat.com Git - summerhack/blob - src/3dengfx/src/n3dmath2/n3dmath2_qdr.cpp

   1 /*
   2 This file is part of the n3dmath2 library.
   3
   4 Copyright (c) 2004, 2005 John Tsiombikas <nuclear@siggraph.org>
   5
   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.
  10
  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.
  15
  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
  19 */
  20
  21 #include "n3dmath2_qdr.hpp"
  22
  23 Surface::Surface(const Vector3 &pos) {
  24         this->pos = pos;
  25 }
  26
  27 Surface::~Surface() {}
  28
  29 void Surface::set_position(const Vector3 &pos) {
  30         this->pos = pos;
  31 }
  32
  33 Vector3 Surface::get_position() const {
  34         return pos;
  35 }
  36
  37 void Surface::set_rotation(const Quaternion &rot) {
  38         this->rot = rot;
  39 }
  40
  41 Quaternion Surface::get_rotation() const {
  42         return rot;
  43 }
  44
  45 //////////////// sphere implementation ///////////////
  46
  47 Sphere::Sphere(const Vector3 &pos, scalar_t rad) : Surface(pos) {
  48         radius = rad;
  49 }
  50
  51 Sphere::~Sphere() {}
  52
  53 void Sphere::set_radius(scalar_t rad) {
  54         radius = rad;
  55 }
  56
  57 scalar_t Sphere::get_radius() const {
  58         return radius;
  59 }
  60
  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);
  65         Vector2 imap;
  66
  67         scalar_t phi = acos(dot_product(normal, pole));
  68         imap.y = phi / pi;
  69
  70         if(imap.y < xsmall_number || 1.0 - imap.y < xsmall_number) {
  71                 imap.x = 0.0;
  72                 return imap;
  73         }
  74
  75         scalar_t theta = acos(dot_product(equator, normal) / sin(phi)) / two_pi;
  76
  77         imap.x = (dot_product(cross_product(pole, equator), normal) < 0.0) ? theta : 1.0 - theta;
  78
  79         return imap;
  80 }
  81
  82 bool Sphere::check_intersection(const Ray &ray) const {
  83         return find_intersection(ray, 0);
  84 }
  85
  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);
  95
  96         // find the discriminant
  97         scalar_t d = SQ(b) - 4.0 * a * c;
  98         if(d < 0.0) return false;
  99
 100         // solve
 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);
 104
 105         if(t1 < error_margin && t2 < error_margin) return false;
 106
 107         if(isect) {
 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;
 112
 113                 isect->normal = (isect->pos - pos) / radius;
 114                 isect->pre_ior = ray.ior;
 115                 //isect->post_ior = mat.ior;
 116         }
 117
 118         return true;
 119 }
 120
 121
 122 Plane::Plane(const Vector3 &pos, const Vector3 &normal) : Surface(pos) {
 123         this->normal = normal;
 124 }
 125
 126 Plane::Plane(const Vector3 &p1, const Vector3 &p2, const Vector3 &p3) : Surface(p1)
 127 {
 128         normal = cross_product(p2 - p1, p3 - p1).normalized();
 129 }
 130
 131 Plane::~Plane() {}
 132
 133 void Plane::set_normal(const Vector3 &normal) {
 134         this->normal = normal;
 135 }
 136
 137 Vector3 Plane::get_normal() const {
 138         return normal;
 139 }
 140
 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);
 145 }
 146
 147 bool Plane::check_intersection(const Ray &ray) const {
 148         return find_intersection(ray, 0);
 149 }
 150
 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;
 154
 155         // TODO: this is only correct if pos is the projection of the origin on the plane
 156         scalar_t d = pos.length();
 157
 158         scalar_t t = -(dot_product(normal, ray.origin) + d) / normal_dot_dir;
 159
 160         if(t < error_margin) return false;
 161
 162         if(isect) {
 163                 isect->pos = ray.origin + ray.dir * t;
 164                 isect->normal = normal;
 165                 isect->t = t;
 166                 isect->pre_ior = ray.ior;
 167                 //isect->post_ior = mat.ior;
 168         }
 169         return true;
 170 }
 171
 172
 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);
 178
 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);
 183 }
 184
 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) {
 187         verts[0] = v0;
 188         verts[1] = v1;
 189         verts[2] = v2;
 190         verts[3] = v3;
 191         verts[4] = v4;
 192         verts[5] = v5;
 193         verts[6] = v6;
 194         verts[7] = v7;
 195 }
 196
 197 Box::Box(const Vector3 *array) {
 198         memcpy(verts, array, sizeof verts);
 199 }
 200
 201 Vector2 Box::inv_map(const Vector3 &pt) const {
 202         return Vector2();       // TODO: implement
 203 }
 204
 205 bool Box::check_intersection(const Ray &ray) const {
 206         return false;   // TODO: implement
 207 }
 208
 209 bool Box::find_intersection(const Ray &ray, SurfPoint *isect) const {
 210         return false;   // TODO: implement
 211 }
 212
 213
 214 /*
 215  * PointOverPlane (MG)
 216  * returns true if the point is in the positive side of the
 217  * plane
 218  */
 219 bool point_over_plane(const Plane &plane, const Vector3 &point)
 220 {
 221         if (dot_product(plane.get_position() - point, plane.get_normal()) < 0)
 222         {
 223                 return true;
 224         }
 225
 226         return false;
 227 }
 228
 229
 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;
 233 }
 234
 235
 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);
 240
 241         if(fabs(dot) < small_number) {
 242                 return false;
 243         }
 244
 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);
 249
 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)) {
 253                         return true;
 254                 }
 255         }
 256
 257         return false;
 258 }