Patch title: Release 90 bulk changes
Abstract:
File: /math/curve.pli
Key:
    Removed line
    Added line
# Copyright  Hubert Tonneau  hubert.tonneau@pliant.cx
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License version 2
# as published by the Free Software Foundation.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# version 2 along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

module "/pliant/language/compiler.pli"
module "/pliant/language/context.pli"
module "/pliant/math/functions.pli"
module "/pliant/math/point.pli"
module "/pliant/math/transform.pli"

constant sqrt2 2^0.5
constant small 1/4
constant big 3/4  



public

  type CurvePoint
    field Float in_x in_y
    field Float x y
    field Float out_x out_y
    field uInt16 flags ; field uInt8 in_mode out_mode
  
  constant curvepoint_automatic 0
  constant curvepoint_directed 1
  constant curvepoint_manual 2
  
  constant curvepoint_selected 1
  constant curvepoint_withinpoint 100h
  constant curvepoint_withoutpoint 200h
  
  constant automatic curvepoint_automatic
  constant directed curvepoint_directed
  constant manual curvepoint_manual
  


module "curve/mode.pli"
constant is_ready (cast 80h Int)

type Curve
  field Array:CurvePoint pts
  field Int status <- standard
 
function 'cast Status' c -> s
  arg Curve c ; arg Status s
  explicit
  s := shunt (c:status .and. is_ready)<>0 success failure


#--------------------------------------------------------------------


function hsin x -> r
  arg Float x r
  r := sin x*pi/2
  
function hcos x -> r
  arg Float x r
  r := cos x*pi/2

function weigh x -> r
  arg Float x r
  r := 1-(cos x*pi/2)

function norme x y -> n
  arg Float x y n
  n := (x*x+y*y)^0.5

function scalarproduct x0 y0 x1 y1 -> p
  arg Float x0 y0 x1 y1 p
  p := x0*x1+y0*y1

function angle x y -> a
  arg Float x y a
  var Float xx yy
  if x>=0
    if y>=0
      xx := x ; yy := y ; a := 0
    else
      xx := -y ; yy := x ; a := 1.5*pi
  else
    if y>=0
      xx := y ; yy := -x ; a := 0.5*pi
    else
      xx := -x ; yy := -y ; a := pi
  check xx>=0 and yy>=0
  if yy<=xx
    if xx=0
      a := 0
    else
      a += atan yy/xx
  else
    a += 0.5*pi-(atan xx/yy)

function min a b c -> r
  arg Float a b c r
  r := min (min a b) c
      
function max a b c -> r
  arg Float a b c r
  r := max (max a b) c
      
function bound a m M -> r
  arg Int a m M r
  r := shunt a<=m m a>=M M a

function distance a b -> d
  arg Point2 a b ; arg Float d
  d := ( (a:x-b:x)*(a:x-b:x) + (a:y-b:y)*(a:y-b:y) )^0.5
  if d=undefined
    d := 0


#--------------------------------------------------------------------


function curve_point x y angle -> p
  arg Float x y ; arg CBool angle ; arg CurvePoint p
  p x := x ; p y := y
  if angle
    p in_x := 0 ; p in_y  := 0
    p out_x := 0 ; p out_y := 0
    p in_mode := manual ; p out_mode := manual
  else
    p in_mode := automatic ; p out_mode := automatic
  p flags := 0

method p tg_in x y
  arg_rw CurvePoint p ; arg Float x y
  p in_x := x ; p in_y := y
  p in_mode := directed

method p tg_out x y
  arg_rw CurvePoint p ; arg Float x y
  p out_x := x ; p out_y := y
  p out_mode := directed

method p in x y
  arg_rw CurvePoint p ; arg Float x y
  p in_x := x ; p in_y := y
  p in_mode := manual

method p out x y
  arg_rw CurvePoint p ; arg Float x y
  p out_x := x ; p out_y := y
  p out_mode := manual

method p is_angle -> c
  arg CurvePoint p ; arg CBool c
  c := p:in_mode<>automatic or p:out_mode<>automatic

export curve_point '. tg_in' '. tg_out' '. in' '. out' '. is_angle'


#--------------------------------------------------------------------


method c reset
  arg_rw Curve c
  c:pts size := 0
  c status := standard
  
method c through x y
  arg_rw Curve c ; arg Float x y
  c pts += curve_point x y false
  if pliant_debugging_level>=2
    c status := standard

