delaunay_sweep.c

A Library of Efficient Data Types and Algorithms,封装了常用的ADT及其相关算法的软件包

    if ( DS.sweep() ) DS.finish(S);
    else  cout << "Fehler beim Sweepen !" << endl;
  else cout << "Initialisierung fehlgeschlagen !" << endl;
} 



GeoBaseScene<list<rat_point> >* shore_sc = 0;
GeoBaseScene<list<rat_point> >* delau_sc = 0;


static GeoWin* Gwp = 0;
static window* Wp = 0;
DelaunaySweep* Dsweep =0;



class WaveFrontRep
{
  public :
  
  point                        p_sweep;
  DelaunaySweep::ystructure*   yst;
  
  WaveFrontRep() : yst(0) {}
  WaveFrontRep(const WaveFrontRep& w) :  p_sweep(w.p_sweep), yst(w.yst) {}
  
  double yparab( double y, point p ) const;
  void xparab( double x, point p, double& y1, double& y2 ) const; 
  void computeX(WAVE_PART* w, double& x0) const;
  void computeY(WAVE_PART* w, double& y0, double& y1) const;
  void plotter(window& w, point p, double y0, double y1, color col) const;
  void show_ystructure() const;
  
  friend unsigned long ID_Number(const WaveFrontRep& x)
  { 
    return ID_Number(x.p_sweep); 
  }
  friend const char* leda_tname(WaveFrontRep* w) {  return "WaveFrontRep"; }
  friend bool identical(const WaveFrontRep& w1, const WaveFrontRep& w2)
  {
    return &w1 == &w2;
  }
};




double WaveFrontRep::yparab( double y, point p ) const
{
  if ( p_sweep.xcoord() - p.xcoord()  == 0 ) return Wp->xmin();
  
  double A = p_sweep.xcoord()*p_sweep.xcoord();
  double B = (y - p.ycoord())*(y - p.ycoord());
  double C = p.xcoord()*p.xcoord();
  double res = (A - B - C)/(2*(p_sweep.xcoord() - p.xcoord()));
  return res;
}

void WaveFrontRep::xparab( double x, point p, double& y1, double& y2 ) const
{
  double A = p_sweep.xcoord()*p_sweep.xcoord();
  double B = 2*x*(p.xcoord() - p_sweep.xcoord());
  double C = p.xcoord()*p.xcoord();
  double D = sqrt(B + A - C);
  
  y1 = p.ycoord() - D;
  y2 = p.ycoord() + D;
}

void  WaveFrontRep::computeX(WAVE_PART* w, double& x0) const
{
  double x;
  point  p;
  
  x0 = Wp->xmin();
  if ( !identical(w->pred, WAVE_PART::INFINITY1) ) 
    {
      p = w->pred.to_point();
      x = yparab(w->p.to_point().ycoord(), p);
      if ( x > x0 ) x0 = x;
    }
  if ( !identical(w->succ, WAVE_PART::INFINITY2 ) )
    {
      p = w->succ.to_point();
      x = yparab(w->p.to_point().ycoord(), p);
      if ( x > x0 ) x0 = x;
    }
}

