  % Title "Intersections Demo"
 '% Animated
 '% ScaleUp
 '% PlaceCentral
 '% AnchorCentral
  uses consoleG
  uses GLO2\GeoPlanar
  uses chaos

  sub drawintersections()
  =======================


  type linex
    float x1,y1 'first point
    float x2,y2 'second point
    float le    'length
    int dl      'deletion flag
    int p1,p2   'point indexes
  end type

  string tab="    "
  point  p
  int    a
  float le

  seed=0x12379823

  pointex pn[100]
  linex dn[10000]
  float r[16]
  indexbase 0

  int i,j,k,q
  '
  '
  /*
  'grid
  q=36
  for i=0 to 5
    for j=0 to 5
      r[0]= i*0.3'+rnd()*.02
      r[1]=-j*0.3'+rnd()*.02
      pn[i*6+j]={ r[0],r[1] }
    next
  next
  '*/
  '
  /*
  'random
  q=24
  for i=0 to q-1
    r[0]=rnd
    r[1]=rnd
    pn[i]={ r[0],r[1]}
  next
  */
  '
  '/*
  'POLYGONS
  % q=31
  point qn[q]
  qn[0]={0,0}
  CreatePolypoints qn[1],0,0,15,1,1,0
  CreatePolypoints qn[16],0.03,0.04,15,0.7,0.75,12.5

  for i=0 to q-1
    pn[i]={ qn[i].x,qn[i].y }
  next
  'perturbPoints @pn[1],q-1,0.05,0.01

  k=0
  float x,y
  int qh=(q-1)/2
  int mi,mj
  j=0
  'RADIALS
  for i=1 to qh
    mi=i+qh
    mj=0
    le=sqdist pn[mi], pn[mj]
    dn[k]={ pn[mi].x,pn[mi].y,pn[mj].x,pn[mj].y,le,0,mi,mj }
    k++
  next
  '
  'CIRCUMS
  '
  for i=1 to qh
    j=i+qh
    le=sqdist pn[i], pn[j]
    dn[k]={ pn[i].x,pn[i].y,pn[j].x,pn[j].y,le,0,i,j }
    k++
    j=i+qh-1 : if i=1 then j=qh+qh 'overlap
    le=sqdist pn[i], pn[j]
    dn[k]={ pn[i].x,pn[i].y,pn[j].x,pn[j].y,le,0,i,j }
    k++
    '
    j=i+1 : if i=qh then j=1 'overlap
    le=sqdist pn[i], pn[j]
    dn[k]={ pn[i].x,pn[i].y,pn[j].x,pn[j].y,le,0,i,j }
    k++
    mi=i+qh
    mj=j+qh
    le=sqdist pn[mi], pn[mj]
    dn[k]={ pn[mi].x,pn[mi].y,pn[mj].x,pn[mj].y,le,0,mi,mj }
    k++
  next
