diff --git a/geometry3Sharp.csproj b/geometry3Sharp.csproj
index fd5761a9..c68413e1 100644
--- a/geometry3Sharp.csproj
+++ b/geometry3Sharp.csproj
@@ -131,9 +131,19 @@
+
+
+
+
+
+
+
+
+
+
diff --git a/intersection/IntrLine2Triangle2.cs b/intersection/IntrLine2Triangle2.cs
new file mode 100644
index 00000000..b1ad8f45
--- /dev/null
+++ b/intersection/IntrLine2Triangle2.cs
@@ -0,0 +1,194 @@
+using System;
+
+namespace g3
+{
+ // ported from WildMagic5
+ public class IntrLine2Triangle2
+ {
+ Line2d line;
+ public Line2d Line
+ {
+ get { return line; }
+ set { line = value; Result = IntersectionResult.NotComputed; }
+ }
+
+ Triangle2d triangle;
+ public Triangle2d Triangle
+ {
+ get { return triangle; }
+ set { triangle = value; Result = IntersectionResult.NotComputed; }
+ }
+
+ public int Quantity = 0;
+ public IntersectionResult Result = IntersectionResult.NotComputed;
+ public IntersectionType Type = IntersectionType.Empty;
+
+ public bool IsSimpleIntersection {
+ get { return Result == IntersectionResult.Intersects && Type == IntersectionType.Point; }
+ }
+
+
+ public Vector2d Point0;
+ public Vector2d Point1;
+ public double Param0;
+ public double Param1;
+
+
+ public IntrLine2Triangle2(Line2d l, Triangle2d t)
+ {
+ line = l; triangle = t;
+ }
+
+
+ public IntrLine2Triangle2 Compute()
+ {
+ Find();
+ return this;
+ }
+
+
+ public bool Find()
+ {
+ if (Result != IntersectionResult.NotComputed)
+ return (Result == IntersectionResult.Intersects);
+
+ // [RMS] if either line direction is not a normalized vector,
+ // results are garbage, so fail query
+ if ( line.Direction.IsNormalized == false ) {
+ Type = IntersectionType.Empty;
+ Result = IntersectionResult.InvalidQuery;
+ return false;
+ }
+
+ Vector3d dist = Vector3d.Zero;
+ Vector3i sign = Vector3i.Zero;
+ int positive = 0, negative = 0, zero = 0;
+ TriangleLineRelations(line.Origin, line.Direction, triangle,
+ ref dist, ref sign, ref positive, ref negative, ref zero);
+
+ if (positive == 3 || negative == 3) {
+ // No intersections.
+ Quantity = 0;
+ Type = IntersectionType.Empty;
+ } else {
+ Vector2d param = Vector2d.Zero;
+ GetInterval(line.Origin, line.Direction, triangle, dist, sign, ref param);
+
+ Intersector1 intr = new Intersector1(param[0], param[1], -double.MaxValue, +double.MaxValue);
+ intr.Find();
+
+ Quantity = intr.NumIntersections;
+ if (Quantity == 2) {
+ // Segment intersection.
+ Type = IntersectionType.Segment;
+ Param0 = intr.GetIntersection(0);
+ Point0 = line.Origin + Param0*line.Direction;
+ Param1 = intr.GetIntersection(1);
+ Point1 = line.Origin + Param1*line.Direction;
+ } else if (Quantity == 1) {
+ // Point intersection.
+ Type = IntersectionType.Point;
+ Param0 = intr.GetIntersection(0);
+ Point0 = line.Origin + Param0*line.Direction;
+ } else {
+ // No intersections.
+ Type = IntersectionType.Empty;
+ }
+ }
+
+ Result = (Type != IntersectionType.Empty) ?
+ IntersectionResult.Intersects : IntersectionResult.NoIntersection;
+ return (Result == IntersectionResult.Intersects);
+ }
+
+
+
+ public static void TriangleLineRelations (
+ Vector2d origin, Vector2d direction,
+ Triangle2d tri, ref Vector3d dist, ref Vector3i sign,
+ ref int positive, ref int negative, ref int zero)
+ {
+ positive = 0;
+ negative = 0;
+ zero = 0;
+ for (int i = 0; i < 3; ++i) {
+ Vector2d diff = tri[i] - origin;
+ dist[i] = diff.DotPerp(direction);
+ if (dist[i] > MathUtil.ZeroTolerance)
+ {
+ sign[i] = 1;
+ ++positive;
+ }
+ else if (dist[i] < -MathUtil.ZeroTolerance)
+ {
+ sign[i] = -1;
+ ++negative;
+ }
+ else
+ {
+ dist[i] = 0.0;
+ sign[i] = 0;
+ ++zero;
+ }
+ }
+ }
+
+
+
+ public static void GetInterval (Vector2d origin, Vector2d direction, Triangle2d tri,
+ Vector3d dist, Vector3i sign, ref Vector2d param)
+ {
+ // Project triangle onto line.
+ Vector3d proj = Vector3d.Zero;
+ int i;
+ for (i = 0; i < 3; ++i) {
+ Vector2d diff = tri[i] - origin;
+ proj[i] = direction.Dot(diff);
+ }
+
+ // Compute transverse intersections of triangle edges with line.
+ double numer, denom;
+ int i0, i1, i2;
+ int quantity = 0;
+ for (i0 = 2, i1 = 0; i1 < 3; i0 = i1++) {
+ if (sign[i0]*sign[i1] < 0) {
+ if ( quantity >= 2 )
+ throw new Exception("IntrLine2Triangle2.GetInterval: too many intersections!");
+ numer = dist[i0]*proj[i1] - dist[i1]*proj[i0];
+ denom = dist[i0] - dist[i1];
+ param[quantity++] = numer/denom;
+ }
+ }
+
+ // Check for grazing contact.
+ if (quantity < 2) {
+ for (i0 = 1, i1 = 2, i2 = 0; i2 < 3; i0 = i1, i1 = i2++) {
+ if (sign[i2] == 0) {
+ if ( quantity >= 2 )
+ throw new Exception("IntrLine2Triangle2.GetInterval: too many intersections!");
+ param[quantity++] = proj[i2];
+ }
+ }
+ }
+
+ // Sort.
+ if ( quantity < 1 )
+ throw new Exception("IntrLine2Triangle2.GetInterval: need at least one intersection");
+ if (quantity == 2)
+ {
+ if (param[0] > param[1])
+ {
+ double save = param[0];
+ param[0] = param[1];
+ param[1] = save;
+ }
+ }
+ else
+ {
+ param[1] = param[0];
+ }
+ }
+
+
+ }
+}
diff --git a/intersection/IntrLine3Box3.cs b/intersection/IntrLine3Box3.cs
new file mode 100644
index 00000000..c59c5e90
--- /dev/null
+++ b/intersection/IntrLine3Box3.cs
@@ -0,0 +1,198 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+
+namespace g3
+{
+ public class IntrLine3Box3
+ {
+ Line3d line;
+ public Line3d Line
+ {
+ get { return line; }
+ set { line = value; Result = IntersectionResult.NotComputed; }
+ }
+
+ Box3d box;
+ public Box3d Box
+ {
+ get { return box; }
+ set { box = value; Result = IntersectionResult.NotComputed; }
+ }
+
+ public int Quantity = 0;
+ public IntersectionResult Result = IntersectionResult.NotComputed;
+ public IntersectionType Type = IntersectionType.Empty;
+
+ public bool IsSimpleIntersection {
+ get { return Result == IntersectionResult.Intersects && Type == IntersectionType.Point; }
+ }
+
+ public double LineParam0, LineParam1;
+ public Vector3d Point0 = Vector3d.Zero;
+ public Vector3d Point1 = Vector3d.Zero;
+
+ public IntrLine3Box3(Line3d l, Box3d b)
+ {
+ line = l; box = b;
+ }
+
+ public IntrLine3Box3 Compute()
+ {
+ Find();
+ return this;
+ }
+
+
+ public bool Find()
+ {
+ if (Result != IntersectionResult.NotComputed)
+ return (Result == IntersectionResult.Intersects);
+
+ // [RMS] if either line direction is not a normalized vector,
+ // results are garbage, so fail query
+ if ( line.Direction.IsNormalized == false ) {
+ Type = IntersectionType.Empty;
+ Result = IntersectionResult.InvalidQuery;
+ return false;
+ }
+
+ LineParam0 = -double.MaxValue;
+ LineParam1 = double.MaxValue;
+ DoClipping(ref LineParam0, ref LineParam1, line.Origin, line.Direction, box,
+ true, ref Quantity, ref Point0, ref Point1, ref Type);
+
+ Result = (Type != IntersectionType.Empty) ?
+ IntersectionResult.Intersects : IntersectionResult.NoIntersection;
+ return (Result == IntersectionResult.Intersects);
+ }
+
+
+
+
+ public bool Test ()
+ {
+ Vector3d AWdU = Vector3d.Zero;
+ Vector3d AWxDdU = Vector3d.Zero;
+ double RHS;
+
+ Vector3d diff = line.Origin - box.Center;
+ Vector3d WxD = line.Direction.Cross(diff);
+
+ AWdU[1] = Math.Abs(line.Direction.Dot(box.AxisY));
+ AWdU[2] = Math.Abs(line.Direction.Dot(box.AxisZ));
+ AWxDdU[0] = Math.Abs(WxD.Dot(box.AxisX));
+ RHS = box.Extent.y*AWdU[2] + box.Extent.z*AWdU[1];
+ if (AWxDdU[0] > RHS) {
+ return false;
+ }
+
+ AWdU[0] = Math.Abs(line.Direction.Dot(box.AxisX));
+ AWxDdU[1] = Math.Abs(WxD.Dot(box.AxisY));
+ RHS = box.Extent.x*AWdU[2] + box.Extent.z*AWdU[0];
+ if (AWxDdU[1] > RHS) {
+ return false;
+ }
+
+ AWxDdU[2] = Math.Abs(WxD.Dot(box.AxisZ));
+ RHS = box.Extent.x*AWdU[1] + box.Extent.y*AWdU[0];
+ if (AWxDdU[2] > RHS) {
+ return false;
+ }
+
+ return true;
+ }
+
+
+
+
+ static public bool DoClipping (ref double t0, ref double t1,
+ Vector3d origin, Vector3d direction,
+ Box3d box, bool solid, ref int quantity,
+ ref Vector3d point0, ref Vector3d point1,
+ ref IntersectionType intrType)
+ {
+ // Convert linear component to box coordinates.
+ Vector3d diff = origin - box.Center;
+ Vector3d BOrigin = new Vector3d(
+ diff.Dot(box.AxisX),
+ diff.Dot(box.AxisY),
+ diff.Dot(box.AxisZ)
+ );
+ Vector3d BDirection = new Vector3d(
+ direction.Dot(box.AxisX),
+ direction.Dot(box.AxisY),
+ direction.Dot(box.AxisZ)
+ );
+
+ double saveT0 = t0, saveT1 = t1;
+ bool notAllClipped =
+ Clip(+BDirection.x, -BOrigin.x-box.Extent.x, ref t0, ref t1) &&
+ Clip(-BDirection.x, +BOrigin.x-box.Extent.x, ref t0, ref t1) &&
+ Clip(+BDirection.y, -BOrigin.y-box.Extent.y, ref t0, ref t1) &&
+ Clip(-BDirection.y, +BOrigin.y-box.Extent.y, ref t0, ref t1) &&
+ Clip(+BDirection.z, -BOrigin.z-box.Extent.z, ref t0, ref t1) &&
+ Clip(-BDirection.z, +BOrigin.z-box.Extent.z, ref t0, ref t1);
+
+ if (notAllClipped && (solid || t0 != saveT0 || t1 != saveT1)) {
+ if (t1 > t0) {
+ intrType = IntersectionType.Segment;
+ quantity = 2;
+ point0 = origin + t0*direction;
+ point1 = origin + t1*direction;
+ } else {
+ intrType = IntersectionType.Point;
+ quantity = 1;
+ point0 = origin + t0*direction;
+ }
+ } else {
+ quantity = 0;
+ intrType = IntersectionType.Empty;
+ }
+
+ return intrType != IntersectionType.Empty;
+ }
+
+
+
+
+ static public bool Clip (double denom, double numer, ref double t0, ref double t1)
+ {
+ // Return value is 'true' if line segment intersects the current test
+ // plane. Otherwise 'false' is returned in which case the line segment
+ // is entirely clipped.
+
+ if (denom > (double)0)
+ {
+ if (numer > denom*t1)
+ {
+ return false;
+ }
+ if (numer > denom*t0)
+ {
+ t0 = numer/denom;
+ }
+ return true;
+ }
+ else if (denom < (double)0)
+ {
+ if (numer > denom*t0)
+ {
+ return false;
+ }
+ if (numer > denom*t1)
+ {
+ t1 = numer/denom;
+ }
+ return true;
+ }
+ else
+ {
+ return numer <= (double)0;
+ }
+ }
+
+
+ }
+}
diff --git a/intersection/IntrRay3Box3.cs b/intersection/IntrRay3Box3.cs
new file mode 100644
index 00000000..03b059d6
--- /dev/null
+++ b/intersection/IntrRay3Box3.cs
@@ -0,0 +1,141 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+
+namespace g3
+{
+ public class IntrRay3Box3
+ {
+ Ray3d ray;
+ public Ray3d Ray
+ {
+ get { return ray; }
+ set { ray = value; Result = IntersectionResult.NotComputed; }
+ }
+
+ Box3d box;
+ public Box3d Box
+ {
+ get { return box; }
+ set { box = value; Result = IntersectionResult.NotComputed; }
+ }
+
+ public int Quantity = 0;
+ public IntersectionResult Result = IntersectionResult.NotComputed;
+ public IntersectionType Type = IntersectionType.Empty;
+
+ public bool IsSimpleIntersection {
+ get { return Result == IntersectionResult.Intersects && Type == IntersectionType.Point; }
+ }
+
+ public double RayParam0, RayParam1;
+ public Vector3d Point0 = Vector3d.Zero;
+ public Vector3d Point1 = Vector3d.Zero;
+
+ public IntrRay3Box3(Ray3d r, Box3d b)
+ {
+ ray = r; box = b;
+ }
+
+ public IntrRay3Box3 Compute()
+ {
+ Find();
+ return this;
+ }
+
+
+ public bool Find()
+ {
+ if (Result != IntersectionResult.NotComputed)
+ return (Result == IntersectionResult.Intersects);
+
+ // [RMS] if either line direction is not a normalized vector,
+ // results are garbage, so fail query
+ if ( ray.Direction.IsNormalized == false ) {
+ Type = IntersectionType.Empty;
+ Result = IntersectionResult.InvalidQuery;
+ return false;
+ }
+
+ RayParam0 = 0.0;
+ RayParam1 = double.MaxValue;
+ IntrLine3Box3.DoClipping(ref RayParam0, ref RayParam1, ray.Origin, ray.Direction, box,
+ true, ref Quantity, ref Point0, ref Point1, ref Type);
+
+ Result = (Type != IntersectionType.Empty) ?
+ IntersectionResult.Intersects : IntersectionResult.NoIntersection;
+ return (Result == IntersectionResult.Intersects);
+ }
+
+
+
+
+ public bool Test ()
+ {
+ Vector3d WdU = Vector3d.Zero;
+ Vector3d AWdU = Vector3d.Zero;
+ Vector3d DdU = Vector3d.Zero;
+ Vector3d ADdU = Vector3d.Zero;
+ Vector3d AWxDdU = Vector3d.Zero;
+ double RHS;
+
+ Vector3d diff = ray.Origin - box.Center;
+
+ WdU[0] = ray.Direction.Dot(box.AxisX);
+ AWdU[0] = Math.Abs(WdU[0]);
+ DdU[0] = diff.Dot(box.AxisX);
+ ADdU[0] = Math.Abs(DdU[0]);
+ if (ADdU[0] > box.Extent.x && DdU[0]*WdU[0] >= (double)0)
+ {
+ return false;
+ }
+
+ WdU[1] = ray.Direction.Dot(box.AxisY);
+ AWdU[1] = Math.Abs(WdU[1]);
+ DdU[1] = diff.Dot(box.AxisY);
+ ADdU[1] = Math.Abs(DdU[1]);
+ if (ADdU[1] > box.Extent.y && DdU[1]*WdU[1] >= (double)0)
+ {
+ return false;
+ }
+
+ WdU[2] = ray.Direction.Dot(box.AxisZ);
+ AWdU[2] = Math.Abs(WdU[2]);
+ DdU[2] = diff.Dot(box.AxisZ);
+ ADdU[2] = Math.Abs(DdU[2]);
+ if (ADdU[2] > box.Extent.z && DdU[2]*WdU[2] >= (double)0)
+ {
+ return false;
+ }
+
+ Vector3d WxD = ray.Direction.Cross(diff);
+
+ AWxDdU[0] = Math.Abs(WxD.Dot(box.AxisX));
+ RHS = box.Extent.y*AWdU[2] + box.Extent.z*AWdU[1];
+ if (AWxDdU[0] > RHS)
+ {
+ return false;
+ }
+
+ AWxDdU[1] = Math.Abs(WxD.Dot(box.AxisY));
+ RHS = box.Extent.x*AWdU[2] + box.Extent.z*AWdU[0];
+ if (AWxDdU[1] > RHS)
+ {
+ return false;
+ }
+
+ AWxDdU[2] = Math.Abs(WxD.Dot(box.AxisZ));
+ RHS = box.Extent.x*AWdU[1] + box.Extent.y*AWdU[0];
+ if (AWxDdU[2] > RHS)
+ {
+ return false;
+ }
+
+ return true;
+ }
+
+
+
+ }
+}
diff --git a/intersection/IntrSegment2Triangle2.cs b/intersection/IntrSegment2Triangle2.cs
new file mode 100644
index 00000000..468685ed
--- /dev/null
+++ b/intersection/IntrSegment2Triangle2.cs
@@ -0,0 +1,108 @@
+using System;
+
+namespace g3
+{
+ // ported from WildMagic5
+ public class IntrSegment2Triangle2
+ {
+ Segment2d segment;
+ public Segment2d Segment
+ {
+ get { return segment; }
+ set { segment = value; Result = IntersectionResult.NotComputed; }
+ }
+
+ Triangle2d triangle;
+ public Triangle2d Triangle
+ {
+ get { return triangle; }
+ set { triangle = value; Result = IntersectionResult.NotComputed; }
+ }
+
+ public int Quantity = 0;
+ public IntersectionResult Result = IntersectionResult.NotComputed;
+ public IntersectionType Type = IntersectionType.Empty;
+
+ public bool IsSimpleIntersection {
+ get { return Result == IntersectionResult.Intersects && Type == IntersectionType.Point; }
+ }
+
+
+ public Vector2d Point0;
+ public Vector2d Point1;
+ public double Param0;
+ public double Param1;
+
+
+ public IntrSegment2Triangle2(Segment2d s, Triangle2d t)
+ {
+ segment = s; triangle = t;
+ }
+
+
+ public IntrSegment2Triangle2 Compute()
+ {
+ Find();
+ return this;
+ }
+
+
+ public bool Find()
+ {
+ if (Result != IntersectionResult.NotComputed)
+ return (Result == IntersectionResult.Intersects);
+
+ // [RMS] if either line direction is not a normalized vector,
+ // results are garbage, so fail query
+ if ( segment.Direction.IsNormalized == false ) {
+ Type = IntersectionType.Empty;
+ Result = IntersectionResult.InvalidQuery;
+ return false;
+ }
+
+ Vector3d dist = Vector3d.Zero;
+ Vector3i sign = Vector3i.Zero;
+ int positive = 0, negative = 0, zero = 0;
+ IntrLine2Triangle2.TriangleLineRelations(segment.Center, segment.Direction, triangle,
+ ref dist, ref sign, ref positive, ref negative, ref zero);
+
+ if (positive == 3 || negative == 3)
+ {
+ // No intersections.
+ Quantity = 0;
+ Type = IntersectionType.Empty;
+ } else {
+ Vector2d param = Vector2d.Zero;
+ IntrLine2Triangle2.GetInterval(segment.Center, segment.Direction, triangle, dist, sign, ref param);
+
+ Intersector1 intr = new Intersector1(param[0], param[1], -segment.Extent, +segment.Extent);
+ intr.Find();
+
+ Quantity = intr.NumIntersections;
+ if (Quantity == 2) {
+ // Segment intersection.
+ Type = IntersectionType.Segment;
+ Param0 = intr.GetIntersection(0);
+ Point0 = segment.Center + Param0*segment.Direction;
+ Param1 = intr.GetIntersection(1);
+ Point1 = segment.Center + Param1*segment.Direction;
+ } else if (Quantity == 1) {
+ // Point intersection.
+ Type = IntersectionType.Point;
+ Param0 = intr.GetIntersection(0);
+ Point0 = segment.Center + Param0*segment.Direction;
+ } else {
+ // No intersections.
+ Type = IntersectionType.Empty;
+ }
+ }
+
+ Result = (Type != IntersectionType.Empty) ?
+ IntersectionResult.Intersects : IntersectionResult.NoIntersection;
+ return (Result == IntersectionResult.Intersects);
+ }
+
+
+
+ }
+}
diff --git a/intersection/IntrSegment3Box3.cs b/intersection/IntrSegment3Box3.cs
new file mode 100644
index 00000000..af560ff2
--- /dev/null
+++ b/intersection/IntrSegment3Box3.cs
@@ -0,0 +1,231 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+
+namespace g3
+{
+ public class IntrSegment3Box3
+ {
+ Segment3d segment;
+ public Segment3d Segment
+ {
+ get { return segment; }
+ set { segment = value; Result = IntersectionResult.NotComputed; }
+ }
+
+ Box3d box;
+ public Box3d Box
+ {
+ get { return box; }
+ set { box = value; Result = IntersectionResult.NotComputed; }
+ }
+
+ bool solid = false;
+ public bool Solid {
+ get { return solid; }
+ set { solid = value; Result = IntersectionResult.NotComputed; }
+ }
+
+ public int Quantity = 0;
+ public IntersectionResult Result = IntersectionResult.NotComputed;
+ public IntersectionType Type = IntersectionType.Empty;
+
+ public bool IsSimpleIntersection {
+ get { return Result == IntersectionResult.Intersects && Type == IntersectionType.Point; }
+ }
+
+ public double SegmentParam0, SegmentParam1;
+ public Vector3d Point0 = Vector3d.Zero;
+ public Vector3d Point1 = Vector3d.Zero;
+
+ // solidBox == false means fully contained segment does not intersect
+ public IntrSegment3Box3(Segment3d s, Box3d b, bool solidBox)
+ {
+ segment = s; box = b; this.solid = solidBox;
+ }
+
+ public IntrSegment3Box3 Compute()
+ {
+ Find();
+ return this;
+ }
+
+
+ public bool Find()
+ {
+ if (Result != IntersectionResult.NotComputed)
+ return (Result == IntersectionResult.Intersects);
+
+ // [RMS] if either line direction is not a normalized vector,
+ // results are garbage, so fail query
+ if ( segment.Direction.IsNormalized == false ) {
+ Type = IntersectionType.Empty;
+ Result = IntersectionResult.InvalidQuery;
+ return false;
+ }
+
+ SegmentParam0 = -segment.Extent;
+ SegmentParam1 = segment.Extent;
+ DoClipping(ref SegmentParam0, ref SegmentParam1, segment.Center, segment.Direction, box,
+ solid, ref Quantity, ref Point0, ref Point1, ref Type);
+
+ Result = (Type != IntersectionType.Empty) ?
+ IntersectionResult.Intersects : IntersectionResult.NoIntersection;
+ return (Result == IntersectionResult.Intersects);
+ }
+
+
+
+
+ public bool Test ()
+ {
+ Vector3d AWdU = Vector3d.Zero;
+ Vector3d ADdU = Vector3d.Zero;
+ Vector3d AWxDdU = Vector3d.Zero;
+ double RHS;
+
+ Vector3d diff = segment.Center - box.Center;
+
+ AWdU[0] = Math.Abs(segment.Direction.Dot(box.AxisX));
+ ADdU[0] = Math.Abs(diff.Dot(box.AxisX));
+ RHS = box.Extent.x + segment.Extent*AWdU[0];
+ if (ADdU[0] > RHS)
+ {
+ return false;
+ }
+
+ AWdU[1] = Math.Abs(segment.Direction.Dot(box.AxisY));
+ ADdU[1] = Math.Abs(diff.Dot(box.AxisY));
+ RHS = box.Extent.y + segment.Extent*AWdU[1];
+ if (ADdU[1] > RHS)
+ {
+ return false;
+ }
+
+ AWdU[2] = Math.Abs(segment.Direction.Dot(box.AxisZ));
+ ADdU[2] = Math.Abs(diff.Dot(box.AxisZ));
+ RHS = box.Extent.z + segment.Extent*AWdU[2];
+ if (ADdU[2] > RHS)
+ {
+ return false;
+ }
+
+ Vector3d WxD = segment.Direction.Cross(diff);
+
+ AWxDdU[0] = Math.Abs(WxD.Dot(box.AxisX));
+ RHS = box.Extent.y*AWdU[2] + box.Extent.z*AWdU[1];
+ if (AWxDdU[0] > RHS)
+ {
+ return false;
+ }
+
+ AWxDdU[1] = Math.Abs(WxD.Dot(box.AxisY));
+ RHS = box.Extent.x*AWdU[2] + box.Extent.z*AWdU[0];
+ if (AWxDdU[1] > RHS)
+ {
+ return false;
+ }
+
+ AWxDdU[2] = Math.Abs(WxD.Dot(box.AxisZ));
+ RHS = box.Extent.x*AWdU[1] + box.Extent.y*AWdU[0];
+ if (AWxDdU[2] > RHS)
+ {
+ return false;
+ }
+
+ return true;
+ }
+
+
+
+
+ static public bool DoClipping (ref double t0, ref double t1,
+ Vector3d origin, Vector3d direction,
+ Box3d box, bool solid, ref int quantity,
+ ref Vector3d point0, ref Vector3d point1,
+ ref IntersectionType intrType)
+ {
+ // Convert linear component to box coordinates.
+ Vector3d diff = origin - box.Center;
+ Vector3d BOrigin = new Vector3d(
+ diff.Dot(box.AxisX),
+ diff.Dot(box.AxisY),
+ diff.Dot(box.AxisZ)
+ );
+ Vector3d BDirection = new Vector3d(
+ direction.Dot(box.AxisX),
+ direction.Dot(box.AxisY),
+ direction.Dot(box.AxisZ)
+ );
+
+ double saveT0 = t0, saveT1 = t1;
+ bool notAllClipped =
+ Clip(+BDirection.x, -BOrigin.x-box.Extent.x, ref t0, ref t1) &&
+ Clip(-BDirection.x, +BOrigin.x-box.Extent.x, ref t0, ref t1) &&
+ Clip(+BDirection.y, -BOrigin.y-box.Extent.y, ref t0, ref t1) &&
+ Clip(-BDirection.y, +BOrigin.y-box.Extent.y, ref t0, ref t1) &&
+ Clip(+BDirection.z, -BOrigin.z-box.Extent.z, ref t0, ref t1) &&
+ Clip(-BDirection.z, +BOrigin.z-box.Extent.z, ref t0, ref t1);
+
+ if (notAllClipped && (solid || t0 != saveT0 || t1 != saveT1)) {
+ if (t1 > t0) {
+ intrType = IntersectionType.Segment;
+ quantity = 2;
+ point0 = origin + t0*direction;
+ point1 = origin + t1*direction;
+ } else {
+ intrType = IntersectionType.Point;
+ quantity = 1;
+ point0 = origin + t0*direction;
+ }
+ } else {
+ quantity = 0;
+ intrType = IntersectionType.Empty;
+ }
+
+ return intrType != IntersectionType.Empty;
+ }
+
+
+
+
+ static public bool Clip (double denom, double numer, ref double t0, ref double t1)
+ {
+ // Return value is 'true' if line segment intersects the current test
+ // plane. Otherwise 'false' is returned in which case the line segment
+ // is entirely clipped.
+
+ if (denom > (double)0)
+ {
+ if (numer > denom*t1)
+ {
+ return false;
+ }
+ if (numer > denom*t0)
+ {
+ t0 = numer/denom;
+ }
+ return true;
+ }
+ else if (denom < (double)0)
+ {
+ if (numer > denom*t0)
+ {
+ return false;
+ }
+ if (numer > denom*t1)
+ {
+ t1 = numer/denom;
+ }
+ return true;
+ }
+ else
+ {
+ return numer <= (double)0;
+ }
+ }
+
+
+ }
+}
diff --git a/intersection/IntrTriangle2Triangle2.cs b/intersection/IntrTriangle2Triangle2.cs
new file mode 100644
index 00000000..7afaf72f
--- /dev/null
+++ b/intersection/IntrTriangle2Triangle2.cs
@@ -0,0 +1,275 @@
+using System;
+
+namespace g3
+{
+ // ported from WildMagic5
+ public class IntrTriangleTriangle2
+ {
+ Triangle2d triangle0;
+ public Triangle2d Triangle0
+ {
+ get { return triangle0; }
+ set { triangle0 = value; Result = IntersectionResult.NotComputed; }
+ }
+
+ Triangle2d triangle1;
+ public Triangle2d Triangle1
+ {
+ get { return triangle1; }
+ set { triangle1 = value; Result = IntersectionResult.NotComputed; }
+ }
+
+
+ public int Quantity = 0;
+ public IntersectionResult Result = IntersectionResult.NotComputed;
+ public IntersectionType Type = IntersectionType.Empty;
+
+ public bool IsSimpleIntersection {
+ get { return Result == IntersectionResult.Intersects && Type == IntersectionType.Point; }
+ }
+
+ // intersection polygon - this array will always be 6 elements long,
+ // however only the first Quantity vertices will be valid
+ public Vector2d[] Points;
+
+ public IntrTriangleTriangle2(Triangle2d t0, Triangle2d t1)
+ {
+ triangle0 = t0;
+ triangle1 = t1;
+ Points = null;
+ }
+
+
+ public bool Test ()
+ {
+ int i0, i1;
+ Vector2d dir = Vector2d.Zero;
+
+ // Test edges of triangle0 for separation.
+ for (i0 = 0, i1 = 2; i0 < 3; i1 = i0++)
+ {
+ // Test axis V0[i1] + t*perp(V0[i0]-V0[i1]), perp(x,y) = (y,-x).
+ dir.x = triangle0[i0].y - triangle0[i1].y;
+ dir.y = triangle0[i1].x - triangle0[i0].x;
+ if (WhichSide(triangle1, triangle0[i1], dir) > 0) {
+ // Triangle1 is entirely on positive side of triangle0 edge.
+ return false;
+ }
+ }
+
+ // Test edges of triangle1 for separation.
+ for (i0 = 0, i1 = 2; i0 < 3; i1 = i0++)
+ {
+ // Test axis V1[i1] + t*perp(V1[i0]-V1[i1]), perp(x,y) = (y,-x).
+ dir.x = triangle1[i0].y - triangle1[i1].y;
+ dir.y = triangle1[i1].x - triangle1[i0].x;
+ if (WhichSide(triangle0, triangle1[i1], dir) > 0) {
+ // Triangle0 is entirely on positive side of triangle1 edge.
+ return false;
+ }
+ }
+
+ return true;
+ }
+
+
+
+
+
+ public IntrTriangleTriangle2 Compute()
+ {
+ Find();
+ return this;
+ }
+
+
+ public bool Find()
+ {
+ if (Result != IntersectionResult.NotComputed)
+ return (Result == IntersectionResult.Intersects);
+
+ // The potential intersection is initialized to triangle1. The set of
+ // vertices is refined based on clipping against each edge of triangle0.
+ Quantity = 3;
+ Points = new Vector2d[6];
+ for (int i = 0; i < 3; ++i) {
+ Points[i] = triangle1[i];
+ }
+
+ for (int i1 = 2, i0 = 0; i0 < 3; i1 = i0++)
+ {
+ // Clip against edge .
+ Vector2d N = new Vector2d(
+ triangle0[i1].y - triangle0[i0].y,
+ triangle0[i0].x - triangle0[i1].x);
+ double c = N.Dot(triangle0[i1]);
+ ClipConvexPolygonAgainstLine(N, c, ref Quantity, ref Points);
+ if (Quantity == 0) {
+ // Triangle completely clipped, no intersection occurs.
+ Type = IntersectionType.Empty;
+ } else if ( Quantity == 1 ) {
+ Type = IntersectionType.Point;
+ } else if ( Quantity == 2 ) {
+ Type = IntersectionType.Segment;
+ } else {
+ Type = IntersectionType.Polygon;
+ }
+ }
+
+ Result = (Type != IntersectionType.Empty) ?
+ IntersectionResult.Intersects : IntersectionResult.NoIntersection;
+ return (Result == IntersectionResult.Intersects);
+ }
+
+
+
+
+ public static int WhichSide (Triangle2d V, Vector2d P, Vector2d D)
+ {
+ // Vertices are projected to the form P+t*D. Return value is +1 if all
+ // t > 0, -1 if all t < 0, 0 otherwise, in which case the line splits the
+ // triangle.
+
+ int positive = 0, negative = 0, zero = 0;
+ for (int i = 0; i < 3; ++i)
+ {
+ double t = D.Dot(V[i] - P);
+ if (t > (double)0)
+ {
+ ++positive;
+ }
+ else if (t < (double)0)
+ {
+ ++negative;
+ }
+ else
+ {
+ ++zero;
+ }
+
+ if (positive > 0 && negative > 0)
+ {
+ return 0;
+ }
+ }
+ return (zero == 0 ? (positive > 0 ? 1 : -1) : 0);
+ }
+
+
+
+ // Vin is input polygon vertices, returns clipped polygon, vertex count of
+ // clipped polygon is returned in quantity
+ // **NOTE** returned array may have more elements than quantity!!
+ public static void ClipConvexPolygonAgainstLine (
+ Vector2d N, double c, ref int quantity, ref Vector2d[] V)
+ {
+ // The input vertices are assumed to be in counterclockwise order. The
+ // ordering is an invariant of this function.
+
+ // Test on which side of line the vertices are.
+ int positive = 0, negative = 0, pIndex = -1;
+ double[] test = new double[6];
+ int i;
+ for (i = 0; i < quantity; ++i)
+ {
+ test[i] = N.Dot(V[i]) - c;
+ if (test[i] > (double)0)
+ {
+ positive++;
+ if (pIndex < 0)
+ {
+ pIndex = i;
+ }
+ }
+ else if (test[i] < (double)0)
+ {
+ negative++;
+ }
+ }
+
+ if (positive > 0)
+ {
+ if (negative > 0)
+ {
+ // Line transversely intersects polygon.
+ Vector2d[] CV = new Vector2d[6];
+ int cQuantity = 0, cur, prv;
+ double t;
+
+ if (pIndex > 0) {
+ // First clip vertex on line.
+ cur = pIndex;
+ prv = cur - 1;
+ t = test[cur]/(test[cur] - test[prv]);
+ CV[cQuantity++] = V[cur] + t*(V[prv] - V[cur]);
+
+ // Vertices on positive side of line.
+ while (cur < quantity && test[cur] > (double)0) {
+ CV[cQuantity++] = V[cur++];
+ }
+
+ // Last clip vertex on line.
+ if (cur < quantity) {
+ prv = cur - 1;
+ } else {
+ cur = 0;
+ prv = quantity - 1;
+ }
+ t = test[cur]/(test[cur] - test[prv]);
+ CV[cQuantity++] = V[cur] + t*(V[prv]-V[cur]);
+ }
+ else // pIndex is 0
+ {
+ // Vertices on positive side of line.
+ cur = 0;
+ while (cur < quantity && test[cur] > (double)0)
+ {
+ CV[cQuantity++] = V[cur++];
+ }
+
+ // Last clip vertex on line.
+ prv = cur - 1;
+ t = test[cur]/(test[cur] - test[prv]);
+ CV[cQuantity++] = V[cur] + t*(V[prv] - V[cur]);
+
+ // Skip vertices on negative side.
+ while (cur < quantity && test[cur] <= (double)0)
+ {
+ ++cur;
+ }
+
+ // First clip vertex on line.
+ if (cur < quantity)
+ {
+ prv = cur - 1;
+ t = test[cur]/(test[cur] - test[prv]);
+ CV[cQuantity++] = V[cur] + t*(V[prv] - V[cur]);
+
+ // Vertices on positive side of line.
+ while (cur < quantity && test[cur] > (double)0)
+ {
+ CV[cQuantity++] = V[cur++];
+ }
+ }
+ else
+ {
+ // cur = 0
+ prv = quantity - 1;
+ t = test[0]/(test[0] - test[prv]);
+ CV[cQuantity++] = V[0] + t*(V[prv] - V[0]);
+ }
+ }
+
+ quantity = cQuantity;
+ Array.Copy(CV, V, cQuantity);
+ }
+ // else polygon fully on positive side of line, nothing to do.
+ } else {
+ // Polygon does not intersect positive side of line, clip all.
+ quantity = 0;
+ }
+ }
+
+
+ }
+}
diff --git a/math/Box2.cs b/math/Box2.cs
new file mode 100644
index 00000000..2127790c
--- /dev/null
+++ b/math/Box2.cs
@@ -0,0 +1,239 @@
+using System;
+
+namespace g3 {
+
+ // partially based on WildMagic5 Box3
+ public struct Box2d
+ {
+ // A box has center C, axis directions U[0] and U[1] (perpendicular and
+ // unit-length vectors), and extents e[0] and e[1] (nonnegative numbers).
+ // A/ point X = C+y[0]*U[0]+y[1]*U[1] is inside or on the box whenever
+ // |y[i]| <= e[i] for all i.
+
+ public Vector2d Center;
+ public Vector2d AxisX;
+ public Vector2d AxisY;
+ public Vector2d Extent;
+
+ public Box2d(Vector2d center) {
+ Center = center;
+ AxisX = Vector2d.AxisX;
+ AxisY = Vector2d.AxisY;
+ Extent = Vector2d.Zero;
+ }
+ public Box2d(Vector2d center, Vector2d x, Vector2d y, Vector2d extent) {
+ Center = center;
+ AxisX = x; AxisY = y;
+ Extent = extent;
+ }
+ public Box2d(AxisAlignedBox2d aaBox) {
+ Extent= 0.5*aaBox.Diagonal;
+ Center = aaBox.Min + Extent;
+ AxisX = Vector2d.AxisX;
+ AxisY = Vector2d.AxisY;
+ }
+
+ public static readonly Box2d Empty = new Box2d(Vector2d.Zero);
+
+
+ public Vector2d Axis(int i)
+ {
+ return (i == 0) ? AxisX : AxisY;
+ }
+
+
+ public Vector2d[] ComputeVertices() {
+ Vector2d[] v = new Vector2d[4];
+ ComputeVertices(v);
+ return v;
+ }
+ public void ComputeVertices (Vector2d[] vertex) {
+ Vector2d extAxis0 = Extent.x*AxisX;
+ Vector2d extAxis1 = Extent.y*AxisY;
+ vertex[0] = Center - extAxis0 - extAxis1;
+ vertex[1] = Center + extAxis0 - extAxis1;
+ vertex[2] = Center + extAxis0 + extAxis1;
+ vertex[3] = Center - extAxis0 + extAxis1;
+ }
+
+
+ // g3 extensions
+ public double MaxExtent {
+ get { return Math.Max(Extent.x, Extent.y); }
+ }
+ public double MinExtent {
+ get { return Math.Min(Extent.x, Extent.y); }
+ }
+ public Vector2d Diagonal {
+ get { return (Extent.x*AxisX + Extent.y*AxisY) -
+ (-Extent.x*AxisX - Extent.y*AxisY); }
+ }
+ public double Area {
+ get { return 2*Extent.x + 2*Extent.y; }
+ }
+
+ public void Contain( Vector2d v) {
+ Vector2d lv = v - Center;
+ for (int k = 0; k < 2; ++k) {
+ double t = lv.Dot(Axis(k));
+ if ( Math.Abs(t) > Extent[k]) {
+ double min = -Extent[k], max = Extent[k];
+ if ( t < min )
+ min = t;
+ else if ( t > max )
+ max = t;
+ Extent[k] = (max-min) * 0.5;
+ Center = Center + ((max+min) * 0.5) * Axis(k);
+ }
+ }
+ }
+
+ // I think this can be more efficient...no? At least could combine
+ // all the axis-interval updates before updating Center...
+ public void Contain( Box2d o ) {
+ Vector2d[] v = o.ComputeVertices();
+ for (int k = 0; k < 4; ++k)
+ Contain(v[k]);
+ }
+
+ public bool Contained( Vector2d v ) {
+ Vector2d lv = v - Center;
+ return (Math.Abs(lv.Dot(AxisX)) <= Extent.x) &&
+ (Math.Abs(lv.Dot(AxisY)) <= Extent.y);
+ }
+
+ public void Expand(double f) {
+ Extent += f;
+ }
+
+ public void Translate( Vector2d v ) {
+ Center += v;
+ }
+
+ }
+
+
+
+
+
+
+
+
+
+
+
+
+ // partially based on WildMagic5 Box3
+ public struct Box2f
+ {
+ // A box has center C, axis directions U[0] and U[1] (perpendicular and
+ // unit-length vectors), and extents e[0] and e[1] (nonnegative numbers).
+ // A/ point X = C+y[0]*U[0]+y[1]*U[1] is inside or on the box whenever
+ // |y[i]| <= e[i] for all i.
+
+ public Vector2f Center;
+ public Vector2f AxisX;
+ public Vector2f AxisY;
+ public Vector2f Extent;
+
+ public Box2f(Vector2f center) {
+ Center = center;
+ AxisX = Vector2f.AxisX;
+ AxisY = Vector2f.AxisY;
+ Extent = Vector2f.Zero;
+ }
+ public Box2f(Vector2f center, Vector2f x, Vector2f y, Vector2f extent) {
+ Center = center;
+ AxisX = x; AxisY = y;
+ Extent = extent;
+ }
+ public Box2f(AxisAlignedBox2f aaBox) {
+ Extent= 0.5f*aaBox.Diagonal;
+ Center = aaBox.Min + Extent;
+ AxisX = Vector2f.AxisX;
+ AxisY = Vector2f.AxisY;
+ }
+
+ public static readonly Box2f Empty = new Box2f(Vector2f.Zero);
+
+
+ public Vector2f Axis(int i)
+ {
+ return (i == 0) ? AxisX : AxisY;
+ }
+
+
+ public Vector2f[] ComputeVertices() {
+ Vector2f[] v = new Vector2f[4];
+ ComputeVertices(v);
+ return v;
+ }
+ public void ComputeVertices (Vector2f[] vertex) {
+ Vector2f extAxis0 = Extent.x*AxisX;
+ Vector2f extAxis1 = Extent.y*AxisY;
+ vertex[0] = Center - extAxis0 - extAxis1;
+ vertex[1] = Center + extAxis0 - extAxis1;
+ vertex[2] = Center + extAxis0 + extAxis1;
+ vertex[3] = Center - extAxis0 + extAxis1;
+ }
+
+
+ // g3 extensions
+ public double MaxExtent {
+ get { return Math.Max(Extent.x, Extent.y); }
+ }
+ public double MinExtent {
+ get { return Math.Min(Extent.x, Extent.y); }
+ }
+ public Vector2f Diagonal {
+ get { return (Extent.x*AxisX + Extent.y*AxisY) -
+ (-Extent.x*AxisX - Extent.y*AxisY); }
+ }
+ public double Area {
+ get { return 2*Extent.x + 2*Extent.y; }
+ }
+
+ public void Contain( Vector2f v) {
+ Vector2f lv = v - Center;
+ for (int k = 0; k < 2; ++k) {
+ double t = lv.Dot(Axis(k));
+ if ( Math.Abs(t) > Extent[k]) {
+ double min = -Extent[k], max = Extent[k];
+ if ( t < min )
+ min = t;
+ else if ( t > max )
+ max = t;
+ Extent[k] = (float)(max-min) * 0.5f;
+ Center = Center + ((float)(max+min) * 0.5f) * Axis(k);
+ }
+ }
+ }
+
+ // I think this can be more efficient...no? At least could combine
+ // all the axis-interval updates before updating Center...
+ public void Contain( Box2f o ) {
+ Vector2f[] v = o.ComputeVertices();
+ for (int k = 0; k < 4; ++k)
+ Contain(v[k]);
+ }
+
+ public bool Contained( Vector2f v ) {
+ Vector2f lv = v - Center;
+ return (Math.Abs(lv.Dot(AxisX)) <= Extent.x) &&
+ (Math.Abs(lv.Dot(AxisY)) <= Extent.y);
+ }
+
+ public void Expand(float f) {
+ Extent += f;
+ }
+
+ public void Translate( Vector2f v ) {
+ Center += v;
+ }
+
+ }
+
+
+
+
+}
diff --git a/math/Box3.cs b/math/Box3.cs
new file mode 100644
index 00000000..65d96358
--- /dev/null
+++ b/math/Box3.cs
@@ -0,0 +1,259 @@
+using System;
+
+namespace g3 {
+
+ // partially based on WildMagic5 Box3
+ public struct Box3d
+ {
+ // A box has center C, axis directions U[0], U[1], and U[2] (mutually
+ // perpendicular unit-length vectors), and extents e[0], e[1], and e[2]
+ // (all nonnegative numbers). A point X = C+y[0]*U[0]+y[1]*U[1]+y[2]*U[2]
+ // is inside or on the box whenever |y[i]| <= e[i] for all i.
+
+ public Vector3d Center;
+ public Vector3d AxisX;
+ public Vector3d AxisY;
+ public Vector3d AxisZ;
+ public Vector3d Extent;
+
+ public Box3d(Vector3d center) {
+ Center = center;
+ AxisX = Vector3d.AxisX;
+ AxisY = Vector3d.AxisY;
+ AxisZ = Vector3d.AxisZ;
+ Extent = Vector3d.Zero;
+ }
+ public Box3d(Vector3d center, Vector3d x, Vector3d y, Vector3d z,
+ Vector3d extent) {
+ Center = center;
+ AxisX = x; AxisY = y; AxisZ = z;
+ Extent = extent;
+ }
+ public Box3d(AxisAlignedBox3d aaBox) {
+ Extent= 0.5*aaBox.Diagonal;
+ Center = aaBox.Min + Extent;
+ AxisX = Vector3d.AxisX;
+ AxisY = Vector3d.AxisY;
+ AxisZ = Vector3d.AxisZ;
+ }
+
+ public static readonly Box3d Empty = new Box3d(Vector3d.Zero);
+
+
+ public Vector3d Axis(int i)
+ {
+ return (i == 0) ? AxisX : (i == 1) ? AxisY : AxisZ;
+ }
+
+
+ public Vector3d[] ComputeVertices() {
+ Vector3d[] v = new Vector3d[8];
+ ComputeVertices(v);
+ return v;
+ }
+ public void ComputeVertices (Vector3d[] vertex) {
+ Vector3d extAxis0 = Extent.x*AxisX;
+ Vector3d extAxis1 = Extent.y*AxisY;
+ Vector3d extAxis2 = Extent.z*AxisZ;
+ vertex[0] = Center - extAxis0 - extAxis1 - extAxis2;
+ vertex[1] = Center + extAxis0 - extAxis1 - extAxis2;
+ vertex[2] = Center + extAxis0 + extAxis1 - extAxis2;
+ vertex[3] = Center - extAxis0 + extAxis1 - extAxis2;
+ vertex[4] = Center - extAxis0 - extAxis1 + extAxis2;
+ vertex[5] = Center + extAxis0 - extAxis1 + extAxis2;
+ vertex[6] = Center + extAxis0 + extAxis1 + extAxis2;
+ vertex[7] = Center - extAxis0 + extAxis1 + extAxis2;
+ }
+
+
+ // g3 extensions
+ public double MaxExtent {
+ get { return Math.Max(Extent.x, Math.Max(Extent.y, Extent.z) ); }
+ }
+ public double MinExtent {
+ get { return Math.Min(Extent.x, Math.Max(Extent.y, Extent.z) ); }
+ }
+ public Vector3d Diagonal {
+ get { return (Extent.x*AxisX + Extent.y*AxisY + Extent.z*AxisZ) -
+ (-Extent.x*AxisX - Extent.y*AxisY - Extent.z*AxisZ); }
+ }
+ public double Volume {
+ get { return 2*Extent.x + 2*Extent.y * 2*Extent.z; }
+ }
+
+ public void Contain( Vector3d v) {
+ Vector3d lv = v - Center;
+ for (int k = 0; k < 3; ++k) {
+ double t = lv.Dot(Axis(k));
+ if ( Math.Abs(t) > Extent[k]) {
+ double min = -Extent[k], max = Extent[k];
+ if ( t < min )
+ min = t;
+ else if ( t > max )
+ max = t;
+ Extent[k] = (max-min) * 0.5;
+ Center = Center + ((max+min) * 0.5) * Axis(k);
+ }
+ }
+ }
+
+ // I think this can be more efficient...no? At least could combine
+ // all the axis-interval updates before updating Center...
+ public void Contain( Box3d o ) {
+ Vector3d[] v = o.ComputeVertices();
+ for (int k = 0; k < 8; ++k)
+ Contain(v[k]);
+ }
+
+ public bool Contained( Vector3d v ) {
+ Vector3d lv = v - Center;
+ return (Math.Abs(lv.Dot(AxisX)) <= Extent.x) &&
+ (Math.Abs(lv.Dot(AxisY)) <= Extent.y) &&
+ (Math.Abs(lv.Dot(AxisZ)) <= Extent.z);
+ }
+
+ public void Expand(double f) {
+ Extent += f;
+ }
+
+ public void Translate( Vector3d v ) {
+ Center += v;
+ }
+
+ }
+
+
+
+
+
+
+
+
+
+
+
+
+ // partially based on WildMagic5 Box3
+ public struct Box3f
+ {
+ // A box has center C, axis directions U[0], U[1], and U[2] (mutually
+ // perpendicular unit-length vectors), and extents e[0], e[1], and e[2]
+ // (all nonnegative numbers). A point X = C+y[0]*U[0]+y[1]*U[1]+y[2]*U[2]
+ // is inside or on the box whenever |y[i]| <= e[i] for all i.
+
+ public Vector3f Center;
+ public Vector3f AxisX;
+ public Vector3f AxisY;
+ public Vector3f AxisZ;
+ public Vector3f Extent;
+
+ public Box3f(Vector3f center) {
+ Center = center;
+ AxisX = Vector3f.AxisX;
+ AxisY = Vector3f.AxisY;
+ AxisZ = Vector3f.AxisZ;
+ Extent = Vector3f.Zero;
+ }
+ public Box3f(Vector3f center, Vector3f x, Vector3f y, Vector3f z,
+ Vector3f extent) {
+ Center = center;
+ AxisX = x; AxisY = y; AxisZ = z;
+ Extent = extent;
+ }
+ public Box3f(AxisAlignedBox3f aaBox) {
+ Extent= 0.5f*aaBox.Diagonal;
+ Center = aaBox.Min + Extent;
+ AxisX = Vector3f.AxisX;
+ AxisY = Vector3f.AxisY;
+ AxisZ = Vector3f.AxisZ;
+ }
+
+ public static readonly Box3f Empty = new Box3f(Vector3f.Zero);
+
+
+ public Vector3f Axis(int i)
+ {
+ return (i == 0) ? AxisX : (i == 1) ? AxisY : AxisZ;
+ }
+
+
+ public Vector3f[] ComputeVertices() {
+ Vector3f[] v = new Vector3f[8];
+ ComputeVertices(v);
+ return v;
+ }
+ public void ComputeVertices (Vector3f[] vertex) {
+ Vector3f extAxis0 = Extent.x*AxisX;
+ Vector3f extAxis1 = Extent.y*AxisY;
+ Vector3f extAxis2 = Extent.z*AxisZ;
+ vertex[0] = Center - extAxis0 - extAxis1 - extAxis2;
+ vertex[1] = Center + extAxis0 - extAxis1 - extAxis2;
+ vertex[2] = Center + extAxis0 + extAxis1 - extAxis2;
+ vertex[3] = Center - extAxis0 + extAxis1 - extAxis2;
+ vertex[4] = Center - extAxis0 - extAxis1 + extAxis2;
+ vertex[5] = Center + extAxis0 - extAxis1 + extAxis2;
+ vertex[6] = Center + extAxis0 + extAxis1 + extAxis2;
+ vertex[7] = Center - extAxis0 + extAxis1 + extAxis2;
+ }
+
+
+ // g3 extensions
+ public double MaxExtent {
+ get { return Math.Max(Extent.x, Math.Max(Extent.y, Extent.z) ); }
+ }
+ public double MinExtent {
+ get { return Math.Min(Extent.x, Math.Max(Extent.y, Extent.z) ); }
+ }
+ public Vector3f Diagonal {
+ get { return (Extent.x*AxisX + Extent.y*AxisY + Extent.z*AxisZ) -
+ (-Extent.x*AxisX - Extent.y*AxisY - Extent.z*AxisZ); }
+ }
+ public double Volume {
+ get { return 2*Extent.x + 2*Extent.y * 2*Extent.z; }
+ }
+
+ public void Contain( Vector3f v) {
+ Vector3f lv = v - Center;
+ for (int k = 0; k < 3; ++k) {
+ double t = lv.Dot(Axis(k));
+ if ( Math.Abs(t) > Extent[k]) {
+ double min = -Extent[k], max = Extent[k];
+ if ( t < min )
+ min = t;
+ else if ( t > max )
+ max = t;
+ Extent[k] = (float)(max-min) * 0.5f;
+ Center = Center + ((float)(max+min) * 0.5f) * Axis(k);
+ }
+ }
+ }
+
+ // I think this can be more efficient...no? At least could combine
+ // all the axis-interval updates before updating Center...
+ public void Contain( Box3f o ) {
+ Vector3f[] v = o.ComputeVertices();
+ for (int k = 0; k < 8; ++k)
+ Contain(v[k]);
+ }
+
+ public bool Contained( Vector3f v ) {
+ Vector3f lv = v - Center;
+ return (Math.Abs(lv.Dot(AxisX)) <= Extent.x) &&
+ (Math.Abs(lv.Dot(AxisY)) <= Extent.y) &&
+ (Math.Abs(lv.Dot(AxisZ)) <= Extent.z);
+ }
+
+ public void Expand(float f) {
+ Extent += f;
+ }
+
+ public void Translate( Vector3f v ) {
+ Center += v;
+ }
+
+ }
+
+
+
+
+}
diff --git a/math/Triangle2.cs b/math/Triangle2.cs
new file mode 100644
index 00000000..0b6870d3
--- /dev/null
+++ b/math/Triangle2.cs
@@ -0,0 +1,71 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+
+namespace g3
+{
+ public struct Triangle2d
+ {
+ public Vector2d V0, V1, V2;
+
+ public Triangle2d(Vector2d v0, Vector2d v1, Vector2d v2)
+ {
+ V0 = v0; V1 = v1; V2 = v2;
+ }
+
+ public Vector2d this[int key]
+ {
+ get { return (key == 0) ? V0 : (key == 1) ? V1 : V2; }
+ set { if (key == 0) V0 = value; else if (key == 1) V1 = value; else V2 = value; }
+ }
+
+ public Vector2d PointAt(double bary0, double bary1, double bary2)
+ {
+ return bary0 * V0 + bary1 * V1 + bary2 * V2;
+ }
+ public Vector2d PointAt(Vector3d bary)
+ {
+ return bary.x* V0 + bary.y* V1 + bary.z* V2;
+ }
+
+ // conversion operators
+ public static implicit operator Triangle2d(Triangle2f v)
+ {
+ return new Triangle2d(v.V0, v.V1, v.V2);
+ }
+ public static explicit operator Triangle2f(Triangle2d v)
+ {
+ return new Triangle2f((Vector2f)v.V0, (Vector2f)v.V1, (Vector2f)v.V2);
+ }
+ }
+
+
+
+ public struct Triangle2f
+ {
+ public Vector2f V0, V1, V2;
+
+ public Triangle2f(Vector2f v0, Vector2f v1, Vector2f v2)
+ {
+ V0 = v0; V1 = v1; V2 = v2;
+ }
+
+ public Vector2f this[int key]
+ {
+ get { return (key == 0) ? V0 : (key == 1) ? V1 : V2; }
+ set { if (key == 0) V0 = value; else if (key == 1) V1 = value; else V2 = value; }
+ }
+
+
+ public Vector2f PointAt(float bary0, float bary1, float bary2)
+ {
+ return bary0 * V0 + bary1 * V1 + bary2 * V2;
+ }
+ public Vector2f PointAt(Vector3f bary)
+ {
+ return bary.x * V0 + bary.y * V1 + bary.z * V2;
+ }
+ }
+
+}
diff --git a/queries/Intersection.cs b/queries/Intersection.cs
index 2dabcd30..90f62191 100644
--- a/queries/Intersection.cs
+++ b/queries/Intersection.cs
@@ -15,7 +15,7 @@ public enum IntersectionResult
public enum IntersectionType
{
- Empty, Point, Segment, Line
+ Empty, Point, Segment, Line, Polygon
}
diff --git a/queries/IntrRay3Triangle3.cs b/queries/IntrRay3Triangle3.cs
index dcd4d033..e303d5b7 100644
--- a/queries/IntrRay3Triangle3.cs
+++ b/queries/IntrRay3Triangle3.cs
@@ -21,8 +21,9 @@ public Triangle3d Triangle
set { triangle = value; Result = IntersectionResult.NotComputed; }
}
- public IntersectionResult Result;
- public IntersectionType Type;
+ public int Quantity = 0;
+ public IntersectionResult Result = IntersectionResult.NotComputed;
+ public IntersectionType Type = IntersectionType.Empty;
public bool IsSimpleIntersection {
get { return Result == IntersectionResult.Intersects && Type == IntersectionType.Point; }
@@ -33,6 +34,12 @@ public bool IsSimpleIntersection {
public Vector3d TriangleBaryCoords;
+ public IntrRay3Triangle3(Ray3d r, Triangle3d t)
+ {
+ ray = r; triangle = t;
+ }
+
+
public IntrRay3Triangle3 Compute()
{
Find();
@@ -85,6 +92,7 @@ public bool Find()
double mTriBary2 = DdE1xQ * inv;
TriangleBaryCoords = new Vector3d(1 - mTriBary1 - mTriBary2, mTriBary1, mTriBary2);
Type = IntersectionType.Point;
+ Quantity = 1;
Result = IntersectionResult.Intersects;
return true;
}