void WaveFrontRep::computeY(WAVE_PART* w, double& y0, double& y1) const
{
  y0 = Wp->ymin();
  y1 = Wp->ymax();
  
  double A, B, C, D, E, F;
  point p, q;
  int   tol;
  
  xparab( Wp->xmin(), w->p.to_point(), E, F );
  
  if( E > y0 ) y0 = E;
  if( F < y1 ) y1 = F;
  
  if( !identical(w->succ, WAVE_PART::INFINITY2 ) )
    {
      p = w->p.to_point();
      q = w->succ.to_point();
      
      A = p.xcoord() - p_sweep.xcoord();
      B = q.xcoord() - p_sweep.xcoord();
      C = q.ycoord() - p.ycoord();
      D = p.xcoord() - q.xcoord();
      
      if ( C == 0 ) 
	{
	  A = p_sweep.xcoord()*p_sweep.xcoord();
	  B = p_sweep.xcoord()*(p.xcoord() + q.xcoord());
	  C = p.xcoord()*q.xcoord();
	  tol = compare(p, q) < 0 ? -1 : 1;
	  F =  p.ycoord() + tol*sqrt(A-B+C); 
	}
      else
	if ( D == 0 )
	  if ( A != 0 && C > 0) F = (q.ycoord() + p.ycoord())/2;
          else; 
	else
	  {
	    if ( C > 0 )
	      F = (A*C)/D + C/D*sqrt(A*(B+D*D/(C*C)*B));
	    else
	      F = (A*C)/D - C/D*sqrt(A*(B+D*D/(C*C)*B));
	    
	    F = F + p.ycoord();
	  }
    }
  
  
  if( !identical(w->pred, WAVE_PART::INFINITY1) )
    {
      p = w->p.to_point();
      q = w->pred.to_point();
      
      A = p.xcoord() - p_sweep.xcoord();
      B = q.xcoord() - p_sweep.xcoord();
      C = q.ycoord() - p.ycoord();
      D = p.xcoord() - q.xcoord();
      
      if ( C == 0 ) 
	{
	  A = p_sweep.xcoord()*p_sweep.xcoord();
	  B = p_sweep.xcoord()*(p.xcoord() + q.xcoord());
	  C = p.xcoord()*q.xcoord();
	  tol = compare(p, q) < 0 ? -1 : 1;
	  E =  p.ycoord() - tol*sqrt(A-B+C); 
	}
      else
	if ( D == 0 ) 
	  if ( A != 0 && C < 0) E = (q.ycoord() + p.ycoord())/2;
	  else;
	else
	  {
	    if ( C > 0 )
	      E = (A*C)/D - C/D*sqrt(A*(B+D*D/(C*C)*B));
	    else
	      E = (A*C)/D + C/D*sqrt(A*(B+D*D/(C*C)*B));
	    
	    E = E + p.ycoord();
	  }
    }
  
  if( E > y0 ) y0 = E;
  if( F < y1 ) y1 = F;
}

void WaveFrontRep::plotter(window& w, point p, double y0, double y1, 
			   color col) const
{
  w.set_color(col);
  
  double y;
  double dy = 1/w.scale();
  double x_new;
  
  if (y1 <= y0) return;
  
  int lw = w.get_line_width();
  w.set_line_width(1);
  
  list<point> L;
  
  x_new = yparab(y0, p);
  L.append(point(x_new, y0));
  
  for(y = y0+dy; y < y1; y+=4*dy)
    {
      x_new = yparab(y, p);
      L.append(point(x_new, y));
    }
  
  x_new = yparab(y1, p);
  L.append(point(x_new, y1));
  
  w.draw_spline( L, 4, col);
  w.set_line_width(lw);  
}


void WaveFrontRep::show_ystructure() const
{
  if (!yst) return;

  if (!shore_sc->is_visible()) return;

  circles.clear();

  DelaunaySweep::yitem   dit;
  forall_items(dit, *yst)
  { WAVE_PART* w = yst->key(dit);
     double y0,y1;
     computeY(w,y0,y1);
     plotter(*Wp,w->p.to_point(),y0,y1,black);
/*
     if( w->describes_circle_event() )
     { CIRCLE ci(w->p, w->pred, w->succ);
       circles.append(ci);
      }
*/
  }
}



// //////////////////////////////////////////////////////////////////////
// routines for using the class in GeoWin

// ps_file output

ps_file& operator<<(ps_file& ps,  const WaveFrontRep& c1)
{
 return ps;
}

// Presentation of the Sweepline :
window& operator<<(window& w, const WaveFrontRep& c1)
{ 
  point r1 = c1.p_sweep;
  point r2 = r1.translate(0, 1);
  w.draw_line(r1, r2, green);
  c1.show_ystructure();
  return w;
}

window& operator >> (window& w, WaveFrontRep& p)
{
  return w;
}

ostream& operator << (ostream& os, const WaveFrontRep& w) { return os; }
istream& operator >> (istream& is, WaveFrontRep& w)      {  return is; }


bool geowin_intersects(const WaveFrontRep& c1, double x1, double y1, 
		       double x2, double y2, bool)
{
  double xc = c1.p_sweep.xcoord();
  return (xc >= x1 && xc <= x2);
}

