sphere_geometry_ogl.h

CGAL is a collaborative effort of several sites in Europe and Israel. The goal is to make the most i

// Copyright (c) 1997-2002  Max-Planck-Institute Saarbruecken (Germany).
// 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.
//
// $Source: /CVSROOT/CGAL/Packages/Nef_S2/include/CGAL/Nef_S2/Sphere_geometry_OGL.h,v $
// $Revision: 1.20.2.3 $ $Date: 2004/12/17 15:03:05 $
// $Name:  $
//
// Author(s)     : Michael Seel  <[e-mail]>

#ifndef CGAL_SPHERE_GEOMETRY_OGL_H
#define CGAL_SPHERE_GEOMETRY_OGL_H

#include <CGAL/Nef_S2/OGL_base_object.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Nef_S2/Sphere_geometry.h>
#include <CGAL/Nef_S2/SM_triangulator.h>
#include <CGAL/IO/Color.h>
#include <qgl.h>
#include <cmath>
#include <cstdlib>
#include <cstdio>
#include <string>
#include <vector>

#undef CGAL_NEF_DEBUG
#define CGAL_NEF_DEBUG 151
#include <CGAL/Nef_2/debug.h>

#define LGREY CGAL::Color(170,170,200)
#define DGREY CGAL::Color(30,30,50)

CGAL_BEGIN_NAMESPACE
namespace OGL {

struct Gen_object {
  Gen_object() {}
  virtual ~Gen_object() {}
  virtual void draw() const {}
  virtual Gen_object* clone() const { return 0; }
  virtual void print() const {}
};

typedef CGAL::Simple_cartesian<double>      VKernel;
typedef VKernel::Vector_3                   VVector;
typedef VKernel::Point_3                    VPoint;
typedef VKernel::Aff_transformation_3       VTrafo;
typedef VKernel::Aff_transformation_3       Affine_3;
typedef std::vector<VPoint>                 VSegment;
typedef VKernel::Triangle_3                 DTriangle;
typedef std::vector<DTriangle>              VTriangle;

const double refinement_angle = 0.1;
const double shrink_fac = 0.995;

template <typename R>
class Approximator {

 public:
  static VPoint approximate(const CGAL::Sphere_point<R>& p) {
    VVector v = convert(p-CGAL::ORIGIN);
    v = v / CGAL_NTS sqrt(v*v) ; // normalize
    return CGAL::ORIGIN+v;
  }
   
  static VSegment approximate(const CGAL::Sphere_segment<R>& s) {
    
    /* we construct the rotation matrix that transfers the x-axis
       into |ps_|, the z-axis into h_.orthogonal_vector() and the
       y-axix into the corresponding crossproduct of the two.*/
    if ( s.is_degenerate() ) {
      VSegment S(1);
      S[0] = approximate(s.source());
      return S;
    }
    
    VVector v0 = convert(s.source()-CGAL::ORIGIN);
    VVector v2 = convert(s.sphere_circle().orthogonal_vector());
    VVector v1(-cross_product(v0,v2));
    VVector v3 = convert(s.target()-CGAL::ORIGIN);
    double v0l = CGAL_NTS sqrt(v0*v0);
    double v1l = CGAL_NTS sqrt(v1*v1);
    double v2l = CGAL_NTS sqrt(v2*v2);
    double v3l = CGAL_NTS sqrt(v3*v3);
    double cosalpha = v0*v3 / v0l / v3l; 
    double alpha = acos(cosalpha);
    const int units_per_halfcircle = 50;
    int units = int(units_per_halfcircle/CGAL_PI * alpha);
    if (units == 0) ++units;
    bool seg_is_short = s.is_short();
    bool seg_is_halfcircle = s.is_halfcircle();
    if ( seg_is_halfcircle ) units = units_per_halfcircle;
    else if ( !seg_is_short ) {
      units = 2*units_per_halfcircle - (units+1);
    } CGAL_NEF_TRACEV(units); CGAL_NEF_TRACEV(cosalpha); CGAL_NEF_TRACEV(alpha);
    
    v0 = v0 / v0l;
    v1 = v1 / v1l;
    v2 = v2 / v2l;
    v3 = v3 / v3l;
    VTrafo T(v0.x(),v1.x(),v2.x(),
	     v0.y(),v1.y(),v2.y(),
	     v0.z(),v1.z(),v2.z());
    VSegment S(units+1);
    for (int i=0; i<units; ++i) 
      S[i] = VPoint(cos(CGAL_PI*i/double(units_per_halfcircle)),
		    sin(CGAL_PI*i/double(units_per_halfcircle)),
		    0.0);
    double sinalpha = 1 - cosalpha*cosalpha;
    if (sinalpha <0) sinalpha = 0; 
    else             sinalpha = CGAL_NTS sqrt(sinalpha);
    if ( seg_is_short ) 
      S[units] = VPoint(cosalpha, sinalpha, 0);
    else
      S[units] = VPoint(cosalpha, -sinalpha, 0);
    VSegment::iterator it;
    for(it = S.begin(); it != S.end(); ++it) { CGAL_NEF_TRACEN(*it<<" "<<T(*it));
    *it = T(*it); 
    } CGAL_NEF_TRACEN("");
    return S;
  }

