/*
   Copyright 2010 Aaron J. Radke

   Licensed under the Apache License, Version 2.0 (the "License");
   you may not use this file except in compliance with the License.
   You may obtain a copy of the License at

       http://www.apache.org/licenses/LICENSE-2.0

   Unless required by applicable law or agreed to in writing, software
   distributed under the License is distributed on an "AS IS" BASIS,
   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   See the License for the pecific language governing permissions and
   limitations under the License.
*/
package cc.drx

import java.awt.image.BufferedImage

//TODO 2017-07-30("Should use an arrow for vector if needed else this is too specific for a general lib", "0.3.0")
case class State(pos:Vec=Vec.zero, vel:Vec=Vec.zero, acc:Vec=Vec.zero, time:Double=0.0){
   def pos(posUpdate:Vec):State      =  copy(pos=posUpdate)
   def vel(velUpdate:Vec):State      =  copy(vel=velUpdate)
   def acc(accUpdate:Vec):State      =  copy(acc=accUpdate)
   def time(timeUpdate:Double):State =  copy(time=timeUpdate)

   def step(vel:Vec, dt:Double):State = {
      val posNext = pos + vel*dt
      State(posNext,vel,acc,time + dt)
   }
}

object Bezier{
   def apply(a:Vec,b:Vec):Bezier = Bezier(a, a, b, b)
   def apply(line:Line):Bezier = Bezier(line.a, line.b)
   def apply(lineA:Line,lineB:Line):Bezier = Bezier(lineA.a,lineA.b,  lineB.b,lineB.a)
}
case class Bezier(override val a:Vec, ca:Vec, cb:Vec, override val b:Vec) extends Segment with Shape.Open{
   def draw(implicit g:DrawContext):Unit = g ! this

   private lazy val c = (a + b)/2 //half point
   def hh = Bezier(a, Vec(c.x,a.y), Vec(c.x,b.y), b) //set control points for horizontally exit to horizontal enter
   def vv = Bezier(a, Vec(a.x,c.y), Vec(b.x,c.y), b) //set control points for vertical exit to vertical enter
   def hv = Bezier(a, Vec(c.x,a.y), Vec(b.x,c.y), b) //set control points for horizontal exit to vertical enter
   def vh = Bezier(a, Vec(a.x,c.y), Vec(c.x,b.y), b) //set control points for vertical exit to horizontal enter

   def last = b

   2017-07-30("since this extends lerp, the apply doesn't do what you may think use applyNonLinear to be more specific","v0.2.15")
   def apply(t:Double):Vec = applyNonLinear(t)
   /**parametric evaluation to get the location point*/
   def applyNonLinear(t:Double):Vec = {
      val q = 1d - t; val q2 = q*q; val q3 = q2*q;
                      val t2 = t*t; val t3 = t2*t;
      a*(q3) + ca*(3*q2*t) + cb*(3*q*t2) + b*t3
   }
   /**parametric evaluation of to get the slope at a point*/
   def gradNonLinear(t:Double):Vec = {
      val q = 1d - t
      (ca-a)*(3*q*q) + (cb-ca)*(6*q*t) + (b-cb)*(3*t*t)
   }

   private val N = 100
   private val dt = 1.0/N
   private lazy val distIndex:Array[Double] = Array.iterate((0.0,0.0,a),N+1){case (d,t,a) =>
       val tb = t + dt
       val b = applyNonLinear(tb)
       (d+(a dist b), tb, b)
   }.map{_._1}
   lazy val dist:Double = distIndex.last
   private def nonLinearLookup(td:Double):Double = {
      if      (td <= 0) 0d
      else if (td >= 1) 1d
      else {
        val d = td*dist
        val i = distIndex.searchBy(d)(identity) - 1
        val t0 = i.toDouble/N
        val d0 = distIndex(i)
        val dd = distIndex(i+1) - d0

        val pd = d - d0
        val t = t0 + pd*(dt/dd) //interpolation between the lookups to prevent jumpy behavior
        t
      }
   }
   def lerp(td:Double):Vec = applyNonLinear(nonLinearLookup(td))
   def grad(td:Double):Vec =  gradNonLinear(nonLinearLookup(td))
}

object Line {
  // def horiz(a:Vec,width:Double):Line = a ~> width.x
  // def vert(a:Vec,height:Double):Line = a ~> height.y
  // def unapply(line:Line):Option[(Vec,Vec)] = Some((line.a, line.b))
  def apply(a:Vec, b:Vec):Line = new Line(a,b)
  def unapply(line:Line):Option[(Vec,Vec)] = Some( line.toTuple )
}
class Line(override val a:Vec, override val b:Vec) extends Segment with Lerp[Vec] with Shape.Open{
   def draw(implicit g:DrawContext):Unit = g ! this
   2017-07-30("use rect for a more direct meaning since bound has other meanings in drx","v0.2.13")
   lazy val bounds:Rect = rect
   lazy val rect:Rect = Rect((a point b){_ min _}, (a point b){_ max _})
   def vertices = List(a,b)
   def toBound:Bound[Vec] = Bound(a,b)
   def toTuple:(Vec,Vec) = (a,b)
   def toArrow:Arrow = Arrow(this)
   def swap = Line(b,a)

   def round = Line(a.round, b.round)

   def extrude(v:Vec):Poly = Poly.convex(List(a,b, a+v, b+v)) //FIXME verify

   lazy val dist:Double = a dist b
   lazy val cm:Vec = lerp(0.5)
   lazy val delta:Vec = b - a

   lazy val slope = (b.y - a.y)/(b.x - a.x) //rise over run
   lazy val yIntercept = a.y - slope*a.x //y = m*x + b    =>  b = y - m*x
   def projectedIntersection(that:Line):Option[Vec]= //TODO maybe generalize this to a 3d case with the z axis and penrose inverse
      if(this.slope == that.slope) None
      else {
         val x = (that.yIntercept - this.yIntercept)/(this.slope - that.slope) // this.slope*x + this.yIntercept == that.slope*x + that.yIntercept
         Some(Vec(x, this.slope*x + this.yIntercept))
      }

   def bezier:Bezier = Bezier(a,b)

   def +(that:Vec):Line = Line(a+that,b+that)
   def -(that:Vec):Line = Line(a-that,b-that)
   def *(x:Double):Line = Line(cm + (a-cm)*x, cm + (b-cm)*x)
   def /(x:Double):Line = *(1.0/x)

   def last = b
   def lerp(t:Double):Vec = (a lerp b)(t)
   private lazy val gradConstant = (b-a).rzRight.unit
   def grad(t:Double):Vec = gradConstant
   //TODO support def ratioOf(v:Vec):Double =   //support LerpInv it may not be used but allows for scale

