SPAR/AIModule/BWTA/vendors/CGAL/CGAL/Arr_geodesic_arc_on_sphere_traits_2.h

Go to the documentation of this file.

// Copyright (c) 2006 Tel-Aviv University (Israel).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you may redistribute it under
// the terms of the Q Public License version 1.0.
// See the file LICENSE.QPL distributed with CGAL.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.5-branch/Arrangement_on_surface_2/include/CGAL/Arr_geodesic_arc_on_sphere_traits_2.h $
// $Id: Arr_geodesic_arc_on_sphere_traits_2.h 49924 2009-06-16 07:14:42Z efif $
// 
//
// Author(s)     : Efi Fogel         <efif@post.tau.ac.il>

#ifndef CGAL_ARR_GEODESIC_ARC_ON_SPHERE_TRAITS_2_H
#define CGAL_ARR_GEODESIC_ARC_ON_SPHERE_TRAITS_2_H

// #define CGAL_FULL_X_MONOTONE_GEODESIC_ARC_ON_SPHERE_IS_SUPPORTED 1

#include <fstream>

#include <CGAL/config.h>
#include <CGAL/tags.h>
#include <CGAL/intersections.h>
#include <CGAL/Arr_tags.h>
#include <CGAL/Arr_enums.h>

CGAL_BEGIN_NAMESPACE

#define CGAL_X_MINUS_1_Y_0      0
#define CGAL_X_MINUS_8_Y_6      1
#define CGAL_X_MINUS_11_Y_7     2

#ifndef CGAL_IDENTIFICATION_XY
#define CGAL_IDENTIFICATION_XY  CGAL_X_MINUS_1_Y_0
#endif

template <typename Kernel> class Arr_x_monotone_geodesic_arc_on_sphere_3;
template <typename Kernel> class Arr_geodesic_arc_on_sphere_3;
template <typename Kernel> class Arr_extended_direction_3;

template <typename T_Kernel>
class Arr_geodesic_arc_on_sphere_traits_2 : public T_Kernel {
  friend class Arr_x_monotone_geodesic_arc_on_sphere_3<T_Kernel>;
  friend class Arr_geodesic_arc_on_sphere_3<T_Kernel>;
  friend class Arr_extended_direction_3<T_Kernel>;

public:
  typedef T_Kernel                              Kernel;
  
  // Category tags:
  typedef Tag_true                              Has_left_category;
  typedef Tag_true                              Has_merge_category;

  typedef Arr_identified_side_tag               Arr_left_side_tag;
  typedef Arr_contracted_side_tag               Arr_bottom_side_tag;
  typedef Arr_contracted_side_tag               Arr_top_side_tag;
  typedef Arr_identified_side_tag               Arr_right_side_tag;

  Arr_geodesic_arc_on_sphere_traits_2(){}

protected:
  typedef typename Kernel::FT                   FT;

  typedef typename Kernel::Direction_3          Direction_3;
  typedef typename Kernel::Vector_3             Vector_3;
  typedef typename Kernel::Direction_2          Direction_2;
  typedef typename Kernel::Vector_2             Vector_2;

  inline static const Direction_3 & pos_pole()
  {
    static const Direction_3 d(0, 0, 1);
    return d;
  }

  inline static const Direction_3 & neg_pole()
  {
    static const Direction_3 d(0, 0, -1);
    return d;
  }

  inline static const Direction_2 & identification_xy()
  {
#if (CGAL_IDENTIFICATION_XY == CGAL_X_MINUS_1_Y_0)
    static const Direction_2 d(-1, 0);
#elif (CGAL_IDENTIFICATION_XY == CGAL_X_MINUS_8_Y_6)    
    static const Direction_2 d(-8, 6);
#elif (CGAL_IDENTIFICATION_XY == CGAL_X_MINUS_11_Y_7)    
    static const Direction_2 d(-11, 7);
#else
#error CGAL_IDENTIFICATION_XY is not defined
#endif
    return d;
  }

  inline static const Direction_3 & identification_normal()
  {
#if (CGAL_IDENTIFICATION_XY == CGAL_X_MINUS_1_Y_0)
    static const Direction_3 d(0, 1, 0);
#elif (CGAL_IDENTIFICATION_XY == CGAL_X_MINUS_8_Y_6)
    static const Direction_3 d(6, 8, 0);
#elif (CGAL_IDENTIFICATION_XY == CGAL_X_MINUS_11_Y_7)
    static const Direction_3 d(7, 11, 0);
#else
#error CGAL_IDENTIFICATION_XY is not defined
#endif
    return d;
  }

  inline static const Direction_2 & neg_x_2()
  {
    static const Direction_2 d(-1, 0);
    return d;
  }

  inline static const Direction_2 & neg_y_2()
  {
    static const Direction_2 d(0, -1);
    return d;
  }

  inline static Sign x_sign(Direction_3 d) { return CGAL::sign(d.dx()); }

  inline static Sign y_sign(Direction_3 d) { return CGAL::sign(d.dy()); }  

  inline static Sign z_sign(Direction_3 d) { return CGAL::sign(d.dz()); }
  
  typedef Direction_2 (*Project)(const Direction_3 & d) ;
  
  inline static Direction_2 project_xy(const Direction_3 & d)
  { return Direction_2(d.dx(), d.dy()); }
  
  inline static Direction_2 project_yz(const Direction_3 & d)
  { return Direction_2(d.dy(), d.dz()); }
  
  inline static Direction_2 project_xz(const Direction_3 & d)
  { return Direction_2(d.dx(), d.dz()); }

  inline Oriented_side oriented_side(const Direction_3 & normal,
                                     const Direction_3 dir) const
  {
    FT dot = normal.vector() * dir.vector();
    return CGAL::sign(dot);
  }
  
  static inline Orientation orientation(const Direction_2 & d1,
                                        const Direction_2 & d2)
  {
    Kernel kernel;
    return kernel.orientation_2_object()(d1.vector(), d2.vector());
  }

  inline bool has_on(const Direction_3 & normal, const Direction_3 & dir) const
  {
    FT dot = normal.vector() * dir.vector();
    return CGAL::sign(dot) == ZERO;
  }
  
public:
  inline Comparison_result compare_y(const Direction_3 & d1,
                                     const Direction_3 & d2) const
  {
    Vector_3 v1 = d1.vector();
    Vector_3 v2 = d2.vector();

    FT dot_p1 = v1.z();
    FT dot_p2 = v2.z();

    Sign s1 = CGAL::sign(dot_p1);
    Sign s2 = CGAL::sign(dot_p2);
    
    if (s1 != s2) {
      if (s1 == NEGATIVE) return SMALLER;
      if (s1 == POSITIVE) return LARGER;
      if (s2 == NEGATIVE) return LARGER;
      if (s2 == POSITIVE) return SMALLER;      
    }
    if (s1 == ZERO) return EQUAL;
    
    FT norm1 = v1 * v1;
    FT norm2 = v2 * v2;

    return (s1 == POSITIVE) ?
      compare(dot_p1 * dot_p1 * norm2, dot_p2 * dot_p2 * norm1) :
      compare(dot_p2 * dot_p2 * norm1, dot_p1 * dot_p1 * norm2);
  }

  inline Comparison_result compare_x(const Direction_2 & d1,
                                     const Direction_2 & d2) const
  {
    const Kernel * kernel = this;
    if (kernel->equal_2_object()(d1, d2)) return EQUAL;
    const Direction_2 & d = identification_xy();
    return (kernel->counterclockwise_in_between_2_object()(d, d1, d2)) ?
      LARGER : SMALLER;
  }
  
  inline Comparison_result compare_x(const Direction_3 & d1,
                                     const Direction_3 & d2) const
  {
    // Compare the projections onto the xy plane:
    Direction_2 d1_2 = project_xy(d1);
    Direction_2 d2_2 = project_xy(d2);
    return compare_x(d1_2, d2_2);
  }

  inline Comparison_result compare_xy(const Direction_3 & d1,
                                      const Direction_3 & d2) const
  {
    Comparison_result res = compare_x(d1, d2);
    if (res == EQUAL) return compare_y(d1, d2);
    return res;
  }
  
public:
  // Traits objects
  typedef Arr_extended_direction_3<Kernel>                Point_2;
  typedef Arr_x_monotone_geodesic_arc_on_sphere_3<Kernel> X_monotone_curve_2;
  typedef Arr_geodesic_arc_on_sphere_3<Kernel>            Curve_2;
  typedef unsigned int                                    Multiplicity;

public:


  class Compare_x_2 {
  protected:
    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;

    const Traits * m_traits;
    
    Compare_x_2(const Traits * traits) : m_traits(traits) {}

    friend class Arr_geodesic_arc_on_sphere_traits_2<Kernel>;
    
  public:    
    Comparison_result operator()(const Point_2 & p1, const Point_2 & p2) const
    {
      CGAL_precondition(p1.is_no_boundary());
      CGAL_precondition(p2.is_no_boundary());

      return m_traits->compare_x(p1, p2);
    }
  };

  Compare_x_2 compare_x_2_object() const { return Compare_x_2(this); }

  class Compare_xy_2 {
  protected:
    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;

    const Traits * m_traits;
    
    Compare_xy_2(const Traits * traits) : m_traits(traits) {}

    friend class Arr_geodesic_arc_on_sphere_traits_2<Kernel>;
    
  public:
    Comparison_result operator()(const Point_2 & p1, const Point_2 & p2) const
    {
      CGAL_precondition(p1.is_no_boundary());
      CGAL_precondition(p2.is_no_boundary());
      
      return m_traits->compare_xy(p1, p2);
    }
  };

  Compare_xy_2 compare_xy_2_object() const { return Compare_xy_2(this); }

  class Construct_min_vertex_2 {
  public:
    const Point_2 & operator()(const X_monotone_curve_2 & xc) const
    { return xc.left(); }
  };

  Construct_min_vertex_2 construct_min_vertex_2_object() const
  { return Construct_min_vertex_2(); }

  class Construct_max_vertex_2 {
  public:
    const Point_2 & operator()(const X_monotone_curve_2 & xc) const
    { return xc.right(); }
  };

  Construct_max_vertex_2 construct_max_vertex_2_object() const
  { return Construct_max_vertex_2(); }

  class Is_vertical_2 {
  public:
    bool operator()(const X_monotone_curve_2 & xc) const
    {
      CGAL_precondition(!xc.is_degenerate());
      return xc.is_vertical();
    }
  };

  Is_vertical_2 is_vertical_2_object() const { return Is_vertical_2(); }

  class Compare_y_at_x_2 {
  protected:
    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;

    const Traits * m_traits;
    
    Compare_y_at_x_2(const Traits * traits) : m_traits(traits) {}

    friend class Arr_geodesic_arc_on_sphere_traits_2<Kernel>;
    