void geowin_fit(const WaveFrontRep& obj, double& x1, double& x2,
		double& y1, double& y2)
{
}

void geowin_move_view(WaveFrontRep& obj, double dx, double dy)
{
  obj.p_sweep = obj.p_sweep.translate(dx,dy);
/*
  point p = obj.p_sweep.translate(dx, dy);
  WaveFrontRep w;
  w.p_sweep = p; 
  w.yst = obj.yst;
  return w;
*/
}


bool sweepline_pre_move_handler(GeoWin&, const WaveFrontRep&, 
                                  double dx, double dy)
{ 
  return dx >0; 
}



//s.n.
//void update_delaunay(const list<WaveFrontRep>& L, GRAPH<rat_point, int>& G)
static void update_delaunay(const list<WaveFrontRep>& L, list<rat_point>&)
{
  circles.clear();

  if( !L.empty() )
  { const WaveFrontRep& wfr = L.head();  

    DelaunaySweep::yitem   dit;
    forall_items(dit, *(wfr.yst))
    { WAVE_PART* w = wfr.yst->key(dit);
      if( w->describes_circle_event() )
        circles.append(CIRCLE(w->p, w->pred, w->succ));
    }

    point pp = wfr.p_sweep;
    rat_point ppp = rat_point(pp);
    rat_circle stop = rat_circle(ppp, ppp, ppp);
    Dsweep->animate(stop, 0);
  }
}


extern bool initialize(DelaunaySweep& DS, const list<rat_point>& L);
extern bool transit(DelaunaySweep& GS );
extern bool finish(DelaunaySweep& GS, const list<rat_point>& L);





void finish_sweep(GeoWin& gw)
{
  if( Dsweep ) 
    { circles.clear();
      Dsweep->sweep();
      Dsweep->finish(Dsweep->get_input());
      geo_scene out = gw.get_scene_with_name("Sweep Output");
      gw.set_visible(gw.get_scene_with_name("Sweepline"),false);
      gw.edit(out);
      gw.activate(gw.get_scene_with_name("Sweepline"));
      delete Dsweep;
      gw.destroy(gw.get_scene_with_name("Sweepline"));
      Dsweep = 0;
      gw.reset_actions();
    }
}

void break_sweep(GeoWin& gw)
{
  if( Dsweep ) 
    {
      Dsweep->finish(Dsweep->get_input());
      geo_scene out = gw.get_scene_with_name("Sweep Output");
      gw.edit(out);
      gw.activate(gw.get_scene_with_name("Sweepline"));
      delete Dsweep;
      gw.destroy(gw.get_scene_with_name("Sweepline"));
      Dsweep = 0;
      gw.reset_actions();
    }
}


// s.n.
static GRAPH<rat_point,int> Gtmp;


class redraw_graph : public geowin_redraw
{
public:
 virtual bool operator()(window& W,color col1, color col2, double,double,double,double)
 { 
  if (!delau_sc->is_visible()) return false; 

  edge_array<bool> drawn(Gtmp,false);
  edge e;
  forall_edges(e,Gtmp) 
  { edge r = Gtmp.reversal(e);
    if (r)
    { if (drawn[e]) continue;
       drawn[r] = drawn[e] = true;
     }
    node v = source(e);
    node w = target(e);
    //W.draw_segment(Gtmp[v].to_point(),Gtmp[w].to_point(),col1);
    W.set_node_width(2);
    W.draw_edge(Gtmp[v].to_point(),Gtmp[w].to_point(),col1);
  }  
  return false;
 }
};

redraw_graph RDG;

geowin_update<list<WaveFrontRep>, list<rat_point> > sweep_update(update_delaunay);


static list<WaveFrontRep> wave_front_list; 
 