   /**point closest on a segment of two points (a to b) to point p*/
   def closest(p:Vec):Vec = {
      val l2 = b dist2 a // {val bma = b-a; bma*bma}
      if(l2 == 0.0) a  // a == b case
      else {
         val s = (p-a) * (b-a) / l2
         if     (s < 0.0)   a            //point a
         else if(s > 1.0)   b            //point b
         else               a + (b-a)*s  //projection onto a & b
      }
   }
   def dist(p:Vec):Double = p dist closest(p)
   def orthoDist(p:Vec):Double = {
      val n = (b-a).unit //unit vector in direction of the line
      val amp = a-p
      (amp - (n*(amp*n))).norm
   }
   def x2d(that:Line):Double = (this.b-that.a) x2d (that.b-that.a)

   def ~(c:Vec):Poly = new Poly(Seq(a,b,c))
   def ^(that:Line):Seq[Vec] = { //intersection points (multiple points of colinear overlay)
      //a1 = p
      //b1 = p + r
      //a2 = q
      //b2 = q + s
      val p = this.a
      val r = this.b - this.a
      val q = that.a
      val s = that.b - that.a

      val qmp = q - p
      val rxs = r x2d s

      val qmpxr = qmp x2d r

      if(rxs == 0.0){//parallel
         if(qmpxr == 0.0){//colinear
            val proj = qmp*r
            if(0.0 <= proj && proj < r*r){ //overlapping
               Seq(
                  qmp + qmp.unit*proj  //TODO: need to verify and calculate other point, realy overlap is a line
               )
            }
            else                          Seq() //disjoint
         }
         else Seq() //parallel not-colinear
      } else {  //not parallel
         val t = qmp x2d s / rxs
         val u = qmp x2d r / rxs
         if(0.0 < t && t < 1.0 && 0.0 < u && u < 1.0) Seq(p + r*t) //intesection
         else Seq() // no intersection
      }

   }
   // def moveStartTo(v:Vec) = Line(0, b-a+v) // this - a + v  // Line(a-a,b-a)+v // = v ~> b

   def right(p:Vec):Boolean = ((b-a) x2d (p-a)) > 0
   def left(p:Vec):Boolean = !right(p)
   def label(text:String,offset:Double=0):Text = Text(text, cm + (b-a).rzRight.mag(offset))
}
object Arrow{
  def apply(line:Line):Arrow = Arrow(line.a, line.b)
  def apply(a:Vec,b:Vec):Arrow = new Arrow(a,b)
}
class Arrow(a:Vec,b:Vec) extends Line(a,b){
   override def draw(implicit g:DrawContext):Unit = g ! this
}


object Img{
  import java.awt
  import javax.swing

  val defaultBox = Rect(100,100)
  def apply(path:String):ImgFile = Img(File(path))
  def apply(file:File):ImgFile = Img(file, defaultBox) //TODO add autodetect size (but for now this feature would depends on the presence of p5 classes)
  def apply(file:File, box:Rect):ImgFile = ImgFile(file, box)
  def apply(img:java.awt.Image):ImgAwt = ImgAwt(toAwt(img))
  def apply(in:Input):ImgAwt = apply(toAwt(in))

  def findInput(file:File,exts:List[String] = List("jpg","png")):Option[Input] = {
    lazy val alts = for(
      dir <- List("img","fig","icon","src/main/resources","../shared/src/main/resources");
      ext <- exts
    ) yield File(dir)/file.companion(ext)
    def findLocal:Option[File] = if(file.isFile) Some(file) else alts.find(_.isFile)
    def findResource:Option[File] = if(Input.isResource(file)) Some(file)
                                    else alts.find(Input.isResource)
    findLocal.map(_.in) orElse findResource.flatMap(Input.resourceOption)
  }
  def loadAwt(f:File):Option[BufferedImage] = Img.findInput(f,List("png","jpg")).map(toAwt)
  def find(f:File):Option[ImgFile] = Img.findInput(f).map{_ => Img(f)}
  //TODO maybe make this depend on an IO src with a name to provide fast caching support without hashing
  //TODO generalize the Img and Svg types so an svg can be an img type


  //---image conversion to java types help
  def toAwt(in:Input):BufferedImage = toAwt(new swing.ImageIcon(in.toByteArray).getImage)
  def toAwt(img:java.awt.Image, bufferedImageType:Int = BufferedImage.TYPE_INT_ARGB):BufferedImage = img match {
    case bi:BufferedImage if bi.getType == bufferedImageType => bi
    case _ =>
      val rgb = new BufferedImage(img.getWidth(null), img.getHeight(null), BufferedImage.TYPE_INT_RGB)
      val g = rgb.createGraphics()
      g.drawImage(img,0,0,null)
      g.dispose()
      rgb
  }
  def size(img:BufferedImage):Vec = Vec( img.getWidth.toDouble,  img.getHeight.toDouble )
  def toARGB(img:java.awt.Image):BufferedImage = toAwt(img, BufferedImage.TYPE_INT_ARGB)
  /**force RGB format of buffered image without an alpha channel*/
  def toRGB( img:java.awt.Image):BufferedImage = toAwt(img, BufferedImage.TYPE_INT_RGB)

  //--- show/fast simple display utilities
  def show(img:BufferedImage):swing.JFrame = show(img,s"Show: $img")
  def show(img:BufferedImage,title:String="Show Image"):swing.JFrame = {
    //--alexander method
    val frame = new swing.JFrame(title)
    // frame.setIconImage(img) //make the popup icon be the view contents
    frame.setResizable(false) //because the image for now is not automatically resizable...
    frame.setDefaultCloseOperation(swing.WindowConstants.DISPOSE_ON_CLOSE)// HIDE_ON_CLOSE) //EXIT_ON_CLOSE)

    val icon = new swing.ImageIcon(img)
    val label = new swing.JLabel()
      label.setIcon(icon)
    frame.getContentPane().add(label, awt.BorderLayout.CENTER)

    frame.pack()
    frame.setLocationRelativeTo(null)
    frame.setVisible(true)
    frame
  }
}
trait Img extends Shape{
  //--required
  def box:Rect
  def moveTo(thatBox:Rect):Img
  def toAwt:BufferedImage

  //--derived
  def fitIn(thatBox:Rect):Img = moveTo(box fitIn thatBox)
  def save(file:File)(implicit g:DrawContext):Try[File] = g.save(this,file)
}
case class ImgFile(file:File,  box:Rect) extends Img{
  def draw(implicit g:DrawContext):Unit = g ! this
  def moveTo(thatBox:Rect) = copy(box = thatBox)
  def toAwt:BufferedImage = Img.loadAwt(file).get //TODO return empty buffer instead of throwing an error
}
object ImgAwt{
  def apply(img:BufferedImage):ImgAwt = {
    val size = Vec( img.getWidth.toDouble,  img.getHeight.toDouble )
    ImgAwt(img, Rect(Vec.zero, size) )
  }
}
case class ImgAwt(img:BufferedImage, box:Rect) extends Img{
  def moveTo(thatBox:Rect):ImgAwt = copy(box = thatBox)
  def draw(implicit g:DrawContext):Unit = g ! this
  def rawSize:Vec = Img.size(img) //TODO require for all Img objects
  def size:Vec = box.size
  def toAwt:BufferedImage = img
}


