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 }