method c angle x y
  arg_rw Curve c ; arg Float x y
  c pts += curve_point x y true
  if pliant_debugging_level>=2
    c status := standard

function '+=' c p
  arg_rw Curve c ; arg CurvePoint p
  c pts += p
  if pliant_debugging_level>=2
    c status := standard


method c compute mode
  arg_rw Curve c ; arg Int mode
  if c:pts:size<2
    return
  for (var Int k) 0 c:pts:size-1
    var Pointer:CurvePoint p0 :> c:pts (shunt k<>0 k-1 c:pts:size-1)
    var Pointer:CurvePoint p :> c:pts k
    var Pointer:CurvePoint p1 :> c:pts (shunt k<>c:pts:size-1 k+1 0)
    if p:in_mode<>manual or p:out_mode<>manual
      var Float a0 := angle p:x-p0:x p:y-p0:y ; check a0=defined
      var Float a1 := angle p1:x-p:x p1:y-p:y ; check a1=defined
      var Float a := min (abs a1-a0-2.0*pi) (abs a1-a0) (abs a1-a0+2.0*pi)
      var Float reduce := shunt a<=pi/2 1 sin:a
      if p:in_mode=automatic or p:out_mode=automatic
        var Float a := (a0+a1)/2
        var Float tx := cos a
        var Float ty := sin a
        if (scalarproduct tx ty p1:x-p0:x p1:y-p0:y)<0 { tx := -tx ; ty := -ty }
        if p:in_mode=automatic { p in_x := -tx ; p in_y := -ty }
        if p:out_mode=automatic { p out_x := tx ; p out_y := ty }
      if p:in_mode<>manual
        var Float a0 := angle p0:x-p:x p0:y-p:y
        var Float a1 := angle p:in_x p:in_y
        var Float a := min (abs a1-a0-2*pi) (abs a1-a0) (abs a1-a0+2*pi)
        var Float n := (1-cos:a)*(norme p0:x-p:x p0:y-p:y)
        var Float d := (sqrt2-1) * sin:a * 2*sin:a * (norme p:in_x p:in_y)
        var Float f := (shunt n<>0 and d<>0 n/d 0)*reduce;
        p in_x *= f ; p in_y *=f
      if p:out_mode<>manual
        var Float a0 := angle p1:x-p:x p1:y-p:y
        var Float a1 := angle p:out_x p:out_y
        var Float a := min (abs a1-a0-2*pi) (abs a1-a0) (abs a1-a0+2*pi)
        var Float n := (1-cos:a) * (norme p1:x-p:x p1:y-p:y)
        var Float d := (sqrt2-1) * sin:a * 2*sin:a * (norme p:out_x p:out_y)
        var Float f := (shunt n<>0 and d<>0 n/d 0)*reduce;
        p out_x *= f ; p out_y *= f
  c status := mode .or. is_ready
  if pliant_debugging_level>=2
    for (var Int k) 0 c:pts:size-1
      var Pointer:CurvePoint p :> c:pts k
      check p:x=defined and p:y=defined
      check p:in_x=defined and p:in_y=defined
      check p:out_x=defined and p:out_y=defined
(the_function '. compute' Curve Int) extra_module :> the_module "/pliant/math/curve/mode.pli"


method c size -> s
  arg Curve c ; arg Int s
  s := c:pts size

method c 'size :=' s
  arg_rw Curve c ; arg Int s
  c:pts size := s
  c status := c:status .and. .not. is_ready

method c point i -> p
  arg Curve c ; arg Int i ; arg_R CurvePoint p
  p :> c:pts i

method c 'point :=' i p
  arg_rw Curve c ; arg Int i ; arg CurvePoint p
  c:pts i := p
  c status := c:status .and. .not. is_ready
  
method c mode -> m
  arg Curve c ; arg Int m
  m := c:status .and. .not. (cast is_ready Int)

function compute_maximum -> m
  arg Float m
  m := 1
  for (var Int step) 0 100
    var Float f := 1+(exp -step)
    while f<>1 and f*m=defined and -f*m=defined
      m *= f
constant maximum compute_maximum