case class Svg(xml:String, box:Rect=Img.defaultBox){ //TODO think about this more. is a filepath, xml string, or bytestring or even byte array better storage?
   def draw(implicit g:DrawContext):Unit = g ! this
}
//TODO a img should have a type as SVG,png,jpg. This should be independent of the source
case class SvgPath(d:String){
  lazy val svg = Svg(s"<svg width='100' height='100'><path d='$d'/></svg>")
  def draw(implicit g:DrawContext):Unit = g ! this
}
trait MediaRemote{
  //--required
  def stop:Unit
  def pause:Unit
  def play:Unit
  def isPlaying:Boolean
  def speedup:Double
  def speedup_=(v:Double):Unit
  def cursor:Time
  def cursor_=(t:Time):Unit
  def visible:Boolean
  def visible_=(v:Boolean):Unit
  def length:Option[Time]
  def size:Option[Vec]
  def snapshot:Img
  def snapshot(size:Vec):Img

  //--derived
  def show:Unit = visible = true
  def hide:Unit = visible = false
  def play(speedupScale:Double):Unit = {speedup = speedupScale; play }
  def play(t:Time,speedupScale:Double):Unit = { speedup = speedupScale; cursor = t; play}
  def toggle:Unit = if(isPlaying) pause else play
}
case class Video(file:File, box:Rect=Img.defaultBox) extends Img{
  // def play:Unit ??? what mechanism will be used here to generically interact with the media playback??
  def remote(implicit g:DrawContext):MediaRemote = g.remote(this)
  def draw(implicit g:DrawContext):Unit = g ! this
  def moveTo(thatBox:Rect):Video = Video(file,thatBox)
  def toAwt:BufferedImage = ???
}
case class Html(file:File, box:Rect=Img.defaultBox, scale:Double = 1d) extends Img{
  // def play:Unit ??? what mechanism will be used here to generically interact with the media playback??
  def draw(implicit g:DrawContext):Unit = g ! this
  def moveTo(thatBox:Rect):Html = copy(box = thatBox) //TODO use generic moveTo and zoomTo for box's
  def zoomTo(region:Rect, viewport:Rect):Html = { //TODO change name from view to zoomTo
    val pos = viewport.a - region.a
    val scale = viewport.width/region.width //FIXME choose width or hight to define zoom
    Html(file, Rect(pos, box.width, box.height), scale)
  }
  def toAwt:BufferedImage = ???
}

object Text{
   def apply(value:String):Text = Text(value, Vec.zero)
}
case class Text(value:String,pos:Vec) extends Shape.Closed{
   //TODO deprecate
   // def string = value
   def -(v:Vec) = Text(value,pos-v)
   def +(v:Vec) = Text(value,pos+v)
   def +(more:String) = Text(value+more,pos)
   def draw(implicit g:DrawContext):Unit = g ! this
   def sizeWith(font:Style.Font)(implicit g:DrawContext):Vec = g.textSize(this,font)
}
object Tri{
   val root3 = math.sqrt(3)
   def iso(c:Vec,r:Double):Tri = {
      val x = Vec(root3*r/2,0) //offset
      val b = Vec(c.x,c.y+r/2)   //down to base
      //Tri(c + r.y.rz(60.deg), center - r.y, center + r.y.rz(-60.deg)
      Tri(c - r.y, b+x, b-x)
   }
   def apply(c:Vec,r:Double):Tri = Tri(c,r*2,r*2)
   def apply(c:Vec,width:Double,height:Double):Tri = {
      val y = (height/2).y
      val x = (width/2).x
      val b = c+y
      Tri(b-x, c-y, b+x)
   }
   def apply(vs:Iterable[Vec]):Tri = {
     val vs3 = vs.take(4).toList  //over-take by one to see if the traversable is only three elements
     require(vs3.size == 3, "A triangle has only have three vertices") //this is messy since we don't know if it will fail until run time
     val List(a,b,c) = vs3
     Tri(a,b,c)
   }
}
case class Tri(a:Vec,b:Vec,c:Vec) extends Shape.Closed{
   def draw(implicit g:DrawContext):Unit = g ! this
   //TODO make it a 3d triangle
   /**circumcircle construction
    * https://en.wikipedia.org/wiki/Circumscribed_circle#Triangle
    * note: the points cannot be co-linear
    *  val r = x*x + y*y  //radius squared relative to the origin
    */
   private def circCenterRelativeToA:Vec = {
      val _b = b - a
      val _c = c - a
      val b2 = _b.normSq
      val c2 = _c.normSq
      val d = 2*(_b.x*_c.y - _b.y*_c.x)
      require(d != 0) //TODO return a line ??
      Vec(_c.y*b2 - _b.y*c2, _b.x*c2 - _c.x*b2)/d
   }
   def circCenter = circCenterRelativeToA + a
   def circ = {
      val c = circCenterRelativeToA //radius is relative to the origin
      Circ(c+a,c.norm)
   }
   /**an optimized faster way of using circ.contains(t) without the sqrt*/
   def circContains(t:Vec):Boolean = {
      val c = circCenterRelativeToA
      (t-a-c).normSq < c.normSq
   }

   private def area2:Double = a.x*(b.y-c.y) + b.x*(c.y-a.y) + c.x*(a.y-b.y)  //TODO test
   def area:Double = area2/2  //TODO test
   def flip = Tri(a,c,b)
   def colinear:Boolean = area == 0
   /**a fast barycentric test with short circuiting*/
   def contains(t:Vec):Boolean = {
      val a2 = area2
      val s = (a.y*c.x - a.x*c.y + (c.y-a.y)*t.x + (a.x-c.x)*t.y)/a2
      s > 0 && s < 1 && {
         val g = (a.x*b.y - a.y*b.x + (a.y-b.y)*t.x + (b.x-a.x)*t.y)/a2
         g > 0 && (1-s-g) > 0
      }
   }
   lazy val toPoly = Poly(List(a, b, c))
   /**point closest including segments*/
   def closest(p:Vec):Vec = toPoly closest p

   def -(r:Vec) = Tri(a-r, b-r, c-r)
   def +(r:Vec) = Tri(a+r, b+r, c+r)
}
object Circ {
   /**origin zero case, with radius*/
   def apply(r:Double):Circ = Circ(Vec.zero, r)
   /**circumcircle construction
    * https://en.wikipedia.org/wiki/Circumscribed_circle#Triangle
    * note: the points cannot be co-linear
    */
   def apply(a:Vec,b:Vec,c:Vec):Circ = Tri(a,b,c).circ
   def apply(tri:Tri):Circ = tri.circ
}
abstract class Containable{ //TODO is this the right scope and names for these kind of things??
     //----Containable attributes
     /**rectangular bounding box*/
     def rect:Rect
     def contains(t:Vec):Boolean
     def area:Double

