home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / diverses / leda / prog / graphics / mulmuley.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-11-15  |  15.6 KB  |  716 lines
  1.  
  2. //------------------------------------------------------------------------------
  3. //
  4. // Mulmuley's randomized algorithm to compute the intersections of a set of 
  5. // line segments 
  6. //
  7. // implemented by Bettina Preis
  8. //
  9. //------------------------------------------------------------------------------
  10.  
  11.  
  12. /*       1.    read segments
  13.          2.    build frame
  14.          3.    draw vertical attachments through all endpoints
  15.          4.    randomly select segment
  16.          5.    start on a endpoint,
  17.                repeat
  18.                   FACE TRAVERSAL
  19.                and
  20.                   FACE TRANSITION,
  21.                until other endpoint reached
  22.          6.    repeat 4. and 5., until all segments are processed
  23. */
  24.  
  25.  
  26.  
  27. #include <LEDA/list.h>
  28. #include <LEDA/graph.h>
  29. #include <LEDA/plane.h>
  30. #include <LEDA/d_array.h>
  31. #include <LEDA/window.h>
  32.  
  33. #define eps 0.00001
  34.  
  35. window W;
  36.  
  37. double xmin, xmax, ymin, ymax;
  38.  
  39. edge frame_e_lower,frame_e_upper;
  40.  
  41.  
  42. int WAIT = 300;
  43.  
  44. void DRAW_EDGE(point a, point b, color c = black) { W.draw_edge(a,b,c); }
  45. void DRAW_EDGE(segment s, color c = black)        { W.draw_edge(s,c); }
  46.  
  47. /*****************************************************************************/
  48.  
  49. /* 
  50.     declared in plane.h:
  51.  
  52. int compare(point& p1,point& p2)
  53. {
  54.   int result = compare(p1.xcoord(),p2.xcoord());
  55.   if (result == 0)
  56.     return compare(p1.ycoord(),p2.ycoord());
  57.   else return result;
  58. }
  59. */
  60.  
  61. /*****************************************************************************/
  62.  
  63. int compare(segment& s1,segment& s2)
  64. {
  65.   int result = compare(s1.start(),s2.start());
  66.   return result;
  67. }
  68.  
  69. /*****************************************************************************/
  70.  
  71. int cmp_points(point& p1,point& p2)
  72. { return compare(p1.xcoord(),p2.xcoord()); }
  73.  
  74. /*****************************************************************************/
  75.  
  76. declare2(GRAPH,point,segment);
  77. GRAPH(point,segment) G;
  78.  
  79. declare2(d_array,point,node);
  80. d_array(point,node) node_of;
  81.  
  82. /*****************************************************************************/
  83.  
  84. struct frame 
  85. {
  86.   point a;              // left lower corner
  87.   point b;              // right upper corner
  88. };
  89.  
  90. frame Fr;
  91.  
  92. /*****************************************************************************/
  93.  
  94. node new_node(point p)
  95. {
  96.   node_of[p] = G.new_node(p);
  97.   return node_of[p];
  98. }
  99.  
  100. /*****************************************************************************/
  101.  
  102. void wait()
  103. {
  104.   for (int i=1;i<1000*WAIT;i++) {}
  105. }
  106.  
  107. /*****************************************************************************/
  108.  
  109. void frame_max(double wert,char ch)
  110. { 
  111.   // adapt frame frontiers to the read segments
  112.  
  113.   if (ch == 'x')
  114.   { 
  115.     if(wert<xmin)      xmin=wert;
  116.     else 
  117.       if(wert>xmax)    xmax=wert;
  118.   }
  119.   else
  120.   {
  121.     if(wert<ymin)      ymin=wert;
  122.     else 
  123.       if(wert>ymax)    ymax=wert;
  124.   }
  125. }
  126.  
  127. /*****************************************************************************/
  128.  
  129. void build_frame(double x1, double x2, double y1, double y2,frame& F)
  130. {
  131.   // build frame from frame frontiers 
  132.   edge edg;
  133.   point p1(x1,y1); 
  134.   point p2(x2,y1);
  135.   point p3(x2,y2);
  136.   point p4(x1,y2);
  137.   node kn1 = G.new_node(p1); W.draw_filled_node(x1,y1);
  138.   node kn2 = G.new_node(p2); W.draw_filled_node(x2,y1);
  139.   node kn3 = G.new_node(p3); W.draw_filled_node(x2,y2);
  140.   node kn4 = G.new_node(p4); W.draw_filled_node(x1,y2);
  141.   segment seg = segment(p1,p2);
  142.   frame_e_lower = G.new_edge(kn1, kn2, seg);
  143.   DRAW_EDGE(p1,p2);
  144.           seg = segment(p2,p3);
  145.           edg = G.new_edge(kn2, kn3, seg);
  146.   DRAW_EDGE(p2,p3);
  147.           seg = segment(p3,p4);
  148.   frame_e_upper = G.new_edge(kn3, kn4, seg);
  149.   DRAW_EDGE(p3,p4);
  150.           seg = segment(p4,p1);
  151.           edg = G.new_edge(kn4, kn1, seg);
  152.   DRAW_EDGE(p4,p1);
  153.   F.a = p1;
  154.   F.b = p3;
  155. }
  156.  
  157. /*****************************************************************************/
  158.  
  159. bool vertical(segment seg)
  160. {
  161.   /* true, if segment vertical
  162.      false, if not              */
  163.  
  164.   double diff = seg.xcoord1() - seg.xcoord2();
  165.  
  166.   return (fabs(diff) < eps);
  167. }
  168.  
  169. /*****************************************************************************/
  170.  
  171. void make_vertical_attachment(point pkt,edge e_lower,edge e_upper)
  172. {
  173.   /* adds the endpoint and the vertical edges from the endpoint to the frame 
  174.      and vice versa to the graph
  175.      furthermore, splits the upper and lower frame edges on the points of
  176.      attachment                                                             
  177.   */
  178.  
  179.   point p1,p2;
  180.   segment seg;
  181.   node kn,kn1,kn2;
  182.   edge edg;
  183.   
  184.   kn = new_node(pkt);
  185.   W.draw_filled_node(pkt);
  186.  
  187.   p1=point(pkt.xcoord(),Fr.b.ycoord());
  188.   kn1 = G.new_node(p1);              // upper node of attachment
  189.   W.draw_filled_node(p1);
  190.  
  191.   p2 = point(pkt.xcoord(),Fr.a.ycoord());
  192.   kn2 = G.new_node(p2);              // lower node of attachment
  193.   W.draw_filled_node(p2);
  194.  
  195.  
  196.   // new vertical edges
  197.  
  198.   seg = segment(p1,pkt);
  199.   edg = G.new_edge(kn1,kn,seg);
  200.   DRAW_EDGE(p1,pkt);
  201.  
  202.   seg = segment(pkt,p2);
  203.   edg = G.new_edge(kn,kn2,seg);
  204.   DRAW_EDGE(pkt,p2);
  205.  
  206.   seg = segment(p2,pkt);
  207.   edg = G.new_edge(kn2,kn,seg);
  208.   DRAW_EDGE(p2,pkt);
  209.   
  210.   seg = segment(pkt,p1);
  211.   edg = G.new_edge(kn,kn1,seg);
  212.   DRAW_EDGE(pkt,p1);
  213.  
  214.   wait();
  215.  
  216.   // new upper frame edges :
  217.  
  218.   DRAW_EDGE(G[e_upper],white);
  219.   seg = segment(p1,G[e_upper].end());
  220.   edg = G.new_edge(kn1,G.target(e_upper),seg);
  221.   DRAW_EDGE(p1,G[G.target(e_upper)]);
  222.   seg = segment(G[e_upper].start(),p1);
  223.   edg = G.new_edge(G.source(e_upper),kn1,seg);
  224.   DRAW_EDGE(G[G.source(e_upper)],p1);
  225.   G.del_edge(e_upper);
  226.   frame_e_upper = edg;
  227.  
  228.   wait();
  229.  
  230.   // new lower frame edges :
  231.  
  232.   DRAW_EDGE(G[e_lower],white);
  233.   seg = segment(G[e_lower].start(),p2);
  234.   edg = G.new_edge(G.source(e_lower),kn2,seg);
  235.   DRAW_EDGE(G[G.source(e_lower)],p2);
  236.   seg = segment(p2,G[e_lower].end());
  237.   edg = G.new_edge(kn2,G.target(e_lower),seg);
  238.   DRAW_EDGE(p2,G[G.target(e_lower)]);
  239.   G.del_edge(e_lower);
  240.   frame_e_lower = edg;
  241.  
  242.   wait();
  243. }
  244.  
  245. /*****************************************************************************/
  246.  
  247. void insert_in_edge(node kn,edge edg,edge& f_edge,edge& s_edge)
  248. {
  249.   /* split edge edg on node kn into two edges f_edge (first edge) 
  250.      and s_edge (second edge)  
  251.   */
  252.  
  253.   segment f_seg,s_seg;
  254.  
  255.   f_seg = segment(G[G.source(edg)],G[kn]);
  256.   f_edge = G.new_edge(G.source(edg),kn,f_seg);
  257.   DRAW_EDGE(f_seg);
  258.   
  259.   wait();
  260.  
  261.   s_seg = segment(G[kn],G[G.target(edg)]);
  262.   s_edge = G.new_edge(kn,G.target(edg),s_seg);
  263.   DRAW_EDGE(s_seg);
  264.  
  265.   wait();
  266.  
  267.   G.del_edge(edg);
  268. }
  269.  
  270. /*****************************************************************************/
  271.  
  272. point face_traversal(segment seg, edge edg, edge& succ_e, edge& pred_e)
  273. {
  274.   /* travelling through a face until the intersection 
  275.      between the segment and a face border is found 
  276.      succ_e is the intersected face border
  277.      pred_e is the predecessor border of succ_edge
  278.   */
  279.  
  280.   point inter;
  281.   double angl;
  282.   double min_angle;
  283.   edge edg2;
  284.  
  285.  
  286.   do 
  287.   {
  288.     min_angle = M_PI;
  289.     forall_adj_edges(edg2,G.target(edg))
  290.     { 
  291.       // angle between edg2 and edg
  292.       segment s = G[edg2];
  293.       segment t = G[edg];
  294.  
  295.       angl = s.angle(t);
  296.  
  297. //      angl = G[edg2].angle(G[edg]);
  298.  
  299.       if ((angl > -M_PI) && (angl <= min_angle))
  300.       {
  301.         min_angle = angl;
  302.         succ_e = edg2;
  303.       }
  304.     }
  305.     pred_e = edg;
  306.     edg = succ_e;
  307.   } 
  308.   while (! seg.intersection(G[succ_e],inter)); 
  309.   return inter;
  310. }
  311.  
  312. /*****************************************************************************/
  313.  
  314. edge face_transition(segment seg, edge edg)
  315. {
  316.   /* find next face that must be traversed
  317.      return first border for next face_traversal   */
  318.   
  319.   double angl;
  320.   point inter;
  321.   edge succ_edge, edg2;
  322.   edge return_edge = nil;
  323.   edge corner_edge = nil;
  324.   edge act_edge = nil;
  325.   
  326.   do 
  327.   { 
  328.     succ_edge = nil;
  329.     edg2 = G.first_adj_edge(G.target(edg));
  330.  
  331.     for (;;)
  332.     {
  333.       angl = G[edg2].angle(G[edg]);
  334.       double diff = fabs(angl) - M_PI;
  335.       if (fabs(diff) < eps)
  336.       { 
  337.         return_edge = edg2;
  338.         break;
  339.       }
  340.       if (fabs(angl) < eps)
  341.         succ_edge = edg2; 
  342.       edg2 = G.adj_succ(edg2);
  343.     }
  344.     edg = succ_edge;
  345.  
  346.   } while (return_edge == nil);
  347.  
  348.   act_edge = return_edge;
  349.  
  350.   while(! seg.intersection(G[act_edge],inter))
  351.   {
  352.     edg2 = G.first_adj_edge(G.target(edg));
  353.     while (edg2 != nil)
  354.     { 
  355.       angl = G[edg2].angle(G[act_edge]);
  356.       if (fabs(angl) < eps )
  357.       {
  358.         act_edge = edg2;
  359.         break;
  360.       }
  361.       edg2 = G.adj_succ(edg2);
  362.     }
  363.   }
  364.  
  365.   return act_edge;
  366. }
  367.         
  368. /*****************************************************************************/
  369.  
  370. void delete_edges(list(edge)& e_list)
  371. {
  372.   // deletes upper vertical attachments
  373.  
  374.   list_item l_it;
  375.   edge edg;
  376.  
  377.   l_it = e_list.first();
  378.  
  379.   while (l_it != nil)
  380.   {
  381.     edg =  e_list.del(l_it);
  382.     G.del_edge(edg);
  383.     l_it = e_list.first();
  384.   }
  385. }
  386.  
  387. /*****************************************************************************/
  388. /*****************************************************************************/
  389.  
  390. list(point) fmul(list(segment)& L)
  391. {
  392.   /* for all segments : xcoord of first endpoint is smaller then
  393.                           xcoord of second endpoint                          */ 
  394.  
  395.   double angl;
  396.   point p,p_pred,p_succ,point_of_inter;
  397.   list_item l_it;
  398.   list(point) p_list,inter_p_list;
  399.   list(edge) edg_list;
  400.   segment seg,seg1,seg2;
  401.   node start_node,kn;
  402.   edge start_edge, succ_edge, pred_edge, first_edge, first2_edge, second_edge, second2_edge, edg1, edg2, edg4, edg8;
  403.  
  404.   /* initialisation of xmin, xmax, ymin, ymax (frame frontiers)
  405.        by the first segment                                                  */
  406.   
  407.   seg = L.head();
  408.  
  409.   xmin = seg.xcoord1();
  410.   xmax = seg.xcoord2();
  411.   
  412.   double r1 = seg.ycoord1();
  413.   double r2 = seg.ycoord2();
  414.  
  415.   if (r1 < r2)
  416.   {
  417.     ymin = r1;
  418.     ymax = r2;
  419.   }
  420.   else
  421.   {
  422.     ymin = r2;
  423.     ymax = r1;
  424.   }
  425.  
  426.   // expand frame frontiers
  427.   forall(seg,L)
  428.   {
  429.     if (seg.vertical()) error_handler(1,"vertical segment");
  430.     
  431.     frame_max(seg.xcoord1(),'x');
  432.     frame_max(seg.ycoord1(),'y');
  433.     frame_max(seg.xcoord2(),'x');
  434.     frame_max(seg.ycoord2(),'y');
  435.   }
  436.  
  437.  
  438.   // initialisation of the frame
  439.   build_frame(xmin-1,xmax+1,ymin-1,ymax+1,Fr);
  440.  
  441.   wait();
  442.  
  443.   /* make a list of all endpoints */
  444.  
  445.   forall(seg,L)
  446.   {
  447.     p = seg.start();
  448.     l_it = p_list.append(p);
  449.     p = seg.end();
  450.     l_it = p_list.append(p);
  451.   }
  452.  
  453.   // sort endpoints from left to right 
  454.   p_list.sort(cmp_points);
  455.  
  456.   // initialisation of the first ("stripped") partition
  457.   forall(p,p_list)
  458.   {
  459.     make_vertical_attachment(p,frame_e_lower,frame_e_upper);
  460.  
  461.   }
  462.  
  463.  
  464.   // algorithm begins
  465.  
  466.   forall(seg,L)
  467.   {
  468.     p = seg.end();                         // endpoint to be inserted
  469.  
  470.     start_node = node_of[p];
  471.     start_edge = nil;
  472.     edg8 = G.first_adj_edge(start_node);
  473.     while (edg8 != nil)
  474.     {
  475.       if (G[edg8].ycoord1() < G[edg8].ycoord2())
  476.       {
  477.         start_edge = edg8;
  478.         break;
  479.       }
  480.       else
  481.         edg8 = G.adj_succ(edg8);
  482.     } 
  483.  
  484.  
  485.     do
  486.     { 
  487.       point_of_inter = face_traversal(seg,start_edge,succ_edge,pred_edge);
  488.  
  489.       if (point_of_inter != seg.start())
  490.       { 
  491.  
  492.         kn = G.new_node(point_of_inter);
  493.         W.draw_filled_node(point_of_inter);
  494.         start_edge = nil;
  495.  
  496.         insert_in_edge(kn,succ_edge,first_edge,second_edge);
  497.    
  498.  
  499.         // new edges through face
  500.     
  501.         seg1 = segment(G[start_node],point_of_inter);
  502.         edg1 = G.new_edge(start_node,kn,seg1);
  503.         DRAW_EDGE(seg1);
  504.  
  505.         wait();
  506.     
  507.         seg2 = segment(point_of_inter,G[start_node]);
  508.         edg2 = G.new_edge(kn,start_node,seg2);
  509.         DRAW_EDGE(seg2);
  510.      
  511.         wait();
  512.       }
  513.  
  514.       else
  515.       {
  516.         kn = node_of[seg.start()];
  517.  
  518.         // new edges through face
  519.     
  520.         seg1 = segment(G[start_node],seg.start());
  521.         edg1 = G.new_edge(start_node,kn,seg1);
  522.         DRAW_EDGE(seg1);
  523.  
  524.         wait();
  525.     
  526.         seg2 = segment(point_of_inter,G[start_node]);
  527.         edg2 = G.new_edge(kn,start_node,seg2);
  528.         DRAW_EDGE(seg2);
  529.      
  530.         wait();
  531.         break;
  532.       }
  533.  
  534.       edg4 = face_transition(seg,second_edge);
  535.  
  536.  
  537.       insert_in_edge(kn,edg4,first2_edge,second2_edge);
  538.  
  539.       start_edge = second2_edge;
  540.  
  541.  
  542.       if (vertical(G[start_edge]))
  543.       { 
  544.         // cut vertical attachment 
  545.  
  546.         angl = G[first_edge].angle(G[pred_edge]);
  547.         if (fabs(angl) > eps)
  548.         {
  549.           if (! node_of.defined(G[G.target(pred_edge)]))
  550.           {
  551.             // delete upper vert_att
  552.  
  553.             edg_list.append(first_edge);
  554.             edg_list.append(second2_edge);
  555.             DRAW_EDGE(G[first_edge],white);
  556.             DRAW_EDGE(G[second2_edge],white);
  557.  
  558.             wait();
  559.           }
  560.           else
  561.           {
  562.             // delete lower vert_att
  563.           
  564.             DRAW_EDGE(G[second_edge],white);
  565.             DRAW_EDGE(G[first2_edge],white);
  566.             G.del_edge(second_edge);
  567.             G.del_edge(first2_edge);
  568.  
  569.             wait();
  570.           }
  571.         }
  572.         else
  573.         {
  574.           // delete lower vert_att
  575.           
  576.           DRAW_EDGE(G[second_edge],white);
  577.           DRAW_EDGE(G[first2_edge],white);
  578.           G.del_edge(second_edge);
  579.           G.del_edge(first2_edge);
  580.           wait();
  581.         }
  582.       }  
  583.       else
  584.  
  585.         inter_p_list.append(point_of_inter);
  586.  
  587.       start_node = kn;
  588.  
  589.     } while (start_edge != nil); 
  590.  
  591.     delete_edges(edg_list);
  592.     edg_list.clear();
  593.   }
  594.  
  595.  
  596.   G.print();
  597.  
  598.   return inter_p_list;
  599. }
  600.  
  601. /*****************************************************************************/
  602. /*****************************************************************************/
  603.  
  604. main()
  605. {
  606.   int    input = 0; 
  607.   int    grid_width = 0;;
  608.   int    N = 10;
  609.   string file_name;
  610.  
  611.   double x1,y1,x2,y2;
  612.   point p,q;
  613.   segment seg;
  614.   list(segment) seg_list;
  615.  
  616.   panel P("Line Segment Intersection (Mulmuley)");
  617.  
  618.   P.choice_item("INPUT", input, "mouse","random","file");
  619.  
  620.   P.int_item("GRID", grid_width, 0,100,20);
  621.   P.int_item("#segments",N);
  622.   P.int_item("WAIT",WAIT,0,300);
  623.  
  624.   P.string_item("FILE",file_name);
  625.  
  626.   P.open();
  627.  
  628.   W.init(-1000,1000,-1000,grid_width);
  629.  
  630.  
  631.   /* make a list of line_segments */
  632.  
  633.   switch(input) {
  634.  
  635.  
  636.   case 0: { // mouse
  637.             W.set_line_style(dashed);
  638.         
  639.             while(W >> seg)
  640.             { if (seg.xcoord1() > seg.xcoord2())
  641.                     seg=segment(seg.end(),seg.start());
  642.               W.draw_segment(seg);
  643.               seg_list.append(seg);
  644.             }
  645.             W.set_line_style(solid);
  646.             break;
  647.           }
  648.  
  649.  
  650.   case 1: { // random
  651.             W.set_line_style(dashed);
  652.             init_random();
  653.             for(int i=0;i<N;i++)
  654.             { x1 = random(-900,-100);
  655.               y1 = random(-900, 900);
  656.               x2 = random( 100, 900);
  657.               y2 = random(-900, 900);
  658.        
  659.               p = point(x1,y1);
  660.               q = point(x2,y2);
  661.               seg = segment(p,q);
  662.               seg_list.append(seg);
  663.       
  664.               W.draw_segment(seg);
  665.             }
  666.             W.set_line_style(solid);
  667.             break;
  668.            }
  669.  
  670.   case 2: { // file 
  671.             file_istream F(file_name);
  672.       
  673.             double x_min=0,x_max=0,y_min=0;
  674.       
  675.             while (F >> x1 >> y1 >> x2 >> y2)
  676.             {
  677.               x_min = Min(x_min,x1);
  678.               x_max = Max(x_max,x2);
  679.               y_min = Min(y_min,y1);
  680.       
  681.               p = point(x1,y1);
  682.               q = point(x2,y2);
  683.               seg = segment(p,q);
  684.               seg_list.append(seg);
  685.               break;
  686.             }
  687.       
  688.             W.init(x_min-1,x_max+1,y_min-1);
  689.           }
  690.  
  691. }//switch
  692.    
  693.  
  694.  
  695.   W.set_node_width(4);
  696.  
  697.  
  698.   list(point) point_list = fmul(seg_list);
  699.  
  700.   point_list.print("point_list = ");
  701.   W.read_mouse();
  702.   
  703.   W.clear();
  704.  
  705.   forall(seg,seg_list) W.draw_segment(seg);
  706.  
  707.   forall(p,point_list) W.draw_filled_node(p);
  708.  
  709.   W.read_mouse();
  710.  
  711. }
  712.  
  713. /*****************************************************************************/
  714. /*****************************************************************************/
  715. /*****************************************************************************/
  716.