  public:
    Comparison_result operator()(const Point_2 & p,
                                 const X_monotone_curve_2 & xc) const
    {
      CGAL_precondition(!p.is_min_boundary() && !p.is_max_boundary());
      CGAL_precondition(xc.is_in_x_range(p));

      if (xc.is_vertical()) {
        // Compare the point with the left endpoint. If smaller, return SMALLER.
        // Otherwise, if EQUAL, return EQUAL.
        // Otherwise, compare with the right endpoint. If larger, return LARGER.
        // Otherwise, return EQUAL:
        if (!xc.left().is_min_boundary()) {
          Comparison_result cr = m_traits->compare_y(p, xc.left());
          if (cr != LARGER) return cr;
        }
        if (xc.right().is_max_boundary()) return EQUAL;
        Comparison_result cr = m_traits->compare_y(p, xc.right());
        return (cr == LARGER) ? LARGER : EQUAL;
      }

      // Compare the point to the underlying plane of xc:
      Oriented_side os = m_traits->oriented_side(xc.normal(), p);
      return (os == ON_ORIENTED_BOUNDARY) ? EQUAL :
        (xc.is_directed_right()) ?
        ((os == ON_NEGATIVE_SIDE) ? SMALLER : LARGER) : 
        ((os == ON_NEGATIVE_SIDE) ? LARGER : SMALLER);
    }
  };

  Compare_y_at_x_2 compare_y_at_x_2_object() const
  { return Compare_y_at_x_2(this); }

  class Compare_y_at_x_left_2 {
  protected:
    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;

    const Traits * m_traits;
    
    Compare_y_at_x_left_2(const Traits * traits) : m_traits(traits) {}

    friend class Arr_geodesic_arc_on_sphere_traits_2<Kernel>;
    
  public:
    Comparison_result operator()(const X_monotone_curve_2 & xc1,
                                 const X_monotone_curve_2 & xc2,
                                 const Point_2 & p) const
    {
      CGAL_precondition(!xc1.is_degenerate());
      CGAL_precondition(!xc2.is_degenerate());
      CGAL_precondition(p == xc1.right());
      CGAL_precondition(p == xc2.right());

      // If Both arcs are vertical, they overlap:
      if (xc1.is_vertical() && xc2.is_vertical()) return EQUAL;
      if (xc1.is_vertical()) return SMALLER;
      if (xc2.is_vertical()) return LARGER;

      // Non of the arc is verticel. Thus, non of the endpoints coincide with
      // a pole.
      // Compare the y-coord. at the x-coord of the most right left-endpoint.
      const Point_2 & l1 = xc1.left();
      const Point_2 & l2 = xc2.left();
      if (!l1.is_no_boundary()) {
        // use l2 and xc1:
        Oriented_side os = m_traits->oriented_side(xc1.normal(), l2);
        return (os == ON_ORIENTED_BOUNDARY) ? EQUAL :
          (xc1.is_directed_right()) ?
          ((os == ON_NEGATIVE_SIDE) ? LARGER : SMALLER) : 
          ((os == ON_NEGATIVE_SIDE) ? SMALLER : LARGER);
      }
      if (!l2.is_no_boundary()) {
        // use l1 and xc2:
        Oriented_side os = m_traits->oriented_side(xc2.normal(), l1);
        return (os == ON_ORIENTED_BOUNDARY) ? EQUAL :
          (xc2.is_directed_right()) ?
          ((os == ON_NEGATIVE_SIDE) ? SMALLER : LARGER) : 
          ((os == ON_NEGATIVE_SIDE) ? LARGER : SMALLER);
      }

      if (m_traits->compare_xy(l1, l2) == SMALLER) {
        // use l2 and xc1:
        Oriented_side os = m_traits->oriented_side(xc1.normal(), l2);
        return (os == ON_ORIENTED_BOUNDARY) ? EQUAL :
          (xc1.is_directed_right()) ?
          ((os == ON_NEGATIVE_SIDE) ? LARGER : SMALLER) : 
          ((os == ON_NEGATIVE_SIDE) ? SMALLER : LARGER);
      }
      // use l1 and xc2: 
      Oriented_side os = m_traits->oriented_side(xc2.normal(), l1);
      return (os == ON_ORIENTED_BOUNDARY) ? EQUAL :
        (xc2.is_directed_right()) ?
        ((os == ON_NEGATIVE_SIDE) ? SMALLER : LARGER) : 
        ((os == ON_NEGATIVE_SIDE) ? LARGER : SMALLER);
    }
  };

  Compare_y_at_x_left_2 compare_y_at_x_left_2_object() const
  { return Compare_y_at_x_left_2(this); }
  
  class Compare_y_at_x_right_2 {
  protected:
    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;

    const Traits * m_traits;
    
    Compare_y_at_x_right_2(const Traits * traits) : m_traits(traits) {}

    friend class Arr_geodesic_arc_on_sphere_traits_2<Kernel>;
    
  public:
    Comparison_result operator()(const X_monotone_curve_2 & xc1,
                                 const X_monotone_curve_2 & xc2,
                                 const Point_2 &) const
    {
      CGAL_precondition(!xc1.is_degenerate());
      CGAL_precondition(!xc2.is_degenerate());

      // CGAL_precondition(p == xc1.left());
      // CGAL_precondition(p == xc2.left());

      // If Both arcs are vertical, they overlap:
      if (xc1.is_vertical() && xc2.is_vertical()) return EQUAL;
      if (xc1.is_vertical()) return LARGER;
      if (xc2.is_vertical()) return SMALLER;

      // Non of the arcs is verticel. Thus, non of the endpoints coincide with
      // a pole.
      // Compare the y-coord. at the x-coord of the most left right-endpoint.
      const Point_2 & r1 = xc1.right();
      const Point_2 & r2 = xc2.right();
      if (!r1.is_no_boundary()) {
        // use r2 and xc1:
        Oriented_side os = m_traits->oriented_side(xc1.normal(), r2);
        return (os == ON_ORIENTED_BOUNDARY) ? EQUAL :
          (xc1.is_directed_right()) ?
          ((os == ON_NEGATIVE_SIDE) ? LARGER : SMALLER) : 
          ((os == ON_NEGATIVE_SIDE) ? SMALLER : LARGER);
      }
      if (!r2.is_no_boundary()) {
        // use r1 and xc2:
        Oriented_side os = m_traits->oriented_side(xc2.normal(), r1);
        return (os == ON_ORIENTED_BOUNDARY) ? EQUAL :
          (xc2.is_directed_right()) ?
          ((os == ON_NEGATIVE_SIDE) ? SMALLER : LARGER) : 
          ((os == ON_NEGATIVE_SIDE) ? LARGER : SMALLER);
      }
      if (m_traits->compare_xy(r1, r2) == LARGER) {
        // use r2 and xc1:
        Oriented_side os = m_traits->oriented_side(xc1.normal(), r2);
        return (os == ON_ORIENTED_BOUNDARY) ? EQUAL :
          (xc1.is_directed_right()) ?
          ((os == ON_NEGATIVE_SIDE) ? LARGER : SMALLER) : 
          ((os == ON_NEGATIVE_SIDE) ? SMALLER : LARGER);
      }
      // use r1 and xc2: 
      Oriented_side os = m_traits->oriented_side(xc2.normal(), r1);
      return (os == ON_ORIENTED_BOUNDARY) ? EQUAL :
        (xc2.is_directed_right()) ?
        ((os == ON_NEGATIVE_SIDE) ? SMALLER : LARGER) : 
        ((os == ON_NEGATIVE_SIDE) ? LARGER : SMALLER);
    }
  };

  Compare_y_at_x_right_2 compare_y_at_x_right_2_object() const
  { return Compare_y_at_x_right_2(this); }

  class Equal_2 {
  protected:
    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;

    const Traits * m_traits;
    
    Equal_2(const Traits * traits) : m_traits(traits) {}

    friend class Arr_geodesic_arc_on_sphere_traits_2<Kernel>;
    
  public:
    bool operator()(const X_monotone_curve_2 & xc1,
                    const X_monotone_curve_2 & xc2) const
    {
      const Kernel * kernel = m_traits;
      typename Kernel::Equal_3 equal_3 = kernel->equal_3_object();
      if (xc1.is_full() || xc2.is_full()) {
        return (xc1.is_full() && xc2.is_full() &&
                equal_3(xc1.left(), xc2.left()));
      }
      
      return (equal_3(xc1.left(), xc2.left()) &&
              equal_3(xc1.right(), xc2.right()));
    }

    bool operator()(const Point_2 & p1, const Point_2 & p2) const
    {
      const Kernel * kernel = m_traits;
      return kernel->equal_3_object()(p1, p2);
    }
  };

  Equal_2 equal_2_object() const { return Equal_2(this); }



  class Parameter_space_in_x_2 {
  public:
    Arr_parameter_space operator()(const X_monotone_curve_2 & xcv,
                                   Arr_curve_end ce) const
    {
      if (xcv.is_vertical()) {
        CGAL_precondition(!xcv.is_on_boundary());
        return ARR_INTERIOR;
      }

      return (ce == ARR_MIN_END) ?
        ((xcv.left().is_no_boundary()) ? ARR_INTERIOR : ARR_LEFT_BOUNDARY) :
        ((xcv.right().is_no_boundary()) ? ARR_INTERIOR : ARR_RIGHT_BOUNDARY);
    }

    Arr_parameter_space operator()(const Point_2 p) const
    {
      return (p.is_mid_boundary()) ? ARR_LEFT_BOUNDARY : ARR_INTERIOR;
    }
  };

  Parameter_space_in_x_2 parameter_space_in_x_2_object() const
  { return Parameter_space_in_x_2(); }
  
  class Parameter_space_in_y_2 {
  public:
    Arr_parameter_space operator()(const X_monotone_curve_2 & xcv,
                                   Arr_curve_end ce) const
    {
      return (ce == ARR_MIN_END) ?
        ((xcv.left().is_min_boundary()) ?  ARR_BOTTOM_BOUNDARY: ARR_INTERIOR) :
        ((xcv.right().is_max_boundary()) ? ARR_TOP_BOUNDARY : ARR_INTERIOR);
    }

    Arr_parameter_space operator()(const Point_2 p) const
    {
      return
        (p.is_min_boundary()) ? ARR_BOTTOM_BOUNDARY :
        (p.is_max_boundary()) ? ARR_TOP_BOUNDARY : ARR_INTERIOR;
    }
  };

  Parameter_space_in_y_2 parameter_space_in_y_2_object() const
  { return Parameter_space_in_y_2(); }

  class Compare_x_near_boundary_2 {
  protected:
    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;

    const Traits * m_traits;
    
    Compare_x_near_boundary_2(const Traits * traits) : m_traits(traits) {}

    friend class Arr_geodesic_arc_on_sphere_traits_2<Kernel>;
    