     //----generic
     /**aspect ratio adjusted version of take*/
     def take(N:Int):Iterable[Vec] = {
       val r = rect
       val extraN = (N*r.area/area).round.toInt
       r take extraN filter contains
     }
     /**take mxn grid points*/
     def take(m:Int,n:Int):Iterable[Vec] = rect.take(m,n) filter contains

     /**sample with area density*/
     def sample(dA:Double):Iterable[Vec] = rect sample dA filter contains
}
case class Circ(c:Vec,r:Double) extends Containable with Shape.Closed{
   def draw(implicit g:DrawContext):Unit = g ! this
   lazy val area:Double = pi*r*r
   lazy val perimeter:Double = 2*pi*r //circumference, length of the path around the circle

   def -(v:Vec) = Circ(c-v,r)
   def +(v:Vec) = Circ(c+v,r)
   def *(s:Double) = Circ(c,r*s)
   def contains(t:Vec):Boolean = (c dist t) <= r

   /**point closest on the circle to a point p*/
   def closest(p:Vec):Vec = c + (p - c).mag(r)

   //TODO line intersects at??

   def rect = Rect(c,r)
}
object Ellipse{
   def apply(r:Vec):Ellipse = Ellipse(Vec.zero, r)
   def apply(center:Vec,width:Double,height:Double):Ellipse = Ellipse(center, Vec(width/2, height/2) )
}
case class Ellipse(c:Vec,r:Vec, rotation:Angle=Angle(0)) extends Shape.Closed{
   def draw(implicit g:DrawContext):Unit = g ! this
   lazy val area:Double = (pi*r.x*r.y).abs
   lazy val perimeter:Double = pi*(3*(r.x+r.y) - math.sqrt((3*r.x+r.y)*(r.x+3*r.y))) //Ramanujan approximation

   lazy val w = r.x.abs*2 //width (without rotation)
   lazy val h = r.y.abs*2 //height

   def major = math.max(r.x.abs,r.y.abs)
   def minor = math.min(r.x.abs,r.y.abs)

   def closest(p:Vec):Vec = {
      val pr = tr(p)
      val s = r.x*r.y/(r.xx*pr.yy + r.yy*pr.xx).sqrt   //TODO math-overflow solution...add reference
      c + (p-c)*s  //not sure why this transforms back correctly
   }
   private def tr(p:Vec) = (p-c) rz -rotation
   def relativeDist(p:Vec):Double = relativeDist2(p).sqrt //a faster version of the more clear: (c dist p) / (c dist closest(p))
   def relativeDist2(p:Vec):Double = {
      val d = tr(p)
      d.xx/r.xx + d.yy/r.yy
   }

   def moveTo(v:Vec) = Ellipse(v , r   , rotation)
   def -(v:Vec)    = Ellipse(c-v , r   , rotation)
   def +(v:Vec)    = Ellipse(c+v , r   , rotation)
   def *(s:Double) = Ellipse(c   , r*s , rotation)
   def /(s:Double) = Ellipse(c   , r/s , rotation)
   def rz(a:Angle) = Ellipse(c   , r   , rotation+a)
   def contains(p:Vec):Boolean = relativeDist2(p) < 1

   def sampleNormal(implicit rand:Rand):Vec = Vec( rand.normal(0,r.x), rand.normal(0,r.y)).rz(rotation) + c

   def pdf(t:Vec):Double = {
     val p = tr(t) //relative version of the lookup point
     val A = 1.0/(tau*r.x*r.y)
     A*math.exp(-(p.x.sq/r.x.sq + p.y.sq/r.y.sq)/2)
   }

   /**represent the ellipse as a path (bezier curves) useful for drawing contexts that don't have circ or arc
    * 2bez: crude approx: https://stackoverflow.com/a/22035723/622016
    * 4bez: https://stackoverflow.com/a/2173084/622016
    * book: interactive bezier book: https://pomax.github.io/bezierinfo/#circles_cubic
    */
   def toPath:Path = {
      val k = 4d/3 //kappa for two beziers
      // val k = 0.551784 //2.sqrt.prev*4/3 //kappa for 4 beziers
      //--shifts
      val d1 = Vec(rotation.sin, rotation.cos)*r.y //TODO change to use a vec rotation
      val d2 = Vec(rotation.cos, rotation.sin)*(r.x*k) //TODO change to use a vec rotation
      //--top
      val tc = Vec(c.x-d1.x, c.y+d1.y) //center
      val tr = c + d2 //right
      val tl = c - d2 //right
      //--bottom
      val bc = Vec(c.x+d1.x, c.y-d1.y) //center
      val br = c + d2 //right
      val bl = c - d2 //right

      Path(List(bc, BezierVertex(br, tr, tc), BezierVertex(tl, bl, bc)))
   }
}

object Star{
   def apply(c:Vec,r:Double,n:Int):Star = Star(c, Bound(r/2,r), n)
   def apply(c:Vec,r:Double):Star = Star(c, r, 5)
}
case class Star(c:Vec, r:Bound[Double], n:Int=5, innerAngle:Angle=Angle(0),outerAngleOffset:Angle=Angle(0)) extends Shape.Closed{
   require(n > 0)
   lazy val vertices = {//FIXME TODO verify this one is working...
      val a = Angle.full/n
      val a2 = a/2
      val initAngle = innerAngle+Angle.quarter
      val finalAngle = Angle.full+innerAngle+Angle.quarter
      @tailrec def vertex(vs:List[Vec],ai:Angle):List[Vec] =
         if(ai < initAngle) vs else vertex(c + Vec.polar(r.max,ai-a2-outerAngleOffset) :: c + Vec.polar(r.min,ai) :: vs, ai-a)
      vertex(Nil, finalAngle)
   }
   def draw(implicit g:DrawContext):Unit = g ! this
}