method c bbox x0 y0 x1 y1
  arg Curve c ; arg_w Float x0 y0 x1 y1
  check (c:status .and. is_ready)<>0
  x0 := maximum ; y0 := maximum
  x1 := -maximum ; y1 := -maximum
  for (var Int i) 0 c:pts:size-1
    var Pointer:CurvePoint p :> c:pts i
    x0 := min x0 (min p:x p:x+p:in_x p:x+p:out_x)
    y0 := min y0 (min p:y p:y+p:in_y p:y+p:out_y)
    x1 := max x1 (max p:x p:x+p:in_x p:x+p:out_x)
    y1 := max y1 (max p:y p:y+p:in_y p:y+p:out_y)
  if x0=maximum
    x0 := undefined ; y0 := undefined ; x1 := undefined ; y1 := undefined


method c pos param0 -> point
  arg Curve c ; arg Float param0 ; arg Point2 point
  check (c:status .and. is_ready)<>0
  var CBool ol := (c:status .and. outline)<>0
  if (c:status .and. bezier)=0
    var Float param := param0 * (shunt ol c:pts:size c:pts:size-1)
    var Int i := bound (cast param-0.5 Int) 0 (shunt ol c:pts:size-1 c:pts:size-2)
    param -= i
    var Float w0 w1
    if c:pts:size=2 and not ol
      w0 := 0 ; w1 := 0
    eif i=0 and not ol and c:pts:0:out_mode=automatic
      w0 := 0 ; w1 := param
    eif i=c:pts:size-2 and not ol and c:pts:(c:pts:size-1):in_mode=automatic
      w0 := 1-param ; w1 := 0
    else
      w0 := weigh 1-param ; w1 := weigh param
    var Pointer:CurvePoint p0 :> c:pts i
    var Pointer:CurvePoint p1 :> c:pts (i+1)%c:pts:size
    point x := w0*(p0:x+param*sqrt2*p0:out_x) + w1*(p1:x+(1-param)*sqrt2*p1:in_x) + (1-w0-w1)*((1-param)*p0:x+param*p1:x)
    point y := w0*(p0:y+param*sqrt2*p0:out_y) + w1*(p1:y+(1-param)*sqrt2*p1:in_y) + (1-w0-w1)*((1-param)*p0:y+param*p1:y)
    check point:x=defined and point:y=defined
  else
    var Float param := param0 * (shunt ol c:pts:size c:pts:size-1)
    var Int i := bound (cast param-0.5 Int) 0 c:pts:size-2
    param -= i
    var Pointer:CurvePoint p0 :> c:pts i
    var Pointer:CurvePoint p1 :> c:pts (i+1)%c:pts:size
    var Float x0 := p0:x
    var Float x1 := p0:out_x+p0:x
    var Float x2 := p1:in_x+p1:x
    var Float x3 := p1:x
    point x := (x3-3*x2+3*x1-x0)*param*param*param+3*(x2-2*x1+x0)*param*param+3*(x1-x0)*param+x0
    var Float y0 := p0:y
    var Float y1 := p0:out_y+p0:y
    var Float y2 := p1:in_y+p1:y
    var Float y3 := p1:y
    point y := (y3-3*y2+3*y1-y0)*param*param*param+3*(y2-2*y1+y0)*param*param+3*(y1-y0)*param+y0;
  check point:x=defined and point:y=defined


#--------------------------------------------------------------------


method c x_param param0 -> x
  arg Curve c ; arg Float param0 x
  check (c:status .and. is_ready)<>0
  var CBool ol := (c:status .and. outline)<>0
  var Float param := param0 * (shunt ol c:pts:size c:pts:size-1)
  var Int i := bound (cast param-0.5 Int) 0 (shunt ol c:pts:size-1 c:pts:size-2)
  param -= i
  var Float w0 w1
  if c:pts:size=2 and not ol
    w0 := 0 ; w1 := 0
  eif i=0 and not ol and c:pts:0:out_mode=automatic
    w0 := 0 ; w1 := param
  eif i=c:pts:size-2 and not ol and c:pts:(c:pts:size-1):in_mode=automatic
    w0 := 1-param ; w1 := 0
  else
    w0 := weigh 1-param ; w1 := weigh param
  var Pointer:CurvePoint p0 :> c:pts i
  var Pointer:CurvePoint p1 :> c:pts (i+1)%c:pts:size
  x := w0*(p0:x+param*sqrt2*p0:out_x) + w1*(p1:x+(1-param)*sqrt2*p1:in_x) + (1-w0-w1)*((1-param)*p0:x+param*p1:x)
  check x=defined

