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 static if( __VERSION__ < 2066 ) private enum nogc = 1; 21 22 /// A Segment is 2 points. 23 /// When considered like a vector, it represents the arrow from a to b. 24 align(1) struct Segment(T, int N) 25 { 26 align(1): 27 public 28 { 29 alias Vector!(T, N) point_t; 30 point_t a, b; 31 } 32 } 33 34 alias Segment!(float, 2) seg2f; /// 2D float segment. 35 alias Segment!(float, 3) seg3f; /// 3D float segment. 36 alias Segment!(double, 2) seg2d; /// 2D double segment. 37 alias Segment!(double, 3) seg3d; /// 3D double segment. 38 39 /// A Triangle is 3 points. 40 align(1) struct Triangle(T, int N) 41 { 42 align(1): 43 public 44 { 45 alias Vector!(T, N) point_t; 46 point_t a, b, c; 47 48 static if (N == 2) 49 { 50 /// Returns: Area of a 2D triangle. 51 @nogc T area() pure const nothrow 52 { 53 return abs(signedArea()); 54 } 55 56 /// Returns: Signed area of a 2D triangle. 57 @nogc T signedArea() pure const nothrow 58 { 59 return ((b.x * a.y - a.x * b.y) 60 + (c.x * b.y - b.x * c.y) 61 + (a.x * c.y - c.x * a.y)) / 2; 62 } 63 } 64 65 static if (N == 3) 66 { 67 /// Returns: Triangle normal. 68 @nogc Vector!(T, 3) computeNormal() pure const nothrow 69 { 70 return cross(b - a, c - a).normalized(); 71 } 72 } 73 } 74 } 75 76 alias Triangle!(float, 2) triangle2f; /// 2D float triangle. 77 alias Triangle!(float, 3) triangle3f; /// 3D float triangle. 78 alias Triangle!(double, 2) triangle2d; /// 2D double triangle. 79 alias Triangle!(double, 3) triangle3d; /// 3D double triangle. 80 81 /// A Sphere is a point + a radius. 82 align(1) struct Sphere(T, int N) 83 { 84 align(1): 85 public nothrow 86 { 87 alias Vector!(T, N) point_t; 88 point_t center; 89 T radius; 90 91 /// Creates a sphere from a point and a radius. 92 @nogc this(in point_t center_, T radius_) pure nothrow 93 { 94 center = center_; 95 radius = radius_; 96 } 97 98 /// Sphere contains point test. 99 /// Returns: true if the point is inside the sphere. 100 @nogc bool contains(in Sphere s) pure const nothrow 101 { 102 if (s.radius > radius) 103 return false; 104 105 T innerRadius = radius - s.radius; 106 return squaredDistanceTo(s.center) < innerRadius * innerRadius; 107 } 108 109 /// Sphere vs point Euclidean distance squared. 110 @nogc T squaredDistanceTo(point_t p) pure const nothrow 111 { 112 return center.squaredDistanceTo(p); 113 } 114 115 /// Sphere vs sphere intersection. 116 /// Returns: true if the spheres intersect. 117 @nogc bool intersects(Sphere s) pure const nothrow 118 { 119 T outerRadius = radius + s.radius; 120 return squaredDistanceTo(s.center) < outerRadius * outerRadius; 121 } 122 123 static if (isFloatingPoint!T) 124 { 125 /// Sphere vs point Euclidean distance. 126 @nogc T distanceTo(point_t p) pure const nothrow 127 { 128 return center.distanceTo(p); 129 } 130 131 static if(N == 2) 132 { 133 /// Returns: Circle area. 134 @nogc T area() pure const nothrow 135 { 136 return PI * (radius * radius); 137 } 138 } 139 } 140 } 141 } 142 143 alias Sphere!(float, 2) sphere2f; /// 2D float sphere (ie. a circle). 144 alias Sphere!(float, 3) sphere3f; /// 3D float sphere. 145 alias Sphere!(double, 2) sphere2d; /// 2D double sphere (ie. a circle). 146 alias Sphere!(double, 3) sphere3d; /// 3D double sphere (ie. a circle). 147 148 149 /// A Ray ir a point + a direction. 150 align(1) struct Ray(T, int N) 151 { 152 align(1): 153 nothrow: 154 public 155 { 156 alias Vector!(T, N) point_t; 157 point_t orig; 158 point_t dir; 159 160 /// Returns: A point further along the ray direction. 161 @nogc point_t progress(T t) pure const nothrow 162 { 163 return orig + dir * t; 164 } 165 166 static if (N == 3) 167 { 168 /// Ray vs triangle intersection. 169 /// See_also: "Fast, Minimum Storage Ray/Triangle intersection", Mommer & Trumbore (1997) 170 /// Returns: Barycentric coordinates, the intersection point is at $(D (1 - u - v) * A + u * B + v * C). 171 @nogc bool intersect(Triangle!(T, 3) triangle, out T t, out T u, out T v) pure const nothrow 172 { 173 point_t edge1 = triangle.b - triangle.a; 174 point_t edge2 = triangle.c - triangle.a; 175 point_t pvec = cross(dir, edge2); 176 T det = dot(edge1, pvec); 177 if (abs(det) < T.epsilon) 178 return false; // no intersection 179 T invDet = 1 / det; 180 181 // calculate distance from triangle.a to ray origin 182 point_t tvec = orig - triangle.a; 183 184 // calculate U parameter and test bounds 185 u = dot(tvec, pvec) * invDet; 186 if (u < 0 || u > 1) 187 return false; 188 189 // prepare to test V parameter 190 point_t qvec = cross(tvec, edge1); 191 192 // calculate V parameter and test bounds 193 v = dot(dir, qvec) * invDet; 194 if (v < 0.0 || u + v > 1.0) 195 return false; 196 197 // calculate t, ray intersects triangle 198 t = dot(edge2, qvec) * invDet; 199 return true; 200 } 201 } 202 } 203 } 204 205 alias Ray!(float, 2) ray2f; /// 2D float ray. 206 alias Ray!(float, 3) ray3f; /// 3D float ray. 207 alias Ray!(double, 2) ray2d; /// 2D double ray. 208 alias Ray!(double, 3) ray3d; /// 3D double ray. 209 210 211 /// 3D plane. 212 /// See_also: Flipcode article by Nate Miller $(WEB www.flipcode.com/archives/Plane_Class.shtml). 213 align(1) struct Plane(T) if (isFloatingPoint!T) 214 { 215 align(1): 216 public 217 { 218 vec3!T n; /// Normal (always stored normalized). 219 T d; 220 221 /// Create from four coordinates. 222 @nogc this(vec4!T abcd) pure nothrow 223 { 224 n = vec3!T(abcd.x, abcd.y, abcd.z).normalized(); 225 d = abcd.z; 226 } 227 228 /// Create from a point and a normal. 229 @nogc this(vec3!T origin, vec3!T normal) pure nothrow 230 { 231 n = normal.normalized(); 232 d = -dot(origin, n); 233 } 234 235 /// Create from 3 non-aligned points. 236 @nogc this(vec3!T A, vec3!T B, vec3!T C) pure nothrow 237 { 238 this(C, cross(B - A, C - A)); 239 } 240 241 /// Assign a plane with another plane. 242 @nogc ref Plane opAssign(Plane other) pure nothrow 243 { 244 n = other.n; 245 d = other.d; 246 return this; 247 } 248 249 /// Returns: signed distance between a point and the plane. 250 @nogc T signedDistanceTo(vec3!T point) pure const nothrow 251 { 252 return dot(n, point) + d; 253 } 254 255 /// Returns: absolute distance between a point and the plane. 256 @nogc T distanceTo(vec3!T point) pure const nothrow 257 { 258 return abs(signedDistanceTo(point)); 259 } 260 261 /// Returns: true if the point is in front of the plane. 262 @nogc bool isFront(vec3!T point) pure const nothrow 263 { 264 return signedDistanceTo(point) >= 0; 265 } 266 267 /// Returns: true if the point is in the back of the plane. 268 @nogc bool isBack(vec3!T point) pure const nothrow 269 { 270 return signedDistanceTo(point) < 0; 271 } 272 273 /// Returns: true if the point is on the plane, with a given epsilon. 274 @nogc bool isOn(vec3!T point, T epsilon) pure const nothrow 275 { 276 T sd = signedDistanceTo(point); 277 return (-epsilon < sd) && (sd < epsilon); 278 } 279 } 280 } 281 282 alias Plane!float planef; /// 3D float plane. 283 alias Plane!double planed; /// 3D double plane. 284 285 unittest 286 { 287 auto p = planed(vec4d(1.0, 2.0, 3.0, 4.0)); 288 auto p2 = planed(vec3d(1.0, 0.0, 0.0), 289 vec3d(0.0, 1.0, 0.0), 290 vec3d(0.0, 0.0, 1.0)); 291 292 assert(p2.isOn(vec3d(1.0, 0.0, 0.0), 1e-7)); 293 assert(p2.isFront(vec3d(1.0, 1.0, 1.0))); 294 assert(p2.isBack(vec3d(0.0, 0.0, 0.0))); 295 } 296 297 298 299 /// 3D frustum. 300 /// See_also: Flipcode article by Dion Picco $(WEB www.flipcode.com/archives/Frustum_Culling.shtml). 301 /// Bugs: verify proper signedness of half-spaces 302 align(1) struct Frustum(T) if (isFloatingPoint!T) 303 { 304 align(1): 305 public 306 { 307 enum int LEFT = 0, 308 RIGHT = 1, 309 TOP = 2, 310 BOTTOM = 3, 311 NEAR = 4, 312 FAR = 5; 313 314 Plane!T[6] planes; 315 316 /// Create a frustum from 6 planes. 317 @nogc this(Plane!T left, Plane!T right, Plane!T top, Plane!T bottom, Plane!T near, Plane!T far) pure nothrow 318 { 319 planes[LEFT] = left; 320 planes[RIGHT] = right; 321 planes[TOP] = top; 322 planes[BOTTOM] = bottom; 323 planes[NEAR] = near; 324 planes[FAR] = far; 325 } 326 327 enum : int 328 { 329 OUTSIDE, /// object is outside the frustum 330 INTERSECT, /// object intersects with the frustum 331 INSIDE /// object is inside the frustum 332 } 333 334 /// Point vs frustum intersection. 335 @nogc bool contains(vec3!T point) pure const nothrow 336 { 337 for(int i = 0; i < 6; ++i) 338 { 339 T distance = planes[i].signedDistanceTo(point); 340 341 if(distance < 0) 342 return false; 343 } 344 return true; 345 } 346 347 /// Sphere vs frustum intersection. 348 /// Returns: Frustum.OUTSIDE, Frustum.INTERSECT or Frustum.INSIDE. 349 @nogc int contains(Sphere!(T, 3) sphere) pure const nothrow 350 { 351 // calculate our distances to each of the planes 352 for(int i = 0; i < 6; ++i) 353 { 354 // find the distance to this plane 355 T distance = planes[i].signedDistanceTo(sphere.center); 356 357 if(distance < -sphere.radius) 358 return OUTSIDE; 359 360 else if (distance < sphere.radius) 361 return INTERSECT; 362 } 363 364 // otherwise we are fully in view 365 return INSIDE; 366 } 367 368 /// AABB vs frustum intersection. 369 /// Returns: Frustum.OUTSIDE, Frustum.INTERSECT or Frustum.INSIDE. 370 @nogc int contains(box3!T box) pure const nothrow 371 { 372 vec3!T corners[8]; 373 int totalIn = 0; 374 375 for (int i = 0; i < 2; ++i) 376 for (int j = 0; j < 2; ++j) 377 for (int k = 0; k < 2; ++k) 378 { 379 auto x = i == 0 ? box.min.x : box.max.x; 380 auto y = i == 0 ? box.min.y : box.max.y; 381 auto z = i == 0 ? box.min.z : box.max.z; 382 corners[i*4 + j*2 + k] = vec3!T(x, y, z); 383 } 384 385 // test all 8 corners against the 6 sides 386 // if all points are behind 1 specific plane, we are out 387 // if we are in with all points, then we are fully in 388 for(int p = 0; p < 6; ++p) 389 { 390 int inCount = 8; 391 int ptIn = 1; 392 393 for(int i = 0; i < 8; ++i) 394 { 395 // test this point against the planes 396 if (planes[p].isBack(corners[i])) 397 { 398 ptIn = 0; 399 --inCount; 400 } 401 } 402 403 // were all the points outside of plane p? 404 if (inCount == 0) 405 return OUTSIDE; 406 407 // check if they were all on the right side of the plane 408 totalIn += ptIn; 409 } 410 411 // so if totalIn is 6, then all are inside the view 412 if(totalIn == 6) 413 return INSIDE; 414 415 // we must be partly in then otherwise 416 return INTERSECT; 417 } 418 419 } 420 }