Sphere-Line Intersection

classic Classic list List threaded Threaded
9 messages Options
Reply | Threaded
Open this post in threaded view
|

Sphere-Line Intersection

Thomas Zangl - Home
 
Dear List,

I need to implement a Sphere-Line intersection which returns the
intersection points.

Is there anything exact already in CGAL hidden, which I am not able to
find? :-)

If no, does somebody have such a construction which does not use any
sqrt or at least computes the sqrt using exact arithmetic?

I already have an implementation. However, look at this fragment:

   if (i>0)
   {
      // first intersection
      mu = (-b + sqrt(to_double(b*b - 4*a*c ))) / (2*a);
      x = x1 + mu*(x2-x1);
      y = y1 + mu*(y2-y1);
      z = z1 + mu*(z2-z1);
      vIntersections.push_back(Point(x,y,z));

      // second intersection
      mu = (-b - sqrt(to_double(b*b - 4*a*c ))) / (2*a);
      x = x1 + mu*(x2-x1);
      y = y1 + mu*(y2-y1);
      z = z1 + mu*(z2-z1);
      vIntersections.push_back(Point(x,y,z));
      return;
   }

As you can see, it requires the sqrt of the term "(b*b - 4*a*c )"
which introduces some inaccuracy. Can I get rid of it? The intersection
routine is based on
http://ozviz.wasp.uwa.edu.au/~pbourke/geometry/sphereline/.

Has anybody something better? TIA!  

Best regards,
--
,yours Thomas Zangl, Bakk.rer.soc.oec. - [hidden email] -
- Freelancer - IT Consulting & Software Development -
- Student of Software Development-Economy (Master) -
            - http://blog.tzis.net -
--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://lists-sop.inria.fr/wws/info/cgal-discuss
Reply | Threaded
Open this post in threaded view
|

Re: Sphere-Line Intersection

Pedro Machado Manhães de Castro
Hi Thomas,

For now, you may want to have a look at CGAL/Root_of_2 class, for example. You will pass 3 parameters: alpha, beta and gamma (alpha + beta sqrt(gamma)) as Field Type (Kernel::FT), and it supports basic arithmetic operations (though it doesn't look for the squarefree part of gamma, but it feets well for sphere intersections application.. do not expect to sum 3 + 2 sqrt(3) with 2 + 5 sqrt(27)).

However, if waiting until next CGAL's release is possible for you, we will have Sphere Line intersection via a 3D version of the Circular_kernel_2.

Cheers,
Pedro Machado

On Thu, Jul 3, 2008 at 8:48 PM, Thomas Zangl - Home <[hidden email]> wrote:

Dear List,

I need to implement a Sphere-Line intersection which returns the
intersection points.

Is there anything exact already in CGAL hidden, which I am not able to
find? :-)

If no, does somebody have such a construction which does not use any
sqrt or at least computes the sqrt using exact arithmetic?

I already have an implementation. However, look at this fragment:

  if (i>0)
  {
     // first intersection
     mu = (-b + sqrt(to_double(b*b - 4*a*c ))) / (2*a);
     x = x1 + mu*(x2-x1);
     y = y1 + mu*(y2-y1);
     z = z1 + mu*(z2-z1);
     vIntersections.push_back(Point(x,y,z));

     // second intersection
     mu = (-b - sqrt(to_double(b*b - 4*a*c ))) / (2*a);
     x = x1 + mu*(x2-x1);
     y = y1 + mu*(y2-y1);
     z = z1 + mu*(z2-z1);
     vIntersections.push_back(Point(x,y,z));
     return;
  }

As you can see, it requires the sqrt of the term "(b*b - 4*a*c )"
which introduces some inaccuracy. Can I get rid of it? The intersection
routine is based on
http://ozviz.wasp.uwa.edu.au/~pbourke/geometry/sphereline/.

Has anybody something better? TIA!

Best regards,
--
,yours Thomas Zangl, Bakk.rer.soc.oec. - [hidden email] -
- Freelancer - IT Consulting & Software Development -
- Student of Software Development-Economy (Master) -
           - http://blog.tzis.net -
--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://lists-sop.inria.fr/wws/info/cgal-discuss

Reply | Threaded
Open this post in threaded view
|

Re:Sphere-Line Intersection

Thomas Zangl - Home
 