  static VSegment approximate(const CGAL::Sphere_circle<R>& s) {
    
    /* we construct the rotation matrix that transfers the x-axis
       into |ps_|, the z-axis into h_.orthogonal_vector() and the
       y-axix into the corresponding crossproduct of the two.*/
    
    VVector v0 = convert(s.base1());
    VVector v1 = convert(s.base2());
    VVector v2 = convert(s.orthogonal_vector());
    double v0l = CGAL_NTS sqrt(v0*v0);
    double v1l = CGAL_NTS sqrt(v1*v1);
    double v2l = CGAL_NTS sqrt(v2*v2);
    const int units = 100;
    v0 = v0 / v0l;
    v1 = v1 / v1l;
    v2 = v2 / v2l;
    VTrafo T(v0.x(),v1.x(),v2.x(),
	     v0.y(),v1.y(),v2.y(),
	     v0.z(),v1.z(),v2.z());
    VSegment S(units);
    for (int i=0; i<units; ++i) {
      S[i] = VPoint(cos(2.0*CGAL_PI*i/double(units)),
		    sin(2.0*CGAL_PI*i/double(units)),
		    0.0);
    }
    VSegment::iterator it;
    for(it = S.begin(); it != S.end(); ++it) *it = T(*it); 
    return S;
  }


/* the following operation refines a sphere triangle as a list of flat
   triangles in 3d. The refinement only works for triangles that are
   contained in a perfect hemisphere (no long sphere segments are
   allowed as triangle segments). We split triangles along at the
   midpoint of their longest side into two. */

  static void refine(const DTriangle& t, VTriangle& T) {
    double angle[3]; int i(0);
    angle[0] = acos((t[0]-CGAL::ORIGIN)*(t[1]-CGAL::ORIGIN));
    angle[1] = acos((t[1]-CGAL::ORIGIN)*(t[2]-CGAL::ORIGIN));
    angle[2] = acos((t[2]-CGAL::ORIGIN)*(t[0]-CGAL::ORIGIN));
    CGAL_NEF_TRACEN("refine "<<angle[0]<<" "<<angle[1]<<" "<<angle[2]);
    if ( angle[1] > angle[0] ) {
      if ( angle[2] > angle[1] ) i=2;
      else                       i=1;
    } else { // angle[0] >= angle[1]
      if ( angle[2] > angle[0] ) i=2;
      else                       i=0;
    }
    // now i references the side of maximal angle
    if ( angle[i] < refinement_angle ) // refinement threshhold
      { T.push_back(t); return; }
    VVector v;
    switch (i) {
    case 0: v = (t[0]-CGAL::ORIGIN)+(t[1]-CGAL::ORIGIN); break;
    case 1: v = (t[1]-CGAL::ORIGIN)+(t[2]-CGAL::ORIGIN); break;
    case 2: v = (t[2]-CGAL::ORIGIN)+(t[0]-CGAL::ORIGIN); break;
    }
    v = v / CGAL_NTS sqrt(v*v) ; // normalize
    VPoint p = CGAL::ORIGIN+v;
    DTriangle t1,t2;
    switch (i) {
    case 0: t1=DTriangle(t[0],p,t[2]); t2=DTriangle(p,t[1],t[2]); break;
    case 1: t1=DTriangle(t[1],p,t[0]); t2=DTriangle(p,t[2],t[0]); break;
    case 2: t1=DTriangle(t[2],p,t[1]); t2=DTriangle(p,t[0],t[1]); break;
    }
    refine(t1,T);
    refine(t2,T);
  }

