CGAL users discussion list

[Sebastien Loriot <>]

Florian Wobbe wrote:
> <div class="moz-text-flowed" style="font-family: -moz-fixed">Hello,
>
> I'm doing my PhD in structural geology as opposed to structural 
> biology. The aim is to do a paleobathymetric reconstruction of the 
> Southern Ocean for circumantarctic paleoclimate modelling. So I have 
> to deal with rotations of plates on the sphere.
>
> I need an algorithm for the calculation of the intersection of a small 
> circle on a sphere (earth) with circular arc segments on the same 
> sphere (e.g. segments of a polygon representing a continent). If you 
> are interested in any details I'd be happy to explain that. However, 
> just for the understanding of the problem one can imagine a latitude 
> (small circle) intersecting the coastline of a continent. The 
> coastline is aproximated by great circle segments given by 
> long/lat-coordinates.
>
> I wonder if I can use your extension to the 3d Spherical Kernel for 
> circles on a reference sphere to approach this problem and not 
> reinvent the wheel. E.g:
>
> 1) construct a reference sphere
> 2) construct a circle on the reference sphere
> 3) iterate through all polygon segments, construct circular arcs for 
> each and calculate intersection with (2)
>
> I'd be much obliged for any feedback and help.
>
> Cheers
> Florian
>
> --
> Florian Wobbe
> Alfred Wegener Institute for Polar and Marine Research
> Division of Geosciences, Section of Geophysics
> Postfach 120161
> 27515 Bremerhaven
>
> </div>
Dear Florian,

From what I understand, you are given  spherical polygons representing 
continents, and you are
interested in splitting the spherical segments of your continents using 
circles in horizontal planes.

Unfortunately, the current public release of CGAL (www.cgal.org) does 
not already contain the extension
for circles on a reference sphere (probably scheduled for the next 
release if accepted).

For the moment, CGAL only contains the Spherical Kernel with general 
objects.
The good point is that you can still use it as there are 3D circles, 
points on spheres and circular arcs on spheres
available. The extension on a reference sphere mainly provides 
additional operations using cylindrical coordinates
and a more appropriate representation of objects.

I attach in this email an example of how to use the CGAL Spherical 
Kernel to compute the intersection of circular arcs
by a circle on a sphere, and split one arc.

You will first need to download and install CGAL at www.cgal.org
Then,  in the directory containing the attach file do
> PATH_TO_CGAL/scripts/cgal_create_cmake_script
> cmake .
> make


The documentation is available at :
http://www.cgal.org/Manual/3.4/doc_html/cgal_manual/Circular_kernel_3/Chapter_main.html 


I am also sending this email to the CGAL mailing list (cgal-discuss) as 
this can be helpful
for other persons.

Best regards,
Sebastien.
#include <CGAL/Exact_spherical_kernel_3.h>
#include <list>

typedef CGAL::Exact_spherical_kernel_3         SK;

//Point in 3D with rational coordinates
typedef CGAL::Point_3<SK>                   Point_3;
//Point in 3D with algebraic numbers of degree 2 as coordinates
typedef CGAL::Circular_arc_point_3<SK>      Circular_arc_point_3;
typedef CGAL::Circular_arc_3<SK>            Circular_arc_3;
typedef CGAL::Circle_3<SK>                  Circle_3;
typedef CGAL::Plane_3<SK>                   Plane_3;
typedef CGAL::Sphere_3<SK>                  Sphere_3;




int main()
{
  //Create a sphere centered at the origin of squared radius 10
  Sphere_3 S(Point_3(0,0,0),10);
  
   
  //create a circle on sphere S on the plane y=0
  Circle_3 C1(S, Plane_3(0,1,0,0) );
  //create a circle on sphere S on the plane x+y=0
  Circle_3 C2(S, Plane_3(1,1,0,0) );
  //create a horizontal circle on sphere S on the plane z=-2
  Circle_3 C3(S, Plane_3(0,0,1,2) );
  
  
  
  //Look at the manual page : http://www.cgal.org/Manual/3.4/doc_html/cgal_manual/Circular_kernel_3_ref/Class_Circular_arc_3.html
  //to see how the source and target of an arc are defined.
  //create arc on C1, endpoints are exactly on the sphere
  Circular_arc_3 arc1 (C1 , Point_3(3,0,-1) , Point_3(-3,0,-1));
  
  
  //create arc on C2, endpoints are (3/sqrt(2),-3/sqrt(2),-1) and (-3/sqrt(2),3/sqrt(2),-1)
  //Point_3 P(3/sqrt(2),-3/sqrt(2),-1)   <------------ this produce an error as the  square root is not allowed using rational numbers
  // Circular_arc_3 arc2 (C2 , Point_3(2.12,-2.12,-1) , Point_3(-2.12,2.12,-1))   <------------ this produce an error as the  points are not exactly on the sphere
  
  //approximate source and target
  Point_3 source_approx(2.12,-2.12,-1);
  Point_3 target_approx(-2.12,2.12,-1);
  //project approximate points onto the sphere S 
  //note that *_approx are in the plane of C2 that passes through the center of S, thus the projections will be on C2
  CGAL::Root_of_traits< SK::FT >::RootOf_2 Root_source=CGAL::make_root_of_2(SK::FT(0),SK::FT(1),S.squared_radius()/(source_approx-S.center()).squared_length());
  CGAL::Root_of_traits< SK::FT >::RootOf_2 Root_target=CGAL::make_root_of_2(SK::FT(0),SK::FT(1),S.squared_radius()/(target_approx-S.center()).squared_length());
  SK::Circular_arc_point_3 source_exact(SK::Root_for_spheres_2_3(source_approx.x()*Root_source,source_approx.y()*Root_source,source_approx.z()*Root_source));
  SK::Circular_arc_point_3 target_exact(SK::Root_for_spheres_2_3(target_approx.x()*Root_target,target_approx.y()*Root_target,target_approx.z()*Root_target));
  
  Circular_arc_3 arc2 (C2 ,source_exact,target_exact);  //<--- this is correct
  
  std::list<CGAL::Object> lst;
  SK::Intersect_3 functor;
  functor(arc1,C3,std::back_inserter(lst));
  functor(arc2,C3,std::back_inserter(lst));
  std::cout << "nb of intersection pts : " << lst.size() << std::endl;
  
  std::list <Circular_arc_point_3>  intersections;
  
  for (std::list<CGAL::Object>::iterator it=lst.begin();it!=lst.end();++it){
    //recover the intersection points
    std::pair <Circular_arc_point_3,unsigned> p=CGAL::object_cast<std::pair <Circular_arc_point_3,unsigned> > (*it);
    intersections.push_back(p.first);
  }
  
  //split arc1 using on of its intersection point with circle C3
  Circular_arc_3 subarc1,subarc2;
  SK::Split_3()(arc1,*intersections.begin(),subarc1,subarc2);
  
  
  return 0;
};