Am Thu, 3 Jul 2008 21:41:14 +0200, schrieb "Pedro_Machado_Manhães_de_Castro" <[hidden email]>:

Hi Pedro!

>For now, you may want to have a look at CGAL/Root_of_2 class, for example. You will pass 3 parameters: alpha, beta and gamma (alpha + beta sqrt(gamma)) as Field Type (Kernel::FT), and it supports basic arithmetic operations (though it doesn&#39;t look for the squarefree part of gamma, but it feets well for sphere intersections application.. do not expect to sum 3 + 2 sqrt(3) with 2 + 5 sqrt(27)).

Can you give me a short example (code) how to use it please?
Especially one, which fits for e.g. the sample code I posted before.
That would help a lot - thank you!

Best regards,
--
,yours Thomas Zangl, Bakk.rer.soc.oec. - [hidden email] -
- Freelancer - IT Consulting & Software Development -
- Student of Software Development-Economy (Master) -
            - http://blog.tzis.net -
--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://lists-sop.inria.fr/wws/info/cgal-discuss
Reply | Threaded
Open this post in threaded view
|

Re: Sphere-Line Intersection

Pedro Machado Manhães de Castro
For example:

typedef typename Kernel::FT                                 FT;
typedef typename Root_of_traits< FT >::RootOf_2 Root_of_2;

const Root_of_2 mu = make_root_of_2(-b/(2*a), FT(1), (b*b - 4*a*c)/(4*a*a));
const Root_of_2 x = x1 + mu*(x2-x1);
const Root_of_2 y = y1 + mu*(y2-y1);
const Root_of_2 z = z1 + mu*(z2-z1);

Note: You can't store them as coordinates for points.
You may want to have a look at CORE::Expr (slow) if you need absolutely to work with CGAL points.
Otherwise, you have to find another way to store them. Just remember that whenever using FT's like Gmpq or MP_Float, you are using something that do not support Root_of_2 (rationals do not include them), and the correspondent Points will not support as well.

Best Regards,
Pedro

On Thu, Jul 3, 2008 at 9:52 PM, Thomas Zangl - Home <[hidden email]> wrote:

Am Thu, 3 Jul 2008 21:41:14 +0200, schrieb "Pedro_Machado_Manhães_de_Castro" <[hidden email]>:

Hi Pedro!

>For now, you may want to have a look at CGAL/Root_of_2 class, for example. You will pass 3 parameters: alpha, beta and gamma (alpha + beta sqrt(gamma)) as Field Type (Kernel::FT), and it supports basic arithmetic operations (though it doesn&#39;t look for the squarefree part of gamma, but it feets well for sphere intersections application.. do not expect to sum 3 + 2 sqrt(3) with 2 + 5 sqrt(27)).

Can you give me a short example (code) how to use it please?
Especially one, which fits for e.g. the sample code I posted before.
That would help a lot - thank you!

Best regards,
--
,yours Thomas Zangl, Bakk.rer.soc.oec. - [hidden email] -
- Freelancer - IT Consulting & Software Development -
- Student of Software Development-Economy (Master) -
           - http://blog.tzis.net -
--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://lists-sop.inria.fr/wws/info/cgal-discuss

Reply | Threaded
Open this post in threaded view
|

Re:Sphere-Line Intersection

Thomas Zangl - Home
 
Am Thu, 3 Jul 2008 22:57:58 +0200, schrieb "Pedro_Machado_Manhães_de_Castro" <[hidden email]>:

Hi!

>typedef typename Kernel::FT                                 FT;
>typedef typename Root_of_traits< FT >::RootOf_2 Root_of_2;
>
>const Root_of_2 mu = make_root_of_2(-b/(2*a), FT(1), (b*b - 4*a*c)/(4*a*a));
>const Root_of_2 x = x1 + mu*(x2-x1);
>const Root_of_2 y = y1 + mu*(y2-y1);
>const Root_of_2 z = z1 + mu*(z2-z1);
>
>Note: You can&#39;t store them as coordinates for points.
>You may want to have a look at CORE::Expr (slow) if you need absolutely to work with CGAL points.
>Otherwise, you have to find another way to store them. Just remember that whenever using FT&#39;s like Gmpq or MP_Float, you are using something that do not support Root_of_2 (rationals do not include them), and the correspondent Points will not support as well.

I already suspected something like this. At least, any ideas how I can
improve accuracy of my intersection code?

The result should be a Point_3 based on Gmpq.

TIA,
--
,yours Thomas Zangl, Bakk.rer.soc.oec. - [hidden email] -
- Freelancer - IT Consulting & Software Development -
- Student of Software Development-Economy (Master) -
            - http://blog.tzis.net -
--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://lists-sop.inria.fr/wws/info/cgal-discuss
Reply | Threaded
Open this post in threaded view
|

Re: Sphere-Line Intersection

Eric Berberich
Thomas Zangl - Home wrote:

>  
> Am Thu, 3 Jul 2008 22:57:58 +0200, schrieb "Pedro_Machado_Manhães_de_Castro" <[hidden email]>:
>
> Hi!
>
>> typedef typename Kernel::FT                                 FT;
>> typedef typename Root_of_traits< FT >::RootOf_2 Root_of_2;
>>
>> const Root_of_2 mu = make_root_of_2(-b/(2*a), FT(1), (b*b - 4*a*c)/(4*a*a));
>> const Root_of_2 x = x1 + mu*(x2-x1);
>> const Root_of_2 y = y1 + mu*(y2-y1);
>> const Root_of_2 z = z1 + mu*(z2-z1);
>>
>> Note: You can&#39;t store them as coordinates for points.
>> You may want to have a look at CORE::Expr (slow) if you need absolutely to work with CGAL points.
>> Otherwise, you have to find another way to store them. Just remember that whenever using FT&#39;s like Gmpq or MP_Float, you are using something that do not support Root_of_2 (rationals do not include them), and the correspondent Points will not support as well.
>
> I already suspected something like this. At least, any ideas how I can
> improve accuracy of my intersection code?
>
> The result should be a Point_3 based on Gmpq.
>
> TIA,

gmpq is not able to exactly represent the square-root that naturally
exists (except for degenerate cases, where the root is a rational) for
the coordinates of the intersection points of a sphere and a line.
Either you compute an interval of two gmpq that contain the exact
coordinates, or you switch to CORE::Expr (or leda::real) or you use the
Root_of number types. I currently don't see another option.

It is the mathematics that result in this "gmpq-is-not-sufficient" case.
Believe us, it would be much nicer, if curved geometry would stay in
rational numbers ;-)

eriC
--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://lists-sop.inria.fr/wws/info/cgal-discuss
Reply | Threaded
Open this post in threaded view
|

consistent cell normals

dfwang
In reply to this post by Thomas Zangl - Home
Dear CGALers

For a given closed triangle mesh, how to make sure that the normal of each
triangle points to the inside of the closed surface in CGAL?

Any comments are appreciated.

Best wishes,
Vincent

--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://lists-sop.inria.fr/wws/info/cgal-discuss
Reply | Threaded
Open this post in threaded view
|

Re: consistent cell normals

andreas.fabri
dfwang wrote:
> Dear CGALers
>
> For a given closed triangle mesh, how to make sure that the normal of
> each triangle points to the inside of the closed surface in CGAL?
>
> Any comments are appreciated.
>
> Best wishes,
> Vincent

Hi Vincent,

if the mesh is in a CGAL::Polyhedron already then the faces are
oriented consistently already, so determining inside/outside for
one face plus propagation should do the job.

If you just have a triangle soup (std::list<CGAL::Triangle_3<K>>)
then you first have to orient your triangles consistently.
There is no out-of the box function in CGAL. It would be a good idea
though.

andreas
--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://lists-sop.inria.fr/wws/info/cgal-discuss
Reply | Threaded
Open this post in threaded view
|

Re: consistent cell normals

Rafael Roberto
Vincent,

just how Andreas said, if your mesh is a CGAL::Polyhedron, then the faces are oriented consistently and, in every case I SAW, they points outside the surface. It make sense to me that the normal points outside, but I don't know if a CGAL::Polyhedron will always put they normal outside.

Just to guarantee that the normal are outside and you have to invert then to point inside, you can take a vector that start at any point of the triangle and point to the center of the polyhedron. If the angle between the normal and the vector that points to the center is bigger then 90° the normal is really point outside the polyhedron and you must invert it. Because all normal are consistent, you have to do this for just one triangle.


Valeu,
Rafael Roberto