  static VTriangle approximate(const CGAL::Sphere_triangle<R>& t) {
    // we subdivide the triangle into a list of triangles until
    // we reach a fine resolution on the surface.
    
    VTriangle T;
    DTriangle td(approximate(t.point(0)), 
		 approximate(t.point(1)),
		 approximate(t.point(2)));
    CGAL_NEF_TRACEN("approximate " << td);
    refine(td,T);
    return T;
  }

};

template <typename R>
VVector convert(const CGAL::Vector_3<R>& v)
{ return VVector(CGAL::to_double(v.x()),
                 CGAL::to_double(v.y()),
                 CGAL::to_double(v.z())); }



template <class R_>
class Sphere_point : public VPoint, public Gen_object {
  typedef R_ R;
  CGAL::Sphere_point<R> p_;
  CGAL::Color           c_;
  unsigned              w_;
public:
  Sphere_point() {}
  Sphere_point(const CGAL::Sphere_point<R>& p,
    CGAL::Color c = CGAL::BLACK, unsigned w = 10) : 
    VPoint(Approximator<R>::approximate(p)), p_(p), c_(c), w_(w) {}
  Sphere_point(const Sphere_point<R>& p) : VPoint(p)
  { p_ = p.p_; c_ = p.c_; w_ = p.w_; }
  Sphere_point<R>& operator=(const Sphere_point<R>& p)
  { VPoint::operator=(p); p_ = p.p_;  c_ = p.c_; w_ = p.w_;
    return *this; }

  virtual ~Sphere_point() {}

  const CGAL::Sphere_point<R>& original() const
  { return p_; }

  virtual Gen_object* clone() const 
  { return new Sphere_point<R>(*this); }
 
  virtual void draw() const { 
    glPointSize(w_);
    glColor3ub(c_.red(),c_.green(),c_.blue());
    glBegin(GL_POINTS);
    glNormal3d(x(),y(),z());
    glVertex3d(x(),y(),z());
    glEnd();
  }

  virtual void print() const 
  { std::cerr << "point " << p_; }

};





template <class R_>
class Sphere_segment : public VSegment, public Gen_object { 
  typedef R_ R;
  CGAL::Sphere_segment<R> s_;
  CGAL::Color             c_;
  unsigned                w_;
public:
  Sphere_segment() {}
  Sphere_segment(const CGAL::Sphere_segment<R>& s,
    CGAL::Color c = CGAL::BLACK, unsigned w = 2) 
    : VSegment(Approximator<R>::approximate(s)), s_(s), c_(c), w_(w) {}
  Sphere_segment(const Sphere_segment<R>& s) : VSegment(s)
  { s_ = s.s_; c_ = s.c_; w_ = s.w_; }
  Sphere_segment<R>& operator=(const Sphere_segment<R>& s)
  { VSegment::operator=(s); s_ = s.s_; c_ = s.c_; w_ = s.w_;
    return *this; }
  virtual ~Sphere_segment() {}

  const CGAL::Sphere_segment<R>& original() const
  { return s_; }

  virtual Gen_object* clone() const 
  { return new Sphere_segment<R>(*this); }

  virtual void draw() const
  { CGAL_NEF_TRACEN("draw "<<s_);
    if ( size() == 1 ) {
      glPointSize(5*w_);
      glColor3ub(c_.red(),c_.green(),c_.blue());
      glBegin(GL_POINTS);
      glNormal3d(begin()->x(),begin()->y(),begin()->z());
      glVertex3d(begin()->x(),begin()->y(),begin()->z());
      glEnd();
    } else {
      glLineWidth(w_);
      glColor3ub(c_.red(),c_.green(),c_.blue()); 
      glBegin(GL_LINE_STRIP);
      VSegment::const_iterator it;
      for(it = begin(); it != end(); ++it) {
        glNormal3d(it->x(),it->y(),it->z());
        glVertex3d(it->x(),it->y(),it->z());
      }
      glEnd();
    }
  }