'*/

  
  /*
  'ALL PERMUTATIONS
  k=0
  float x,y
  for i=0 to q-2
    for j=i+1 to q-1
      le=sqdist pn[i], pn[j]
      dn[k]={ pn[i].x,pn[i].y,pn[j].x,pn[j].y,le,0,i,j }
      k++
    next
  next
  */
  '
  macro compare(c,a,b)
    if dn(a).le>dn(b).le
      c=1
    endif
  end macro
  '
  NewSortIndex(idx, k, kf, filter, compare)
  '
  int mi,mj
  int ct
  line *dni,*dnj 'pointers for intersection test
  for i=0 to k-2
    mi=idx[i]
    @dni=@dn[mi]
    if dn[mi].dl=0
      for j=i+1 to k-1
        mj=idx[j]
        @dnj=@dn[mj]
        if dn[mj].dl=0
          'mark longer line for deletion
          a=intersected dni,dnj,p,1 'insiderange
          if a=16 'excluding ends <16
            ct++
            'dn[mj].dl=1
          endif
        endif
      next
    endif
  next
  '

  'shading
  'flat 'default

  '
  color 0,1,1
  printl "intersect" tab k tab ct

  scale 10,10
  move 1,-1
  '
  color 0,1,1
  int e
  e=k-1
  linex dln[1000]
  int dc
  int p1,p2,n
  pointex*ph
  for i=0 to e
    mi=idx[i]
    if not dn[mi].dl
      dln[dc]=dn[mi]
      p1=dln[dc].p1 'first point index
      p2=dln[dc].p2 'second point index
      @ph=@pn[p1] 'first hub
      n=ph.n
      ph.pp[n]=p2 'store other hub index
      ph.n++
      @ph=@pn[p2] 'second hub
      n=ph.n
      ph.pp[n]=p1 'store other hub index
      ph.n++
      'color 1.0,1.0,0
      color 0.7,0,0
      drawline @dln[dc] 'DELAUNAY LINES
      dc++ 'index next dln
    endif
  next
  '
  'SORT EACH HUB LINES IN ASCENDING ANGLE ORDER
  pointex*pt
  float ang[32],an,xe,ye
  int m
  string s,t
  'HUBS
  for i=1 to q-1
    @pt=@pn[i]
    n=pt.n
    x=pt.x
    y=pt.y
    'COLLECT LINE ANGLES
    t=""
    for j=0 to n-1
      m=pt.pp[j]
      xe=pn[m].x
      ye=pn[m].y
      an=deg(atan(ye-y, xe-x))
      if an<0 then an+=360
      ang[j]=an
    next
    'SORT INTO ASCENDING ORDER
    for j=0 to n-2
      for k=j+1 to n-1
        if ang[j]>ang[k]
          swap ang[k],ang[j]
          swap pt.pp[k],pt.pp[j] 
        endif
      next
    next
    't="hub  " i "  " n " lines" cr
    for j=0 to n-1
      m=pt.pp[j]
      xe=pn[m].x
      ye=pn[m].y
      line l1={x,y,xe,ye}
      'drawline @l1
      't+=" " j "  " int(ang(j)) cr
    next
    s+=t
    'GET MIDPOINTS OF EACH ADJACENT LINE
    int m1,m2
    point qc1,qc2,qi
    point q1,q2,q3,q4,q5
    q1={x,y}
    pt.pp[n]=pt.pp[0] 'overlap
    ang[n]=ang[0]+360 'overlap
    for j=0 to n-1
      k=j+1
      m1=pt.pp[j]
      m2=pt.pp[k]
      q2={pn[m1].x, pn[m1].y}
      q3={pn[m2].x, pn[m2].y}
      'line l1={x,y,q2.x,q2.y}
      'drawline @l1
      midpoint q1, q2, qc1
      midpoint q1, q3, qc2
      pointsize 3
      drawpoint @qc1
      drawpoint @qc2
      '
      'EXCLUDE
      if ang[k]-ang[j]>=180
        exit for 'j
      end if
      'GET INTERSECTION POINT OF EACH LINE PAIR  PERPENDICULARS
      rot90  qc1,q1,q4
      rotm90 qc2,q1,q5
      line l1,l2,l3,l4
      l1={qc1.x,qc1.y,q4.x,q4.y}
      l2={qc2.x,qc2.y,q5.x,q5.y}
      a=intersect l1,l2,qi
      if a'=16
        if inrange -1,1,qi.x
          if inrange -1,1,qi.y
            pointsize 4
            color 1,0,0
            drawpoint @qi
            'VORONOI LINES
            color 0,1,1
            'color 0.5,0,0
            l3={qc1.x,qc1.y,qi.x,qi.y}
            drawline @l3 
            l4={qc2.x,qc2.y,qi.x,qi.y}
            drawline @l4
          endif
        endif  
      endif
    next 'j
  next 'i
'putfile "t.txt",s
s=""
 
          
    

  pointsize 7
  color 1,1,0
  '
  for i=0 to q-1
    drawpoint @pn[i]
  next
  del idx
  end sub


  sub main()
  ==========
  string s
  cls 0,0,0
  pushstate
  scale 2,2
  'typeface=1
  printl "Intersection Points"
  popstate
  '
  move 0,-2
  '
  PushState
  move 1,-1
  DrawIntersections
  PopState
  '
  
  end sub

  EndScript