  public:
    Comparison_result operator()(const Point_2 & point,
                                 const X_monotone_curve_2 & xcv,
                                 Arr_curve_end ce) const
    {
      CGAL_precondition(point.is_no_boundary());
      CGAL_precondition_code
        (const Point_2 & p2 = (ce == ARR_MIN_END) ? xcv.left() : xcv.right(););
      CGAL_precondition(!p2.is_no_boundary());

      if (xcv.is_vertical()) {
        CGAL_precondition(!xcv.is_on_boundary());

        /* The following is replaced by the precondition above.
         * If xc coincides with the discontinuity arc, all its points are
         * assumed to be smaller than non-boundary points:
         * if (xcv.is_on_boundary()) return LARGER;
         */
        
        // xcv is vertical, but does not coincide with the discontinuity arc.
        // Obtain the direction contained in the underlying plane, which is
        // also on the xy-plane:
        Direction_3 normal = xcv.normal();
        Direction_2 q = (xcv.is_directed_right()) ?
          Direction_2(-(normal.dy()), normal.dx()) :
          Direction_2(normal.dy(), -(normal.dx()));
        Direction_2 p = Traits::project_xy(point);
        return m_traits->compare_x(p, q);
      }

      // xcv is not a vertical sphercial_arc:
      return (ce == ARR_MIN_END) ? LARGER : SMALLER;
    }

    Comparison_result operator()(const X_monotone_curve_2 & xcv1,
                                 Arr_curve_end ce1,
                                 const X_monotone_curve_2 & xcv2,
                                 Arr_curve_end ce2) const
    {
      CGAL_precondition_code
        (const Point_2 & p1 = (ce1 == ARR_MIN_END) ? xcv1.left() : xcv1.right(););
      CGAL_precondition(!p1.is_no_boundary());
      CGAL_precondition_code
        (const Point_2 & p2 = (ce2 == ARR_MIN_END) ? xcv2.left() : xcv2.right(););
      CGAL_precondition(!p2.is_no_boundary());

      if (xcv1.is_vertical() && xcv2.is_vertical()) {
        CGAL_precondition(!xcv1.is_on_boundary());
        CGAL_precondition(!xcv2.is_on_boundary());

        /* The following is replaced by the code above
         * if (xcv1.is_on_boundary() && xcv2.is_on_boundary()) return EQUAL;
         * if (xcv1.is_on_boundary()) return SMALLER;
         * if (xcv2.is_on_boundary()) return LARGER;
         */
        
        // Non of the arcs coincide with the discontinuity arc:
        // Obtain the directions contained in the underlying planes, which are
        // also on the xy-plane:
        Direction_3 normal1 = xcv1.normal();
        Direction_2 p = (xcv1.is_directed_right()) ?
          Direction_2(-(normal1.dy()), normal1.dx()) :
          Direction_2(normal1.dy(), -(normal1.dx()));
        Direction_3 normal2 = xcv2.normal();
        Direction_2 q = (xcv2.is_directed_right()) ?
          Direction_2(-(normal2.dy()), normal2.dx()) :
          Direction_2(normal2.dy(), -(normal2.dx()));
        return m_traits->compare_x(p, q);
      }
      if (xcv1.is_vertical()) {
        CGAL_precondition(!xcv1.is_on_boundary());
        /* The following is replaced by the code above
         * if (xcv1.is_on_boundary()) return SMALLER;
         */
        return (ce2 == ARR_MAX_END) ? SMALLER : LARGER;
      }
      if (xcv2.is_vertical()) {
        CGAL_precondition(!xcv2.is_on_boundary());
        /* The following is replaced by the code above
         * if (xcv2.is_on_boundary()) return LARGER;
         */
        return (ce1 == ARR_MAX_END) ? LARGER : SMALLER;
      }
      // Non of the arcs are vertical:
      if (ce1 == ce2) return EQUAL;
      if (ce1 == ARR_MIN_END) return SMALLER;
      return LARGER;
    }
  };

  Compare_x_near_boundary_2 compare_x_near_boundary_2_object() const
  { return Compare_x_near_boundary_2(this); }
    

  class Compare_y_near_boundary_2 {
  protected:
    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;

    const Traits * m_traits;
    
    Compare_y_near_boundary_2(const Traits * traits) : m_traits(traits) {}

    friend class Arr_geodesic_arc_on_sphere_traits_2<Kernel>;
    
  public:
    Comparison_result operator()(const X_monotone_curve_2 & xcv1,
                                 const X_monotone_curve_2 & xcv2,
                                 Arr_curve_end ce) const
    {
      CGAL_precondition(!xcv1.is_degenerate());
      CGAL_precondition(!xcv2.is_degenerate());

      const Point_2 & l1 = xcv1.left();
      const Point_2 & r1 = xcv1.right();
      const Point_2 & l2 = xcv2.left();
      const Point_2 & r2 = xcv2.right();

      // If xcv1 is vertical, xcv1 coincides with the discontinuity arc:
      if (xcv1.is_vertical()) {
        CGAL_precondition(!l1.is_no_boundary());
        CGAL_precondition(!r1.is_no_boundary());
      }

      // If xcv2 is vertical, xcv2 coincides with the discontinuity arc:
      if (xcv2.is_vertical()) {
        CGAL_precondition(!l2.is_no_boundary());
        CGAL_precondition(!r2.is_no_boundary());
      }

      if (ce == ARR_MIN_END) {
        // Handle the south pole. It has the smallest y coords:
        if (l1.is_min_boundary())
          return (l2.is_min_boundary()) ? EQUAL : SMALLER;
        if (l2.is_min_boundary()) return LARGER;

        // None of xcv1 and xcv2 endpoints coincide with a pole:
        Comparison_result cr = m_traits->compare_y(l1, l2);
        if (cr != EQUAL) return cr;

        // If Both arcs are vertical, they overlap:
        if (xcv1.is_vertical() && xcv2.is_vertical()) return EQUAL;
        if (xcv1.is_vertical()) return LARGER;
        if (xcv2.is_vertical()) return SMALLER;

        // Non of the arcs is verticel. Thus, non of the endpoints coincide
        // with a pole.
        // Compare the y-coord. at the x-coord of the most left right-endpoint.
        CGAL_assertion(r1.is_no_boundary());
        CGAL_assertion(r2.is_no_boundary());
        
        if (m_traits->compare_xy(r1, r2) == LARGER) {
          // use r2 and xcv1:
          Oriented_side os =
            m_traits->oriented_side(xcv1.normal(), r2);
          return (os == ON_ORIENTED_BOUNDARY) ? EQUAL :
            (xcv1.is_directed_right()) ?
            ((os == ON_NEGATIVE_SIDE) ? LARGER : SMALLER) : 
            ((os == ON_NEGATIVE_SIDE) ? SMALLER : LARGER);
        }
        // use r1 and xcv2: 
        Oriented_side os = m_traits->oriented_side(xcv2.normal(), r1);
        return (os == ON_ORIENTED_BOUNDARY) ? EQUAL :
          (xcv2.is_directed_right()) ?
          ((os == ON_NEGATIVE_SIDE) ? SMALLER : LARGER) : 
          ((os == ON_NEGATIVE_SIDE) ? LARGER : SMALLER);
      }

      // ce == ARR_MAX_END
      
      // Handle the north pole. It has the largest y coords:
      if (r1.is_max_boundary()) return (r2.is_max_boundary()) ? EQUAL : LARGER;
      if (r2.is_max_boundary()) return SMALLER;

      // None of xcv1 and xcv2 endpoints coincide with a pole:
      Direction_2 r1_xy = Traits::project_xy(r1);
      Comparison_result cr = m_traits->compare_y(r1, r2);
      if (cr != EQUAL) return cr;

      // If Both arcs are vertical, they overlap:
      if (xcv1.is_vertical() && xcv2.is_vertical()) return EQUAL;
      if (xcv1.is_vertical()) return LARGER;
      if (xcv2.is_vertical()) return SMALLER;
        
      // Compare to the left:
      Direction_2 p_r1 = Traits::project_xy(r1);
      cr = m_traits->compare_y(r1, r2);
      if (cr != EQUAL) return cr;

      // Non of the arcs is verticel. Thus, non of the endpoints coincide with
      // a pole.
      // Compare the y-coord. at the x-coord of the most right left-endpoint.
      CGAL_assertion(l1.is_no_boundary());
      CGAL_assertion(l2.is_no_boundary());

      if (m_traits->compare_xy(l1, l2) == SMALLER) {
        // use l2 and xcv1:
        Oriented_side os = m_traits->oriented_side(xcv1.normal(), l2);
        return (os == ON_ORIENTED_BOUNDARY) ? EQUAL :
          (xcv1.is_directed_right()) ?
          ((os == ON_NEGATIVE_SIDE) ? LARGER : SMALLER) : 
          ((os == ON_NEGATIVE_SIDE) ? SMALLER : LARGER);
      }
      // use l1 and xcv2: 
      Oriented_side os = m_traits->oriented_side(xcv2.normal(), l1);
      return (os == ON_ORIENTED_BOUNDARY) ? EQUAL :
        (xcv2.is_directed_right()) ?
        ((os == ON_NEGATIVE_SIDE) ? SMALLER : LARGER) : 
        ((os == ON_NEGATIVE_SIDE) ? LARGER : SMALLER);
    }
  };

  Compare_y_near_boundary_2 compare_y_near_boundary_2_object() const
  { return Compare_y_near_boundary_2(this); }
  
  class Is_on_x_identification_2 {
  protected:
    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;
    
  public:
    bool operator()(const Point_2 & p) const
    {
      CGAL_error_msg("There is no horizontal identification arc!");
      return false;
    }

    bool operator()(const X_monotone_curve_2 & xcv) const
    {
      CGAL_error_msg("There is no horizontal identification arc!");
      return false;
    }
  };
  
  Is_on_x_identification_2 is_on_x_identification_2_object() const
  { return Is_on_x_identification_2(); }
  
  class Is_on_y_identification_2 {
  protected:
    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;

    const Traits * m_traits;
    
    Is_on_y_identification_2(const Traits * traits) : m_traits(traits) {}

    friend class Arr_geodesic_arc_on_sphere_traits_2<Kernel>;    
    
  public:
    bool operator()(const Point_2 & p) const { return p.is_mid_boundary(); }

    bool operator()(const X_monotone_curve_2 & xcv) const
    { return xcv.is_on_boundary(); }
  };
  
  Is_on_y_identification_2 is_on_y_identification_2_object() const
  { return Is_on_y_identification_2(this); }
  
  class Compare_x_on_boundary_2 {
  public:
    Comparison_result operator()(const Point_2 & p1, const Point_2 & p2) const
    {
      CGAL_error_msg("There is no horizontal identification arc!");
      return SMALLER;
    }
  };

  Compare_x_on_boundary_2 compare_x_on_boundary_2_object() const
  { return Compare_x_on_boundary_2(); }

  class Compare_y_on_boundary_2 {
  protected:
    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;

    const Traits * m_traits;
    
    Compare_y_on_boundary_2(const Traits * traits) : m_traits(traits) {}

    friend class Arr_geodesic_arc_on_sphere_traits_2<Kernel>;
    