method c y_param param0 -> y
  arg Curve c ; arg Float param0 y
  check (c:status .and. is_ready)<>0
  var CBool ol := (c:status .and. outline)<>0
  var Float param := param0 * (shunt ol c:pts:size c:pts:size-1)
  var Int i := bound (cast param-0.5 Int) 0 (shunt ol c:pts:size-1 c:pts:size-2)
  param -= i
  var Float w0 w1
  if c:pts:size=2 and not ol
    w0 := 0 ; w1 := 0
  eif i=0 and not ol and c:pts:0:out_mode=automatic
    w0 := 0 ; w1 := param
  eif i=c:pts:size-2 and not ol and c:pts:(c:pts:size-1):in_mode=automatic
    w0 := 1-param ; w1 := 0
  else
    w0 := weigh 1-param ; w1 := weigh param
  var Pointer:CurvePoint p0 :> c:pts i
  var Pointer:CurvePoint p1 :> c:pts (i+1)%c:pts:size
  y := w0*(p0:y+param*sqrt2*p0:out_y) + w1*(p1:y+(1-param)*sqrt2*p1:in_y) + (1-w0-w1)*((1-param)*p0:y+param*p1:y)
  check y=defined


method c param_x x deltax -> param
  arg Curve c ; arg Float x deltax param
  check (c:status .and. is_ready)<>0
  var Float p0 := 0 ; var Float x0 := c x_param p0 ; if x0=x return:p0
  var Float p1 := 1 ; var Float x1 := c x_param p1 ; if x1=x return:p1
  if x0>x1  { swap x0 x1 ; swap p0 p1 }
  if x<x0 return:p0
  if x>x1 return:p1
  var Int mode := 0
  var Float fm pm xm
  while true
    if mode=0
      fm := (x-x0)/(x1-x0)
      pm := p0+fm*(p1-p0)
    eif  mode=1
      fm := (x-x0)/(x1-x0)
      fm := shunt fm>=0.5 0.5+(fm-0.5)*(fm-0.5)*2 0.5-(0.5-fm)*(0.5-fm)*2
      pm := p0+fm*(p1-p0)
    else
      pm := (p0+p1)/2
    check pm=defined
    if pm=p0 or pm=p1
      return pm
    check (p1>p0 and pm>p0 and pm<p1) or (p1<p0 and pm>p1 and pm<p0)
    xm := c x_param pm
    if (abs xm-x)<=deltax
      return pm
    eif xm<x
      x0 := xm ; p0 := pm
      mode := shunt fm<=small (min mode+1 2) (max mode-1 0)
    eif xm>x
      x1 := xm ; p1 := pm
      mode := shunt fm>big (min mode+1 2) (max mode-1 0) 
    else
      return pm

method c param_y y deltay -> param
  arg Curve c ; arg Float y deltay param
  check (c:status .and. is_ready)<>0
  var Float p0 := 0 ; var Float y0 := c y_param p0 ; if y0=y return:p0
  var Float p1 := 1 ; var Float y1 := c y_param p1 ; if y1=y return:p1
  if y0>y1  { swap y0 y1 ; swap p0 p1 }
  if y<y0 return:p0
  if y>y1 return:p1
  var Int mode := 0
  var Float fm pm ym
  while true
    if mode=0
      fm := (y-y0)/(y1-y0)
      pm := p0+fm*(p1-p0)
    eif  mode=1
      fm := (y-y0)/(y1-y0)
      fm := shunt fm>=0.5 0.5+(fm-0.5)*(fm-0.5)*2 0.5-(0.5-fm)*(0.5-fm)*2
      pm := p0+fm*(p1-p0)
    else
      pm := (p0+p1)/2
    check pm=defined
    if pm=p0 or pm=p1
      return pm
    check (p1>p0 and pm>p0 and pm<p1) or (p1<p0 and pm>p1 and pm<p0)
    ym := c y_param pm
    if (abs ym-y)<=deltay
      return pm
    eif ym<y
      y0 := ym ; p0 := pm
      mode := shunt fm<=small (min mode+1 2) (max mode-1 0)
    eif ym>y
      y1 := ym ; p1 := pm
      mode := shunt fm>big (min mode+1 2) (max mode-1 0) 
    else
      return pm