object Arc{
   def apply(angle:Bound[Angle]):Arc = Arc(Vec.zero, Bound(0.0,1.0), angle)
   def apply(radius:Bound[Double],angle:Bound[Angle]):Arc = Arc(Vec.zero, radius, angle)
   def apply(center:Vec,radius:Bound[Double],angle:Bound[Angle]):Arc = {
      val a = Vec.polar(radius.min, angle.min)
      val b = Vec.polar(radius.max, angle.max)
      Arc(center,a,b)
   }
}
case class Arc(c:Vec, a:Vec, b:Vec){
   def draw(implicit g:DrawContext):Unit = g ! this
   lazy val r = (a.norm + b.norm)/2
   lazy val w = (a.norm - b.norm).abs
   lazy val (inner,outer) = if(a.norm < b.norm) (a,b) else (b,a)
   lazy val innerRadius = inner.norm
   lazy val outerRadius = outer.norm
   //lazy val startAngle = a.heading
   //lazy val endAngle = startAngle + (b.heading rzDelta a.heading)  //this keeps the end angle always larger than the start angle
   //lazy val angleBound = Bound(startAngle, endAngle)
   lazy val angle = Bound(a.heading.norm, b.heading.norm)
   def rz(theta:Angle):Arc = Arc(c, a rz theta, b rz theta)
   def +(t:Vec) = Arc(c+t, a, b)

   def contains(theta:Angle):Boolean = angle contains theta.norm //contains checks both directions for containment
   def contains(t:Vec):Boolean = {
      val d = c dist t
      val h = (t-c).heading
      innerRadius <= d && d <= outerRadius  && contains(h)
   }
}

/* trivial
class Axis[A,B](domain:Bound[A], range:Segment)(implicit ta:Tickable[A]) extends Scale[A,Vec]{
   val scale = Scale(domain, range)
   def apply(x:A):Vec = scale(x)
   def inv(rv:Vec):A = scale.inv(rv)
   def ticks:Iterable[Text] = domain.ticksOn(range)
}
*/

object Rect{
   /** make a rect around a center point with a given radius (circle like constructor)
    * TODO the center based constructor trips me up sometimes after some time when I am expecting a corner based constructor
    * for the corner based constructor try using a line with the Vect offset constructor
    * */
   def apply(center:Vec,r:Double):Rect = Rect(center, 2*r,2*r)
   /** make a rect around a center point with a given width, height (radius like constructor)*/
   def apply(center:Vec,width:Double,height:Double):Rect = Rect(Vec(center.x-width/2, center.y-height/2), Vec(center.x+width/2, center.y+height/2))
   def radius(center:Vec,r:Vec):Rect = Rect(center, r.x*2, r.y*2)
   //TODO should this really be created at cm???
   def apply(width:Double,height:Double):Rect = Rect(Vec.zero, width, height)
   def apply(b:Vec):Rect = Rect(Vec.zero, b)
   def apply(line:Line):Rect = line.rect //TODO make this the case class constructor so the a,b apply based a,b constructor is center and r

   private val goldenRatio = 5.sqrt.next/2 // (1 + sqrt(5)) / 2
   def golden(width:Double):Rect = Rect(Vec.zero, Vec(1,goldenRatio.inv)*width)
}

