

Dear List,
I need to implement a SphereLine 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*(x2x1);
y = y1 + mu*(y2y1);
z = z1 + mu*(z2z1);
vIntersections.push_back(Point(x,y,z));
// second intersection
mu = (b  sqrt(to_double(b*b  4*a*c ))) / (2*a);
x = x1 + mu*(x2x1);
y = y1 + mu*(y2y1);
z = z1 + mu*(z2z1);
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 DevelopmentEconomy (Master) 
 http://blog.tzis.net 

You are currently subscribed to cgaldiscuss.
To unsubscribe or access the archives, go to
https://listssop.inria.fr/wws/info/cgaldiscuss


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 SphereLine 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*(x2x1);
y = y1 + mu*(y2y1);
z = z1 + mu*(z2z1);
vIntersections.push_back(Point(x,y,z));
// second intersection
mu = (b  sqrt(to_double(b*b  4*a*c ))) / (2*a);
x = x1 + mu*(x2x1);
y = y1 + mu*(y2y1);
z = z1 + mu*(z2z1);
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 DevelopmentEconomy (Master) 
 http://blog.tzis.net 

You are currently subscribed to cgaldiscuss.
To unsubscribe or access the archives, go to
https://listssop.inria.fr/wws/info/cgaldiscuss


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'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 DevelopmentEconomy (Master) 
 http://blog.tzis.net 

You are currently subscribed to cgaldiscuss.
To unsubscribe or access the archives, go to
https://listssop.inria.fr/wws/info/cgaldiscuss


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*(x2x1); const Root_of_2 y = y1 + mu*(y2y1); const Root_of_2 z = z1 + mu*(z2z1); 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'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!


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*(x2x1);
>const Root_of_2 y = y1 + mu*(y2y1);
>const Root_of_2 z = z1 + mu*(z2z1);
>
>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.
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 DevelopmentEconomy (Master) 
 http://blog.tzis.net 

You are currently subscribed to cgaldiscuss.
To unsubscribe or access the archives, go to
https://listssop.inria.fr/wws/info/cgaldiscuss


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*(x2x1);
>> const Root_of_2 y = y1 + mu*(y2y1);
>> const Root_of_2 z = z1 + mu*(z2z1);
>>
>> 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.
>
> 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 squareroot 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 "gmpqisnotsufficient" case.
Believe us, it would be much nicer, if curved geometry would stay in
rational numbers ;)
eriC

You are currently subscribed to cgaldiscuss.
To unsubscribe or access the archives, go to
https://listssop.inria.fr/wws/info/cgaldiscuss


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 cgaldiscuss.
To unsubscribe or access the archives, go to
https://listssop.inria.fr/wws/info/cgaldiscuss


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 outof the box function in CGAL. It would be a good idea
though.
andreas

You are currently subscribed to cgaldiscuss.
To unsubscribe or access the archives, go to
https://listssop.inria.fr/wws/info/cgaldiscuss


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