method c y x deltax -> y
  arg Curve c ; arg Float x deltax y
  check (c:status .and. is_ready)<>0 and (c:status .and. y_from_x)<>0
  var Float p := c param_x x deltax
  y := c y_param p

method c x y deltay -> x
  arg Curve c ; arg Float y deltay x
  check (c:status .and. is_ready)<>0 and (c:status .and. x_from_y)<>0
  var Float p := c param_y y deltay
  x := c x_param p


#--------------------------------------------------------------------


method c polyline table t0 p0 t1 p1 m epsilon
  arg Curve c ; arg_rw Array:Point2 table ; arg Float t0 t1 epsilon ; arg Transform2 m ; arg Point2 p0 p1
  var Float t05 := (t0+t1)/2.0
  var Point2 p05 := m (c pos t05)
  if (distance (point (p0:x+p1:x)*0.5 (p0:y+p1:y)*0.5) p05)<=epsilon
    table += p05 ; table += p1
  else
    c polyline table t0 p0 t05 p05 m epsilon
    c polyline table t05 p05 t1 p1 m epsilon

method c polyline m epsilon -> table
  arg Curve c ; arg Transform2 m ; arg Float epsilon ; arg Array:Point2 table
  check (c:status .and. is_ready)<>0
  check epsilon>0
  table size := 1 ; table 0 := m (c pos 0)
  var Int n := shunt (c:status .and. outline)<>0 c:pts:size c:pts:size-1
  for (var Int i) 0 n-1
    var Pointer:CurvePoint cp0 :> c:pts i
    var Pointer:CurvePoint cp1 :> c:pts (i+1)%c:pts:size
    if cp0:out_x=0 and cp0:out_y=0 and cp1:in_x=0 and cp1:in_y=0
      table += m (point cp1:x cp1:y)
    else
      var Float t0 := i/n ; var Point2 p0 := table table:size-1
      var Float t05 := (i+0.5)/n ; var Point2 p05 := m (c pos t05)
      var Float t1 := (i+1)/n ; var Point2 p1 := m (c pos t1)
      c polyline table  t0 p0  t05 p05 m epsilon
      c polyline table  t05 p05 t1 p1 m epsilon


method c polyline table p0 p1 m epsilon
  arg Curve c ; arg_rw Array:Point3 table ; arg Transform2 m ; arg Float epsilon ; arg Point3 p0 p1
  var Point3 p05
  p05 z := (p0:z+p1:z)/2.0
  addressof:p05 map Point2 := m (c pos p05:z)
  if (distance (point (p0:x+p1:x)*0.5 (p0:y+p1:y)*0.5) (addressof:p05 map Point2))<=epsilon
    table += p05 ; table+=p1
  else
    c polyline table p0 p05 m epsilon
    c polyline table p05 p1 m epsilon

method c polyline3 m epsilon -> table
  arg Curve c ; arg Transform2 m ; arg Float epsilon ; arg Array:Point3 table
  check (c:status .and. is_ready)<>0
  check epsilon>0
  table size := 1 ; table:0 z := 0 ; (addressof table:0) map Point2 := m (c pos 0)
  var Int n := shunt (c:status .and. outline)<>0 c:pts:size c:pts:size-1
  for (var Int i) 0 n-1
    var Pointer:CurvePoint cp0 :> c:pts i
    var Pointer:CurvePoint cp1 :> c:pts (i+1)%c:pts:size
    if cp0:out_x=0 and cp0:out_y=0 and cp1:in_x=0 and cp1:in_y=0
      var Point3 p1 ; p1 z := (i+1)/n ; addressof:p1 map Point2 := m (point cp1:x cp1:y)
      table += p1
    else
      var Point3 p0 := table table:size-1
      var Point3 p05 ; p05 z := (i+0.5)/n ; addressof:p05 map Point2 := m (c pos p05:z)
      var Point3 p1 ; p1 z := (i+1)/n ; addressof:p1 map Point2 := m (c pos p1:z)
      c polyline table p0 p05 m epsilon
      c polyline table p05 p1 m epsilon


export Curve '. reset' '. through' '. angle' '+=' '. compute'
export '. size' '. size :=' '. point' '. point :=' '. mode' '. bbox' '. pos'
export '. x_param' '. y_param' '. param_x' '. param_y' '. x' '. y'
export '. polyline' '. polyline3'
export 'cast Status'