  public:
    Comparison_result operator()(const Point_2 & p1, const Point_2 & p2) const
    {
      CGAL_precondition(!p1.is_no_boundary());
      CGAL_precondition(!p2.is_no_boundary());
      return m_traits->compare_y(p1, p2);
    }
  };

  Compare_y_on_boundary_2 compare_y_on_boundary_2_object() const
  { return Compare_y_on_boundary_2(this); }



  class Make_x_monotone_2 {
  protected:
    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;

    const Traits * m_traits;

    Make_x_monotone_2(const Traits * traits) : m_traits(traits) {}

    friend class Arr_geodesic_arc_on_sphere_traits_2<Kernel>;
    
  public:    
    template<typename OutputIterator>
    OutputIterator operator()(const Curve_2 & c, OutputIterator oi) const
    {
      if (c.is_degenerate()) {
        // The spherical_arc is a degenerate point - wrap it with an object:
        *oi++ = make_object(c.right());
        return oi;
      }

      if (c.is_x_monotone()) {
        // The spherical arc is monotone - wrap it with an object:
        // *oi++ = make_object(X_monotone_curve_2(c));
        const X_monotone_curve_2 * xc = &c;
        *oi++ = make_object(*xc);
        return oi;
      }

      if (c.is_full()) {
        // The spherical arc is full
        if (c.is_vertical()) {
          // The arc is vertical => divide it into 2 half arcs;
          const Direction_3 & np = m_traits->neg_pole();
          const Direction_3 & pp = m_traits->pos_pole();
          X_monotone_curve_2 xc1(np, pp, c.normal(), true, true);
          X_monotone_curve_2 xc2(pp, np, c.normal(), true, false);
          *oi++ = make_object(xc1);
          *oi++ = make_object(xc2);
          return oi;
        }
#if defined(CGAL_FULL_X_MONOTONE_GEODESIC_ARC_ON_SPHERE_IS_SUPPORTED)
        // The arc is not vertical => break it at the discontinuity arc:
        const X_monotone_curve_2 xc(c.normal());
        *oi++ = make_object(xc);
#else
        // Full x-monotone arcs are not supported!
        // Split the arc at the intersection point with the complement of the
        // discontinuity arc:
        Direction_3 normal = c.normal();
        bool directed_right = Traits::x_sign(normal) == POSITIVE;
        Direction_3 d1(-(normal.dz()), 0, normal.dx());
        Direction_3 d2(normal.dz(), 0, -(normal.dx()));
        X_monotone_curve_2 xc1(d1, d2, normal, false, directed_right);
        X_monotone_curve_2 xc2(d2, d1, normal, false, directed_right);
        *oi++ = make_object(xc1);
        *oi++ = make_object(xc2);
#endif
        return oi;
      }
      
      const Point_2 & source = c.source();
      const Point_2 & target = c.target();
      const Direction_3 & normal = c.normal();

      if (c.is_vertical()) {
        /* If one of the endpoints coincide with a pole, divide the arc at
         * the opposite pole:
         */
        const Direction_3 & np = m_traits->neg_pole();
        const Direction_3 & pp = m_traits->pos_pole();
        if (source.is_min_boundary() || target.is_min_boundary()) {
          X_monotone_curve_2 xc1(source, pp, normal, true, true);
          X_monotone_curve_2 xc2(pp, target, normal, true, false);
          *oi++ = make_object(xc1);
          *oi++ = make_object(xc2);
          return oi;
        }

        if (source.is_max_boundary() || target.is_max_boundary()) {
          X_monotone_curve_2 xc1(source, np, normal, true, false);
          X_monotone_curve_2 xc2(np, target, normal, true, true);
          *oi++ = make_object(xc1);
          *oi++ = make_object(xc2);
          return oi;
        }

        // None of the enpoints coincide with a pole.
        bool s_is_positive, t_is_positive, plane_is_positive;
        CGAL::Sign xsign = Traits::x_sign(normal);
        if (xsign == ZERO) {
          s_is_positive = Traits::x_sign(source) == POSITIVE;
          t_is_positive = Traits::x_sign(target) == POSITIVE;
          plane_is_positive = Traits::y_sign(normal) == NEGATIVE;
        } else {
          s_is_positive = Traits::y_sign(source) == POSITIVE;
          t_is_positive = Traits::y_sign(target) == POSITIVE;
          plane_is_positive = xsign == POSITIVE;
        }
        bool ccw = ((plane_is_positive && s_is_positive) ||
                    (!plane_is_positive && !s_is_positive));
        const Point_2 & pole1 = (ccw) ? pp : np;
        X_monotone_curve_2 xc1(source, pole1, normal, true, ccw);
        *oi++ = make_object(xc1);
        if (s_is_positive != t_is_positive) {
          // Construct 1 more arc:
          X_monotone_curve_2 xc2(pole1, target, normal, true, !ccw);
          *oi++ = make_object(xc2);
          return oi;
        }
        // Construct 2 more arcs:
        const Point_2 & pole2 = (ccw) ? np : pp;
        X_monotone_curve_2 xc2(pole1, pole2, normal, true, !ccw);
        *oi++ = make_object(xc2);
        X_monotone_curve_2 xc3(pole2, target, normal, true, ccw);
        *oi++ = make_object(xc3);
        return oi;
      }

      // The curve is not vertical, (none of the enpoints coincide with a pole)
#if (CGAL_IDENTIFICATION_XY == CGAL_X_MINUS_1_Y_0)
      Direction_3 dp = (CGAL::sign(normal.dz()) == POSITIVE) ?
        Direction_3(-(normal.dz()), 0, normal.dx()) :
        Direction_3(normal.dz(), 0, -(normal.dx()));
#else
      const Direction_2 & xy = Traits::identification_xy();
      FT x = xy.dx();
      FT y = xy.dy();
      FT z((xy.dx() * normal.dx() + xy.dy() * normal.dy()) /  -(normal.dz()));
      Direction_3 dp(x, y, z);
#endif
      Point_2 p(dp, Point_2::MID_BOUNDARY_LOC);
      Direction_2 s = Traits::project_xy(source);
      Direction_2 t = Traits::project_xy(target);
      const Direction_2 & d = Traits::identification_xy();
      const Kernel * kernel = m_traits;
      bool directed_right =
        kernel->counterclockwise_in_between_2_object()(d, s, t);
      
      X_monotone_curve_2 xc1(source, p, normal, false, directed_right);
      X_monotone_curve_2 xc2(p, target, normal, false, directed_right);
      *oi++ = make_object(xc1);
      *oi++ = make_object(xc2);
      return oi;
    }
  };

  Make_x_monotone_2 make_x_monotone_2_object() const
  { return Make_x_monotone_2(this); }

  class Split_2 {
  protected:
    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;

    const Traits * m_traits;
    
    Split_2(const Traits * traits) : m_traits(traits) {}

    friend class Arr_geodesic_arc_on_sphere_traits_2<Kernel>;
    
  public:
    void operator()(const X_monotone_curve_2 & xc, const Point_2 & p,
                    X_monotone_curve_2 & xc1, X_monotone_curve_2 & xc2) const
    {
      CGAL_precondition(!xc.is_degenerate());
      const Point_2 & source = xc.source();
      const Point_2 & target = xc.target();
      CGAL_precondition_code(const Kernel * kernel = m_traits);
      CGAL_precondition_code
        (typename Kernel::Equal_3 equal_3 = kernel->equal_3_object());
      CGAL_precondition(!equal_3(p, source));
      CGAL_precondition(!equal_3(p, target));

      xc1.normal(xc.normal());
      xc1.set_is_vertical(xc.is_vertical());
      xc1.set_is_degenerate(false);
      xc1.set_is_empty(false);

      xc2.normal(xc.normal());
      xc2.set_is_vertical(xc.is_vertical());
      xc2.set_is_empty(false);
      
      if (xc.is_directed_right()) {
        xc1.set_source(source);
        xc1.set_target(p);
        xc1.set_is_directed_right(true);
        xc2.set_source(p);
        xc2.set_target(target);
        xc2.set_is_directed_right(true);
      } else {
        xc1.set_source(p);
        xc1.set_target(target);
        xc1.set_is_directed_right(false);
        xc2.set_source(source);
        xc2.set_target(p);
        xc2.set_is_directed_right(false);
      }
    }
  };

  Split_2 split_2_object() const { return Split_2(this); }

  class Clockwise_in_between_2 {
  protected:
    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;

    const Traits * m_traits;
    
    Clockwise_in_between_2(const Traits * traits) : m_traits(traits) {}

    friend class Arr_geodesic_arc_on_sphere_traits_2<Kernel>;

  public:
    bool operator()(const Direction_2 & d,
                    const Direction_2 & d1, const Direction_2 & d2) const
    {
      const Kernel * kernel = m_traits;
      return kernel->counterclockwise_in_between_2_object()(d, d2, d1);
    }
  };

  Clockwise_in_between_2 clockwise_in_between_2_object() const
  { return Clockwise_in_between_2(this); }
  
  class Intersect_2 {
  private:
    
    template <typename In_between, typename OutputIterator>
    OutputIterator compute_intersection(const Point_2 & l1_3,
                                        const Point_2 r1_3,
                                        const Point_2 & l2_3,
                                        const Point_2 r2_3,
                                        const Direction_3 & normal,
                                        bool vertical,
                                        const Direction_2 & start,
                                        const In_between & in_between,
                                        Project project,
                                        OutputIterator oi) const
    {
      typedef std::pair<Point_2,unsigned int>                   Point_2_pair;
      const Kernel * kernel = m_traits;
      typename Kernel::Equal_2 equal = kernel->equal_2_object();

      Direction_2 l1 = project(l1_3);
      Direction_2 r1 = project(r1_3);
      Direction_2 l2 = project(l2_3);
      Direction_2 r2 = project(r2_3);

      if (equal(l1, l2)) {
        const Point_2 & trg = (in_between(r1, l2, r2)) ? r1_3 : r2_3;
        X_monotone_curve_2 xc(l1_3, trg, normal, vertical, true);
        *oi++ = make_object(xc);
        return oi;
      }

      bool l1_eq_start = equal(l1, start);
      bool l2_eq_start = equal(l2, start);

      if (l1_eq_start || (!l2_eq_start && in_between(l1, start, l2))) {
        // The following applies only to full circles:
        if (l1_eq_start && equal(r2, start))
          *oi++ = make_object(Point_2_pair(r2_3, 1));
        if (in_between(r1, l1, l2)) return oi;      // no intersection
        if (equal(r1, l2)) {
          *oi++ = make_object(Point_2_pair(r1_3, 1));
          return oi;
        }
        const Point_2 & trg = (in_between(r1, l2, r2)) ? r1_3 : r2_3;
        X_monotone_curve_2 xc(l2_3, trg, normal, vertical, true);
        *oi++ = make_object(xc);
        return oi;
      }
      CGAL_assertion(l2_eq_start || in_between(l2, start, l1));
      // The following applies only to full circles:
      if (l2_eq_start && equal(r1, start))
        *oi++ = make_object(Point_2_pair(r1_3, 1));
      if (in_between(r2, l2, l1)) return oi;      // no intersection
      if (equal(r2, l1)) {
        *oi++ = make_object(Point_2_pair(r2_3, 1));
        return oi;
      }
      const Point_2 & trg = (in_between(r1, l2, r2)) ? r1_3 : r2_3;
      X_monotone_curve_2 xc(l1_3, trg, normal, vertical, true);
      *oi++ = make_object(xc);
      return oi;
    }
    
