axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / CLIP.spad

)abbrev package CLIP TwoDimensionalPlotClipping
++ Author: Clifton J. Williamson
++ Date Created: 22 December 1989
++ Date Last Updated: 10 July 1990
++ Description:
++ Automatic clipping for 2-dimensional plots
++ The purpose of this package is to provide reasonable plots of
++ functions with singularities.
 
TwoDimensionalPlotClipping() : SIG == CODE where

  B      ==> Boolean
  L      ==> List
  SEG    ==> Segment
  RN     ==> Fraction Integer
  SF     ==> DoubleFloat
  Pt     ==> Point DoubleFloat
  PLOT   ==> Plot
  CLIPPED ==> Record(brans: L L Pt,xValues: SEG SF,yValues: SEG SF)
 
  SIG ==> with

    clip : PLOT -> CLIPPED
      ++ clip(p) performs two-dimensional clipping on a plot, p, from
      ++ the domain \spadtype{Plot} for the graph of one variable,
      ++ \spad{y = f(x)}; the default parameters \spad{1/4} for the fraction 
      ++ and \spad{5/1} for the scale are used in the \spadfun{clip} function.

    clip : (PLOT,RN,RN) -> CLIPPED
      ++ clip(p,frac,sc) performs two-dimensional clipping on a plot, p, 
      ++ from the domain \spadtype{Plot} for the graph of one variable
      ++ \spad{y = f(x)}; the fraction parameter is specified by \spad{frac} 
      ++ and the scale parameter is specified by \spad{sc} for use in the 
      ++ \spadfun{clip} function.

    clipParametric : PLOT -> CLIPPED
      ++ clipParametric(p) performs two-dimensional clipping on a plot,
      ++ p, from the domain \spadtype{Plot} for the parametric curve
      ++ \spad{x = f(t)}, \spad{y = g(t)}; the default parameters \spad{1/2} 
      ++ for the fraction and \spad{5/1} for the scale are used in the 
      ++ \fakeAxiomFun{iClipParametric} subroutine, which is called by this
      ++ function.

    clipParametric : (PLOT,RN,RN) -> CLIPPED
      ++ clipParametric(p,frac,sc) performs two-dimensional clipping on a 
      ++ plot, p, from the domain \spadtype{Plot} for the parametric curve
      ++ \spad{x = f(t)}, \spad{y = g(t)}; the fraction parameter is 
      ++ specified by \spad{frac} and the scale parameter is specified 
      ++ by \spad{sc} for use in the \fakeAxiomFun{iClipParametric} subroutine,
      ++ which is called by this function.

    clipWithRanges : (L L Pt,SF,SF,SF,SF) -> CLIPPED
      ++ clipWithRanges(pointLists,xMin,xMax,yMin,yMax) performs clipping
      ++ on a list of lists of points, \spad{pointLists}.  Clipping is 
      ++ done within the specified ranges of \spad{xMin}, \spad{xMax} and 
      ++ \spad{yMin}, \spad{yMax}.  This function is used internally by 
      ++ the \fakeAxiomFun{iClipParametric} subroutine in this package.

    clip : L Pt -> CLIPPED
      ++ clip(l) performs two-dimensional clipping on a curve l, which is
      ++ a list of points; the default parameters \spad{1/2} for the 
      ++ fraction and \spad{5/1} for the scale are used in the 
      ++ \fakeAxiomFun{iClipParametric} subroutine, which is called by this 
      ++ function.

    clip : L L Pt -> CLIPPED
      ++ clip(ll) performs two-dimensional clipping on a list of lists 
      ++ of points, \spad{ll}; the default parameters \spad{1/2} for
      ++ the fraction and \spad{5/1} for the scale are used in the 
      ++ \fakeAxiomFun{iClipParametric} subroutine, which is called by this 
      ++ function.
 
  CODE ==> add

    import PointPackage(DoubleFloat)
    import ListFunctions2(Point DoubleFloat,DoubleFloat)
 
    point:(SF,SF) -> Pt
    intersectWithHorizLine:(SF,SF,SF,SF,SF) -> Pt
    intersectWithVertLine:(SF,SF,SF,SF,SF) -> Pt
    intersectWithBdry:(SF,SF,SF,SF,Pt,Pt) -> Pt
    discardAndSplit: (L Pt,Pt -> B,SF,SF,SF,SF) -> L L Pt
    norm: Pt -> SF
    iClipParametric: (L L Pt,RN,RN) -> CLIPPED
    findPt: L L Pt -> Union(Pt,"failed")
    Fnan?: SF ->Boolean
    Pnan?:Pt ->Boolean

    Fnan? x == x~=x

    Pnan? p == any?(Fnan?,p)
   
    iClipParametric(pointLists,fraction,scale) ==
      -- error checks and special cases
      (fraction < 0) or (fraction > 1) =>
        error "clipDraw: fraction should be between 0 and 1"
      empty? pointLists => [nil(),segment(0,0),segment(0,0)]
      -- put all points together , sort them according to norm
      sortedList := sort((x:Pt,y:Pt):Boolean +-> norm(x) < norm(y),
                     select((z:Pt):Boolean +-> not Pnan? z,concat pointLists))
      empty? sortedList => [nil(),segment(0,0),segment(0,0)]
      n := # sortedList 
      num := numer fraction
      den := denom fraction
      clipNum := (n * num) quo den
      lastN := n - 1 - clipNum
      firstPt := first sortedList
      xMin : SF := xCoord firstPt
      xMax : SF := xCoord firstPt
      yMin : SF := yCoord firstPt 
      yMax : SF := yCoord firstPt
      -- calculate min/max for the first (1-fraction)*N points
      -- this contracts the range
      -- this unnecessarily clips monotonic functions 
      -- (step-function, x^(high power),etc.)
      for k in 0..lastN  for pt in rest sortedList repeat
        xMin := min(xMin,xCoord pt)
        xMax := max(xMax,xCoord pt)
        yMin := min(yMin,yCoord pt)
        yMax := max(yMax,yCoord pt)
      xDiff := xMax - xMin; yDiff := yMax - yMin
      xDiff = 0 =>
        yDiff = 0 =>
          [pointLists,segment(xMin-1,xMax+1),segment(yMin-1,yMax+1)]
        [pointLists,segment(xMin-1,xMax+1),segment(yMin,yMax)]
      yDiff = 0 =>
        [pointLists,segment(xMin,xMax),segment(yMin-1,yMax+1)]
      numm := numer scale; denn := denom scale
      -- now expand the range by scale
      xMin := xMin - (numm :: SF) * xDiff / (denn :: SF)
      xMax := xMax + (numm :: SF) * xDiff / (denn :: SF)
      yMin := yMin - (numm :: SF) * yDiff / (denn :: SF)
      yMax := yMax + (numm :: SF) * yDiff / (denn :: SF)
      -- clip with the calculated range
      newclip:=clipWithRanges(pointLists,xMin,xMax,yMin,yMax)
      -- if we split the lists use the new clip
      # (newclip.brans) > # pointLists   => newclip
      -- calculate extents
      xs :L SF:= map (xCoord,sortedList)
      ys :L SF:= map (yCoord,sortedList)
      xMin :SF :=reduce (min,xs)
      yMin :SF :=reduce (min,ys)
      xMax :SF :=reduce (max,xs)
      yMax :SF :=reduce (max,ys) 
      xseg:SEG SF :=xMin..xMax
      yseg:SEG SF :=yMin..yMax
      -- return original
      [pointLists,xseg,yseg]@CLIPPED
 
    point(xx,yy) == point(l : L SF := [xx,yy])
 
    intersectWithHorizLine(x1,y1,x2,y2,yy) ==
      x1 = x2 => point(x1,yy)
      point(x1 + (x2 - x1)*(yy - y1)/(y2 - y1),yy)
 
    intersectWithVertLine(x1,y1,x2,y2,xx) ==
      y1 = y2 => point(xx,y1)
      point(xx,y1 + (y2 - y1)*(xx - x1)/(x2 - x1))
 
    intersectWithBdry(xMin,xMax,yMin,yMax,pt1,pt2) ==
      -- pt1 is in rectangle, pt2 is not
      x1 := xCoord pt1; y1 := yCoord pt1
      x2 := xCoord pt2; y2 := yCoord pt2
      if y2 > yMax then
        pt2 := intersectWithHorizLine(x1,y1,x2,y2,yMax)
        x2 := xCoord pt2; y2 := yCoord pt2
      if y2 < yMin then
        pt2 := intersectWithHorizLine(x1,y1,x2,y2,yMin)
        x2 := xCoord pt2; y2 := yCoord pt2
      if x2 > xMax then
        pt2 := intersectWithVertLine(x1,y1,x2,y2,xMax)
        x2 := xCoord pt2; y2 := yCoord pt2
      if x2 < xMin then
        pt2 := intersectWithVertLine(x1,y1,x2,y2,xMin)
      pt2
 
    discardAndSplit(pointList,pred,xMin,xMax,yMin,yMax) ==
      ans : L L Pt := nil()
      list : L Pt  := nil()
      lastPt? : B  := false
      lastPt : Pt  := point(0,0)
      while not empty? pointList repeat
        pt := first pointList
        pointList := rest pointList
        pred(pt) =>
          if (empty? list) and lastPt? then
            bdryPt := intersectWithBdry(xMin,xMax,yMin,yMax,pt,lastPt)
            -- print bracket [ coerce bdryPt ,coerce pt ]
            --list := cons(bdryPt,list)
          list := cons(pt,list)
        if not empty? list then
          bdryPt := intersectWithBdry(xMin,xMax,yMin,yMax,first list,pt)
          -- print bracket [ coerce bdryPt,coerce first list]
          --list := cons(bdryPt,list)
          ans := cons( list,ans)
        lastPt := pt 
        lastPt? := true
        list := nil()
      empty? list => ans
      reverse_! cons(reverse_! list,ans)
 
    clip(plot,fraction,scale) ==
      (fraction < 0) or (fraction > 1/2) =>
        error "clipDraw: fraction should be between 0 and 1/2"
      xVals := xRange plot
      empty?(pointLists := listBranches plot) =>
        [nil(),xVals,segment(0,0)]
      more?(pointLists := listBranches plot,1) =>
        error "clipDraw: plot has more than one branch"
      empty?(pointList := first pointLists) =>
        [nil(),xVals,segment(0,0)]
      sortedList := sort((x,y)+->yCoord(x) < yCoord(y),pointList)
      n := # sortedList; num := numer fraction; den := denom fraction
      clipNum := (n * num) quo den
      -- throw out points with large and small y-coordinates
      yMin := yCoord(sortedList.clipNum)
      yMax := yCoord(sortedList.(n - 1 - clipNum))
      if Fnan? yMin then yMin : SF := 0
      if Fnan? yMax then yMax : SF := 0
      (yDiff := yMax - yMin) = 0 =>
        [pointLists,xRange plot,segment(yMin - 1,yMax + 1)]
      numm := numer scale; denn := denom scale
      xMin := lo xVals; xMax := hi xVals
      yMin := yMin - (numm :: SF) * yDiff / (denn :: SF)
      yMax := yMax + (numm :: SF) * yDiff / (denn :: SF)
      lists := discardAndSplit(pointList,_
         x +-> (yCoord(x) < yMax) and (yCoord(x) > yMin),
           xMin,xMax,yMin,yMax)
      yMin := yCoord(sortedList.clipNum)
      yMax := yCoord(sortedList.(n - 1 - clipNum))
      if Fnan? yMin then yMin : SF := 0
      if Fnan? yMax then yMax : SF := 0
      for list in lists repeat
        for pt in list repeat
          if not Fnan?(yCoord pt) then
            yMin := min(yMin,yCoord pt)
            yMax := max(yMax,yCoord pt)
      [lists,xVals,segment(yMin,yMax)]
 
    clip(plot:PLOT) == clip(plot,1/4,5/1)
 
    norm(pt) == 
      x := xCoord(pt); y := yCoord(pt)
      if Fnan? x then
        if Fnan? y then
          r:SF := 0
        else
          r:SF := y**2
      else
        if Fnan? y then
          r:SF := x**2
        else
          r:SF := x**2 + y**2
      r
 
    findPt lists ==
      for list in lists repeat
        not empty? list => 
             for p in list repeat 
               not Pnan? p => return p
      "failed"

    clipWithRanges(pointLists,xMin,xMax,yMin,yMax) ==
      lists : L L Pt := nil()
      for pointList in pointLists repeat
        lists := concat(lists,discardAndSplit(pointList,
          (x:Pt):Boolean +-> (xCoord(x) <= xMax) and (xCoord(x) >= xMin) and 
            (yCoord(x) <= yMax) and (yCoord(x) >= yMin), 
             xMin,xMax,yMin,yMax))
      (pt := findPt lists) case "failed" =>
        [nil(),segment(0,0),segment(0,0)]
      firstPt := pt :: Pt
      xMin : SF := xCoord firstPt; xMax : SF := xCoord firstPt
      yMin : SF := yCoord firstPt; yMax : SF := yCoord firstPt
      for list in lists repeat
        for pt in list repeat
          if not Pnan? pt then
            xMin := min(xMin,xCoord pt)
            xMax := max(xMax,xCoord pt)
            yMin := min(yMin,yCoord pt)
            yMax := max(yMax,yCoord pt)
      [lists,segment(xMin,xMax),segment(yMin,yMax)]
 
    clipParametric(plot,fraction,scale) ==
      iClipParametric(listBranches plot,fraction,scale)
 
    clipParametric plot == clipParametric(plot,1/2,5/1)
 
    clip(l: L Pt)   == iClipParametric(list l,1/2,5/1)

    clip(l: L L Pt) == iClipParametric(l,1/2,5/1)