open-axiom repository from github
\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra clip.spad}
\author{Clifton J. Williamson}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package CLIP TwoDimensionalPlotClipping}
<<package CLIP TwoDimensionalPlotClipping>>=
)abbrev package CLIP TwoDimensionalPlotClipping
++ Automatic clipping for 2-dimensional plots
++ Author: Clifton J. Williamson
++ Date Created: 22 December 1989
++ Date Last Updated: 10 July 1990
++ Keywords: plot, singularity
++ Examples:
++ References:
TwoDimensionalPlotClipping(): Exports == Implementation where
++ The purpose of this package is to provide reasonable plots of
++ functions with singularities.
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)
Exports ==> 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.
Implementation ==> 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")
Pnan?:Pt ->Boolean
Pnan? p == any?(nan?,p)
iClipParametric(pointLists,fraction,scale) ==
import Point SF
import List Point SF
-- error checks and special cases
negative? fraction 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(norm(#1) < norm(#2),select(not Pnan? #1,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) ==
-- sayBrightly([" clip: "::OutputForm]$List(OutputForm))$Lisp
negative? fraction 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)]
#(pointLists := listBranches plot) > 1 =>
error "clipDraw: plot has more than one branch"
empty?(pointList := first pointLists) =>
[nil(),xVals,segment(0,0)]
sortedList := sort(yCoord(#1) < yCoord(#2),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 nan? yMin then yMin : SF := 0
if nan? 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,_
(yCoord(#1) < yMax) and (yCoord(#1) > yMin),xMin,xMax,yMin,yMax)
yMin := yCoord(sortedList.clipNum)
yMax := yCoord(sortedList.(n - 1 - clipNum))
if nan? yMin then yMin : SF := 0
if nan? yMax then yMax : SF := 0
for list in lists repeat
for pt in list repeat
if not nan?(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 nan? x then
if nan? y then
r:SF := 0
else
r:SF := y**2
else
if nan? 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,_
(xCoord(#1) <= xMax) and (xCoord(#1) >= xMin) and _
(yCoord(#1) <= yMax) and (yCoord(#1) >= 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: local 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)
@
\section{License}
<<license>>=
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
-- - Redistributions of source code must retain the above copyright
-- notice, this list of conditions and the following disclaimer.
--
-- - Redistributions in binary form must reproduce the above copyright
-- notice, this list of conditions and the following disclaimer in
-- the documentation and/or other materials provided with the
-- distribution.
--
-- - Neither the name of The Numerical ALgorithms Group Ltd. nor the
-- names of its contributors may be used to endorse or promote products
-- derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
@
<<*>>=
<<license>>
<<package CLIP TwoDimensionalPlotClipping>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}