    bool is_in_between(const Point_2 & point,
                       const X_monotone_curve_2 & xc) const
    {
      const Kernel * kernel = m_traits;
      CGAL_precondition(m_traits->has_on(xc.normal(), point));
      
      const Point_2 & left = xc.left();
      const Point_2 & right = xc.right();

      // Handle the poles:
      if (point.is_max_boundary()) return (right.is_max_boundary());
      if (point.is_min_boundary()) return (left.is_min_boundary());

      if (xc.is_vertical()) {
        // Compare the x coordinates. If they are not equal, return false:
        Direction_3 normal = xc.normal();
        bool plane_is_positive, p_is_positive;
        CGAL::Sign xsign = Traits::x_sign(normal);
        if (xsign == ZERO) {
          plane_is_positive = Traits::y_sign(normal) == NEGATIVE;
          p_is_positive = Traits::x_sign(point) == POSITIVE;
        } else {
          plane_is_positive = xsign == POSITIVE;
          p_is_positive = Traits::y_sign(point) == POSITIVE;
        }
        
        bool xc_is_positive = ((plane_is_positive && xc.is_directed_right()) ||
                               (!plane_is_positive && !xc.is_directed_right()));

        if ((xc_is_positive && !p_is_positive) ||
            (!xc_is_positive && p_is_positive))
          return false;

        // Compare the y-coords:
        return (((left.is_min_boundary()) ||
                 (m_traits->compare_y(point, left) != SMALLER)) &&
                ((right.is_max_boundary()) ||
                 (m_traits->compare_y(point, right) != LARGER)));
      }

      // The arc is not vertical. Compare the projections onto the xy-plane:
      typename Kernel::Equal_2 equal_2 = kernel->equal_2_object(); 
      Direction_2 p = Traits::project_xy(point);
      Direction_2 r = Traits::project_xy(right);
      if (equal_2(p, r)) return true;
      Direction_2 l = Traits::project_xy(left);
      if (equal_2(p, l)) return true;
      return kernel->counterclockwise_in_between_2_object()(p, l, r);
    }

  protected:
    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;

    const Traits * m_traits;
    
    Intersect_2(const Traits * traits) : m_traits(traits) {}

    friend class Arr_geodesic_arc_on_sphere_traits_2<Kernel>;
    
  public:
    template<typename OutputIterator>
    OutputIterator operator()(const X_monotone_curve_2 & xc1,
                              const X_monotone_curve_2 & xc2,
                              OutputIterator oi) const
    {
      CGAL_precondition(!xc1.is_degenerate());
      CGAL_precondition(!xc2.is_degenerate());

      typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;
      typedef typename Kernel::Equal_2                          Equal_2;
      typedef typename Kernel::Counterclockwise_in_between_2
        Counterclockwise_in_between_2;
      typedef typename Traits::Clockwise_in_between_2
        Clockwise_in_between_2;
      typedef typename Kernel::Equal_3                          Equal_3;
        
      typedef std::pair<Point_2,unsigned int>         Point_2_pair;
      const Kernel * kernel = m_traits;

      Equal_3 equal_3 = kernel->equal_3_object();
      const Direction_3 & normal1 = xc1.normal();
      const Direction_3 & normal2 = xc2.normal();

      Direction_3 opposite_normal1 = 
        kernel->construct_opposite_direction_3_object()(normal1);

      if (equal_3(normal1, normal2) || equal_3(opposite_normal1, normal2)) {
        // The underlying planes are the same
        Equal_2 equal = kernel->equal_2_object();
        Counterclockwise_in_between_2 ccib =
          kernel->counterclockwise_in_between_2_object();
        typename Traits::Clockwise_in_between_2 cib =
          m_traits->clockwise_in_between_2_object();

        if (xc1.is_vertical()) {
          // Both arcs are vertical
          bool res = kernel->equal_3_object()(normal1, normal2);
          if ((!res && (xc1.is_directed_right() == xc2.is_directed_right())) ||
              (res && (xc1.is_directed_right() != xc2.is_directed_right())))
          {
            if (xc1.left().is_min_boundary() && xc2.left().is_min_boundary())
              *oi++ = make_object(Point_2_pair(xc1.left(), 1));
            if (xc1.right().is_max_boundary() && xc2.right().is_max_boundary())
              *oi++ = make_object(Point_2_pair(xc1.right(), 1));
            return oi;
          }
            
          if (xc1.left().is_min_boundary() && xc1.right().is_max_boundary()) {
            *oi++ = make_object(xc2);
            return oi;
          }
          if (xc2.left().is_min_boundary() && xc2.right().is_max_boundary()) {
            *oi++ = make_object(xc1);
            return oi;
          }
          const Point_2 & point =
            xc1.left().is_min_boundary() ? xc1.right() : xc1.left();

          CGAL::Sign xsign = Traits::x_sign(normal1);
          bool xz_plane = xsign == ZERO;
          Project project =
            (xz_plane) ? Traits::project_xz : Traits::project_yz;

          Direction_3 normal = (xc1.is_directed_right()) ?
            normal1 : opposite_normal1;
          
          bool p_x_is_positive = Traits::x_sign(point) == POSITIVE;
          bool p_y_is_positive = Traits::y_sign(point) == POSITIVE;

          if ((xz_plane && p_x_is_positive) || (!xz_plane && p_y_is_positive)) {
            // The endpoints reside in the positive x-halfspace:
            return compute_intersection(xc1.left(), xc1.right(),
                                        xc2.left(), xc2.right(),
                                        normal, true, Traits::neg_y_2(),
                                        ccib, project, oi);
          }
          // The endpoints reside in the negative x-halfspace:
          return compute_intersection(xc1.left(), xc1.right(),
                                      xc2.left(), xc2.right(),
                                      normal, true, Traits::neg_y_2(),
                                      cib, project, oi);
        }

        // The arcs are not vertical:
        bool plane_is_positive = (Traits::z_sign(normal1) == POSITIVE);
        Direction_3 normal =
          (plane_is_positive) ? normal1 : opposite_normal1;
        return compute_intersection(xc1.left(), xc1.right(),
                                    xc2.left(), xc2.right(),
                                    normal, false, Traits::identification_xy(),
                                    ccib, Traits::project_xy, oi);
      }

      Vector_3 v =
        kernel->
        construct_cross_product_vector_3_object()(xc1.normal().vector(),
                                                  xc2.normal().vector());

      // Determine which one of the two directions:
      Point_2 ed(v.direction());
      if (is_in_between(ed, xc1) && is_in_between(ed, xc2)) {
        *oi++ = make_object(Point_2_pair(ed, 1));
        return oi;
      }
      
      Vector_3 vo(kernel->construct_opposite_vector_3_object()(v));
      Point_2 edo(vo.direction());
      if (is_in_between(edo, xc1) && is_in_between(edo, xc2)) {
        *oi++ = make_object(Point_2_pair(edo, 1));
        return oi;
      }
      return oi;
    }
  };

  Intersect_2 intersect_2_object() const { return Intersect_2(this); }

  class Are_mergeable_2 {
    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;

    const Traits * m_traits;
    
    Are_mergeable_2(const Traits * traits) : m_traits(traits) {}

    friend class Arr_geodesic_arc_on_sphere_traits_2<Kernel>;

  public:
    bool operator()(const X_monotone_curve_2 & xc1,
                    const X_monotone_curve_2 & xc2) const
    {
      return false;
      
      if (xc1.is_empty() || xc2.is_empty()) return true;
      if (xc1.is_full() && xc2.is_full()) return false;

      const Kernel * kernel = m_traits;
      typename Kernel::Equal_3 equal = kernel->equal_3_object();

      if (xc1.is_degenerate() && xc2.is_degenerate())
        return equal(xc1.left(), xc2.left());
      if (xc1.is_full() && xc2.is_degenerate())
        return xc1.has_on(xc2.left());
      if (xc2.is_full() && xc1.is_degenerate())
        return xc2.has_on(xc1.left());

      const Direction_3 & normal1 = xc1.normal();
      const Direction_3 & normal2 = xc2.normal();
      Direction_3 opposite_normal1 = 
        kernel->construct_opposite_direction_3_object()(normal1);
      if (!equal(normal1, normal2) && !equal(opposite_normal1, normal2))
        return false;

      bool eq1 = equal(xc1.right(), xc2.left());
      bool eq2 = equal(xc1.left(), xc2.right());

#if defined(CGAL_FULL_X_MONOTONE_GEODESIC_ARC_ON_SPHERE_IS_SUPPORTED)
      if (eq1 && eq2) return true;
#else
      if (eq1 && eq2) return false;
#endif

      if (eq1 && xc2.left().is_no_boundary()) return true;
      if (eq2 && xc1.left().is_no_boundary()) return true;
      return false;
    }
  };

  Are_mergeable_2 are_mergeable_2_object() const
  { return Are_mergeable_2(this); }

  class Merge_2 {
  protected:
    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;

    const Traits * m_traits;

    Merge_2(const Traits * traits) : m_traits(traits) {}

    friend class Arr_geodesic_arc_on_sphere_traits_2<Kernel>;
    
  public:
    void operator()(const X_monotone_curve_2 & xc1,
                    const X_monotone_curve_2 & xc2,
                    X_monotone_curve_2 & xc) const
    {
      CGAL_precondition (m_traits->are_mergeable_2_object()(xc1, xc2) == true);

      if (xc1.is_degenerate() || xc1.is_empty()) {
        xc = xc2;
        return;
      }

      if (xc2.is_degenerate() || xc2.is_empty()) {
        xc = xc1;
        return;
      }
      
      const Kernel * kernel = m_traits;
      typename Kernel::Equal_3 equal = kernel->equal_3_object();

      xc.set_is_degenerate(false);
      xc.set_is_empty(false);
      xc.set_is_vertical(xc1.is_vertical());
      
      bool eq1 = equal(xc1.right(), xc2.left());

#if defined(CGAL_FULL_X_MONOTONE_GEODESIC_ARC_ON_SPHERE_IS_SUPPORTED)
      bool eq2 = equal(xc1.left(), xc2.right());
      if (eq1 && eq2) {
        const Point_2 & p =
          xc1.source().is_mid_boundary() ? xc1.source() : xc1.target();
        xc.set_source(p);
        xc.set_target(p);
        xc.normal(xc1.normal());
        xc.set_is_full(true);
      }
#endif
      
      if (xc1.is_directed_right() || xc2.is_directed_right()) {
        xc.normal(xc1.is_directed_right() ?
                               xc1.normal() : xc2.normal());
        xc.set_is_directed_right(true);

        if (eq1) {
          xc.set_source(xc1.left());
          xc.set_target(xc2.right());
        } else {
          CGAL_assertion(equal(xc1.left(), xc2.right()));
          xc.set_source(xc2.left());
          xc.set_target(xc1.right());
        }
      } else {
        xc.normal(xc1.normal());
        xc.set_is_directed_right(false);

        if (eq1) {
          xc.set_source(xc2.right());
          xc.set_target(xc1.left());
        } else {
          CGAL_assertion(equal(xc1.left(), xc2.right()));
          xc.set_source(xc1.right());
          xc.set_target(xc2.left());
        }
      }
    }      
  };

