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 }