case class Rect(a:Vec,b:Vec) extends Containable with Shape.Closed{
   def toAwt:java.awt.Rectangle = new java.awt.Rectangle(a.x.toInt, a.y.toInt,  size.x.toInt, size.y.toInt)
   def draw(implicit g:DrawContext):Unit = g ! this
   lazy val width = b.x - a.x //TODO maybe change these to width: dx, and height: dy
   lazy val height = b.y - a.y
   lazy val size = Vec(width, height)
   /**aspect ration w:h => width/height*/
   lazy val aspect = width/height
   lazy val cm:Vec = Vec((a.x+b.x)/2, (a.y+b.y)/2)
   lazy val area:Double = (b.x-a.x)*(b.y-a.y)
   lazy val diag:Vec = b - a

   def axes[A,B](xDomain:Bound[A], yDomain:Bound[B])(implicit ta:Tickable[A],tb:Tickable[B]) = new Axes(xDomain,yDomain,this)

   //--corners
   def nw = a
   def ne = Vec(b.x,a.y)
   def se = b
   def sw = Vec(a.x,b.y)
   def center = cm


   //--get sides
   def top = Line(nw,ne)
   def bottom = Line(sw,se)
   def right = Line(ne,se)
   def left = Line(nw,sw)

   def side(prop:Style.Property):Line = prop match {
     case Style.Top => top
     case Style.Bottom => bottom
     case Style.Left => left
     case Style.Right=> right
     case _ => ??? //FIXME use a proper property trait space the spans the Properties.Side cases so this case does not happen
   }

   //--get other type of lines
   def hLine(t:Double):Line = top  + (t*height).y
   def vLine(t:Double):Line = left + (t*width).x

   def midLine:Line    = hLine(0.5) //cm +- width.x
   def centerLine:Line = vLine(0.5) //cm +- height.y

   //-- move{Side}To
   def topTo(y:Double):Rect    = Rect( Vec(a.x,       y)   , b)
   def bottomTo(y:Double):Rect = Rect( a                   , Vec(b.x    , y) )
   def leftTo(x:Double):Rect   = Rect( Vec(x,       a.y)   , b)
   def rightTo(x:Double):Rect  = Rect( a                   , Vec(x      , b.y) )
   //move{Side}By
   def topBy(dy:Double):Rect    = Rect( Vec(a.x,    a.y+dy), b)
   def bottomBy(dy:Double):Rect = Rect( a                  , Vec(b.x    , b.y+dy) )
   def leftBy(dx:Double):Rect   = Rect( Vec(a.x+dx, a.y)   , b)
   def rightBy(dx:Double):Rect  = Rect( a                  , Vec(b.x+dx , b.y) )

   //--combinator  //TODO think about the relevance of this?  Or is a layout manager beter?
   // def below_:(rect:Rect)//TODO allow all shape types instead of just rect

   /**expand/shrink a shape by a level. +offset: makes a larger shape -offset: makes a smaller rect
    * keywords: pad, margin, padding, border, expand, reduce*/
   def pad(v:Double):Rect   = {val vm = Vec(v,v); Rect( a-vm, b+vm)} //TODO check for under padded size to return a zero
   import Style._
   def pad(v:Double, align:AlignHorizontal):Rect = align match {
     case Left   => Rect( a-v.x, b)
     case Right  => Rect( a, b+v.x)
     case Center => Rect( a-v.x, b+v.x)
   }
   def pad(v:Double, align:AlignVertical):Rect = align match {
     case Top    => Rect( a-v.y, b)
     case Bottom => Rect( a, b+v.y)
     case Midline=> Rect( a-v.y, b+v.y)
   }

   2017-07-30("use top","v0.2.15")
   def n = Line(nw,ne)
   2017-07-30("use bottom","v0.2.15")
   def s = Line(sw,se)
   2017-07-30("use right","v0.2.15")
   def e = Line(ne,se)
   2017-07-30("use left","v0.2.15")
   def w = Line(nw,sw)

   lazy val vertices = List(nw, ne, se, sw)

   def rect = this
   /**take points with m-rows and n-columns*/
   override def take(m:Int,n:Int):Iterable[Vec] = for(x <- a.x ~ b.x take n; y <- a.y ~ b.y take m) yield Vec(x,y)
   /**take N points from the rectangle area where N = n*m and n = m*aspect */
   override def take(N:Int):Iterable[Vec] = {
      val m = math.sqrt(N/aspect) //number of columns
      val n = aspect*m //number of rows
      take(m.toInt,n.toInt)
   }
   /**sample a dA (sample area size) */
   override def sample(dA:Double):Iterable[Vec] = this take (area/dA).round.toInt

   private def scaleToAt(that:Rect,newCenter:Vec):Rect =
      //use the height as the limit
      if(that.aspect > this.aspect)   Rect(newCenter,              this.width*that.height/this.height, that.height)
      //use the width as the limit
      else                            Rect(newCenter,  that.width, this.height*that.width/this.width)
   def scaleTo(that:Rect):Rect = scaleToAt(that,this.cm)
   def fitIn(that:Rect):Rect = scaleToAt(that, that.cm)
   def moveTo(pos:Vec):Rect = Rect(this.cm+pos, this.width, this.height)

   def round = Rect(a.round, b.round)
   def floor = Rect(a.floor, b.floor)
   def ceil  = Rect(a.ceil,  b.ceil)

   def -(r:Vec) = Rect(a-r, b-r)
   def +(r:Vec) = Rect(a+r, b+r)
   def *(s:Double) = Rect(cm,width*s,height*s)
   def /(s:Double) = Rect(cm,width/s,height/s)
   def contains(t:Vec):Boolean = !(t.x < a.x || t.x > b.x || t.y < a.y  || t.y > b.y)
   def poly = Poly(Seq(a, Vec(a.x,b.y), b, Vec(b.x,a.y)))

   lazy val bound  = Bound(a, b)
   lazy val xBound = Bound(a.x, b.x)
   lazy val yBound = Bound(a.y, b.y)
   def mobius(v:Vec):Vec = Vec( xBound mobius v.x, yBound mobius v.y)

   lazy val toPoly = Poly(vertices)
   /**point closest including segments*/
   def closest(p:Vec):Vec = toPoly closest p

   //--- split features

   /**landscape => vsplit, else portrait => hsplit*/
   def split(n:Int):Vector[Rect] = splitToRect(n,Rect.goldenRatio)

   def hsplit(n:Int):Vector[Rect] = {
     val dy = (b.y-a.y)/n
     Vector.tabulate(n){i => Rect(Vec(a.x, a.y + dy*i), Vec(b.x, a.y + dy*(i+1) ))}
   }
   def vsplit(n:Int):Vector[Rect] =  {
     val dx = (b.x-a.x)/n
     Vector.tabulate(n){i => Rect(Vec(a.x + dx*i, a.y), Vec(a.x + dx*(i+1), b.y))}
   }

   def hsplit():Vector[Rect] = hsplit(2)
   def vsplit():Vector[Rect] = vsplit(2)
   def split():Vector[Rect] = split(2)

   /**hsplit at a ratio*/ //TODO add an option if the ratio is greater than 1 to return toInt rects
   def hsplit(p:Double):Vector[Rect] = {
     val y = a.y + (b.y-a.y)*p  //where to draw the split
     Vector(
       Rect(a, Vec(b.x, y)    ),
       Rect(   Vec(a.x, y),  b)
     )
   }
   /**vsplit at a ratio*/
   def vsplit(p:Double):Vector[Rect] = {
     val x = a.x + (b.x-a.x)*p //where to draw the split
     Vector(
       Rect(a, Vec(x, b.y)    ),
       Rect(   Vec(x, b.y),  b)
     )
   }

   /**auto split and zip plus round since most often used for size by size pixel viz (note: requires a .size call)*/
   def vsplit[A](vs:Iterable[A]):Iterable[(Rect,A)] = vsplit(vs.size) zip vs
   def hsplit[A](vs:Iterable[A]):Iterable[(Rect,A)] = hsplit(vs.size) zip vs
   def split[A](vs:Iterable[A]):Iterable[(Rect,A)] = split(vs.size) zip vs

   /**split m rows an n cols*/
   def split(m:Int,n:Int):Vector[Rect] = this hsplit m flatMap {_ vsplit n}
   def splitToSquares(N:Int):Vector[Rect] = splitToRect(N,1d)
   def split(N:Int,aspectRatioGoal:Ratio):Vector[Rect] = splitToRect(N, aspectRatioGoal.toDouble)
   //TODO 2017-07-30("use split(N,Ratio(1,2)) instead")
   def splitToRect(N:Int,aspectRatio:Double = Rect.goldenRatio):Vector[Rect] = {
      if(      N <= 0) Vector()
      else if (N == 1) Vector(this)
      else {
        val s = (area/N).sqrt // s*s = w*h
        val h = s/aspectRatio  //target height per rect
        val m = (height/h).toInt  //rows

        // val w = s*aspectRatio //optimal width and height
        // val n = (width/w).round.toInt
        // println(s"s:$s h:$h N:$N m:$m    w:$w n:$n")

        if     (m >= N) hsplit(N) //desired aspect does fill column then force an single column hsplit
        else if(m <= 1) vsplit(N) //ensure at least 1 row and a fast soln
        else {
          val n = N / m  //cols
          val e = N % m  //extra
          val e_m = e / m  //extra to add per row
          val e_e = e % m  //extra overflow up to row
          // assert(   N == m*(n + e_m) + e_e)  //???
          hsplit(m).zipWithIndex.flatMap{ case (r,i) =>
            r.vsplit(
              n + e_m + (if(i < e_e) 1 else 0)
            )
          }
        }
      }
   }

}
object Poly{
   // def apply(vs:Vec*):Poly = new Poly(vs)
   // def apply(vs:Seq[Vec]):Poly = new Poly(vs)

   def points(vertices:(Double,Double)*):Poly = new Poly(vertices map {p => Vec(p)} )
   /**andrew monotone chain  algorithm for convex hull generation
    * reference: http://en.wikibooks.org/wiki/Algorithm_implmentation/Geometry/Convex_hull/Monotone_chain#Haskell */
   def convex(points:Seq[Vec]):Poly = {

      if(points.size < 3) new Poly(Seq()) else {
         def clockwise(o:Vec,a:Vec,b:Vec):Boolean = ((a - o) x2d (b - o)) <= 0
         def chain(acc:List[Vec],ps:List[Vec]):List[Vec] = (acc,ps) match {
            case (r1::r2::rs, x::xs) => if( clockwise(r2,r1,x)) chain(r2::rs, x::xs)  //clockwise:     remove recent chain
                                               else          chain(x::acc, xs)     //anticlockwise: append to chain
            case (acc, x::xs)        => chain(x::acc, xs)  //only one point: add next vÑ–sited
            case (acc, Nil)          => acc.tail //no more points: all done (to make it clockwise add a reverse and use lower++upper instead of upper++lower)
         }

         val sorted = points.sortBy(p => p.x -> p.y).toList
         val lower =  chain(Nil, sorted)
         val upper =  chain(Nil, sorted.reverse)
         new Poly(upper ++ lower)
      }
   }
}