  Merge_2 merge_2_object() const { return Merge_2(this); }


  typedef double                          Approximate_number_type;

  class Approximate_2 {
  public:

    Approximate_number_type operator()(const Point_2 & p, int i) const
    {
      CGAL_precondition(i == 0 || i == 1);
      return (i == 0) ? CGAL::to_double(p.x()) : CGAL::to_double(p.y());
    }
  };

  Approximate_2 approximate_2_object() const { return Approximate_2(); }

  class Construct_x_monotone_curve_2 {
  public:

    X_monotone_curve_2 operator()(const Point_2 & p, const Point_2 & q) const
    { return X_monotone_curve_2(p, q); }
  };

  Construct_x_monotone_curve_2 construct_x_monotone_curve_2_object() const
  { return Construct_x_monotone_curve_2(); }



 
  class Compare_endpoints_xy_2 {
  public:

    Comparison_result operator()(const X_monotone_curve_2 & xc)
    { return (xc.is_directed_right()) ? SMALLER : LARGER; }
  };

  Compare_endpoints_xy_2 compare_endpoints_xy_2_object() const
  { return Compare_endpoints_xy_2(); }

  class Construct_opposite_2 {
  public:
    X_monotone_curve_2 operator()(const X_monotone_curve_2 & xc)
    { return xc.opposite(); }
  };

  Construct_opposite_2 construct_opposite_2_object() const
  { return Construct_opposite_2(); }

#if 0

  template <typename OutputStream>
  friend OutputStream & operator<<(OutputStream & os, const Point_2 & p)
  {
    CGAL::To_double<typename Kernel::FT> todouble;
    os << static_cast<float>(todouble(p.dx())) << ", "
       << static_cast<float>(todouble(p.dy())) << ", "
       << static_cast<float>(todouble(p.dz()));
    return os;
  }

  template <typename OutputStream>
  friend OutputStream & operator<<(OutputStream & os,
                                   const X_monotone_curve_2 & xc)
  {
    os << "(" << xc.left() << "), (" << xc.right() << ")";
    return os;
  }

  template <typename InputStream>
  friend InputStream & operator>>(InputStream & is, X_monotone_curve_2 & arc)
  {
    std::cerr << "Not implemented yet!" << std::endl;
    return is;
  }  
#endif
};

template <typename Kernel>
class Arr_extended_direction_3 : public Kernel::Direction_3 {
public:
  typedef typename Kernel::FT                           FT;
  typedef typename Kernel::Direction_3                  Direction_3;

  enum Location_type {
    NO_BOUNDARY_LOC,
    MIN_BOUNDARY_LOC,
    MID_BOUNDARY_LOC,
    MAX_BOUNDARY_LOC
  };

private:
  typedef typename Kernel::Direction_2                  Direction_2;

  Location_type m_location;

  inline Sign x_sign(Direction_3 d) const { return CGAL::sign(d.dx()); }

  inline Sign y_sign(Direction_3 d) const { return CGAL::sign(d.dy()); }  

  inline Sign z_sign(Direction_3 d) const { return CGAL::sign(d.dz()); }
  
public:
  Arr_extended_direction_3() : 
    Direction_3(0, 0, 1),
    m_location(MAX_BOUNDARY_LOC)
  {}
      
  Arr_extended_direction_3(const Direction_3 & dir, Location_type location) :
    Direction_3(dir),
    m_location(location)
  {}

  Arr_extended_direction_3(const FT & x, const FT & y, const FT & z) :
    Direction_3(x, y, z)
  { init(Direction_3(x, y, z)); }
  
  Arr_extended_direction_3(const Direction_3 & dir) : Direction_3(dir)
  { init(dir); }

  void init(const Direction_3 & dir)
  {
#if (CGAL_IDENTIFICATION_XY == CGAL_X_MINUS_1_Y_0)
    if (y_sign(dir) != ZERO) {
      m_location = NO_BOUNDARY_LOC;
      return;
    }
    CGAL::Sign signx = x_sign(dir);
    m_location =
      (signx == POSITIVE) ? NO_BOUNDARY_LOC :
      ((signx == NEGATIVE) ? MID_BOUNDARY_LOC :
       ((z_sign(dir) == NEGATIVE) ? MIN_BOUNDARY_LOC : MAX_BOUNDARY_LOC));
#else
    if ((x_sign(dir) == ZERO) && (y_sign(dir) == ZERO)) {
      m_location =
        (z_sign(dir) == NEGATIVE) ? MIN_BOUNDARY_LOC : MAX_BOUNDARY_LOC;
      return;
    }

    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;
    Direction_2 dir_xy = Traits::project_xy(dir);
    Kernel kernel;
    typename Kernel::Equal_2 equal_2 = kernel.equal_2_object();
    const Direction_2 & xy = Traits::identification_xy();
    m_location = equal_2(dir_xy, xy) ? MID_BOUNDARY_LOC : NO_BOUNDARY_LOC;
#endif
  }
  
  Arr_extended_direction_3(const Arr_extended_direction_3 & ed) :
    Direction_3(static_cast<const Direction_3&>(ed))
  {
    m_location = ed.discontinuity_type();
  }

  Arr_extended_direction_3 & operator=(const Arr_extended_direction_3 & ed)
  {
    *(static_cast<Direction_3*>(this)) = static_cast<const Direction_3&>(ed);
    m_location = ed.discontinuity_type();
    return (*this);
  }

  Location_type discontinuity_type() const
  { return m_location; }
  
  bool is_no_boundary() const { return (m_location == NO_BOUNDARY_LOC); }
  
  bool is_min_boundary() const { return (m_location == MIN_BOUNDARY_LOC); }
  
  bool is_mid_boundary() const { return (m_location == MID_BOUNDARY_LOC); }
  
  bool is_max_boundary() const { return (m_location == MAX_BOUNDARY_LOC); }
};

template <typename T_Kernel>
class Arr_x_monotone_geodesic_arc_on_sphere_3 {
protected:
  typedef T_Kernel                                    Kernel;
  
  typedef typename Kernel::Plane_3                    Plane_3;
  typedef typename Kernel::Vector_3                   Vector_3;
  typedef typename Kernel::Direction_2                Direction_2;
  typedef typename Kernel::Direction_3                Direction_3;

  // For some reason compilation under Windows fails without the qualifier
  typedef CGAL::Arr_extended_direction_3<Kernel>      Arr_extended_direction_3;

  Arr_extended_direction_3 m_source;

  Arr_extended_direction_3 m_target;

  Direction_3 m_normal;

  bool m_is_vertical;

  bool m_is_directed_right;

  bool m_is_full;

  /* The arc is degenerate - it consists of a single point */
  bool m_is_degenerate;

  bool m_is_empty;
  
  inline Sign x_sign(Direction_3 d) const { return CGAL::sign(d.dx()); }

  inline Sign y_sign(Direction_3 d) const { return CGAL::sign(d.dy()); }  

  inline Sign z_sign(Direction_3 d) const { return CGAL::sign(d.dz()); }

  inline Direction_3 construct_normal_3(const Direction_3 & d1,
                                        const Direction_3 & d2) const
  {
    Kernel kernel;
    Vector_3 v = kernel.construct_cross_product_vector_3_object()(d1.vector(),
                                                                  d2.vector());
    return v.direction();
  }

public:
  Arr_x_monotone_geodesic_arc_on_sphere_3() :
    m_is_vertical(false),
    m_is_directed_right(false),
    m_is_full(false),
    m_is_degenerate(false),
    m_is_empty(true)
  {}
  
  Arr_x_monotone_geodesic_arc_on_sphere_3
  (const Arr_extended_direction_3 & src,
   const Arr_extended_direction_3 & trg,
   const Direction_3 & normal,
   bool is_vertical, bool is_directed_right,
   bool is_full = false, bool is_degenerate = false, bool is_empty = false) :
    m_source(src),
    m_target(trg),
    m_normal(normal),
    m_is_vertical(is_vertical),
    m_is_directed_right(is_directed_right),
    m_is_full(is_full),
    m_is_degenerate(is_degenerate),
    m_is_empty(is_empty)
  {}

  Arr_x_monotone_geodesic_arc_on_sphere_3
  (const Arr_x_monotone_geodesic_arc_on_sphere_3 & other)
  {
    m_source = other.m_source;
    m_target = other.m_target;
    m_normal = other.m_normal;
    m_is_vertical = other.m_is_vertical;
    m_is_directed_right = other.m_is_directed_right;
    m_is_full = other.m_is_full;
    m_is_degenerate = other.m_is_degenerate;
    m_is_empty = other.m_is_empty;
  }

  Arr_x_monotone_geodesic_arc_on_sphere_3 & operator=
  (const Arr_x_monotone_geodesic_arc_on_sphere_3 & other)
  {
    m_source = other.m_source;
    m_target = other.m_target;
    m_normal = other.m_normal;
    m_is_vertical = other.m_is_vertical;
    m_is_directed_right = other.m_is_directed_right;
    m_is_full = other.m_is_full;
    m_is_degenerate = other.m_is_degenerate;
    m_is_empty = other.m_is_empty;
    return (*this);
  }
  
