BoundingSphere ,RayIntersection implementation ...
RuiTheAlmeida wrote at 2012-11-29 11:48:
Hi Objo and everybody,
First of all I want to congratulate you for this great Helix 3D tookit
I had to use one of the methods of class BoundingSphere and as it was not implemented so I did it. Feel free to test and include the code I'm sending here.
Best Regards,
Rui
public bool RayIntersection(Ray3D ray, out Point3D[] result)
{
result = null;
bool inter = false;
double h = this.center.X ;
double j = this.center.Y;
double k = this.center.Z;
double ra = this.radius;
double x1 = ray.Origin.X ;
double y1 = ray.Origin.Y;
double z1 = ray.Origin.Z;
double a = x1 - ray.Origin.X * ray.Direction.X;
double b = y1 - ray.Origin.Y * ray.Direction.Y;
double c = z1 - ray.Origin.Z * ray.Direction.Z;
//Quadratic solving
double aQ = a * a + b * b + c * c;
double bQ = 2 * a * (x1 - h) + 2 * b * (y1 - j) + 2 * c * (z1 - k);
double cQ = x1 * x1 + y1 * y1 + z1 * z1 + h * h + k * k + j * j - 2 * (j * y1 + k * z1 + h * x1) - ra * ra;
//Dscriminant
double Q = bQ * bQ - 4 * aQ * cQ;
if (Q >= 0)
{
double rx1 = (-bQ + Math.Sqrt(bQ * bQ - 4 * aQ * cQ)) / (2 * aQ); //First Root
double rx2 = (-bQ - Math.Sqrt(bQ * bQ - 4 * aQ * cQ)) / (2 * aQ); //Second Root
//first Intersection
Point3D I1 = new Point3D(x1 + a * rx1,y1 + b * rx1 , z1 + c * rx1);
//Second Intersection
Point3D I2 = new Point3D(x1 + a * rx2, y1 + b * rx2, z1 + c * rx2);
//ray inside sphere
if ((ray.Origin - this.center).Length <= this.radius)
{
result = new Point3D[1];
//Check if the vectors have the same direction
if (Vector3D.DotProduct(ray.Direction, I1 - ray.Origin) >= 0)
{
result[0] = I1;
}
else
{
result[0] = I2;
}
}
else
{
result = new Point3D[2];
result[0] = I1;
result[1] = I2;
}
inter = true;
}
else
{
inter = false;
}
return inter;
}
objo wrote at 2012-12-04 19:45:
thanks for the code! Did you make this yourself, or do you have a reference?
I created some unit tests for this code, but could not get my last case (no intersection) to work:
[Test, Ignore]
public void RayInterSection_NoIntersections()
{
var ray = new Ray3D(new Point3D(0.2, 0.3, 1), new Vector3D(1, 0.1, 0.1));
var sphere = new BoundingSphere(new Point3D(0.2, 0.3, 0), 1.0);
Point3D[] result;
Assert.IsFalse(sphere.RayIntersection(ray, out result));
Assert.AreEqual(sphere.Radius, sphere.Center.DistanceTo(result[0]), 1e-6, "Point is not on sphere.");
Assert.AreEqual(0, ray.GetNearest(result[0]).DistanceTo(result[0]), 1e-6, "Point is not on ray.");
}
RuiTheAlmeida wrote at 2012-12-07 05:12:
Hi objo,
I did the all code myself, and yes I have some references from Wikipedia and other common links that I will give.
I'm also finishing a detailed presentation with the all theoretical stuff, and I will post it all here.
Regards,
Rui
objo wrote at 2012-12-07 08:01:
Good! Can you also look at the "no intersections" case? Is there an error in my test?
I found some efficient implementations at
- http://wiki.cgsociety.org/index.php/Ray_Sphere_Intersection
- http://www.siggraph.org/education/materials/HyperGraph/raytrace/rtinter1.htm
- http://paulbourke.net/geometry/circlesphere/
- http://www.lighthouse3d.com/tutorials/maths/ray-sphere-intersection/
- http://www.scratchapixel.com/lessons/3d-basic-lessons/lesson-7-intersecting-simple-shapes/ray-sphere-intersection/
- http://www.ccs.neu.edu/home/fell/CSU540/programs/RayTracingFormulas.htm
RuiTheAlmeida wrote at 2012-12-09 19:11:
Hi Objo,
I've made some improvements, and you can find the new code below. I'm also posting some test cases. (I've changed your "no intersections" case )
You can also find my detailed explanation and references here
thanks,
Rui
public bool RayIntersection(Ray3D ray, out Point3D[] result) { result = new Point3D[2]; result[0] = new Point3D(double.PositiveInfinity, double.PositiveInfinity, double.PositiveInfinity); result[1] = new Point3D(double.PositiveInfinity, double.PositiveInfinity, double.PositiveInfinity); bool inter = false; double h = this.center.X; double j = this.center.Y; double k = this.center.Z; double ra = this.radius; double x1 = ray.Origin.X; double y1 = ray.Origin.Y; double z1 = ray.Origin.Z; //[a,b,c] = Ray.direction double a = ray.Direction.X; double b = ray.Direction.Y; double c = ray.Direction.Z; //Quadratic solving double aQ = a * a + b * b + c * c; double bQ = 2 * a * (x1 - h) + 2 * b * (y1 - j) + 2 * c * (z1 - k); double cQ = x1 * x1 + y1 * y1 + z1 * z1 + h * h + k * k + j * j - 2 * (j * y1 + k * z1 + h * x1) - ra * ra; //Dscriminant double Q = bQ * bQ - 4 * aQ * cQ; if (Q >= 0) //We have at least one possible Intersection { double rx1 = (-bQ + Math.Sqrt(bQ * bQ - 4 * aQ * cQ)) / (2 * aQ); //First Root double rx2 = (-bQ - Math.Sqrt(bQ * bQ - 4 * aQ * cQ)) / (2 * aQ); //Second Root //first Possible Intersection Point3D I1 = new Point3D(x1 + a * rx1, y1 + b * rx1, z1 + c * rx1); //Second Possible Intersection Point3D I2 = new Point3D(x1 + a * rx2, y1 + b * rx2, z1 + c * rx2); //Check if I1 is intersection by laying in the ray if (Vector3D.DotProduct(ray.Direction, I1 - ray.Origin) >= 0) { result[0] = I1; inter = true; } //Check if I2 is intersection by laying in the ray if (Vector3D.DotProduct(ray.Direction, I2 - ray.Origin) >= 0) { result[1] = I2; inter = true; } } return inter; }
TESTS:
[Test]
public void RayInterSection_HasIntersections()
{
var ray = new Ray3D(new Point3D(0.2, 0.3, 0), new Vector3D(1, 0, 0));
var sphere = new BoundingSphere(new Point3D(0.2, 0.3, 0), 1.0);
Point3D[] result;
Assert.IsTrue(sphere.RayIntersection(ray, out result));
Assert.AreEqual(sphere.Radius, sphere.Center.DistanceTo(result[0]), 1e-6, "Point is not on sphere.");
Assert.AreEqual(0, ray.GetNearest(result[0]).DistanceTo(result[0]), 1e-6, "Point "+ result[0] +" is not on ray.");
}
[Test]
public void RayInterSection_NoIntersections()
{
var ray = new Ray3D(new Point3D(0.2, 0.3, 1), new Vector3D(1, 0.1, 0.1));
var sphere = new BoundingSphere(new Point3D(0.2, 0.3, 0), 1.0);
Point3D[] result;
Assert.IsFalse(sphere.RayIntersection(ray, out result));
Assert.IsTrue( double.IsInfinity(sphere.Center.DistanceTo(result[0])), "There is a finite solution 1.");
Assert.AreNotEqual(0, ray.GetNearest(result[0]).DistanceTo(result[0]), "Point 1 is on ray.");
Assert.IsTrue(double.IsInfinity(sphere.Center.DistanceTo(result[1])), "There is a finite solution 2.");
Assert.AreNotEqual(0, ray.GetNearest(result[0]).DistanceTo(result[1]), "Point 2 is on ray.");
}
[Test]
public void RayInterSection_TwoIntersections()
{
var sphere = new BoundingSphere(new Point3D(0.2, 0.3, 0), 1.0);
var ray = new Ray3D(new Point3D(12, 23, 32), sphere.Center - new Point3D(12, 23, 32)); // vector away and pointing to the center of sphere
Point3D[] result;
Assert.IsTrue(sphere.RayIntersection(ray, out result));
Assert.AreEqual(sphere.Radius, sphere.Center.DistanceTo(result[0]), 1e-6, "Point 1 is not on sphere.");
Assert.AreEqual(0, ray.GetNearest(result[0]).DistanceTo(result[0]), 1e-6, "Point 1" + result[0] + " is not on ray.");
Assert.AreEqual(sphere.Radius, sphere.Center.DistanceTo(result[1]), 1e-6, "Point 2 is not on sphere.");
Assert.AreEqual(0, ray.GetNearest(result[1]).DistanceTo(result[1]), 1e-6, "Point 2 " + result[0] + " is not on ray.");
}
[Test]
public void RayInterSection_OnlyOneIntersection()
{
var sphere = new BoundingSphere(new Point3D(0.2, 0.3, 0), 1.0);
var ray = new Ray3D(sphere.Center, new Vector3D(1, 0.1, 0.1)); // vector origin in the center of sphere pointing away
Point3D[] result;
Assert.IsTrue(sphere.RayIntersection(ray, out result), "No intersection");
Assert.IsTrue(sphere.Center.DistanceTo(ray.Origin ) <= sphere.Radius,"Ray origin outside sphere");
Assert.IsFalse( double.IsInfinity(sphere.Center.DistanceTo(result[0])) && double.IsInfinity(sphere.Center.DistanceTo(result[1])), "No Intersections");
Assert.IsFalse(!double.IsInfinity(sphere.Center.DistanceTo(result[0])) && !double.IsInfinity(sphere.Center.DistanceTo(result[1])), "Two Intersections");
}
[Test]
public void RayInterSection_TwoTangentIntersections()
{
var sphere = new BoundingSphere(new Point3D(0.2, 0.3, 0), 1.0);
var ray = new Ray3D(new Point3D(sphere.Center.X + sphere.Radius, sphere.Center.Y + sphere.Radius, sphere.Center.Z + 50), new Point3D(sphere.Center.X + sphere.Radius, sphere.Center.Y, sphere.Center.Z) - new Point3D(sphere.Center.X + sphere.Radius, sphere.Center.Y + sphere.Radius, sphere.Center.Z + 50)); // ray tangent to sphere
Point3D[] result;
Assert.IsTrue(sphere.RayIntersection(ray, out result), "No intersection");
Assert.IsTrue(!double.IsInfinity(sphere.Center.DistanceTo(result[0])) && !double.IsInfinity(sphere.Center.DistanceTo(result[1])), "Less then two Intersections");
Assert.IsTrue(result[0]==result[1], "Solutions are not equal");
}
objo wrote at 2012-12-13 23:08:
thanks for the unit tests. I have done some modifications and submitted the changes!
Customer support service by UserEcho