cubemapper

annotate src/geom.cc @ 0:8fc9e1d3aad2

find changesets by author, revision, files, or words in the commit message
initial commit
author John Tsiombikas <nuclear@member.fsf.org>
date Thu, 27 Jul 2017 20:36:12 +0300
parents
children 2bfafdced01a
rev   line source
nuclear@0 1 #include <assert.h>
nuclear@0 2 #include <float.h>
nuclear@0 3 #include <algorithm>
nuclear@0 4 #include "geom.h"
nuclear@0 5
nuclear@0 6 GeomObject::~GeomObject()
nuclear@0 7 {
nuclear@0 8 }
nuclear@0 9
nuclear@0 10
nuclear@0 11 Sphere::Sphere()
nuclear@0 12 {
nuclear@0 13 radius = 1.0;
nuclear@0 14 }
nuclear@0 15
nuclear@0 16 Sphere::Sphere(const Vec3 &cent, float radius)
nuclear@0 17 : center(cent)
nuclear@0 18 {
nuclear@0 19 this->radius = radius;
nuclear@0 20 }
nuclear@0 21
nuclear@0 22 void Sphere::set_union(const GeomObject *obj1, const GeomObject *obj2)
nuclear@0 23 {
nuclear@0 24 const Sphere *sph1 = dynamic_cast<const Sphere*>(obj1);
nuclear@0 25 const Sphere *sph2 = dynamic_cast<const Sphere*>(obj2);
nuclear@0 26
nuclear@0 27 if(!sph1 || !sph2) {
nuclear@0 28 fprintf(stderr, "Sphere::set_union: arguments must be spheres");
nuclear@0 29 return;
nuclear@0 30 }
nuclear@0 31
nuclear@0 32 float dist = length(sph1->center - sph2->center);
nuclear@0 33 float surf_dist = dist - (sph1->radius + sph2->radius);
nuclear@0 34 float d1 = sph1->radius + surf_dist / 2.0;
nuclear@0 35 float d2 = sph2->radius + surf_dist / 2.0;
nuclear@0 36 float t = d1 / (d1 + d2);
nuclear@0 37
nuclear@0 38 if(t < 0.0) t = 0.0;
nuclear@0 39 if(t > 1.0) t = 1.0;
nuclear@0 40
nuclear@0 41 center = sph1->center * t + sph2->center * (1.0 - t);
nuclear@0 42 radius = std::max(dist * t + sph2->radius, dist * (1.0f - t) + sph1->radius);
nuclear@0 43 }
nuclear@0 44
nuclear@0 45 void Sphere::set_intersection(const GeomObject *obj1, const GeomObject *obj2)
nuclear@0 46 {
nuclear@0 47 fprintf(stderr, "Sphere::intersection undefined\n");
nuclear@0 48 }
nuclear@0 49
nuclear@0 50 bool Sphere::intersect(const Ray &ray, HitPoint *hit) const
nuclear@0 51 {
nuclear@0 52 float a = dot(ray.dir, ray.dir);
nuclear@0 53 float b = 2.0 * ray.dir.x * (ray.origin.x - center.x) +
nuclear@0 54 2.0 * ray.dir.y * (ray.origin.y - center.y) +
nuclear@0 55 2.0 * ray.dir.z * (ray.origin.z - center.z);
nuclear@0 56 float c = dot(ray.origin, ray.origin) + dot(center, center) -
nuclear@0 57 2.0 * dot(ray.origin, center) - radius * radius;
nuclear@0 58
nuclear@0 59 float discr = b * b - 4.0 * a * c;
nuclear@0 60 if(discr < 1e-4) {
nuclear@0 61 return false;
nuclear@0 62 }
nuclear@0 63
nuclear@0 64 float sqrt_discr = sqrt(discr);
nuclear@0 65 float t0 = (-b + sqrt_discr) / (2.0 * a);
nuclear@0 66 float t1 = (-b - sqrt_discr) / (2.0 * a);
nuclear@0 67
nuclear@0 68 if(t0 < 1e-4)
nuclear@0 69 t0 = t1;
nuclear@0 70 if(t1 < 1e-4)
nuclear@0 71 t1 = t0;
nuclear@0 72
nuclear@0 73 float t = t0 < t1 ? t0 : t1;
nuclear@0 74 if(t < 1e-4) {
nuclear@0 75 return false;
nuclear@0 76 }
nuclear@0 77
nuclear@0 78 // fill the HitPoint structure
nuclear@0 79 if(hit) {
nuclear@0 80 hit->obj = this;
nuclear@0 81 hit->dist = t;
nuclear@0 82 hit->pos = ray.origin + ray.dir * t;
nuclear@0 83 hit->normal = (hit->pos - center) / radius;
nuclear@0 84 }
nuclear@0 85 return true;
nuclear@0 86 }
nuclear@0 87
nuclear@0 88
nuclear@0 89 AABox::AABox()
nuclear@0 90 {
nuclear@0 91 }
nuclear@0 92
nuclear@0 93 AABox::AABox(const Vec3 &vmin, const Vec3 &vmax)
nuclear@0 94 : min(vmin), max(vmax)
nuclear@0 95 {
nuclear@0 96 }
nuclear@0 97
nuclear@0 98 void AABox::set_union(const GeomObject *obj1, const GeomObject *obj2)
nuclear@0 99 {
nuclear@0 100 const AABox *box1 = dynamic_cast<const AABox*>(obj1);
nuclear@0 101 const AABox *box2 = dynamic_cast<const AABox*>(obj2);
nuclear@0 102
nuclear@0 103 if(!box1 || !box2) {
nuclear@0 104 fprintf(stderr, "AABox::set_union: arguments must be AABoxes too\n");
nuclear@0 105 return;
nuclear@0 106 }
nuclear@0 107
nuclear@0 108 min.x = std::min(box1->min.x, box2->min.x);
nuclear@0 109 min.y = std::min(box1->min.y, box2->min.y);
nuclear@0 110 min.z = std::min(box1->min.z, box2->min.z);
nuclear@0 111
nuclear@0 112 max.x = std::max(box1->max.x, box2->max.x);
nuclear@0 113 max.y = std::max(box1->max.y, box2->max.y);
nuclear@0 114 max.z = std::max(box1->max.z, box2->max.z);
nuclear@0 115 }
nuclear@0 116
nuclear@0 117 void AABox::set_intersection(const GeomObject *obj1, const GeomObject *obj2)
nuclear@0 118 {
nuclear@0 119 const AABox *box1 = dynamic_cast<const AABox*>(obj1);
nuclear@0 120 const AABox *box2 = dynamic_cast<const AABox*>(obj2);
nuclear@0 121
nuclear@0 122 if(!box1 || !box2) {
nuclear@0 123 fprintf(stderr, "AABox::set_intersection: arguments must be AABoxes too\n");
nuclear@0 124 return;
nuclear@0 125 }
nuclear@0 126
nuclear@0 127 for(int i=0; i<3; i++) {
nuclear@0 128 min[i] = std::max(box1->min[i], box2->min[i]);
nuclear@0 129 max[i] = std::min(box1->max[i], box2->max[i]);
nuclear@0 130
nuclear@0 131 if(max[i] < min[i]) {
nuclear@0 132 max[i] = min[i];
nuclear@0 133 }
nuclear@0 134 }
nuclear@0 135 }
nuclear@0 136
nuclear@0 137 bool AABox::intersect(const Ray &ray, HitPoint *hit) const
nuclear@0 138 {
nuclear@0 139 Vec3 param[2] = {min, max};
nuclear@0 140 Vec3 inv_dir(1.0 / ray.dir.x, 1.0 / ray.dir.y, 1.0 / ray.dir.z);
nuclear@0 141 int sign[3] = {inv_dir.x < 0, inv_dir.y < 0, inv_dir.z < 0};
nuclear@0 142
nuclear@0 143 float tmin = (param[sign[0]].x - ray.origin.x) * inv_dir.x;
nuclear@0 144 float tmax = (param[1 - sign[0]].x - ray.origin.x) * inv_dir.x;
nuclear@0 145 float tymin = (param[sign[1]].y - ray.origin.y) * inv_dir.y;
nuclear@0 146 float tymax = (param[1 - sign[1]].y - ray.origin.y) * inv_dir.y;
nuclear@0 147
nuclear@0 148 if(tmin > tymax || tymin > tmax) {
nuclear@0 149 return false;
nuclear@0 150 }
nuclear@0 151 if(tymin > tmin) {
nuclear@0 152 tmin = tymin;
nuclear@0 153 }
nuclear@0 154 if(tymax < tmax) {
nuclear@0 155 tmax = tymax;
nuclear@0 156 }
nuclear@0 157
nuclear@0 158 float tzmin = (param[sign[2]].z - ray.origin.z) * inv_dir.z;
nuclear@0 159 float tzmax = (param[1 - sign[2]].z - ray.origin.z) * inv_dir.z;
nuclear@0 160
nuclear@0 161 if(tmin > tzmax || tzmin > tmax) {
nuclear@0 162 return false;
nuclear@0 163 }
nuclear@0 164 if(tzmin > tmin) {
nuclear@0 165 tmin = tzmin;
nuclear@0 166 }
nuclear@0 167 if(tzmax < tmax) {
nuclear@0 168 tmax = tzmax;
nuclear@0 169 }
nuclear@0 170
nuclear@0 171 float t = tmin < 1e-4 ? tmax : tmin;
nuclear@0 172 if(t >= 1e-4) {
nuclear@0 173
nuclear@0 174 if(hit) {
nuclear@0 175 hit->obj = this;
nuclear@0 176 hit->dist = t;
nuclear@0 177 hit->pos = ray.origin + ray.dir * t;
nuclear@0 178
nuclear@0 179 float min_dist = FLT_MAX;
nuclear@0 180 Vec3 offs = min + (max - min) / 2.0;
nuclear@0 181 Vec3 local_hit = hit->pos - offs;
nuclear@0 182
nuclear@0 183 static const Vec3 axis[] = {
nuclear@0 184 Vec3(1, 0, 0), Vec3(0, 1, 0), Vec3(0, 0, 1)
nuclear@0 185 };
nuclear@0 186 //int tcidx[][2] = {{2, 1}, {0, 2}, {0, 1}};
nuclear@0 187
nuclear@0 188 for(int i=0; i<3; i++) {
nuclear@0 189 float dist = fabs((max[i] - offs[i]) - fabs(local_hit[i]));
nuclear@0 190 if(dist < min_dist) {
nuclear@0 191 min_dist = dist;
nuclear@0 192 hit->normal = axis[i] * (local_hit[i] < 0.0 ? 1.0 : -1.0);
nuclear@0 193 //hit->texcoord = Vec2(hit->pos[tcidx[i][0]], hit->pos[tcidx[i][1]]);
nuclear@0 194 }
nuclear@0 195 }
nuclear@0 196 }
nuclear@0 197 return true;
nuclear@0 198 }
nuclear@0 199 return false;
nuclear@0 200
nuclear@0 201 }
nuclear@0 202
nuclear@0 203 Plane::Plane()
nuclear@0 204 : normal(0.0, 1.0, 0.0)
nuclear@0 205 {
nuclear@0 206 }
nuclear@0 207
nuclear@0 208 Plane::Plane(const Vec3 &p, const Vec3 &norm)
nuclear@0 209 : pt(p)
nuclear@0 210 {
nuclear@0 211 normal = normalize(norm);
nuclear@0 212 }
nuclear@0 213
nuclear@0 214 Plane::Plane(const Vec3 &p1, const Vec3 &p2, const Vec3 &p3)
nuclear@0 215 : pt(p1)
nuclear@0 216 {
nuclear@0 217 normal = normalize(cross(p2 - p1, p3 - p1));
nuclear@0 218 }
nuclear@0 219
nuclear@0 220 Plane::Plane(const Vec3 &normal, float dist)
nuclear@0 221 {
nuclear@0 222 this->normal = normalize(normal);
nuclear@0 223 pt = this->normal * dist;
nuclear@0 224 }
nuclear@0 225
nuclear@0 226 void Plane::set_union(const GeomObject *obj1, const GeomObject *obj2)
nuclear@0 227 {
nuclear@0 228 fprintf(stderr, "Plane::set_union undefined\n");
nuclear@0 229 }
nuclear@0 230
nuclear@0 231 void Plane::set_intersection(const GeomObject *obj1, const GeomObject *obj2)
nuclear@0 232 {
nuclear@0 233 fprintf(stderr, "Plane::set_intersection undefined\n");
nuclear@0 234 }
nuclear@0 235
nuclear@0 236 bool Plane::intersect(const Ray &ray, HitPoint *hit) const
nuclear@0 237 {
nuclear@0 238 float ndotdir = dot(normal, ray.dir);
nuclear@0 239 if(fabs(ndotdir) < 1e-4) {
nuclear@0 240 return false;
nuclear@0 241 }
nuclear@0 242
nuclear@0 243 if(hit) {
nuclear@0 244 Vec3 ptdir = pt - ray.origin;
nuclear@0 245 float t = dot(normal, ptdir) / ndotdir;
nuclear@0 246
nuclear@0 247 hit->pos = ray.origin + ray.dir * t;
nuclear@0 248 hit->normal = normal;
nuclear@0 249 hit->obj = this;
nuclear@0 250 }
nuclear@0 251 return true;
nuclear@0 252 }
nuclear@0 253
nuclear@0 254 float sphere_distance(const Vec3 &cent, float rad, const Vec3 &pt)
nuclear@0 255 {
nuclear@0 256 return length(pt - cent) - rad;
nuclear@0 257 }
nuclear@0 258
nuclear@0 259 // TODO version which takes both radii into account
nuclear@0 260 float capsule_distance(const Vec3 &a, float ra, const Vec3 &b, float rb, const Vec3 &pt)
nuclear@0 261 {
nuclear@0 262 Vec3 ab_dir = b - a;
nuclear@0 263 float ab_len_sq = length_sq(ab_dir);
nuclear@0 264
nuclear@0 265 if(fabs(ab_len_sq) < 1e-5) {
nuclear@0 266 // if a == b, the capsule is a sphere with radius the maximum of the capsule radii
nuclear@0 267 return sphere_distance(a, std::max(ra, rb), pt);
nuclear@0 268 }
nuclear@0 269 float ab_len = sqrt(ab_len_sq);
nuclear@0 270
nuclear@0 271 Vec3 ap_dir = pt - a;
nuclear@0 272
nuclear@0 273 float t = dot(ap_dir, ab_dir / ab_len) / ab_len;
nuclear@0 274 if(t < 0.0) {
nuclear@0 275 return sphere_distance(a, ra, pt);
nuclear@0 276 }
nuclear@0 277 if(t >= 1.0) {
nuclear@0 278 return sphere_distance(b, rb, pt);
nuclear@0 279 }
nuclear@0 280
nuclear@0 281 Vec3 pproj = a + ab_dir * t;
nuclear@0 282 return length(pproj - pt) - ra;
nuclear@0 283 }
nuclear@0 284
nuclear@0 285 #if 0
nuclear@0 286 float capsule_distance(const Vec3 &a, float ra, const Vec3 &b, float rb, const Vec3 &pt)
nuclear@0 287 {
nuclear@0 288 Vec3 ab_dir = b - a;
nuclear@0 289
nuclear@0 290 if(fabs(length_sq(ab_dir)) < 1e-5) {
nuclear@0 291 // if a == b, the capsule is a sphere with radius the maximum of the capsule radii
nuclear@0 292 return sphere_distance(a, std::max(ra, rb), pt);
nuclear@0 293 }
nuclear@0 294 float ab_len = length(ab_dir);
nuclear@0 295
nuclear@0 296 Vec3 ap_dir = pt - a;
nuclear@0 297 Vec3 rotaxis = normalize(cross(ab_dir, ap_dir));
nuclear@0 298
nuclear@0 299 Mat4 rmat;
nuclear@0 300 rmat.set_rotation(rotaxis, M_PI / 2.0);
nuclear@0 301 Vec3 right = rmat * ab_dir / ab_len;
nuclear@0 302
nuclear@0 303 // XXX I think this check is redundant, always false, due to the cross product order
nuclear@0 304 //assert(dot(right, ab_dir) >= 0.0);
nuclear@0 305 if(dot(right, ab_dir) < 0.0) {
nuclear@0 306 right = -right;
nuclear@0 307 }
nuclear@0 308 Vec3 aa = a + right * ra;
nuclear@0 309 Vec3 bb = b + right * rb;
nuclear@0 310
nuclear@0 311 // project pt to the line segment bb-aa, see if the projection lies within the interval [0, 1)
nuclear@0 312 Vec3 aabb_dir = bb - aa;
nuclear@0 313 float aabb_len = length(aabb_dir);
nuclear@0 314 Vec3 aap_dir = pt - aa;
nuclear@0 315
nuclear@0 316 float t = dot(aap_dir, aabb_dir / aabb_len) / aabb_len;
nuclear@0 317 if(t < 0.0) {
nuclear@0 318 return sphere_distance(a, ra, pt);
nuclear@0 319 }
nuclear@0 320 if(t >= 1.0) {
nuclear@0 321 return sphere_distance(b, rb, pt);
nuclear@0 322 }
nuclear@0 323
nuclear@0 324 Vec3 ppt = aa + aabb_dir * t;
nuclear@0 325 Vec3 norm = ppt - pt;
nuclear@0 326 float dist = length(norm);
nuclear@0 327
nuclear@0 328 if(dot(norm, right) < 0.0) {
nuclear@0 329 // inside the cone
nuclear@0 330 dist = -dist;
nuclear@0 331 }
nuclear@0 332 return dist;
nuclear@0 333 }
nuclear@0 334 #endif