  Arr_x_monotone_geodesic_arc_on_sphere_3
  (const Arr_extended_direction_3 & source,
   const Arr_extended_direction_3 & target) :
    m_source(source),
    m_target(target),
    m_is_full(false),
    m_is_degenerate(false),
    m_is_empty(false)
  {
    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;

    Kernel kernel;
    CGAL_precondition(!kernel.equal_3_object()(source, target));
    CGAL_precondition(!kernel.equal_3_object()
                      (kernel.construct_opposite_direction_3_object()(source),
                       target));
    m_normal = construct_normal_3(source, target);
      
    // Check whether any one of the endpoint coincide with a pole:
    if (source.is_max_boundary()) {
      set_is_vertical(true);
      set_is_directed_right(false);
      return;
    }
    if (source.is_min_boundary()) {
      set_is_vertical(true);
      set_is_directed_right(true);
      return;
    }
    if (target.is_max_boundary()) {
      set_is_vertical(true);
      set_is_directed_right(true);
      return;
    }
    if (target.is_min_boundary()) {
      set_is_vertical(true);
      set_is_directed_right(false);
      return;
    }

    // None of the enpoints coincide with a pole:
    typename Kernel::Equal_2 equal_2 = kernel.equal_2_object();
    Direction_2 s = Traits::project_xy(source);
    Direction_2 t = Traits::project_xy(target);

    Orientation orient = Traits::orientation(s, t);
    if (orient == COLLINEAR) {
      set_is_vertical(true);
      const Direction_2 & nx = Traits::neg_x_2();
      if (Traits::orientation(nx, s) == COLLINEAR) {
        // Project onto xz plane:
        s = Traits::project_xz(source);
        t = Traits::project_xz(target);
        const Direction_2 & ny = Traits::neg_y_2();
        Orientation orient1 = Traits::orientation(ny, s);
        CGAL_assertion_code(Orientation orient2 = Traits::orientation(ny, t));
        CGAL_assertion(orient1 == orient2);
        orient = Traits::orientation(s, t);
        CGAL_assertion(orient != COLLINEAR);
        if (orient1 == LEFT_TURN) {
          set_is_directed_right(orient == LEFT_TURN);
          return;
        }
        set_is_directed_right(orient == RIGHT_TURN);
        return;
      }
      // Project onto yz plane:
      s = Traits::project_yz(source);
      t = Traits::project_yz(target);
      const Direction_2 & ny = Traits::neg_y_2();
      Orientation orient1 = Traits::orientation(ny, s);
      CGAL_assertion_code(Orientation orient2 = Traits::orientation(ny, t));
      CGAL_assertion(orient1 == orient2);
      if (orient1 == LEFT_TURN) {
        orient = Traits::orientation(s, t);
        CGAL_assertion(orient != COLLINEAR);
        set_is_directed_right(orient == LEFT_TURN);
        return;
      }
      orient = Traits::orientation(s, t);
      CGAL_assertion(orient != COLLINEAR);
      set_is_directed_right(orient == RIGHT_TURN);
      return;
    }
    
    // The arc is not vertical!
    set_is_vertical(false);
    set_is_directed_right(orient == LEFT_TURN);
    set_is_full(kernel.equal_3_object()(source, target));
  }

  Arr_x_monotone_geodesic_arc_on_sphere_3(const Direction_3 & normal) :
    m_normal(normal),
    m_is_vertical(false),
    m_is_directed_right(z_sign(normal) == POSITIVE),
    m_is_full(true),
    m_is_degenerate(false),
    m_is_empty(false)
  {
    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;
    typedef typename Kernel::FT                         FT;
    
    CGAL_precondition(z_sign(normal) != ZERO);

#if (CGAL_IDENTIFICATION_XY == CGAL_X_MINUS_1_Y_0)    
    Direction_3 d = (CGAL::sign(normal.dz()) == POSITIVE) ?
      Direction_3(-(normal.dz()), 0, normal.dx()) :
      Direction_3(normal.dz(), 0, -(normal.dx()));
#else
    const Direction_2 & xy = Traits::identification_xy();
    FT x = xy.dx();
    FT y = xy.dy();
    FT z((xy.dx() * normal.dx() + xy.dy() * normal.dy()) / -(normal.dz()));
    Direction_3 d(x, y, z);
#endif
    m_source = m_target =
      Arr_extended_direction_3(d, Arr_extended_direction_3::MID_BOUNDARY_LOC);
  }

  Arr_x_monotone_geodesic_arc_on_sphere_3
  (const Arr_extended_direction_3 & point,
   const Direction_3 & normal) :
    m_source(point),
    m_target(point),
    m_normal(normal),
    m_is_vertical(false),
    m_is_directed_right(z_sign(normal) == POSITIVE),
    m_is_full(true),
    m_is_degenerate(false),
    m_is_empty(false)
  {
    CGAL_precondition(has_on(point));
    CGAL_precondition(z_sign(normal) != ZERO);
#if !defined(CGAL_FULL_X_MONOTONE_GEODESIC_ARC_ON_SPHERE_IS_SUPPORTED)
    CGAL_error_msg( "Full x-monotone arcs are not supported!");
#endif
  }
  
  Arr_x_monotone_geodesic_arc_on_sphere_3
  (const Arr_extended_direction_3 & source,
   const Arr_extended_direction_3 & target,
   const Direction_3 & normal) :
    m_source(source),
    m_target(target),
    m_normal(normal),
    m_is_full(false),
    m_is_degenerate(false),
    m_is_empty(false)
  {
    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;

    CGAL_precondition(has_on(source));
    CGAL_precondition(has_on(target));

    // Check whether any one of the endpoint coincide with a pole:
    if (source.is_max_boundary()) {
      set_is_vertical(true);
      set_is_directed_right(false);
      return;
    }
    if (source.is_min_boundary()) {
      set_is_vertical(true);
      set_is_directed_right(true);
      return;
    }
    if (target.is_max_boundary()) {
      set_is_vertical(true);
      set_is_directed_right(true);
      return;
    }
    if (target.is_min_boundary()) {
      set_is_vertical(true);
      set_is_directed_right(false);
      return;
    }

    if (z_sign(normal) == ZERO) {
      set_is_vertical(true);
      bool s_is_positive, plane_is_positive;
      CGAL::Sign xsign = x_sign(normal);
      if (xsign == ZERO) {
        s_is_positive = x_sign(source) == POSITIVE;
        plane_is_positive = y_sign(normal) == NEGATIVE;
      } else {
        s_is_positive = y_sign(source) == POSITIVE;
        plane_is_positive = xsign == POSITIVE;
      }
      bool ccw = ((plane_is_positive && s_is_positive) ||
                  (!plane_is_positive && !s_is_positive));
      set_is_directed_right(ccw);
      return;
    }
      
    // The arc is not vertical!
    set_is_vertical(false);
    set_is_directed_right(z_sign(normal) == POSITIVE);
  }

  void set_source(const Direction_3 & p) { m_source = p; }

  void set_target(const Direction_3 & p) { m_target = p; }

  void normal(const Direction_3 & normal) { m_normal = normal; }

  void set_is_vertical(bool flag) { m_is_vertical = flag; }
  void set_is_directed_right(bool flag) { m_is_directed_right = flag; }
  void set_is_full(bool flag) { m_is_full = flag; }
  void set_is_degenerate(bool flag) { m_is_degenerate = flag; }
  void set_is_empty(bool flag) { m_is_empty = flag; }

  const Arr_extended_direction_3 & source() const { return m_source; }

  const Arr_extended_direction_3 & target() const { return m_target; }
    
  const Direction_3 & normal() const { return m_normal; }

  const Arr_extended_direction_3 & left() const
  { return (m_is_directed_right ? m_source : m_target); }

  const Arr_extended_direction_3 & right() const
  { return (m_is_directed_right ? m_target : m_source); }

  bool is_vertical() const { return m_is_vertical; }

  bool is_directed_right() const { return m_is_directed_right; }

  bool is_full() const { return m_is_full; }

  bool is_degenerate() const { return m_is_degenerate; }

  bool is_empty() const { return m_is_empty; }
  
  bool is_on_boundary() const
  {
    /* If the curve is not vertical and non of its endpoints lie on the
     * boundary, the arc itself cannot lie on the identification arc.
     */
    if (m_source.is_no_boundary() || m_target.is_no_boundary() ||
        !is_vertical())
      return false;

    if (m_source.is_mid_boundary() || m_target.is_mid_boundary())
      return true;

    /* Both endpoints lie on opposite poles respectively. If the normal
     * coincides with the normal of the plane that contains the identification
     * arc, the arc lies on the identification arc.
     */
#if (CGAL_IDENTIFICATION_XY == CGAL_X_MINUS_1_Y_0)
    return ((x_sign(m_normal) == ZERO) &&
            (((y_sign(m_normal) == NEGATIVE) && !is_directed_right()) ||
             ((y_sign(m_normal) == POSITIVE) && is_directed_right())));
#else
    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;

    const Direction_3 & iden_normal = Traits::identification_normal();
    const Direction_2 iden_normal_xy = Traits::project_xy(iden_normal);
    Direction_2 normal_xy = Traits::project_xy(m_normal);
    Kernel kernel;
    if (is_directed_right()) {
      return kernel.equal_2_object()(normal_xy, iden_normal_xy);
    } else {
      Direction_2 opposite_normal_xy = 
        kernel.construct_opposite_direction_2_object()(normal_xy);
      return kernel.equal_2_object()(opposite_normal_xy, iden_normal_xy);
    }
#endif
  }
  
  bool is_in_x_range(const Arr_extended_direction_3 & point) const
  {
    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;

    CGAL_precondition(!point.is_min_boundary());
    CGAL_precondition(!point.is_max_boundary());
    
    Direction_2 p = Traits::project_xy(point);
    if (is_vertical()) {
      Direction_2 q = (is_directed_right()) ?
        Direction_2(-(m_normal.dy()), m_normal.dx()) :
        Direction_2(m_normal.dy(), -(m_normal.dx()));
      Kernel kernel;
      return kernel.equal_2_object()(p, q);
    }

    // The curve is not vertical:
    Direction_2 r = Traits::project_xy(right());
    Kernel kernel;
    if (kernel.equal_2_object()(p, r)) return true;
    Direction_2 l = Traits::project_xy(left());
    if (kernel.equal_2_object()(p, l)) return true;
    return kernel.counterclockwise_in_between_2_object()(p, l, r);
  }

#if 0

  Bbox_2 bbox() const
  {
    Kernel kernel;
    Segment_2 seg = kernel.construct_spherical_arc_2_object()(this->m_source,
                                                              this->m_target);
    return kernel.construct_bbox_2_object()(seg);
  }
#endif
  
  Arr_x_monotone_geodesic_arc_on_sphere_3 opposite() const
  {
    Arr_x_monotone_geodesic_arc_on_sphere_3 opp;
    opp.m_source = this->m_target;
    opp.m_target = this->m_source;
    opp.m_normal = this->m_normal;
    opp.m_is_directed_right = !(this->is_directed_right());
    opp.m_is_vertical = this->is_vertical();
    opp.m_is_full = this->is_full();
    opp.m_is_degenerate = this->is_degenerate();
    opp.m_is_empty = this->is_empty();
    return opp;
  }

  inline bool has_on(const Direction_3 & dir) const
  {
    typename Kernel::FT dot = normal().vector() * dir.vector();
    return CGAL::sign(dot) == ZERO;
  }
};