  virtual void print() const 
  { std::cerr << "segment " << s_; }

};



template <class R_>
class Sphere_circle : public VSegment, public Gen_object { 
  typedef R_ R;
  CGAL::Sphere_circle<R> s_;
  CGAL::Color            c_;
  unsigned               w_;
public:
  Sphere_circle() {}
  Sphere_circle(const CGAL::Sphere_circle<R>& s,
    CGAL::Color c = CGAL::BLACK, unsigned w = 2) 
    : VSegment(Approximator<R>::approximate(s)), s_(s), c_(c), w_(w) {}
  Sphere_circle(const Sphere_circle<R>& s) : VSegment(s)
  { s_ = s.s_; c_ = s.c_; w_ = s.w_; }
  Sphere_circle<R>& operator=(const Sphere_circle<R>& s)
  { VSegment::operator=(s); s_ = s.s_; c_ = s.c_; w_ = s.w_;
    return *this; }
  virtual ~Sphere_circle() {}

  const CGAL::Sphere_circle<R>& original() const
  { return s_; }

  virtual Gen_object* clone() const 
  { return new Sphere_circle<R>(*this); }

  virtual void draw() const
  { CGAL_NEF_TRACEN("draw "<<s_);
    glLineWidth(w_);
    glColor3ub(c_.red(),c_.green(),c_.blue()); 
    glBegin(GL_LINE_LOOP);
    VSegment::const_iterator it;
    for(it = begin(); it != end(); ++it) {
      glNormal3d(it->x(),it->y(),it->z());
      glVertex3d(it->x(),it->y(),it->z());
    }
    glEnd();
  }

  virtual void print() const 
  { std::cerr << "circle " << s_; }

};


/* The following class approximates a spherical triangle by a list
   of flat triangles */

template <class R_>
class Sphere_triangle : public VTriangle, public Gen_object { 
  typedef R_ R;
  CGAL::Sphere_triangle<R> t_;
  CGAL::Color              c_;
public:
  Sphere_triangle() {}

  Sphere_triangle(const CGAL::Sphere_triangle<R>& t,
    CGAL::Color c = CGAL::Color(100,100,120)) 
    : VTriangle(Approximator<R>::approximate(t)), t_(t), c_(c) {}

  Sphere_triangle(const Sphere_triangle<R>& t) : VTriangle(t)
  { t_ = t.t_; c_ = t.c_; }

  Sphere_triangle<R>& operator=(const Sphere_triangle<R>& t)
  { VTriangle::operator=(t); t_ = t.t_; c_ = t.c_; return *this; }

  virtual ~Sphere_triangle() {}

  const CGAL::Sphere_triangle<R>& original() const
  { return t_; }

  virtual Gen_object* clone() const 
  { return new Sphere_triangle<R>(*this); }

  virtual void draw() const { 
    CGAL_NEF_TRACEN("draw "<<t_ << " in " << c_);
    VTriangle::const_iterator it;
    VPoint p;
    glColorMaterial(GL_FRONT, GL_DIFFUSE);
    glEnable(GL_COLOR_MATERIAL);
    glColor3ub(c_.red(),c_.green(),c_.blue());
    glBegin(GL_TRIANGLES);
    for(it = begin(); it != end(); ++it) {
      p = it->vertex(0); 
      glNormal3d(p.x(),p.y(),p.z()); 	//glVertex3d(p.x(),p.y(),p.z());
      glVertex3d(shrink_fac*p.x(),shrink_fac*p.y(),shrink_fac*p.z());
      p = it->vertex(1); 
      glNormal3d(p.x(),p.y(),p.z()); 
      glVertex3d(shrink_fac*p.x(),shrink_fac*p.y(),shrink_fac*p.z());
      p = it->vertex(2); 
      glNormal3d(p.x(),p.y(),p.z()); 
      glVertex3d(shrink_fac*p.x(),shrink_fac*p.y(),shrink_fac*p.z());
    }