trait Segment extends Lerp[Vec]{
   def lerp(t:Double):Vec
   def grad(t:Double):Vec

   val a:Vec = lerp(0d)
   val b:Vec = lerp(1d)

   // def tangent(t:Double):Vec
   def arrowHead(size:Double):Tri = {
     val ortho = grad(1.0).mag(size/2)
     val b2 = b - ortho.rzLeft*2
     Tri(b2-ortho, b, b2+ortho)
   }

   def dist:Double
   def last:Vec

   def step(stepDist:Double) = take( tickCount(stepDist) )

   def tickCount(stepDist:Double):Int = {
     val n0 = (dist/stepDist).round.toInt  //auto detect ~ tick count from given spacing
     if(n0 < 2) 2 else n0
   }
}
/**Vertex append to points to make segments(paths in a more drawing context language)
 *   Vertex        context-api
 *   ------------  ---------
 *   PathStart     moveTo
 *   Vec           lineTo
 *   BezierVertex  bezierCurveTo
 *   PathClose     close (closed area type closed))
 *   PathEnd       end (line type)
 * */
trait Vertex{
   //note that a Vec extends the Vertex class
   def last:Vec
   def segmentFrom(a:Vec):Segment
}
case class BezierVertex(ca:Vec,cb:Vec,b:Vec) extends Vertex{
  def last = b
  def segmentFrom(a:Vec) = Bezier(a, ca, cb, b)
}

object Path{
  /**Visvalingam & Whyatt triangular area algorithm to simplify and preserve
   * This makes use of a mutable minHeap-priority queue and hashMap for updating areas
   *https://bost.ocks.org/mike/simplify/
   *https://hydra.hull.ac.uk/resources/hull:8338
   */
  private def simplifyIndices(vs:Array[Vec], maxPoints:Int, minArea:Double):IndexedSeq[Int] = {
    //Note: the mutable mess but the logic is contained inside this function
    val indices = 0 until vs.size

    //--mutable stores
    val removed = scala.collection.mutable.Set.empty[Int]
    val area = scala.collection.mutable.Map.empty[Int,Double]

    //--move around the points
    @tailrec def next(i:Int):Int = if(removed(i+1)) next(i+1) else i+1
    @tailrec def last(i:Int):Int = if(removed(i-1)) last(i-1) else i-1

    //--Point calculations
    def isEndPoint(i:Int):Boolean = i == 0 || i == indices.last

    //--Priority queue data and features load up the priority queue (minHeap)
    object MinArea extends Ordering[Int]{
      def compare(a:Int, b:Int) = area(b) compare area(a)  //here area is a closure
    }
    val minHeap = scala.collection.mutable.PriorityQueue.empty(MinArea)

    //--update the new area value and re-add it to the priority queue for sorting
    def update(i:Int):Unit = if(!removed(i) && !isEndPoint(i) ){
      area(i) = Tri(vs(last(i)), vs(i), vs(next(i))).area.abs
      minHeap += i
      () //make sure it returns a unit function
    }

    //--add to the removed set and recompute areas for the new left and right points
    //TODO be careful if the next calculated area is even less than the current triangle
    def remove(i:Int):Unit = {removed += i; update(last(i)); update(next(i)) }

    //---initialize data structures
    indices foreach update //calculate the area and add each point to the minHeap priority queue

    //---do a simplification step
    def step:Unit = remove( minHeap.dequeue() )
    // Log(size, removed, minHeap)

    implicit val logLevel = LogLevel.Info
    //--remove points until meeting a condition //TODO add some smarter stopping condition
    def size:Int = vs.size - removed.size
    if(size > 2)
      while(size > maxPoints || area(minHeap.head) < minArea ) step //short circuit on easier

    //--return the remaining points as a Path with the (z-axis set to the area)
    // def vecWithArea(i:Int):Vec = vs(i) mz area.getOrElse(i,Double.MaxValue)
    // indices filterNot removed map vecWithArea

    //--return the remaining indices
    indices filterNot removed
  }
  def simplify(vs:Array[Vec], maxPoints:Int, minArea:Double=0.0):IndexedSeq[Vec] = {
    simplifyIndices(vs, maxPoints, minArea) map vs
  }
  def simplifyBy[A](xs:IndexedSeq[A], maxPoints:Int, minArea:Double=0.0)(f: A=>Vec):IndexedSeq[A] = {
    //Note: the mutable mess but the logic is contained inside this function
    //--make an array of Vec (the index represents the point
    val vs = (xs map f).toArray  //here vs are now Vec, a concrete type of Array[Vec]
    //--return the remaining indices of the src collection
    simplifyIndices(vs, maxPoints, minArea) map xs
  }
  def apply(vs:Vertex*):Path = Path(vs)
}
/**Path collection of vertex drawing points*/
case class Path(vertices:Iterable[Vertex], isClosed:Boolean=false) extends Segment with Shape.Open{
  def draw(implicit g:DrawContext):Unit = g ! this
  private lazy val segmentIndex:Array[(Double,Segment)] = {
    val vs = vertices.toArray
    val segments = vs zip vs.tail map {case (a,b) => b segmentFrom a.last}
    val si = segments.tail.scanLeft((segments.head.dist,segments.head)){case ((d,a),b) => (d + b.dist, b)}
    if(isClosed){
      val close = Line(vs.last.last, vs.head.last)
      si :+ ( (close.dist, close) )
    }
    else si
  }
  lazy val dist = segmentIndex.last._1
  lazy val last = vertices.last.last //TODO speed up find of last
  lazy val head = vertices.head.last //TODO speed up find of head since this is across an iterable
  lazy val size = vertices.size
  private def findSubSegmentLerp(td:Double):(Segment,Double) = {
      if      (td <= 0) (segmentIndex.head._2,0d)
      else if (td >= 1) (segmentIndex.last._2,1d)
      else {
        val d = td*dist
        val i = segmentIndex.searchBy(d)(_._1)
        val (db,sb) = segmentIndex(i)
        val t = if(sb.dist > 0) (d - (db - sb.dist))/sb.dist else 0
        (sb,t)
      }
  }
  def lerp(td:Double):Vec = findSubSegmentLerp(td) match {case (segment,t) => segment.lerp(t)}
  def grad(td:Double):Vec = findSubSegmentLerp(td) match {case (segment,t) => segment.grad(t)}

  def simplify(maxPoints:Int, minArea:Double=0.0):Path = {
    val vs:Array[Vec] = vertices.map{_.last}.toArray
    val vsSimple = Path.simplify(vs, maxPoints, minArea)
    Path(vsSimple)
  }
}