template <typename T_Kernel>
class Arr_geodesic_arc_on_sphere_3 :
  public Arr_x_monotone_geodesic_arc_on_sphere_3<T_Kernel> {
protected:
  typedef T_Kernel                                              Kernel;
  typedef Arr_x_monotone_geodesic_arc_on_sphere_3<Kernel> Base;

  typedef typename Base::Plane_3                                Plane_3;
  typedef typename Base::Direction_3                            Direction_3;
  typedef typename Base::Direction_2                            Direction_2;
  
  // For some reason compilation under Windows fails without the qualifier
  typedef CGAL::Arr_extended_direction_3<Kernel>    Arr_extended_direction_3;

  bool m_is_x_monotone;
  
public:
  Arr_geodesic_arc_on_sphere_3() : Base(), m_is_x_monotone(true) {}
  
  Arr_geodesic_arc_on_sphere_3
  (const Arr_geodesic_arc_on_sphere_3 & other) : Base(other)
  {
    m_is_x_monotone = other.m_is_x_monotone;
  }
  
  Arr_geodesic_arc_on_sphere_3(const Arr_extended_direction_3 & src,
                                     const Arr_extended_direction_3 & trg,
                                     const Direction_3 & normal,
                                     bool is_x_monotone, bool is_vertical,
                                     bool is_directed_right,
                                     bool is_full = false,
                                     bool is_degenerate = false,
                                     bool is_empty = false) :
    Base(src, trg, normal,
         is_vertical, is_directed_right, is_full, is_degenerate, is_empty),
    m_is_x_monotone(is_x_monotone)
  {
    CGAL_precondition(this->has_on(src));
    CGAL_precondition(this->has_on(trg));
  }

  Arr_geodesic_arc_on_sphere_3(const Arr_extended_direction_3 & source,
                               const Arr_extended_direction_3 & target)
  {
    this->set_source(source);
    this->set_target(target);
    this->set_is_full(false);
    this->set_is_degenerate(false);
    this->set_is_empty(false);

    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel>         Traits;
    typedef typename Kernel::Direction_2                        Direction_2;
    typedef typename Kernel::Direction_3                        Direction_3;

    Kernel kernel;
    CGAL_precondition(!kernel.equal_3_object()(source, target));
    CGAL_precondition(!kernel.equal_3_object()
                      (kernel.construct_opposite_direction_3_object()(source),
                       target));
    this->m_normal = this->construct_normal_3(source, target);

    // Check whether one of the endpoints coincides with a pole: */
    if (source.is_max_boundary()) {
      this->set_is_vertical(true);
      this->set_is_directed_right(false);
      set_is_x_monotone(true);
      return;
    }
    if (source.is_min_boundary()) {
      this->set_is_vertical(true);
      this->set_is_directed_right(true);
      set_is_x_monotone(true);
      return;
    }
    if (target.is_max_boundary()) {
      this->set_is_vertical(true);
      this->set_is_directed_right(true);
      set_is_x_monotone(true);
      return;
    }
    if (target.is_min_boundary()) {
      this->set_is_vertical(true);
      this->set_is_directed_right(false);
      set_is_x_monotone(true);
      return;
    }

    // None of the enpoints coincide with a pole:
    Direction_3 normal = this->m_normal;
    if (z_sign(normal) == ZERO) {
      // The arc is vertical
      this->set_is_vertical(true);
      bool s_is_positive, t_is_positive, plane_is_positive;
      CGAL::Sign xsign = x_sign(normal);
      if (xsign == ZERO) {
        s_is_positive = x_sign(source) == POSITIVE;
        t_is_positive = x_sign(target) == POSITIVE;
        plane_is_positive = y_sign(normal) == NEGATIVE;
      } else {
        s_is_positive = y_sign(source) == POSITIVE;
        t_is_positive = y_sign(target) == POSITIVE;
        plane_is_positive = xsign == POSITIVE;
      }
      set_is_x_monotone(s_is_positive == t_is_positive);
      bool ccw = ((plane_is_positive && s_is_positive) ||
                  (!plane_is_positive && !s_is_positive));
      this->set_is_directed_right(ccw);
      return;
    }
    
    // The arc is not vertical!
    this->set_is_vertical(false);
    Direction_2 s = Traits::project_xy(source);
    Direction_2 t = Traits::project_xy(target);
    Orientation orient = Traits::orientation(s, t);
    
    const Direction_2 & d = Traits::identification_xy();
    if (orient == LEFT_TURN) {
      this->set_is_directed_right(true);
      set_is_x_monotone(!kernel.counterclockwise_in_between_2_object()(d, s, t));
      return;
    }        

    // (orient == RIGHT_TURN)
    this->set_is_directed_right(false);
    set_is_x_monotone(!kernel.counterclockwise_in_between_2_object()(d, t, s));
    return;
  }

  Arr_geodesic_arc_on_sphere_3(const Arr_extended_direction_3 & source,
                               const Arr_extended_direction_3 & target,
                               const Direction_3 & normal)
  {
    Kernel kernel;

    this->set_source(source);
    this->set_target(target);
    this->normal(normal);
    this->set_is_degenerate(false);
    this->set_is_empty(false);

    typedef Arr_geodesic_arc_on_sphere_traits_2<Kernel> Traits;

    CGAL_precondition(this->has_on(source));
    CGAL_precondition(this->has_on(target));

    if (z_sign(normal) == ZERO) {
      this->set_is_vertical(true);
    
      // Check whether both endpoint coincide with the poles:
      if (source.is_min_boundary() && target.is_max_boundary()) {
        // Both endpoints coincide with the 2 poles respectively.
        this->set_is_directed_right(true);
        this->set_is_full(false);
        set_is_x_monotone(true);
        return;
      }
      
      if (source.is_max_boundary() && target.is_min_boundary()) {
        // Both endpoints coincide with the 2 poles respectively.
        this->set_is_directed_right(false);
        this->set_is_full(false);
        set_is_x_monotone(true);
        return;
      }

      CGAL::Sign xsign = x_sign(normal);
      bool xz_plane = xsign == ZERO;
      bool s_is_positive, t_is_positive, plane_is_positive;
      if (xz_plane) {
        s_is_positive = x_sign(source) == POSITIVE;
        t_is_positive = x_sign(target) == POSITIVE;
        plane_is_positive = y_sign(normal) == NEGATIVE;
      } else {
        s_is_positive = y_sign(source) == POSITIVE;
        t_is_positive = y_sign(target) == POSITIVE;
        plane_is_positive = xsign == POSITIVE;
      }

      // Process degenerate cases:
      if (source.is_min_boundary()) {
        this->set_is_directed_right(true);
        set_is_x_monotone((plane_is_positive && t_is_positive) ||
                          (!plane_is_positive && !t_is_positive));
        return;
      }
      if (source.is_max_boundary()) {
        this->set_is_directed_right(false);
        set_is_x_monotone((plane_is_positive && !t_is_positive) ||
                          (!plane_is_positive && t_is_positive));
        return;
      }
      if (target.is_min_boundary()) {
        this->set_is_directed_right(false);
        set_is_x_monotone((plane_is_positive && !s_is_positive) ||
                          (!plane_is_positive && s_is_positive));
        return;
      }
      if (target.is_max_boundary()) {
        this->set_is_directed_right(true);
        set_is_x_monotone((plane_is_positive && s_is_positive) ||
                          (!plane_is_positive && !s_is_positive));
        return;
      }
      if (s_is_positive != t_is_positive) {
        set_is_x_monotone(false);
        return;
      }

      /* Non of the endpoints coincide with a pole.
       * The projections of both endpoints lie on the same hemi-circle.
       * Thus, either the arc is x-monotone, or it includes both poles.
       * This means that it is sufficient to check whether one pole lies
       * on the arc in order to determine x-monotonicity
       */
      
      typename Traits::Project project =
        (xz_plane) ? Traits::project_xz : Traits::project_yz;
      Direction_2 s = project(source);
      Direction_2 t = project(target);
      const Direction_2 & ny = Traits::neg_y_2();
      typename Kernel::Counterclockwise_in_between_2 ccib =
        kernel.counterclockwise_in_between_2_object();
      set_is_x_monotone((plane_is_positive && !ccib(ny, s, t)) ||
                        (!plane_is_positive && !ccib(ny, t, s)));

      bool ccw = ((plane_is_positive && s_is_positive) ||
                  (!plane_is_positive && !s_is_positive));
      this->set_is_directed_right(ccw);
      return;
    }
      
    // The arc is not vertical!
    this->set_is_vertical(false);
    this->set_is_directed_right(z_sign(normal) == POSITIVE);
    const Direction_2 & d = Traits::identification_xy();
    Direction_2 s = Traits::project_xy(source);
    Direction_2 t = Traits::project_xy(target);    
    typename Kernel::Counterclockwise_in_between_2 ccib =
      kernel.counterclockwise_in_between_2_object();
    bool plane_is_positive = (z_sign(normal) == POSITIVE);
    set_is_x_monotone((plane_is_positive && !ccib(d, s, t)) ||
                      (!plane_is_positive && !ccib(d, t, s)));
  }

  Arr_geodesic_arc_on_sphere_3(const Direction_3 & normal)
  {
    this->normal(normal);
    this->set_is_vertical(CGAL::sign(normal.dz()) == ZERO);
    this->set_is_directed_right(true);
    this->set_is_full(true);
    this->set_is_degenerate(false);
    this->set_is_empty(false);
    set_is_x_monotone(false);
  }

  bool is_x_monotone() const { return m_is_x_monotone; }

  void set_is_x_monotone(bool flag) { m_is_x_monotone = flag; }
};

template <typename Kernel, typename OutputStream>
OutputStream & operator<<(OutputStream & os,
                          const Arr_extended_direction_3<Kernel> & ed)
{
  CGAL::To_double<typename Kernel::FT> todouble;
#if defined(CGAL_ARR_GEODESIC_ARC_ON_SPHERE_DETAILS)
  os << "(";
#endif
#if 0
  os << ed.dx() << ", " << ed.dy() << ", " << ed.dz();
#else
  os << static_cast<float>(todouble(ed.dx())) << ", "
     << static_cast<float>(todouble(ed.dy())) << ", "
     << static_cast<float>(todouble(ed.dz()));
#endif
#if defined(CGAL_ARR_GEODESIC_ARC_ON_SPHERE_DETAILS)
  os << ")"
     << ", "
     <<
    (ed.is_min_boundary() ? "min" :
     ed.is_max_boundary() ? "max" :
     ed.is_mid_boundary() ? "dis" : "reg");
#endif
  return os;
}

template <typename Kernel, typename OutputStream>
OutputStream &
operator<<(OutputStream & os,
           const Arr_x_monotone_geodesic_arc_on_sphere_3<Kernel> & arc)
{
#if defined(CGAL_ARR_GEODESIC_ARC_ON_SPHERE_DETAILS)
  os << "(";
#endif
  os << "(" << arc.source() << "), (" << arc.target() << ")";
#if defined(CGAL_ARR_GEODESIC_ARC_ON_SPHERE_DETAILS)
  os << "("
     << ", " << (arc.is_vertical() ? " |" : "!|")
     << ", " << (arc.is_directed_right() ? "=>" : "<=")
     << ", " << (arc.is_full() ? "o" : "/");
#endif
  return os;
}

template <typename Kernel, typename InputStream>
InputStream &
operator>>(InputStream & is,
           const Arr_x_monotone_geodesic_arc_on_sphere_3<Kernel> & arc)
{
  std::cerr << "Not implemented yet!" << std::endl;
  return is;
}  

CGAL_END_NAMESPACE

#endif