0

BoundingSphere ,RayIntersection implementation ...

Anonymous 10 years ago 0
This discussion was imported from CodePlex

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



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!