float2.cpp

1 /* Copyright Jukka Jyl�nki
2
3    Licensed under the Apache License, Version 2.0 (the "License");
4    you may not use this file except in compliance with the License.
5    You may obtain a copy of the License at
6
7        http://www.apache.org/licenses/LICENSE-2.0
8
9    Unless required by applicable law or agreed to in writing, software
10    distributed under the License is distributed on an "AS IS" BASIS,
11    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12    See the License for the specific language governing permissions and
13    limitations under the License. */
14
15 /** @file float2.cpp
16         @author Jukka Jyl�nki
17         @brief */
18 #include "float2.h"
19 #include "float3.h"
20 #include "float4.h"
21 #include "MathFunc.h"
22 #include "../Algorithm/Random/LCG.h"
23 #include "assume.h"
24 #include <string.h>
25 #include <stdlib.h>
26 #include <locale.h>
27
28 #ifdef MATH_ENABLE_STL_SUPPORT
29 #include "myassert.h"
30 #include <iostream>
31 #include <utility>
32 #include <algorithm>
33 #endif
34
35 MATH_BEGIN_NAMESPACE
36
37 using namespace std;
38
39 float2::float2(float x_, float y_)
40 :x(x_), y(y_)
41 {
42 }
43
44 float2::float2(float scalar)
45 :x(scalar), y(scalar)
46 {
47 }
48
49 float2::float2(const float *data)
50 {
51         assume(data);
52 #ifndef MATH_ENABLE_INSECURE_OPTIMIZATIONS
53         if (!data)
54                 return;
55 #endif
56         x = data[0];
57         y = data[1];
58 }
59
60 CONST_WIN32 float float2::At(int index) const
61 {
62         assume(index >= 0);
63         assume(index < Size);
64 #ifndef MATH_ENABLE_INSECURE_OPTIMIZATIONS
65         if (index < 0 || index >= Size)
66                 return FLOAT_NAN;
67 #endif
68         return ptr()[index];
69 }
70
71 float &float2::At(int index)
72 {
73         assume(index >= 0);
74         assume(index < Size);
75 #ifndef MATH_ENABLE_INSECURE_OPTIMIZATIONS
76         if (index < 0 || index >= Size)
77                 return ptr()[0];
78 #endif
79         return ptr()[index];
80 }
81
82 float2 float2::Swizzled(int i, int j) const
83 {
84         return float2(At(i), At(j));
85 }
86
87 float3 float2::Swizzled(int i, int j, int k) const
88 {
89         return float3(At(i), At(j), At(k));
90 }
91
92 float4 float2::Swizzled(int i, int j, int k, int l) const
93 {
94         return float4(At(i), At(j), At(k), At(l));
95 }
96
97 float float2::LengthSq() const
98 {
99         return x*x + y*y;
100 }
101
102 float float2::Length() const
103 {
104         return Sqrt(LengthSq());
105 }
106
107 void float2::SetFromPolarCoordinates(float theta, float length)
108 {
109         float sin, cos;
110         SinCos(theta, sin, cos);
111         x = cos * length;
112         y = sin * length;
113 }
114
115 float2 float2::FromPolarCoordinates(float theta, float length)
116 {
117         float2 euclidean;
118         euclidean.SetFromPolarCoordinates(theta, length);
119         return euclidean;
120 }
121
122 float2 float2::ToPolarCoordinates() const
123 {
124         float radius = Length();
125         if (radius > 1e-4f)
126                 return float2(atan2(y, x), radius);
127         else
128                 return float2::zero;
129 }
130
131 float float2::AimedAngle() const
132 {
133         assume(!IsZero());
134         return atan2(y, x);
135 }
136
137 float float2::Normalize()
138 {
139         assume(IsFinite());
140         float lengthSq = LengthSq();
141         if (lengthSq > 1e-6f)
142         {
143                 float length = Sqrt(lengthSq);
144                 *this *= 1.f / length;
145                 return length;
146         }
147         else
148         {
149                 Set(1.f, 0.f); // We will always produce a normalized vector.
150                 return 0; // But signal failure, so user knows we have generated an arbitrary normalization.
151         }
152 }
153
154 float2 float2::Normalized() const
155 {
156         float2 copy = *this;
157         float oldLength = copy.Normalize();
158         assume(oldLength > 0.f && "float2::Normalized() failed!");
159         MARK_UNUSED(oldLength);
160         return copy;
161 }
162
163 float float2::ScaleToLength(float newLength)
164 {
165         float length = LengthSq();
166         if (length < 1e-6f)
167                 return 0.f;
168
169         length = Sqrt(length);
170         float scalar = newLength / length;
171         x *= scalar;
172         y *= scalar;
173         return length;
174 }
175
176 float2 float2::ScaledToLength(float newLength) const
177 {
178         assume(!IsZero());
179
180         float2 v = *this;
181         v.ScaleToLength(newLength);
182         return v;
183 }
184
185 bool float2::IsNormalized(float epsilonSq) const
186 {
187         return MATH_NS::Abs(LengthSq()-1.f) <= epsilonSq;
188 }
189
190 bool float2::IsZero(float epsilonSq) const
191 {
192         return LengthSq() <= epsilonSq;
193 }
194
195 bool float2::IsFinite() const
196 {
197         return MATH_NS::IsFinite(x) && MATH_NS::IsFinite(y);
198 }
199
200 bool float2::IsPerpendicular(const float2 &other, float epsilonSq) const
201 {
202         float dot = Dot(other);
203         return dot*dot <= epsilonSq * LengthSq() * other.LengthSq();
204 }
205
206 bool float2::Equals(const float2 &rhs, float epsilon) const
207 {
208         return EqualAbs(x, rhs.x, epsilon) && EqualAbs(y, rhs.y, epsilon);
209 }
210
211 bool float2::Equals(float x_, float y_, float epsilon) const
212 {
213         return EqualAbs(x, x_, epsilon) && EqualAbs(y, y_, epsilon);
214 }
215
216 bool float2::BitEquals(const float2 &other) const
217 {
218         return ReinterpretAsU32(x) == ReinterpretAsU32(other.x) &&
219                 ReinterpretAsU32(y) == ReinterpretAsU32(other.y);
220 }
221
222 /// It is too performance-heavy to set the locale in each serialization and deserialization function call.
223 /// Therefore expect the user to has a proper locale set up for the application at startup. This is assert()ed
224 /// at debug runs.
225 bool IsNeutralCLocale()
226 {
227         // Android NDK locale.h does not have struct lconv or localeconv() implemented, and only contains stub
228         // symbols with a comment 'MISSING FROM BIONIC - DEFINED TO MAKE libstdc++-v3 happy'
229 #ifndef ANDROID
230         lconv *lc = localeconv();
231         if (strcmp(lc->decimal_point, "."))
232                 return false;
233 #endif
234         return true;
235 }
236
237 #ifdef MATH_ENABLE_STL_SUPPORT
238 std::string float2::ToString() const
239 {
240         char str[256];
241         sprintf(str, "(%f, %f)", x, y);
242         return std::string(str);
243 }
244
245 std::string float2::SerializeToString() const
246 {
247         char str[256];
248         char *s = SerializeFloat(x, str); *s = ','; ++s;
249         s = SerializeFloat(y, s);
250         assert(s+1 - str < 256);
251         MARK_UNUSED(s);
252         return str;
253 }
254
255 std::string float2::SerializeToCodeString() const
256 {
257         return "float2(" + SerializeToString() + ")";
258 }
259 #endif
260
261 float2 float2::FromString(const char *str, const char **outEndStr)
262 {
263         assert(IsNeutralCLocale());
264         assume(str);
265         if (!str)
266                 return float2::nan;
267         MATH_SKIP_WORD(str, "float2");
268         MATH_SKIP_WORD(str, "(");
269         float2 f;
270         f.x = DeserializeFloat(str, &str);
271         f.y = DeserializeFloat(str, &str);
272         if (*str == ')')
273                 ++str;
274         if (*str == ',')
275                 ++str;
276         if (outEndStr)
277                 *outEndStr = str;
278         return f;
279 }
280
281 float float2::SumOfElements() const
282 {
283         return x + y;
284 }
285
286 float float2::ProductOfElements() const
287 {
288         return x * y;
289 }
290
291 float float2::AverageOfElements() const
292 {
293         return (x + y) * 0.5f;
294 }
295
296 float float2::MinElement() const
297 {
298         return MATH_NS::Min(x, y);
299 }
300
301 int float2::MinElementIndex() const
302 {
303         return (x <= y) ? 0 : 1;
304 }
305
306 float float2::MaxElement() const
307 {
308         return  MATH_NS::Max(x, y);
309 }
310
311 int float2::MaxElementIndex() const
312 {
313         return (x > y) ? 0 : 1;
314 }
315
316 float2 float2::Abs() const
317 {
318         return float2(MATH_NS::Abs(x), MATH_NS::Abs(y));
319 }
320
321 float2 float2::Neg() const
322 {
323         return float2(-x, -y);
324 }
325
326 float2 float2::Recip() const
327 {
328         return float2(1.f/x, 1.f/y);
329 }
330
331 float2 float2::Min(float floor) const
332 {
333         return float2(MATH_NS::Min(x, floor),  MATH_NS::Min(x, floor));
334 }
335
336 float2 float2::Min(const float2 &floor) const
337 {
338         return float2(MATH_NS::Min(x, floor.x),  MATH_NS::Min(x, floor.x));
339 }
340
341 float2 float2::Max(float ceil) const
342 {
343         return float2(MATH_NS::Max(x, ceil),  MATH_NS::Max(x, ceil));
344 }
345
346 float2 float2::Max(const float2 &ceil) const
347 {
348         return float2(MATH_NS::Max(x, ceil.x),  MATH_NS::Max(x, ceil.x));
349 }
350
351 float2 float2::Clamp(const float2 &floor, const float2 &ceil) const
352 {
353         return float2(MATH_NS::Clamp(x, floor.x, ceil.x),  MATH_NS::Clamp(y, floor.y, ceil.y));
354 }
355
356 float2 float2::Clamp(float floor, float ceil) const
357 {
358         return float2(MATH_NS::Clamp(x, floor, ceil),  MATH_NS::Clamp(y, floor, ceil));
359 }
360
361 float2 float2::Clamp01() const
362 {
363         return Clamp(0.f, 1.f);
364 }
365
366 float float2::DistanceSq(const float2 &rhs) const
367 {
368         float dx = x - rhs.x;
369         float dy = y - rhs.y;
370         return dx*dx + dy*dy;
371 }
372
373 float float2::Distance(const float2 &rhs) const
374 {
375         return Sqrt(DistanceSq(rhs));
376 }
377
378 float float2::Dot(const float2 &rhs) const
379 {
380         return x * rhs.x + y * rhs.y;
381 }
382
383 float2 float2::Perp() const
384 {
385         return float2(-y, x);
386 }
387
388 float float2::PerpDot(const float2 &rhs) const
389 {
390         return x * rhs.y - y * rhs.x;
391 }
392
393 float2 float2::Reflect(const float2 &normal) const
394 {
395         assume2(normal.IsNormalized(), normal.SerializeToCodeString(), normal.Length());
396         return 2.f * this->ProjectToNorm(normal) - *this;
397 }
398
399 /// Implementation from http://www.flipcode.com/archives/reflection_transmission.pdf .
400 float2 float2::Refract(const float2 &normal, float negativeSideRefractionIndex, float positiveSideRefractionIndex) const
401 {
402         // This code is duplicated in float3::Refract.
403         float n = negativeSideRefractionIndex / positiveSideRefractionIndex;
404         float cosI = this->Dot(normal);
405         float sinT2 = n*n*(1.f - cosI*cosI);
406         if (sinT2 > 1.f) // Total internal reflection occurs?
407                 return (-*this).Reflect(normal);
408         return n * *this - (n + Sqrt(1.f - sinT2)) * normal;
409 }
410
411 float2 float2::ProjectTo(const float2 &direction) const
412 {
413         assume(!direction.IsZero());
414         return direction * this->Dot(direction) / direction.LengthSq();
415 }
416
417 float2 float2::ProjectToNorm(const float2 &direction) const
418 {
419         assume(direction.IsNormalized());
420         return direction * this->Dot(direction);
421 }
422
423 float float2::AngleBetween(const float2 &other) const
424 {
425         return acos(Dot(other)) / Sqrt(LengthSq() * other.LengthSq());
426 }
427
428 float float2::AngleBetweenNorm(const float2 &other) const
429 {
430         assume(this->IsNormalized());
431         assume(other.IsNormalized());
432         return acos(Dot(other));
433 }
434
435 float2 float2::Lerp(const float2 &b, float t) const
436 {
437         assume(0.f <= t && t <= 1.f);
438         return (1.f - t) * *this + t * b;
439 }
440
441 float2 float2::Lerp(const float2 &a, const float2 &b, float t)
442 {
443         return a.Lerp(b, t);
444 }
445
446 void float2::Decompose(const float2 &direction, float2 &outParallel, float2 &outPerpendicular) const
447 {
448         assume(direction.IsNormalized());
449         outParallel = this->Dot(direction) * direction;
450         outPerpendicular = *this - outParallel;
451 }
452
453 void float2::Orthogonalize(const float2 &a, float2 &b)
454 {
455         assume(!a.IsZero());
456         b -= a.Dot(b) / a.Length() * a;
457 }
458
459 bool float2::AreOrthogonal(const float2 &a, const float2 &b, float epsilon)
460 {
461         return a.IsPerpendicular(b, epsilon);
462 }
463
464
465 void float2::Orthonormalize(float2 &a, float2 &b)
466 {
467         assume(!a.IsZero());
468         a.Normalize();
469         b -= a.Dot(b) * a;
470 }
471
472 float2 float2::FromScalar(float scalar)
473 {
474         return float2(scalar, scalar);
475 }
476
477 void float2::SetFromScalar(float scalar)
478 {
479         x = scalar;
480         y = scalar;
481 }
482
483 void float2::Set(float x_, float y_)
484 {
485         x = x_;
486         y = y_;
487 }
488
489 void float2::Rotate90CW()
490 {
491         float oldX = x;
492         x = y;
493         y = -oldX;
494 }
495
496 float2 float2::Rotated90CW() const
497 {
498         return float2(y, -x);
499 }
500
501 void float2::Rotate90CCW()
502 {
503         float oldX = x;
504         x = -y;
505         y = oldX;
506 }
507
508 float2 float2::Rotated90CCW() const
509 {
510         return float2(-y, x);
511 }
512
513 bool float2::OrientedCCW(const float2 &a, const float2 &b, const float2 &c)
514 {
515         // Compute the determinant
516         // | ax ay 1 |
517         // | bx by 1 |
518         // | cx cy 1 |
519         // See Christer Ericson, Real-Time Collision Detection, p.32.
520         return (a.x-c.x)*(b.y-c.y) - (a.y-c.y)*(b.x-c.x) >= 0.f;
521 }
522
523 class SortByPolarAngle
524 {
525 public:
526         float2 perspective;
527
528         bool operator()(const float2 &a, const float2 &b) const
529         {
530                 float2 A = a - perspective;
531                 float2 B = b - perspective;
532                 return B.x*A.y < A.x*B.y;
533         }
534 };
535
536 #ifdef MATH_ENABLE_STL_SUPPORT
537 void float2::ConvexHull(const float2 *pointArray, int numPoints, std::vector<float2> &outConvexHull)
538 {
539         outConvexHull.clear();
540         if (numPoints == 0)
541                 return;
542         outConvexHull.insert(outConvexHull.end(), pointArray, pointArray + numPoints);
543         int convexHullSize = ConvexHullInPlace(&outConvexHull[0], (int)outConvexHull.size());
544         outConvexHull.resize(convexHullSize);
545 }
546 #endif
547
548 #ifdef MATH_ENABLE_STL_SUPPORT
549 /** This function implements the Graham's Scan algorithm for finding the convex hull of
550         a 2D point set. The running time is O(nlogn). For details, see
551         "Introduction to Algorithms, 2nd ed.", by Cormen, Leiserson, Rivest, p.824, or
552         a lecture by Shai Simonson: http://www.aduni.org/courses/algorithms/index.php?view=cw , lecture 02-13-01. */
553 int float2::ConvexHullInPlace(float2 *p, int n)
554 {
555         if (n <= 2)
556                 return n;
557
558         if (n >= 50)
559         {
560                 /* Perform Akl�Toussaint heuristic. The limit n=50 is arbitrary and based on quick profiling:
561                  Without heuristic:
562                    n=10: 1143 ticks
563                    n=50: 8657 ticks
564                    n=100: 19533 ticks
565                  With heuristic:
566                    n=10: 1322 ticks
567                    n=50: 6759 ticks
568                    n=100: 14448 ticks
569                 */
570                 int minX = 0, minY = 0, maxX = 0, maxY = 0;
571                 for(int i = 1; i < n; ++i)
572                 {
573                         if (p[i].x < p[minX].x) minX = i;
574                         else if (p[i].x > p[maxX].x) maxX = i;
575                         if (p[i].y < p[minY].y) minY = i;
576                         else if (p[i].y > p[maxY].y) maxY = i;
577                 }
578                 // Direction vectors which point inside the convex hull.
579                 float2 e0 = (p[maxX] - p[minY]).Rotated90CCW();
580                 float2 e1 = (p[maxY] - p[maxX]).Rotated90CCW();
581                 float2 e2 = (p[minX] - p[maxY]).Rotated90CCW();
582                 float2 e3 = (p[minY] - p[minX]).Rotated90CCW();
583
584                 // Add a small epsilon so that the four extreme points on the convex hull will not get pruned
585                 // due to floating point imprecision.
586                 const float eps = 1e-6f;
587                 float e0_d = e0.Dot(p[minY]) + eps;
588                 float e1_d = e1.Dot(p[maxX]) + eps;
589                 float e2_d = e2.Dot(p[maxY]) + eps;
590                 float e3_d = e3.Dot(p[minX]) + eps;
591
592                 for(int i = 0; i < n; ++i)
593                         if (e0.Dot(p[i]) > e0_d && e1.Dot(p[i]) > e1_d && e2.Dot(p[i]) > e2_d && e3.Dot(p[i]) > e3_d)
594                                 Swap(p[i--], p[--n]);
595         }
596
597         // Find the lowest point of the set.
598         SortByPolarAngle pred;
599         pred.perspective = p[0];
600         int smallestY = 0;
601         for(int i = 1; i < n; ++i)
602                 if (p[i].y < pred.perspective.y || (p[i].y == pred.perspective.y && p[i].x < pred.perspective.x))
603                 {
604                         pred.perspective = p[i];
605                         smallestY = i;
606                 }
607         Swap(p[0], p[smallestY]);
608
609         // For robustness, remove duplicates of the perspective pivot points.
610         // This is because duplicates on that element will cause the sorting to be nontransitive and break
611         // the whole sort.
612         int d = 0;
613         for(int i = 1; i < n; ++i)
614                 if (!p[i].Equals(p[0]))
615                         p[++d] = p[i];
616         n = d+1;
617
618         std::sort(&p[1], &p[n], pred);
619
620         // For robustness, remove duplicate input values.
621         d = 0;
622         for(int i = 1; i < n; ++i)
623                 if (!p[i].Equals(p[d]))
624                         p[++d] = p[i];
625         n = d+1;
626
627         int h = 1; // Points to the index of the last point added to the hull so far. The first two points are in the hull to start.
628
629         float2 a = p[h] - p[h-1];
630         const float epsilon = 1e-5f;
631         for(int i = 2; i < n; ++i)
632         {
633                 // The last two added points determine a line, check which side of that line the next point to be added lies in.
634                 float2 d = p[i] - p[h-1];
635                 float dir = d.x*a.y - d.y*a.x;
636                 // Remove previous points from the convex hull until we have a left turn. Also for numerical stability,
637                 // in the case of three collinear points, remove the middle point.
638                 while(dir > epsilon || (dir > -epsilon && d.Dot(d) >= a.Dot(a)))
639                 {
640                         --h;
641                         if (h >= 1)
642                         {
643                                 a = p[h] - p[h-1];
644                                 d = p[i] - p[h-1];
645                                 dir = d.x*a.y - d.y*a.x;
646                         }
647                         else
648                                 break;
649                 }
650                 p[++h] = p[i];
651                 a = p[i] - p[h-1];
652         }
653
654         // Return the number of points on the new hull.
655         return h+1;
656 }
657
658 bool float2::ConvexHullContains(const float2 *convexHull, int numPointsInConvexHull, const float2 &point)
659 {
660         int j = numPointsInConvexHull-1;
661         for(int i = 0; i < numPointsInConvexHull; ++i)
662         {
663                 float2 d = (convexHull[i] - convexHull[j]).Rotated90CCW(); // Points inwards the convex hull.
664                 float2 n = point - convexHull[j];
665                 if (n.IsZero()) return true;
666                 if (n.Dot(d) < 0.f)
667                         return false;
668                 j = i;
669         }
670         return true;
671 }
672
673 #endif
674
675 #define NEXT_P(ptr) ((ptr)+1 < (pEnd) ? (ptr)+1 : (p))
676
677 float float2::MinAreaRectInPlace(float2 *p, int n, float2 &center, float2 &uDir, float2 &vDir, float &minU, float &maxU, float &minV, float &maxV)
678 {
679         assume(p || n == 0);
680         if (!p)
681                 return 0.f;
682
683         // As a preparation, need to compute the convex hull so that points are CCW-oriented,
684         // and this also greatly reduces the number of points for performance.
685         n = float2::ConvexHullInPlace(p, n);
686
687         // e[i] point to the antipodal point pairs: e[0] and e[2] are pairs, so are e[1] and e[3].
688         float2 *e[4] = { p, p, p, p };
689
690         // Compute the initial AABB rectangle antipodal points for the rotating calipers method.
691         // Order the initial vertices minX -> minY -> maxX -> maxY to establish
692         // a counter-clockwise orientation.
693         for(int i = 1; i < n; ++i)
694         {
695                 if (p[i].x < e[0]->x) e[0] = &p[i];
696                 else if (p[i].x > e[2]->x) e[2] = &p[i];
697                 if (p[i].y < e[1]->y) e[1] = &p[i];
698                 else if (p[i].y > e[3]->y) e[3] = &p[i];
699         }
700
701         // Direction vector of the edge that the currently tested rectangle is in contact with.
702         // This specifies the reference frame for the rectangle, and this is the direction the
703         // convex hull points toward at the antipodal point e[0].
704         float2 ed = -float2::unitY;
705         float minArea = FLOAT_INF; // Track the area of the best rectangle seen so far.
706         const float2 * const pEnd = p + n; // For wraparound testing in NEXT_P().
707
708         // These track directions the convex hull is pointing towards at each antipodal point.
709         float2 d[4];
710         d[0] = (*NEXT_P(e[0]) - *e[0]).Normalized();
711         d[1] = (*NEXT_P(e[1]) - *e[1]).Normalized();
712         d[2] = (*NEXT_P(e[2]) - *e[2]).Normalized();
713         d[3] = (*NEXT_P(e[3]) - *e[3]).Normalized();
714
715         // Rotate the calipers 90 degrees to see through each possible edge that might support
716         // the bounding rectangle.
717         while(ed.y <= 0.f)
718         {
719                 // Compute how much each edge can at most rotate before hitting the next vertex in the convex hull.
720                 float cosA0 =  ed.Dot(d[0]);
721                 float cosA1 =  ed.PerpDot(d[1]);
722                 float cosA2 = -ed.Dot(d[2]);
723                 float cosA3 = -ed.PerpDot(d[3]);
724
725                 float maxCos = MATH_NS::Max(MATH_NS::Max(cosA0, cosA1), MATH_NS::Max(cosA2, cosA3));
726                 // Pick the smallest angle (largest cosine of that angle) and increment the antipodal point index to travel the edge.
727                 if (cosA0 >= maxCos)      { ed = d[0];                e[0] = NEXT_P(e[0]); d[0] = (*NEXT_P(e[0]) - *e[0]).Normalized(); }
728                 else if (cosA1 >= maxCos) { ed = d[1].Rotated90CW();  e[1] = NEXT_P(e[1]); d[1] = (*NEXT_P(e[1]) - *e[1]).Normalized(); }
729                 else if (cosA2 >= maxCos) { ed = -d[2];               e[2] = NEXT_P(e[2]); d[2] = (*NEXT_P(e[2]) - *e[2]).Normalized(); }
730                 else                      { ed = d[3].Rotated90CCW(); e[3] = NEXT_P(e[3]); d[3] = (*NEXT_P(e[3]) - *e[3]).Normalized(); }
731
732                 // Check if the area of the new rectangle is smaller than anything seen so far.
733                 float minu = ed.PerpDot(*e[0]);
734                 float maxu = ed.PerpDot(*e[2]);
735                 float minv = ed.Dot(*e[1]);
736                 float maxv = ed.Dot(*e[3]);
737
738                 float area = MATH_NS::Abs((maxu-minu) * (maxv-minv));
739                 if (area < minArea)
740                 {
741                         vDir = ed;
742                         minArea = area;
743                         minU = MATH_NS::Min(minu, maxu);
744                         maxU = MATH_NS::Max(minu, maxu);
745                         minV = MATH_NS::Min(minv, maxv);
746                         maxV = MATH_NS::Max(minv, maxv);
747                 }
748         }
749         uDir = vDir.Rotated90CCW();
750         center = 0.5f * (uDir * (minU+maxU) + vDir * (minV+maxV));
751
752         return minArea;
753 }
754
755 float2 float2::RandomDir(LCG &lcg, float r)
756 {
757         assume(r > 1e-3f);
758         for(int i = 0; i < 1000; ++i)
759         {
760                 float x = lcg.Float(-r, r);
761                 float y = lcg.Float(-r, r);
762                 float lenSq = x*x + y*y;
763                 if (lenSq >= 1e-6f && lenSq <= r*r)
764                         return r / Sqrt(lenSq) * float2(x,y);
765         }
766         assume(false && "Failed to generate a random float2 direction vector!");
767         return float2(r, 0);
768 }
769
770 MUST_USE_RESULT float2 float2::RandomBox(LCG &lcg, float minElem, float maxElem)
771 {
772         float x = lcg.Float(minElem, maxElem);
773         float y = lcg.Float(minElem, maxElem);
774         return float2(x, y);
775 }
776
777 float2 float2::operator +(const float2 &rhs) const
778 {
779         return float2(x + rhs.x, y + rhs.y);
780 }
781
782 float2 float2::operator -(const float2 &rhs) const
783 {
784         return float2(x - rhs.x, y - rhs.y);
785 }
786
787 float2 float2::operator -() const
788 {
789         return float2(-x, -y);
790 }
791
792 float2 float2::operator *(float scalar) const
793 {
794         return float2(x * scalar, y * scalar);
795 }
796
797 float2 operator *(float scalar, const float2 &rhs)
798 {
799         return float2(scalar * rhs.x, scalar * rhs.y);
800 }
801
802 float2 float2::operator /(float scalar) const
803 {
804         float invScalar = 1.f / scalar;
805         return float2(x * invScalar, y * invScalar);
806 }
807
808 float2 &float2::operator +=(const float2 &rhs)
809 {
810         x += rhs.x;
811         y += rhs.y;
812
813         return *this;
814 }
815
816 float2 &float2::operator -=(const float2 &rhs)
817 {
818         x -= rhs.x;
819         y -= rhs.y;
820
821         return *this;
822 }
823
824 float2 &float2::operator *=(float scalar)
825 {
826         x *= scalar;
827         y *= scalar;
828
829         return *this;
830 }
831
832 float2 float2::Add(float s) const
833 {
834         return float2(x + s, y + s);
835 }
836
837 float2 float2::Sub(float s) const
838 {
839         return float2(x - s, y - s);
840 }
841
842 float2 float2::SubLeft(float s) const
843 {
844         return float2(s - x, s - y);
845 }
846
847 float2 float2::DivLeft(float s) const
848 {
849         return float2(s / x, s / y);
850 }
851
852 float2 float2::Mul(const float2 &rhs) const
853 {
854         return float2(x * rhs.x, y * rhs.y);
855 }
856
857 float2 float2::Div(const float2 &rhs) const
858 {
859         return float2(x / rhs.x, y / rhs.y);
860 }
861
862 float2 &float2::operator /=(float scalar)
863 {
864         float invScalar = 1.f / scalar;
865         x *= invScalar;
866         y *= invScalar;
867
868         return *this;
869 }
870
871 #ifdef MATH_ENABLE_STL_SUPPORT
872 std::ostream &operator <<(std::ostream &out, const float2 &rhs)
873 {
874         std::string str = rhs.ToString();
875         out << str;
876         return out;
877 }
878 #endif
879
880 const float2 float2::zero = float2(0, 0);
881 const float2 float2::one = float2(1, 1);
882 const float2 float2::unitX = float2(1, 0);
883 const float2 float2::unitY = float2(0, 1);
884 const float2 float2::nan = float2(FLOAT_NAN, FLOAT_NAN);
885 const float2 float2::inf = float2(FLOAT_INF, FLOAT_INF);
886
887 MATH_END_NAMESPACE