1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/src/geom.cc	Sat Feb 14 10:08:00 2015 +0200
     1.3 @@ -0,0 +1,240 @@
     1.4 +#include <float.h>
     1.5 +#include <algorithm>
     1.6 +#include "geom.h"
     1.7 +
     1.8 +static bool col_sphere_sphere(const Sphere *a, const Sphere *b, HitPoint *hit);
     1.9 +static bool col_plane_plane(const Plane *a, const Plane *b, HitPoint *hit);
    1.10 +static bool col_sphere_plane(const Sphere *sph, const Plane *plane, HitPoint *hit);
    1.11 +
    1.12 +GeomShape::GeomShape()
    1.13 +{
    1.14 +	type = GEOM_UNKNOWN;
    1.15 +}
    1.16 +
    1.17 +GeomShape::~GeomShape()
    1.18 +{
    1.19 +}
    1.20 +
    1.21 +GeomShape::Type GeomShape::get_type() const
    1.22 +{
    1.23 +	return type;
    1.24 +}
    1.25 +
    1.26 +// ---- Sphere class ----
    1.27 +
    1.28 +Sphere::Sphere()
    1.29 +{
    1.30 +	type = GEOM_SPHERE;
    1.31 +	radius = 1.0f;
    1.32 +}
    1.33 +
    1.34 +Sphere::Sphere(const Vector3 &c, float rad)
    1.35 +	: center(c)
    1.36 +{
    1.37 +	type = GEOM_SPHERE;
    1.38 +	radius = rad;
    1.39 +}
    1.40 +
    1.41 +bool Sphere::intersect(const Ray &ray, HitPoint *hit) const
    1.42 +{
    1.43 +	float a = dot_product(ray.dir, ray.dir);
    1.44 +	float b = 2.0 * ray.dir.x * (ray.origin.x - center.x) +
    1.45 +		2.0 * ray.dir.y * (ray.origin.y - center.y) +
    1.46 +		2.0 * ray.dir.z * (ray.origin.z - center.z);
    1.47 +	float c = dot_product(ray.origin, ray.origin) + dot_product(center, center) -
    1.48 +		2.0 * dot_product(ray.origin, center) - radius * radius;
    1.49 +
    1.50 +	float discr = b * b - 4.0 * a * c;
    1.51 +	if(discr < 1e-4) {
    1.52 +		return false;
    1.53 +	}
    1.54 +
    1.55 +	float sqrt_discr = sqrt(discr);
    1.56 +	float t0 = (-b + sqrt_discr) / (2.0 * a);
    1.57 +	float t1 = (-b - sqrt_discr) / (2.0 * a);
    1.58 +
    1.59 +	if(t0 < 1e-4)
    1.60 +		t0 = t1;
    1.61 +	if(t1 < 1e-4)
    1.62 +		t1 = t0;
    1.63 +
    1.64 +	float t = t0 < t1 ? t0 : t1;
    1.65 +	if(t < 1e-4) {
    1.66 +		return false;
    1.67 +	}
    1.68 +
    1.69 +	// fill the HitPoint structure
    1.70 +	if(hit) {
    1.71 +		hit->shape = this;
    1.72 +		hit->data = 0;
    1.73 +		hit->dist = t;
    1.74 +		hit->pos = ray.origin + ray.dir * t;
    1.75 +		hit->normal = (hit->pos - center) / radius;
    1.76 +	}
    1.77 +	return true;
    1.78 +}
    1.79 +
    1.80 +bool Sphere::collide(const GeomShape *geom, HitPoint *hit) const
    1.81 +{
    1.82 +	switch(geom->get_type()) {
    1.83 +	case GEOM_SPHERE:
    1.84 +		return col_sphere_sphere(this, (const Sphere*)geom, hit);
    1.85 +	case GEOM_PLANE:
    1.86 +		return col_sphere_plane(this, (const Plane*)geom, hit);
    1.87 +	default:
    1.88 +		break;
    1.89 +	}
    1.90 +	return false;
    1.91 +}
    1.92 +
    1.93 +float Sphere::distance(const Vector3 &pt) const
    1.94 +{
    1.95 +	return sqrt(distance_sq(pt));
    1.96 +}
    1.97 +
    1.98 +float Sphere::distance_sq(const Vector3 &pt) const
    1.99 +{
   1.100 +	return (pt - center).length_sq();
   1.101 +}
   1.102 +
   1.103 +// ---- Plane class ----
   1.104 +
   1.105 +Plane::Plane()
   1.106 +	: normal(0, 1, 0)
   1.107 +{
   1.108 +	type = GEOM_PLANE;
   1.109 +	dist = 0.0f;
   1.110 +}
   1.111 +
   1.112 +Plane::Plane(const Vector3 &n, float d)
   1.113 +	: normal(n)
   1.114 +{
   1.115 +	type = GEOM_PLANE;
   1.116 +	dist = d;
   1.117 +
   1.118 +	normal.normalize();
   1.119 +}
   1.120 +
   1.121 +Plane::Plane(float a, float b, float c, float d)
   1.122 +	: normal(a, b, c)
   1.123 +{
   1.124 +	type = GEOM_PLANE;
   1.125 +	dist = d;
   1.126 +
   1.127 +	normal.normalize();
   1.128 +}
   1.129 +
   1.130 +Plane::Plane(const Vector3 &pos, const Vector3 &norm)
   1.131 +	: normal(norm)
   1.132 +{
   1.133 +	type = GEOM_PLANE;
   1.134 +
   1.135 +	dist = 0.0f;
   1.136 +	dist = distance(pos);
   1.137 +}
   1.138 +
   1.139 +Plane::Plane(const Vector3 &p1, const Vector3 &p2, const Vector3 &p3)
   1.140 +{
   1.141 +	type = GEOM_PLANE;
   1.142 +
   1.143 +	normal = cross_product(p2 - p1, p3 - p1).normalized();
   1.144 +	dist = 0.0f;
   1.145 +	dist = distance(p1);
   1.146 +}
   1.147 +
   1.148 +
   1.149 +bool Plane::intersect(const Ray &ray, HitPoint *hit) const
   1.150 +{
   1.151 +	Vector3 pt = Vector3(0, 0, 0) + normal * dist;
   1.152 +
   1.153 +	float ndotdir = dot_product(normal, ray.dir);
   1.154 +	if(fabs(ndotdir) < 1e-4) {
   1.155 +		return false;
   1.156 +	}
   1.157 +
   1.158 +	if(hit) {
   1.159 +		Vector3 ptdir = pt - ray.origin;
   1.160 +		float t = dot_product(normal, ptdir) / ndotdir;
   1.161 +
   1.162 +		hit->pos = ray.origin + ray.dir * t;
   1.163 +		hit->normal = normal;
   1.164 +		hit->shape = this;
   1.165 +	}
   1.166 +	return true;
   1.167 +}
   1.168 +
   1.169 +bool Plane::collide(const GeomShape *geom, HitPoint *hit) const
   1.170 +{
   1.171 +	switch(geom->get_type()) {
   1.172 +	case GEOM_SPHERE:
   1.173 +		{
   1.174 +			bool res = col_sphere_plane((const Sphere*)geom, this, hit);
   1.175 +			if(hit) {
   1.176 +				hit->normal = -hit->normal;
   1.177 +			}
   1.178 +			return res;
   1.179 +		}
   1.180 +
   1.181 +	case GEOM_PLANE:
   1.182 +		return col_plane_plane(this, (const Plane*)geom, hit);
   1.183 +
   1.184 +	default:
   1.185 +		break;
   1.186 +	}
   1.187 +	return false;
   1.188 +}
   1.189 +
   1.190 +float Plane::distance(const Vector3 &v) const
   1.191 +{
   1.192 +	Vector3 pt = normal * dist;
   1.193 +	return dot_product(v - pt, normal);
   1.194 +}
   1.195 +
   1.196 +float Plane::distance_sq(const Vector3 &v) const
   1.197 +{
   1.198 +	float d = distance(v);
   1.199 +	return d * d;
   1.200 +}
   1.201 +
   1.202 +
   1.203 +// collision detection functions
   1.204 +
   1.205 +static bool col_sphere_sphere(const Sphere *a, const Sphere *b, HitPoint *hit)
   1.206 +{
   1.207 +	float sum_rad = a->radius + b->radius;
   1.208 +	Vector3 dir = b->center - a->center;
   1.209 +	float dist_sq = dir.length_sq();
   1.210 +	bool res = dist_sq <= sum_rad * sum_rad;
   1.211 +
   1.212 +	if(res && hit) {
   1.213 +		hit->pos = a->center + dir * 0.5;
   1.214 +		hit->shape = b;
   1.215 +		hit->normal = -dir;
   1.216 +	}
   1.217 +	return res;
   1.218 +}
   1.219 +
   1.220 +static bool col_plane_plane(const Plane *a, const Plane *b, HitPoint *hit)
   1.221 +{
   1.222 +	fprintf(stderr, "%s: not implemented\n", __FUNCTION__);
   1.223 +	return false;
   1.224 +}
   1.225 +
   1.226 +static bool col_sphere_plane(const Sphere *sph, const Plane *plane, HitPoint *hit)
   1.227 +{
   1.228 +	float dist = plane->distance(sph->center);
   1.229 +	bool res = fabs(dist) <= sph->radius;
   1.230 +
   1.231 +	if(res && hit) {
   1.232 +		Vector3 planept = plane->normal * plane->dist;
   1.233 +		Vector3 sphdir = sph->center - planept;
   1.234 +		if(dot_product(sphdir, plane->normal) >= 0.0f) {
   1.235 +			hit->normal = plane->normal;
   1.236 +		} else {
   1.237 +			hit->normal = -plane->normal;
   1.238 +		}
   1.239 +		hit->pos = sph->center - plane->normal * fabs(dist);
   1.240 +		hit->shape = plane;
   1.241 +	}
   1.242 +	return res;
   1.243 +}