namespace eval bezierCrossing {
  namespace import ::tcl::mathop::* ::tcl::mathfunc::* 
  variable nearby 0.01
  namespace export bezXbez
}

proc ::bezierCrossing::bez1stHalf {x0 y0 x1 y1 x2 y2 x3 y3} {
  # return first half of bezier
  set x01 [/ [+ $x0 $x1] 2.0]
  set x12 [/ [+ $x1 $x2] 2.0]
  set x23 [/ [+ $x2 $x3] 2.0]
  set x012 [/ [+ $x01 $x12] 2.0]
  set x123 [/ [+ $x12 $x23] 2.0]
  set x0123 [/ [+ $x012 $x123] 2.0]
  #
  set y01 [/ [+ $y0 $y1] 2.0]
  set y12 [/ [+ $y1 $y2] 2.0]
  set y23 [/ [+ $y2 $y3] 2.0]
  set y012 [/ [+ $y01 $y12] 2.0]
  set y123 [/ [+ $y12 $y23] 2.0]
  set y0123 [/ [+ $y012 $y123] 2.0]
  #
  list $x0 $y0 $x01 $y01 $x012 $y012 $x0123 $y0123
}

proc ::bezierCrossing::bez2ndHalf args {
  # return second half of bezier
  coordsReverse [bez1stHalf {*}[coordsReverse $args]]
}

proc ::bezierCrossing::bezAt {bez f} {
  # return coordinates of bezier on position f where 0 <= f <= 1
  # https://de.wikipedia.org/wiki/B%C3%A9zierkurve
  lassign $bez x0 y0 x1 y1 x2 y2 x3 y3
  list\
    [expr {(-$x0 + 3*$x1 - 3*$x2 + $x3) * $f**3 +
           (3*$x0 - 6*$x1 + 3*$x2) * $f**2 +
           (-3*$x0 + 3*$x1) * $f +
           $x0}]\
    [expr {(-$y0 + 3*$y1 - 3*$y2 + $y3) * $f**3 +
           (3*$y0 - 6*$y1 + 3*$y2) * $f**2 +
           (-3*$y0 + 3*$y1) * $f +
           $y0}]
}

proc ::bezierCrossing::bezCenter {x0 y0 x1 y1 x2 y2 x3 y3} {
  # return center coords of bezier
  set x01 [/ [+ $x0 $x1] 2.0]
  set x12 [/ [+ $x1 $x2] 2.0]
  set x23 [/ [+ $x2 $x3] 2.0]
  set x012 [/ [+ $x01 $x12] 2.0]
  set x123 [/ [+ $x12 $x23] 2.0]
  #
  set y01 [/ [+ $y0 $y1] 2.0]
  set y12 [/ [+ $y1 $y2] 2.0]
  set y23 [/ [+ $y2 $y3] 2.0]
  set y012 [/ [+ $y01 $y12] 2.0]
  set y123 [/ [+ $y12 $y23] 2.0]
  #
  list [/ [+ $x012 $x123] 2.0] [/ [+ $y012 $y123] 2.0]
}

proc ::bezierCrossing::coordsReverse coords {
  # revert $coords {x0 y0 ... xn yn} -> {xn yn ... x0 y0}
  concat {*}[lmap {a b} [lreverse $coords] {list $b $a}]
}

proc ::bezierCrossing::distance {x0 y0 x1 y1} {
  # return distance of coord pairs
  hypot [- $x1 $x0] [- $y1 $y0]
}

proc ::bezierCrossing::calcXY {func args} {
  # apply function on pairwise x, y arguments, return results as list
  foreach {x y} $args {
    lappend xx $x
    lappend yy $y
  }
  list [{*}$func {*}$xx] [{*}$func {*}$yy] 
}

proc ::bezierCrossing::disiunct {b0 b1} {
  # test if beziers don't touch each other
  lassign [calcXY min {*}$b0] x0min y0min
  lassign [calcXY max {*}$b0] x0max y0max
  lassign [calcXY min {*}$b1] x1min y1min
  lassign [calcXY max {*}$b1] x1max y1max
  expr {($x0min > $x1max) || ($x1min > $x0max) ||
        ($y0min > $y1max) || ($y1min > $y0max)}
}

proc ::bezierCrossing::bezSize bez {
  # return size of bounding box of bezier
  lassign [calcXY min {*}$bez] xmin ymin
  lassign [calcXY max {*}$bez] xmax ymax
  max [- $xmax $xmin] [- $ymax $ymin]
}

proc ::bezierCrossing::bezXbezRaw {b0 b1} {
  # return list of crossing fractions {f0 g0 f1 g1 ...} of beziers
  # hurry and don't worry about nearby-doublettes
  if {[disiunct $b0 $b1]} then return
  variable nearby
  set result ""
  lassign [bezCenter {*}$b0] x0 y0
  lassign [bezCenter {*}$b1] x1 y1
  if {[distance $x0 $y0 $x1 $y1] < $nearby} then {
    lappend result 0.5 0.5
  } 
  if {[bezSize "$b0 $b1"] > 2*$nearby} then {
    foreach {f0 f1} [bezXbezRaw [bez1stHalf {*}$b0] [bez1stHalf {*}$b1]] {
      lappend result [* 0.5 $f0] [* 0.5 $f1]
    }
    foreach {f0 f1} [bezXbezRaw [bez1stHalf {*}$b0] [bez2ndHalf {*}$b1]] {
      lappend result [* 0.5 $f0] [+ 0.5 [* 0.5 $f1]]
    }
    foreach {f0 f1} [bezXbezRaw [bez2ndHalf {*}$b0] [bez1stHalf {*}$b1]] {
      lappend result [+ 0.5 [* 0.5 $f0]] [* 0.5 $f1]
    }
    foreach {f0 f1} [bezXbezRaw [bez2ndHalf {*}$b0] [bez2ndHalf {*}$b1]] {
      lappend result [+ 0.5 [* 0.5 $f0]] [+ 0.5 [* 0.5 $f1]]
    }
  }
  set result
}

proc ::bezierCrossing::bezXbez {b0 b1} {
  # list crossings of beziers without nearby-doublettes
  variable nearby
  set result {}
  set pairs [lmap {f0 f1} [bezXbezRaw $b0 $b1] {list $f0 $f1}]
  if {$pairs eq ""} then return
  set sorted [lsort -unique -real -index 0 $pairs]
  lassign $sorted testDot
  foreach pair [lrange $sorted 1 end] {
    lassign $testDot f0
    set p0 [bezAt $b0 $f0]
    lassign $pair f1
    set p1 [bezAt $b0 $f1]
    if {[distance {*}$p0 {*}$p1] > $nearby*2} then {
      lappend result {*}$testDot
      set testDot $pair
    }
  }
  lappend result {*}$testDot
}

# example:
# set blue "85.0 282.0 136.0 202.0 226.0 206.0 263.0 288.0"
# set red "104.0 221.192152 155.0 301.192152 245.0 297.192152 282.0 215.192152"
# bezierCrossing::bezXbez $blue $red
# --> 0.20660400390625 0.10235595703125 0.8758544921875 0.7569580078125