开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 查看源码
seg_int.c
资源名称：leda.tar.gz　[点击查看]
上传用户：gzelex
上传日期：2007-01-07
资源大小：707k
文件大小：14k
源码类别：

数值算法/人工智能

开发平台：

MultiPlatform

seg_int.c：源码内容
  2. #include <LEDA/segment.h>
  3. #include <LEDA/graph.h>
  4. #include <LEDA/prio.h>
  6. #include <LEDA/_sortseq.h>
  7. #include <LEDA/impl/skiplist.h>
  8. #include <LEDA/impl/ab_tree.h>
  10. #include <LEDA/plane_alg.h>
  12. #include <math.h>
  15. #define EPS  0.00001
  16. #define EPS2 0.0000000001
  18. class sw_POINT;
  19. class sw_SEGMENT;
  20. typedef sw_POINT* sw_point;
  21. typedef sw_SEGMENT* sw_segment;
  23. enum sw_point_type {Cross=0,Rightend=1,Leftend=2}; 
  25. class sw_POINT
  26. {
  27.   friend class sw_SEGMENT;
  29.   sw_segment seg;
  30.   int     kind;
  31.   double  x;
  32.   double  y;
  34.   public:
  36.   sw_POINT(double a,double b)  
  37.   { 
  38.     x=a; y=b; seg=0; kind=Cross;
  39.    }
  43.   sw_POINT(point p)         
  44.   { 
  45.     x=p.xcoord();y=p.ycoord();seg=0;kind=Cross;
  46.    }
  49.   LEDA_MEMORY(sw_POINT);
  51.   friend double     get_x(sw_point p);
  52.   friend double     get_y(sw_point p);
  53.   friend int        get_kind(sw_point p);
  54.   friend sw_segment get_seg(sw_point p);
  56.   friend bool intersection(sw_segment, sw_segment, sw_point&);
  57. };
  59. inline double     get_x(sw_point p)    { return p->x; }
  60. inline double     get_y(sw_point p)    { return p->y; }
  61. inline int        get_kind(sw_point p) { return p->kind; }
  62. inline sw_segment get_seg(sw_point p)  { return p->seg; }   
  66. static int compare(const sw_point& p1, const sw_point& p2)
  67. { if (p1==p2) return 0;
  69.   double diffx = get_x(p1) - get_x(p2);
  70.   if (fabs(diffx) > EPS2 ) return (diffx > 0.0) ? 1 : -1;
  72.   int    diffk = get_kind(p1)-get_kind(p2);
  73.   if (diffk != 0) return diffk;
  75.   double diffy = get_y(p1) - get_y(p2);
  76.   if (fabs(diffy) > EPS2 ) return (diffy > 0.0) ? 1 : -1;
  78.   return 0;
  80. }
  83. class sw_SEGMENT
  84. {
  85.   sw_point startpoint;
  86.   sw_point endpoint;
  87.   double  slope;
  88.   double  yshift;
  89.   node  left_node;
  90.   int   orient;
  91.   int   color;
  92.   int   name;
  95.   public:
  97.   sw_SEGMENT(sw_point, sw_point,int,int);     
  99.  ~sw_SEGMENT() { delete startpoint; delete endpoint; }     
  101.   LEDA_MEMORY(sw_SEGMENT);
  103.   friend sw_point get_startpoint(sw_segment seg);
  104.   friend sw_point get_endpoint(sw_segment seg);
  105.   friend double get_slope(sw_segment seg);
  106.   friend double get_yshift(sw_segment seg);
  107.   friend int  get_orient(sw_segment seg);
  108.   friend int  get_color(sw_segment seg);
  109.   friend int  get_name(sw_segment seg);
  110.   friend node get_left_node(sw_segment seg);
  111.   friend void set_left_node(sw_segment seg, node v);
  113.   friend bool intersection(sw_segment, sw_segment, sw_point&);
  114. };
  116. inline sw_point get_startpoint(sw_segment seg)    { return seg->startpoint; }
  117. inline sw_point get_endpoint(sw_segment seg)      { return seg->endpoint; }
  118. inline double get_slope(sw_segment seg)           { return seg->slope; }
  119. inline double get_yshift(sw_segment seg)          { return seg->yshift; }
  120. inline int  get_orient(sw_segment seg)            { return seg->orient; }
  121. inline int  get_color(sw_segment seg)             { return seg->color; }
  122. inline int  get_name(sw_segment seg)              { return seg->name; }
  123. inline node get_left_node(sw_segment seg)         { return seg->left_node; }
  124. inline void set_left_node(sw_segment seg, node v) { seg->left_node = v; }
  127. sw_SEGMENT::sw_SEGMENT(sw_point p1,sw_point p2,int c, int n)    
  128.   {
  129.     left_node  = nil;
  130.     color      = c;
  131.     name       = n;
  133.     if (compare(p1,p2) < 0)
  134.      { startpoint = p1; 
  135.        endpoint = p2; 
  136.        orient = 0;
  137.       }
  138.     else
  139.      { startpoint = p2; 
  140.        endpoint = p1; 
  141.        orient = 1;
  142.       }
  144.     startpoint->kind = Leftend; 
  145.     endpoint->kind = Rightend; 
  146.     startpoint->seg = this; 
  147.     endpoint->seg = this;
  149.     if (endpoint->x != startpoint->x)
  150.     {
  151.       slope = (endpoint->y - startpoint->y)/(endpoint->x - startpoint->x);
  152.       yshift = startpoint->y - slope * startpoint->x;
  154.       startpoint->x -= EPS;
  155.       startpoint->y -= EPS * slope;
  156.       endpoint->x += EPS;
  157.       endpoint->y += EPS * slope;
  158.     }
  159.     else //vertical segment
  160.     { startpoint->y -= EPS;
  161.       endpoint->y   += EPS;
  162.      }
  163.   }
  166. static double x_sweep;
  167. static double y_sweep;
  170. static int compare(const sw_segment& s1, const sw_segment& s2)
  171. {
  172.   double y1 = get_slope(s1)*x_sweep+get_yshift(s1);
  173.   double y2 = get_slope(s2)*x_sweep+get_yshift(s2);
  175.   double diff = y1-y2;
  177.   if (fabs(diff) > EPS2) return (diff > 0.0) ? 1 : -1;
  179.   if (get_slope(s1) == get_slope(s2)) 
  180.         return compare(get_x(get_startpoint(s1)), get_x(get_startpoint(s2)));
  182.   if (y1 <= y_sweep+EPS2)
  183.         return compare(get_slope(s1),get_slope(s2));
  184.   else
  185.         return compare(get_slope(s2),get_slope(s1));
  187. }
  190. void Print(sw_segment& x) 
  191. { sw_point s = get_startpoint(x);
  192.   sw_point e = get_endpoint(x);
  193.   cout << 
  194.     string("[(%f,%f)----(%f,%f)]",get_x(s),get_y(s), get_x(e),get_y(e));
  195. }
  198. static priority_queue<seq_item,sw_point>* X_structure;
  200. static sortseq<sw_segment,pq_item>* Y_structure;
  204. bool intersection(sw_segment seg1,sw_segment seg2, sw_point& inter)
  205. {
  206.   if (seg1->slope == seg2->slope)
  207.     return false;
  208.   else
  209.   {
  210.     //x-coordinate  of the intersection
  211.     double Cross_x = (seg2->yshift - seg1->yshift) / (seg1->slope - seg2->slope);
  212.  
  213.     if (Cross_x <= x_sweep) return false;
  215.     double s1 = seg1->startpoint->x;
  216.     double s2 = seg2->startpoint->x;
  217.     double e1 = seg1->endpoint->x;
  218.     double e2 = seg2->endpoint->x;
  220.     if (s2>Cross_x || s1>Cross_x || Cross_x>e1 || Cross_x>e2) return false;
  222.     //y-coordinate of the intersection
  223.     double Cross_y = (seg1->slope * Cross_x + seg1->yshift);
  224.     inter =  new sw_POINT(Cross_x,Cross_y);
  226.     return true;
  227.   }
  228. }
  231. static pq_item Xinsert(seq_item i, sw_point p) 
  232. { 
  233.   return X_structure->insert(i,p);
  234. }
  236. static sw_point Xdelete(pq_item i) 
  237. {
  238.   sw_point p = X_structure->inf(i);
  239.   X_structure->del_item(i);
  240.   return p;
  241. }
  244. static node New_Node(GRAPH<point,int>& G,double x, double y )
  245. { return G.new_node(point(x,y)); }
  247. static void New_Edge(GRAPH<point,int>& G,node v, node w, sw_segment l )
  248. { if (get_orient(l)==0)
  249.        G.new_edge(v,w,get_color(l));
  250.   else G.new_edge(w,v,get_color(l));
  251. }
  254. void vertical_segment(GRAPH<point,int>& SUB, sw_segment l)
  255. { 
  256.   sw_point p = new sw_POINT(get_x(get_startpoint(l)),get_y(get_startpoint(l)));
  257.   sw_point q = new sw_POINT(get_x(get_endpoint(l)),get_y(get_endpoint(l)));
  259.   sw_point r = new sw_POINT(get_x(p)+1,get_y(p));
  260.   sw_point s = new sw_POINT(get_x(q)+1,get_y(q));
  262.   sw_segment bot = new sw_SEGMENT(p,r,0,0);
  263.   sw_segment top = new sw_SEGMENT(q,s,0,0);
  265.   seq_item bot_it = Y_structure->insert(bot,nil);
  266.   seq_item top_it = Y_structure->insert(top,nil);
  267.   seq_item sit;
  269.   node u,v,w;
  270.   sw_segment seg;
  271.   
  273.   for(sit=Y_structure->succ(bot_it); sit != top_it; sit=Y_structure->succ(sit))
  274.   { seg = Y_structure->key(sit);
  276.     double cross_y = (get_slope(seg) * get_x(p) + get_yshift(seg));
  278.     v = get_left_node(seg);
  279.     if (v==nil)
  280.     { w = New_Node(SUB,get_x(p),cross_y);
  281.       set_left_node(seg,w);
  282.      }
  283.     else
  284.     { double vx = SUB[v].xcoord();
  285.       if ( vx < get_x(p)-EPS2) 
  286.       { w = New_Node(SUB,get_x(p),cross_y);
  287.         New_Edge(SUB,v,w,seg);
  288.         set_left_node(seg,w);
  289.        }
  290.       else w = v;
  291.      }
  293.     u = get_left_node(l);
  294.     if (u!=nil && u!=w) New_Edge(SUB,u,w,l);
  295.     set_left_node(l,w);
  297.    }
  298.   
  299.   delete l;
  300.   delete top;
  301.   delete bot;
  302.     
  303.   Y_structure->del_item(bot_it);
  304.   Y_structure->del_item(top_it);
  306.  }
  309. void Sweep_Segments(const list<segment>& L1, const list<segment>& L2, 
  310.                     GRAPH<point,int>& SUB,
  311.                     priority_queue<seq_item,sw_point>& Q,
  312.                     sortseq<sw_segment,pq_item>& S)
  313. {
  314.   sw_point    p,inter;
  315.   sw_segment  seg, l,lsit,lpred,lsucc,lpredpred;
  316.   pq_item  pqit,pxmin;
  317.   seq_item sitmin,sit,sitpred,sitsucc,sitpredpred;
  319.   X_structure = &Q;
  320.   Y_structure = &S;
  322.   int count=1;
  323.  
  324.   //initialization of the X-structure
  326.   segment s;
  328.   forall(s,L1) 
  329.    { sw_point p = new sw_POINT(s.start());
  330.      sw_point q = new sw_POINT(s.end());
  331.      seg = new sw_SEGMENT(p,q,0,count++);
  332.      Xinsert(nil,get_startpoint(seg));
  333.    }
  335.   count = -1;
  337.   forall(s,L2) 
  338.    { sw_point p = new sw_POINT(s.start());
  339.      sw_point q = new sw_POINT(s.end());
  340.      seg = new sw_SEGMENT(p,q,1,count--);
  341.      Xinsert(nil,get_startpoint(seg));
  342.    }
  345.   count = 0;
  347.   x_sweep = -MAXINT;
  348.   y_sweep = -MAXINT;
  351.   while(!X_structure->empty())
  352.   {
  353.     pxmin = X_structure->find_min();
  354.     p = X_structure->inf(pxmin);
  356.     sitmin = X_structure->key(pxmin);
  358.     Xdelete(pxmin);
  362.     if (get_kind(p) == Leftend)
  364.     //left endpoint
  365.     { 
  367.       l = get_seg(p); 
  369.       x_sweep = get_x(p);
  370.       y_sweep = get_y(p);
  373.       if (get_x(p) == get_x(get_endpoint(l)))
  374.         vertical_segment(SUB,l);
  375.       else
  376.       {
  378. /*
  379.       sit = Y_structure->lookup(l);
  380.       if (sit!=nil)  
  381.            error_handler(1,"plane sweep: sorry, overlapping segments");
  382. */
  385.       sit = Y_structure->insert(l,nil);
  387.       Xinsert(sit,get_endpoint(l));
  389.       sitpred = Y_structure->pred(sit);
  390.       sitsucc = Y_structure->succ(sit);
  392.       if (sitpred != nil) 
  393.       { if ((pqit = Y_structure->inf(sitpred)) != nil)
  394.           delete Xdelete(pqit);
  396.         lpred = Y_structure->key(sitpred);
  398.         Y_structure->change_inf(sitpred,nil);
  400.         if (intersection(lpred,l,inter))
  401.             Y_structure->change_inf(sitpred,Xinsert(sitpred,inter));
  402.       }
  405.       if (sitsucc != nil)
  406.       { lsucc = Y_structure->key(sitsucc);
  407.         if (intersection(lsucc,l,inter))
  408.            Y_structure->change_inf(sit,Xinsert(sit,inter));
  409.       }
  410.      } /* else if vertical */
  412.     }
  413.     else if (get_kind(p) == Rightend)
  414.          //right endpoint
  415.          { 
  416.            x_sweep = get_x(p);
  417.            y_sweep = get_y(p);
  419.            sit = sitmin;
  421.            sitpred = Y_structure->pred(sit);
  422.            sitsucc = Y_structure->succ(sit);
  424.            sw_segment seg = Y_structure->key(sit);
  426.            Y_structure->del_item(sit);
  428.            delete seg;
  430.            if((sitpred != nil)&&(sitsucc != nil))
  431.            {
  432.              lpred = Y_structure->key(sitpred);
  433.              lsucc = Y_structure->key(sitsucc);
  434.              if (intersection(lsucc,lpred,inter))
  435.                 Y_structure->change_inf(sitpred,Xinsert(sitpred,inter));
  436.            }
  437.          }
  438.          else
  439.          /*point of intersection*/
  440.          { 
  441.            node w = New_Node(SUB,get_x(p),get_y(p));
  443.            count++;
  445.            /* Let L = list of all lines intersecting in p 
  446.  
  447.               we compute sit     = L.head();
  448.               and        sitpred = L.tail();
  450.               by scanning the Y_structure in both directions 
  451.               starting at sitmin;
  453.            */
  455.            /* search for sitpred upwards from sitmin: */
  457.            Y_structure->change_inf(sitmin,nil);
  459.            sitpred = Y_structure->succ(sitmin);
  462.            while ((pqit=Y_structure->inf(sitpred)) != nil)
  463.            { sw_point q = X_structure->inf(pqit);
  464.              if (compare(p,q) != 0) break; 
  465.              X_structure->del_item(pqit);
  466.              Y_structure->change_inf(sitpred,nil);
  467.              sitpred = Y_structure->succ(sitpred);
  468.             }
  471.            /* search for sit downwards from sitmin: */
  473.            sit = sitmin;
  475.            seq_item sit1;
  476.            
  477.            while ((sit1=Y_structure->pred(sit)) != nil)
  478.            { pqit = Y_structure->inf(sit1);
  479.              if (pqit == nil) break;
  480.              sw_point q = X_structure->inf(pqit);
  481.              if (compare(p,q) != 0) break; 
  482.              X_structure->del_item(pqit);
  483.              Y_structure->change_inf(sit1,nil);
  484.              sit = sit1;
  485.             }
  489.            // insert edges to p for all sw_segments in sit, ..., sitpred into SUB
  490.            // and set left node to w 
  492.            lsit = Y_structure->key(sitpred);
  494.            node v = get_left_node(lsit);
  495.            if (v!=nil) New_Edge(SUB,v,w,lsit);
  496.            set_left_node(lsit,w);
  498.            for(sit1=sit; sit1!=sitpred; sit1 = Y_structure->succ(sit1))
  499.            { lsit = Y_structure->key(sit1);
  501.              v = get_left_node(lsit);
  502.              if (v!=nil) New_Edge(SUB,v,w,lsit);
  503.              set_left_node(lsit,w);
  504.             }
  506.            lsit = Y_structure->key(sit);
  507.            lpred=Y_structure->key(sitpred);
  508.            sitpredpred = Y_structure->pred(sit);
  509.            sitsucc=Y_structure->succ(sitpred);
  512.            if (sitpredpred != nil)
  513.             { 
  514.               lpredpred=Y_structure->key(sitpredpred);
  516.               if ((pqit = Y_structure->inf(sitpredpred)) != nil)
  517.                 delete Xdelete(pqit);
  519.               Y_structure->change_inf(sitpredpred,nil);
  522.               if (intersection(lpred,lpredpred,inter))
  523.                 Y_structure->change_inf(sitpredpred,
  524.                                        Xinsert(sitpredpred,inter));
  525.              }
  528.            if (sitsucc != nil)
  529.             {
  530.               lsucc=Y_structure->key(sitsucc);
  532.               if ((pqit = Y_structure->inf(sitpred)) != nil)
  533.                 delete Xdelete(pqit);
  534.                  
  535.               Y_structure->change_inf(sitpred,nil);
  537.               if (intersection(lsucc,lsit,inter))
  538.                   Y_structure->change_inf(sit,Xinsert(sit,inter));
  539.              }
  542. // reverse the subsequence sit, ... ,sitpred  in the Y_structure
  544.            x_sweep = get_x(p);
  545.            y_sweep = get_y(p);
  547.            Y_structure->reverse_items(sit,sitpred);
  549.           delete p;
  551.          } // intersection
  553.   }
  555.   X_structure->clear();
  557. }
  561. main()
  562. { 
  563.   int N = read_int("N = ");
  565.   list<segment> seglist1,seglist2;
  566.   
  567.   while (N--) 
  568.   { double x1 = rand_int(-1000,-100);
  569.     double y1 = rand_int(-1000,1000);
  570.     double x2 = rand_int(100,1000);
  571.     double y2 = rand_int(-1000,1000);
  572.     seglist1.append(segment(x1,y1,x2,y2));
  573.    }
  575.   GRAPH<point,segment> SUB0;
  576.   GRAPH<point,int>    SUB;
  578.   sortseq<sw_segment,pq_item>           rst_seq;
  580.  _sortseq<sw_segment,pq_item,ab_tree>   abt_seq;
  581.  _sortseq<sw_segment,pq_item,skiplist>  skip_seq;
  583.   priority_queue<seq_item,sw_point>     Q;
  586.   float T;
  588.   T = used_time();
  589.   cout << "SWEEP_SEGMENTS: ";
  590.   cout.flush(); 
  591.   SWEEP_SEGMENTS(seglist1,SUB0);
  592.   cout<< string(" # = %d time = %6.2f sec",SUB0.number_of_nodes(), used_time(T));
  593.   newline;
  595.   T = used_time();
  596.   cout << "rs_tree:        ";
  597.   cout.flush(); 
  598.   Sweep_Segments(seglist1,seglist2,SUB,Q,rst_seq);
  599.   cout<< string(" # = %d time = %6.2f sec",SUB.number_of_nodes(), used_time(T));
  600.   newline;
  602.   SUB.clear();
  604.   T = used_time();
  605.   cout << "ab_tree:        ";
  606.   cout.flush(); 
  607.   Sweep_Segments(seglist1,seglist2,SUB,Q,abt_seq);
  608.   cout<< string(" # = %d time = %6.2f sec",SUB.number_of_nodes(), used_time(T));
  609.   newline;
  611.   SUB.clear();
  613.   T = used_time();
  614.   cout << "skiplist:       ";
  615.   cout.flush(); 
  616.   Sweep_Segments(seglist1,seglist2,SUB,Q,skip_seq);
  617.   cout<< string(" # = %d time = %6.2f sec",SUB.number_of_nodes(), used_time(T));
  618.   newline;
  620. return 0;
  622. }