1 /** 2 This module implements some abstract geometric shapes: 3 $(UL 4 $(LI Line segments.) 5 $(LI Triangle.) 6 $(LI Circles/spheres.) 7 $(LI Rays) 8 $(LI Planes) 9 $(LI Frustum) 10 ) 11 */ 12 module gfm.math.shapes; 13 14 import std.math, 15 std.traits; 16 17 import gfm.math.vector, 18 gfm.math.box; 19 20 /// A Segment is 2 points. 21 /// When considered like a vector, it represents the arrow from a to b. 22 struct Segment(T, int N) 23 { 24 public 25 { 26 alias Vector!(T, N) point_t; 27 point_t a, b; 28 } 29 } 30 31 alias Segment!(float, 2) seg2f; /// 2D float segment. 32 alias Segment!(float, 3) seg3f; /// 3D float segment. 33 alias Segment!(double, 2) seg2d; /// 2D double segment. 34 alias Segment!(double, 3) seg3d; /// 3D double segment. 35 alias Segment!(int, 2) seg2i; /// 2D integer segment. 36 alias Segment!(int, 3) seg3i; /// 3D integer segment. 37 38 /// A Triangle is 3 points. 39 struct Triangle(T, int N) 40 { 41 public 42 { 43 alias Vector!(T, N) point_t; 44 point_t a, b, c; 45 46 static if (N == 2) 47 { 48 /// Returns: Area of a 2D triangle. 49 @nogc T area() pure const nothrow 50 { 51 return abs(signedArea()); 52 } 53 54 /// Returns: Signed area of a 2D triangle. 55 @nogc T signedArea() pure const nothrow 56 { 57 return ((b.x * a.y - a.x * b.y) 58 + (c.x * b.y - b.x * c.y) 59 + (a.x * c.y - c.x * a.y)) / 2; 60 } 61 } 62 63 static if (N == 3 && isFloatingPoint!T) 64 { 65 /// Returns: Triangle normal. 66 @nogc Vector!(T, 3) computeNormal() pure const nothrow 67 { 68 return cross(b - a, c - a).normalized(); 69 } 70 } 71 } 72 } 73 74 alias Triangle!(float, 2) triangle2f; /// 2D float triangle. 75 alias Triangle!(float, 3) triangle3f; /// 3D float triangle. 76 alias Triangle!(double, 2) triangle2d; /// 2D double triangle. 77 alias Triangle!(double, 3) triangle3d; /// 3D double triangle. 78 79 /// A Sphere is a point + a radius. 80 struct Sphere(T, int N) 81 { 82 public nothrow 83 { 84 alias Vector!(T, N) point_t; 85 point_t center; 86 T radius; 87 88 /// Creates a sphere from a point and a radius. 89 @nogc this(in point_t center_, T radius_) pure nothrow 90 { 91 center = center_; 92 radius = radius_; 93 } 94 95 /// Sphere contains point test. 96 /// Returns: true if the point is inside the sphere. 97 @nogc bool contains(in Sphere s) pure const nothrow 98 { 99 if (s.radius > radius) 100 return false; 101 102 T innerRadius = radius - s.radius; 103 return squaredDistanceTo(s.center) < innerRadius * innerRadius; 104 } 105 106 /// Sphere vs point Euclidean distance squared. 107 @nogc T squaredDistanceTo(point_t p) pure const nothrow 108 { 109 return center.squaredDistanceTo(p); 110 } 111 112 /// Sphere vs sphere intersection. 113 /// Returns: true if the spheres intersect. 114 @nogc bool intersects(Sphere s) pure const nothrow 115 { 116 T outerRadius = radius + s.radius; 117 return squaredDistanceTo(s.center) < outerRadius * outerRadius; 118 } 119 120 static if (isFloatingPoint!T) 121 { 122 /// Sphere vs point Euclidean distance. 123 @nogc T distanceTo(point_t p) pure const nothrow 124 { 125 return center.distanceTo(p); 126 } 127 128 static if(N == 2) 129 { 130 /// Returns: Circle area. 131 @nogc T area() pure const nothrow 132 { 133 return PI * (radius * radius); 134 } 135 } 136 } 137 } 138 } 139 140 alias Sphere!(float, 2) sphere2f; /// 2D float sphere (ie. a circle). 141 alias Sphere!(float, 3) sphere3f; /// 3D float sphere. 142 alias Sphere!(double, 2) sphere2d; /// 2D double sphere (ie. a circle). 143 alias Sphere!(double, 3) sphere3d; /// 3D double sphere (ie. a circle). 144 145 146 /// A Ray ir a point + a direction. 147 struct Ray(T, int N) 148 { 149 nothrow: 150 public 151 { 152 alias Vector!(T, N) point_t; 153 point_t orig; 154 point_t dir; 155 156 /// Returns: A point further along the ray direction. 157 @nogc point_t progress(T t) pure const nothrow 158 { 159 return orig + dir * t; 160 } 161 162 static if (N == 3 && isFloatingPoint!T) 163 { 164 /// Ray vs triangle intersection. 165 /// See_also: "Fast, Minimum Storage Ray/Triangle intersection", Mommer & Trumbore (1997) 166 /// Returns: Barycentric coordinates, the intersection point is at $(D (1 - u - v) * A + u * B + v * C). 167 @nogc bool intersect(Triangle!(T, 3) triangle, out T t, out T u, out T v) pure const nothrow 168 { 169 point_t edge1 = triangle.b - triangle.a; 170 point_t edge2 = triangle.c - triangle.a; 171 point_t pvec = cross(dir, edge2); 172 T det = dot(edge1, pvec); 173 if (abs(det) < T.epsilon) 174 return false; // no intersection 175 T invDet = 1 / det; 176 177 // calculate distance from triangle.a to ray origin 178 point_t tvec = orig - triangle.a; 179 180 // calculate U parameter and test bounds 181 u = dot(tvec, pvec) * invDet; 182 if (u < 0 || u > 1) 183 return false; 184 185 // prepare to test V parameter 186 point_t qvec = cross(tvec, edge1); 187 188 // calculate V parameter and test bounds 189 v = dot(dir, qvec) * invDet; 190 if (v < 0.0 || u + v > 1.0) 191 return false; 192 193 // calculate t, ray intersects triangle 194 t = dot(edge2, qvec) * invDet; 195 return true; 196 } 197 } 198 } 199 } 200 201 alias Ray!(float, 2) ray2f; /// 2D float ray. 202 alias Ray!(float, 3) ray3f; /// 3D float ray. 203 alias Ray!(double, 2) ray2d; /// 2D double ray. 204 alias Ray!(double, 3) ray3d; /// 3D double ray. 205 206 207 /// 3D plane. 208 /// See_also: Flipcode article by Nate Miller $(WEB www.flipcode.com/archives/Plane_Class.shtml). 209 struct Plane(T) if (isFloatingPoint!T) 210 { 211 public 212 { 213 vec3!T n; /// Normal (always stored normalized). 214 T d; 215 216 /// Create from four coordinates. 217 @nogc this(vec4!T abcd) pure nothrow 218 { 219 n = vec3!T(abcd.x, abcd.y, abcd.z).normalized(); 220 d = abcd.z; 221 } 222 223 /// Create from a point and a normal. 224 @nogc this(vec3!T origin, vec3!T normal) pure nothrow 225 { 226 n = normal.normalized(); 227 d = -dot(origin, n); 228 } 229 230 /// Create from 3 non-aligned points. 231 @nogc this(vec3!T A, vec3!T B, vec3!T C) pure nothrow 232 { 233 this(C, cross(B - A, C - A)); 234 } 235 236 /// Assign a plane with another plane. 237 @nogc ref Plane opAssign(Plane other) pure nothrow 238 { 239 n = other.n; 240 d = other.d; 241 return this; 242 } 243 244 /// Returns: signed distance between a point and the plane. 245 @nogc T signedDistanceTo(vec3!T point) pure const nothrow 246 { 247 return dot(n, point) + d; 248 } 249 250 /// Returns: absolute distance between a point and the plane. 251 @nogc T distanceTo(vec3!T point) pure const nothrow 252 { 253 return abs(signedDistanceTo(point)); 254 } 255 256 /// Returns: true if the point is in front of the plane. 257 @nogc bool isFront(vec3!T point) pure const nothrow 258 { 259 return signedDistanceTo(point) >= 0; 260 } 261 262 /// Returns: true if the point is in the back of the plane. 263 @nogc bool isBack(vec3!T point) pure const nothrow 264 { 265 return signedDistanceTo(point) < 0; 266 } 267 268 /// Returns: true if the point is on the plane, with a given epsilon. 269 @nogc bool isOn(vec3!T point, T epsilon) pure const nothrow 270 { 271 T sd = signedDistanceTo(point); 272 return (-epsilon < sd) && (sd < epsilon); 273 } 274 } 275 } 276 277 alias Plane!float planef; /// 3D float plane. 278 alias Plane!double planed; /// 3D double plane. 279 280 unittest 281 { 282 auto p = planed(vec4d(1.0, 2.0, 3.0, 4.0)); 283 auto p2 = planed(vec3d(1.0, 0.0, 0.0), 284 vec3d(0.0, 1.0, 0.0), 285 vec3d(0.0, 0.0, 1.0)); 286 287 assert(p2.isOn(vec3d(1.0, 0.0, 0.0), 1e-7)); 288 assert(p2.isFront(vec3d(1.0, 1.0, 1.0))); 289 assert(p2.isBack(vec3d(0.0, 0.0, 0.0))); 290 } 291 292 293 294 /// 3D frustum. 295 /// See_also: Flipcode article by Dion Picco $(WEB www.flipcode.com/archives/Frustum_Culling.shtml). 296 /// Bugs: verify proper signedness of half-spaces 297 struct Frustum(T) if (isFloatingPoint!T) 298 { 299 public 300 { 301 enum int LEFT = 0, 302 RIGHT = 1, 303 TOP = 2, 304 BOTTOM = 3, 305 NEAR = 4, 306 FAR = 5; 307 308 Plane!T[6] planes; 309 310 /// Create a frustum from 6 planes. 311 @nogc this(Plane!T left, Plane!T right, Plane!T top, Plane!T bottom, Plane!T near, Plane!T far) pure nothrow 312 { 313 planes[LEFT] = left; 314 planes[RIGHT] = right; 315 planes[TOP] = top; 316 planes[BOTTOM] = bottom; 317 planes[NEAR] = near; 318 planes[FAR] = far; 319 } 320 321 enum : int 322 { 323 OUTSIDE, /// object is outside the frustum 324 INTERSECT, /// object intersects with the frustum 325 INSIDE /// object is inside the frustum 326 } 327 328 /// Point vs frustum intersection. 329 @nogc bool contains(vec3!T point) pure const nothrow 330 { 331 for(int i = 0; i < 6; ++i) 332 { 333 T distance = planes[i].signedDistanceTo(point); 334 335 if(distance < 0) 336 return false; 337 } 338 return true; 339 } 340 341 /// Sphere vs frustum intersection. 342 /// Returns: Frustum.OUTSIDE, Frustum.INTERSECT or Frustum.INSIDE. 343 @nogc int contains(Sphere!(T, 3) sphere) pure const nothrow 344 { 345 // calculate our distances to each of the planes 346 for(int i = 0; i < 6; ++i) 347 { 348 // find the distance to this plane 349 T distance = planes[i].signedDistanceTo(sphere.center); 350 351 if(distance < -sphere.radius) 352 return OUTSIDE; 353 354 else if (distance < sphere.radius) 355 return INTERSECT; 356 } 357 358 // otherwise we are fully in view 359 return INSIDE; 360 } 361 362 /// AABB vs frustum intersection. 363 /// Returns: Frustum.OUTSIDE, Frustum.INTERSECT or Frustum.INSIDE. 364 @nogc int contains(box3!T box) pure const nothrow 365 { 366 vec3!T[8] corners; 367 int totalIn = 0; 368 369 for (int i = 0; i < 2; ++i) 370 for (int j = 0; j < 2; ++j) 371 for (int k = 0; k < 2; ++k) 372 { 373 auto x = i == 0 ? box.min.x : box.max.x; 374 auto y = i == 0 ? box.min.y : box.max.y; 375 auto z = i == 0 ? box.min.z : box.max.z; 376 corners[i*4 + j*2 + k] = vec3!T(x, y, z); 377 } 378 379 // test all 8 corners against the 6 sides 380 // if all points are behind 1 specific plane, we are out 381 // if we are in with all points, then we are fully in 382 for(int p = 0; p < 6; ++p) 383 { 384 int inCount = 8; 385 int ptIn = 1; 386 387 for(int i = 0; i < 8; ++i) 388 { 389 // test this point against the planes 390 if (planes[p].isBack(corners[i])) 391 { 392 ptIn = 0; 393 --inCount; 394 } 395 } 396 397 // were all the points outside of plane p? 398 if (inCount == 0) 399 return OUTSIDE; 400 401 // check if they were all on the right side of the plane 402 totalIn += ptIn; 403 } 404 405 // so if totalIn is 6, then all are inside the view 406 if(totalIn == 6) 407 return INSIDE; 408 409 // we must be partly in then otherwise 410 return INTERSECT; 411 } 412 413 } 414 } 415 416 unittest 417 { 418 seg2f se; 419 triangle3f tr; 420 Frustum!double frust; 421 planed pl; 422 } 423 424 /// True if `T` is a kind of Segment 425 enum isSegment(T) = is(T : Segment!U, U...); 426 427 /// True if `T` is a kind of Triangle 428 enum isTriangle(T) = is(T : Triangle!U, U...); 429 430 /// True if `T` is a kind of Sphere 431 enum isSphere(T) = is(T : Sphere!U, U...); 432 433 /// True if `T` is a kind of Ray 434 enum isRay(T) = is(T : Ray!U, U...); 435 436 /// True if `T` is a kind of Plane 437 enum isPlane(T) = is(T : Plane!U, U); 438 439 /// True if `T` is a kind of Frustum 440 enum isFrustum(T) = is(T : Frustum!U, U); 441 442 /// True if `T` is a kind of 2 dimensional Segment 443 enum isSegment2D(T) = is(T : Segment!(U, 2), U); 444 445 /// True if `T` is a kind of 2 dimensional Triangle 446 enum isTriangle2D(T) = is(T : Triangle!(U, 2), U); 447 448 /// True if `T` is a kind of 2 dimensional Sphere 449 enum isSphere2D(T) = is(T : Sphere!(U, 2), U); 450 451 /// True if `T` is a kind of 2 dimensional Ray 452 enum isRay2D(T) = is(T : Ray!(U, 2), U); 453 454 /// True if `T` is a kind of 3 dimensional Segment 455 enum isSegment3D(T) = is(T : Segment!(U, 3), U); 456 457 /// True if `T` is a kind of 3 dimensional Triangle 458 enum isTriangle3D(T) = is(T : Triangle!(U, 3), U); 459 460 /// True if `T` is a kind of 3 dimensional Sphere 461 enum isSphere3D(T) = is(T : Sphere!(U, 3), U); 462 463 /// True if `T` is a kind of 3 dimensional Ray 464 enum isRay3D(T) = is(T : Ray!(U, 3), U); 465 466 unittest 467 { 468 enum ShapeType 469 { 470 segment, 471 triangle, 472 sphere, 473 ray, 474 plane, 475 frustum 476 } 477 478 void test(T, ShapeType type, int dim)() 479 { 480 with (ShapeType) 481 { 482 static assert(isSegment!T == (type == segment )); 483 static assert(isTriangle!T == (type == triangle)); 484 static assert(isSphere!T == (type == sphere )); 485 static assert(isRay!T == (type == ray )); 486 static assert(isPlane!T == (type == plane )); 487 static assert(isFrustum!T == (type == frustum )); 488 489 static assert(isSegment2D!T == (type == segment && dim == 2)); 490 static assert(isTriangle2D!T == (type == triangle && dim == 2)); 491 static assert(isSphere2D!T == (type == sphere && dim == 2)); 492 static assert(isRay2D!T == (type == ray && dim == 2)); 493 494 static assert(isSegment3D!T == (type == segment && dim == 3)); 495 static assert(isTriangle3D!T == (type == triangle && dim == 3)); 496 static assert(isSphere3D!T == (type == sphere && dim == 3)); 497 static assert(isRay3D!T == (type == ray && dim == 3)); 498 } 499 } 500 501 with (ShapeType) 502 { 503 // test case type #dimensions 504 test!(seg2f , segment , 2); 505 test!(seg3d , segment , 3); 506 test!(triangle2d , triangle, 2); 507 test!(triangle3f , triangle, 3); 508 test!(sphere2d , sphere , 2); 509 test!(Sphere!(uint, 3), sphere , 3); 510 test!(ray2f , ray , 2); 511 test!(Ray!(real, 3) , ray , 3); 512 test!(planed , plane , 0); // ignore dimension (always 3D) 513 test!(Plane!float , plane , 0); 514 test!(Frustum!double , frustum , 0); 515 } 516 } 517 518 /// Get the numeric type used to measure a shape's dimensions. 519 alias DimensionType(T : Segment!U, U...) = U[0]; 520 /// ditto 521 alias DimensionType(T : Triangle!U, U...) = U[0]; 522 /// ditto 523 alias DimensionType(T : Sphere!U, U...) = U[0]; 524 /// ditto 525 alias DimensionType(T : Ray!U, U...) = U[0]; 526 /// ditto 527 alias DimensionType(T : Plane!U, U) = U; 528 /// ditto 529 alias DimensionType(T : Frustum!U, U) = U; 530 531 /// 532 unittest 533 { 534 static assert(is(DimensionType!seg2i == int)); 535 static assert(is(DimensionType!triangle3d == double)); 536 static assert(is(DimensionType!sphere2d == double)); 537 static assert(is(DimensionType!ray3f == float)); 538 static assert(is(DimensionType!planed == double)); 539 static assert(is(DimensionType!(Frustum!real) == real)); 540 }