case class Clip(path:Path)


case class Poly(val vertices:Seq[Vec]) extends Shape.Closed{
   def draw(implicit g:DrawContext):Unit = g ! this
   val valid = vertices.take(3).size > 2
   lazy val segments:Seq[Line] = vertices.zip(vertices.tail :+ vertices.head).map{case (a,b) => Line(a,b)}
   /*
   lazy val cm:Vec = Vec(vertices.map(_.x).sum/vertices.size,
                    vertices.map(_.y).sum/vertices.size)
  */

   lazy val convex:Poly = Poly.convex(vertices)

   //http://en.wikipedia.org/wiki/Centroid  for ccw computed area
   lazy val area:Double = segments.map{case Line(a,b) => b.x*a.y - a.x*b.y}.sum/2
   lazy val cm:Vec = segments.foldLeft(Vec(0,0)){
     case (acc, Line(a,b)) =>
         val m = b.x*a.y - a.x*b.y
         acc + Vec((a.x+b.x)*m, (a.y+b.y)*m)
     case (acc, _) => acc
   }/6/area

   def +(r:Vec) = Poly(vertices.map(_ + r))
   def *(s:Double) = Poly(vertices.map(v => cm + (v - cm)*s))

   2017-07-30("use rect for a more direct meaning since bound has other meanings in drx","v0.2.13")
   lazy val bounds:Rect = rect
   lazy val rect:Rect = vertices.tail.foldLeft(Rect(vertices.head, vertices.head)){case (Rect(min,max),v) => Rect((min point v){_ min _}, (max point v){_ max _}) }

   def contains(t:Vec):Boolean = {
      (rect contains t) &&
      {
         var c = false
         segments.foreach{case Line(a,b) =>
            val ytest = (a.y > t.y) != (b.y > t.y)
            if( ytest && ( t.x    <    (b.x-a.x) * (t.y-a.y) / (b.y-a.y) + a.x)) c = !c
         }
         c
      }
   }

   def closest(p:Vec):Vec = segments.map{_ closest p}.sortBy(_ dist p).head
   def dist(p:Vec):Double = closest(p) dist p


   //Weiler-Atherton clipping algorithm (returns separate polygons
   //innefficient but simple implementation
   // 1) could find intersections once instead of twice for both lists
   // 2) could use linked lists rather than vector lookups

   2017-07-30("use clip instead","0.1")
   def ^(that:Poly):Seq[Poly] = clip(that)
   def clip(that:Poly):Seq[Poly] = weiler_atherton(that, false)._1
   def merge(that:Poly):Seq[Poly] = weiler_atherton(that, true)._1


   //---weiler atherton clipping and merging algorithm, exports additional information for better vizualization of the algorithm
   sealed trait Kind{def color:Color}
   case object Enter extends Kind{def color = Green}
   case object Exit extends Kind{def color = Red}
   case object Clip extends Kind{def color = Pink}
   case object Subj extends Kind{def color = Blue}

   //--Vertex
   case class V(v:Vec,k:Kind){ def color:Color = k.color }

   def weiler_atherton(that:Poly, merge:Boolean):(Seq[Poly],Seq[V],Seq[V]) = {
      case class I(v:Vec,a:Line,b:Line,k:Kind){
         def toV:V = V(v,k)
      }

      def kind(a:Line,b:Line):Kind = if((a x2d b) < 0) Enter else Exit

      //--segement lists
      val as = this.segments
      val bs = that.segments

      //--intersections
      val is = as.flatMap{a =>  bs.flatMap{b => (a ^ b) map {i => I(i, a,b, kind(a,b))} } }

      if(is.isEmpty){
         //--clip a entirely encompases subject b; return the subject polygon
         if(this contains that.vertices.head) (Seq(that),Seq(),Seq())
         //--subect b is entirely outside of clip a; return no polygons
         else (Seq(),Seq(),Seq())
      //--compute intersection polys; return a sequence of the polygons
      }else {
         //--new Vertex lists
         val A = as.flatMap{a => V(a.a,Clip) +: is.filter(_.a == a).sortBy{i => a.a dist i.v}.map{_.toV}}
         val B = bs.flatMap{b => V(b.a,Subj) +: is.filter(_.b == b).sortBy{i => b.a dist i.v}.map{_.toV}}

         //val startKind = if(merge) Exit else Enter
         val starts = A.filter(_.k == Enter)

         def nextV(vs:Seq[V],v:V):V = {
            val i = vs.indexOf(v)
            if(i == -1) println(s"Error could not find $v")
            val iNext = if(i==vs.size-1) 0 else i+1
            vs(iNext)
         }
         def collect(vs:Seq[V]):Seq[V] = {
            val head = vs.head
            val last = vs.last
            val next = last.k match {
               case Enter => nextV(if(merge) A else B,last)
               case Exit =>  nextV(if(merge) B else A,last)
               case Clip =>  nextV(A,last)
               case Subj =>  nextV(B,last)
            }
            //println(f"${next.k}%5s$i%2d")
            if(vs.size > (this.vertices.size + that.vertices.size)*2) {println(s"Error endless loop: $last to $next");vs}
            else if(head == next) vs
            else collect(vs :+ next)
         }

         val polys = starts.foldLeft(Seq[Poly]()){case (ps, start) =>
            if(ps.find(_.vertices contains start.v).isDefined) ps
            else ps :+ Poly(collect(Seq(start)).map(_.v))
         }
         (polys,A,B)
      }
   }

   //Radke quick overlap detection algorithm (faster unordered polygon) with short circuit
   def overlaps(that:Poly):Boolean = {
      this.vertices.exists(that contains _) ||
      that.vertices.exists(this contains _) ||
      this.segments.exists{lineA => that.segments.exists{lineB => (lineA ^ lineB).nonEmpty}}
   }

   def ~(v:Vec):Poly = Poly(vertices :+ v)
   def ~(that:Poly):Polys = Polys(Seq(this,that))

}
case class Polys(polys:Seq[Poly]=Seq()) extends Shape.Closed{
   def draw(implicit g:DrawContext):Unit = g ! this
   def findContaining(t:Vec):Option[Poly] = polys find {_ contains t}
   def contains(t:Vec):Boolean = findContaining(t).isDefined
   def ~(that:Poly):Polys = Polys(polys :+ that)
   def ~(that:Polys):Polys = Polys(polys ++ that.polys)
   def closest(t:Vec):Vec = polys.map{_ closest t}.sortBy{_ dist t}.head

   def foreach(func:Poly => Unit):Unit = polys.foreach(func)

   def size = polys.size
   lazy val area:Double = polys.map(_.area).sum
}