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