void start_sweep( GeoWin& gw)
{
  // is a sweep running ?
  if( Dsweep ) return;
  
  // mouse actions:
  gw.clear_actions();
  gw.set_action(A_LEFT| A_DRAG | A_OBJECT, geo_object_dragging);
  gw.del_pin_point();
  
  // get input
  list<rat_point> L;
  //gw.get_objects(L);
  geowin_get_objects(gw, gw.get_active_scene(), L);
   

  // make a scene for the sweepline
  wave_front_list.clear();

  GeoEditScene< list<WaveFrontRep> >*
  SweeplineScene = (GeoEditScene< list<WaveFrontRep> >*)
      //geowin_new_scene(gw, (list<WaveFrontRep>*)0); 
      geowin_new_scene(gw, wave_front_list); 

  gw.set_name(SweeplineScene, "Sweepline");
  
  geowin_set_pre_move_handler(SweeplineScene, sweepline_pre_move_handler);
  SweeplineScene->set_box_intersection_fcn(geowin_intersects);
  SweeplineScene->set_get_bounding_box_fcn(geowin_fit);
  SweeplineScene->set_move_fcn(geowin_move_view);
  
  // scene for graph :
  // link it to the  SweeplineScene 
  // s.n.
  GeoBaseScene<list<rat_point> >* tmp;
  
  tmp = geowin_new_scene(gw, sweep_update, SweeplineScene, "Sweep Output");  

  gw.set_visible(tmp, true);
  gw.set_color(tmp, blue);
  gw.set_color2(tmp, blue);

  gw.set_redraw_pt((geo_scene)tmp, &RDG);
  
  // create a DelaunaySweepObjekt
  // s.n.
  //Dsweep = new DelaunaySweep(&L, &(tmp->get_objref()));   
  Dsweep = new DelaunaySweep(&L, &Gtmp);

  Dsweep->set_init_handler(initialize);
  Dsweep->set_event_handler(transit);
  Dsweep->set_finish_handler(finish);
  
  if( !Dsweep->init(L) )
    {
      gw.message("Initialisierung fehlgeschlagen !");
      finish_sweep(gw);
      return;
    }
  
  DelaunaySweep::xkey swp = Dsweep->sweep_position();
  list<WaveFrontRep> LL;
  WaveFrontRep startfront;
  circle ccc = swp.to_circle();
  point ppp  = ccc.center();
  startfront.p_sweep = ppp.translate(ccc.radius(), 0);
  startfront.yst =&(Dsweep->get_ystructure());
  LL.append(startfront);
  
  SweeplineScene->set_objects(LL);
  gw.set_visible(SweeplineScene, true); 
  gw.activate(SweeplineScene);
  gw.redraw();
}


#include <LEDA/rat_geo_alg.h>

geowin_graph_update<list<rat_point>,GRAPH<rat_circle,rat_point> > voro(VORONOI);

int main()
{
  rand_int.set_seed(1234567);

  GeoWin GW("DELAUNAY SWEEP DEMO");
  Gwp = &GW;

  list<rat_point> L;
  geo_scene sc = geowin_new_scene(GW,L);
  GW.set_point_style(sc,circle_point);
  GW.set_visible(sc,true);
  GW.set_name(sc,"sweep_points");

  // make a (dummy) scene for the shore line 
  // s.n.
  list<rat_point> dummy_L;
  shore_sc = geowin_new_scene(GW,dummy_L);  
  GW.set_name(shore_sc,"shoreline");  
  //GW.set_visible(shore_sc, true);

  // and another one for the delaunay triangulation
  // s.n.
  list<rat_point> dummy_L1;
  delau_sc = geowin_new_scene(GW,dummy_L1);  
  GW.set_name(delau_sc,"delaunay");  
  GW.set_visible(delau_sc, true);



  geo_scene circle_sc = geowin_new_scene(GW,circles);
  GW.set_color(circle_sc,red);
  GW.set_fill_color(circle_sc,invisible);
  GW.set_name(circle_sc,"Circle Events");

  geo_scene sc_voro = geowin_new_scene(GW,voro,sc, "Voronoi Diagram");
  GW.set_color(sc_voro,brown);

  Wp = &GW.get_window();

  GW.display(window::center,window::center);

  for(;;)
  { GW.set_fill_color(sc,yellow);
    GW.edit(sc);
    GW.set_fill_color(sc,ivory);
    start_sweep(GW);
    GW.edit(GW.get_scene_with_name("Sweepline"));
    finish_sweep(GW);
    GW.edit(GW.get_scene_with_name("Sweep Output"));
   }

  break_sweep(GW);
 
  return 0;
}