tube.spad

)abbrev domain TUBE TubePlot
++ Author: Clifton J. Williamson
++ Date Created: Bastille Day 1989
++ Keywords:
++ Examples:
++ Description:
++   Package for constructing tubes around 3-dimensional parametric curves.
++ Domain of tubes around 3-dimensional parametric curves.
TubePlot(Curve) : Exports == Implementation where
  Curve : PlottableSpaceCurveCategory
  B   ==> Boolean
  L   ==> List
  Pt  ==> Point DoubleFloat

  Exports ==> with
    getCurve : % -> Curve
      ++ getCurve(t) returns the \spadtype{PlottableSpaceCurveCategory}
      ++ representing the parametric curve of the given tube plot t.
    listLoops : % -> L L Pt
      ++ listLoops(t) returns the list of lists of points, or the 'loops',
      ++ of the given tube plot t.
    closed? : % -> B
      ++ closed?(t) tests whether the given tube plot t is closed.
    open? : % -> B
      ++ open?(t) tests whether the given tube plot t is open.
    setClosed : (%, B) -> B
      ++ setClosed(t, b) declares the given tube plot t to be closed if
      ++ b is true, or if b is false, t is set to be open.
    tube : (Curve, L L Pt, B) -> %
      ++ tube(c, ll, b) creates a tube of the domain \spadtype{TubePlot} from a
      ++ space curve c of the category \spadtype{PlottableSpaceCurveCategory},
      ++ a list of lists of points (loops) ll and a boolean b which if
      ++ true indicates a closed tube, or if false an open tube.

  Implementation ==> add

--% representation

    Rep := Record(parCurve : Curve, loops : L L Pt, closedTube? : B)

    getCurve plot == plot.parCurve

    listLoops plot == plot.loops

    closed? plot == plot.closedTube?
    open? plot   == not plot.closedTube?

    setClosed(plot, flag) == plot.closedTube? := flag

    tube(curve, ll, b) == [curve, ll, b]

)abbrev package TUBETOOL TubePlotTools
++ Author: Clifton J. Williamson
++ Date Created: Bastille Day 1989
++ Keywords:
++ Examples:
++ Description:
++   Tools for constructing tubes around 3-dimensional parametric curves.
TubePlotTools() : Exports == Implementation where
  I   ==> Integer
  SF  ==> DoubleFloat
  L   ==> List
  Pt  ==> Point SF

  Exports ==> with
    point : (SF, SF, SF, SF) -> Pt
      ++ point(x1, x2, x3, c) creates and returns a point from the three
      ++ specified coordinates \spad{x1}, \spad{x2}, \spad{x3}, and also
      ++ a fourth coordinate, c, which is generally used to specify the
      ++ color of the point.
    "*" : (SF,Pt) -> Pt
      ++ s * p returns a point whose coordinates are the scalar multiple
      ++ of the point p by the scalar s, preserving the color, or fourth
      ++ coordinate, of p.
    "+" : (Pt,Pt) -> Pt
      ++ p + q computes and returns a point whose coordinates are the sums
      ++ of the coordinates of the two points \spad{p} and \spad{q}, using
      ++ the color, or fourth coordinate, of the first point \spad{p}
      ++ as the color also of the point \spad{q}.
    "-" : (Pt,Pt) -> Pt
      ++ p - q computes and returns a point whose coordinates are the
      ++ differences of the coordinates of two points \spad{p} and \spad{q},
      ++ using the color, or fourth coordinate, of the first point \spad{p}
      ++ as the color also of the point \spad{q}.
    dot : (Pt, Pt) -> SF
      ++ dot(p, q) computes the dot product of the two points \spad{p}
      ++ and \spad{q} using only the first three coordinates, and returns
      ++ the resulting \spadtype{DoubleFloat}.
    cross : (Pt, Pt) -> Pt
      ++ cross(p, q) computes the cross product of the two points \spad{p}
      ++ and \spad{q} using only the first three coordinates, and keeping
      ++ the color of the first point p.  The result is returned as a point.
    unitVector : Pt -> Pt
      ++ unitVector(p) creates the unit vector of the point p and returns
      ++ the result as a point.  Note: \spad{unitVector(p) = p/|p|}.
    cosSinInfo : I -> L L SF
      ++ cosSinInfo(n) returns the list of lists of values for n, in the
      ++ form:  \spad{[[cos(n - 1) a, sin(n - 1) a], ..., [cos 2 a, sin 2 a], [cos a, sin a]]}
      ++ where \spad{a = 2 pi/n}.  Note: n should be greater than 2.
    loopPoints : (Pt, Pt, Pt, SF, L L SF) -> L Pt
      ++ loopPoints(p, n, b, r, lls) creates and returns a list of points
      ++ which form the loop with radius r, around the center point
      ++ indicated by the point p, with the principal normal vector of
      ++ the space curve at point p given by the point(vector) n, and the
      ++ binormal vector given by the point(vector) b, and a list of lists,
      ++ lls, which is the \spadfun{cosSinInfo} of the number of points
      ++ defining the loop.

  Implementation ==> add
    import from PointPackage(SF)

    point(x, y, z, c) == point(l : L SF := [x, y, z, c])

    getColor : Pt -> SF
    getColor pt == (maxIndex pt > 3 => color pt; 0)

    getColor2 : (Pt, Pt) -> SF
    getColor2(p0, p1) ==
      maxIndex p0 > 3 => color p0
      maxIndex p1 > 3 => color p1
      0

    a * p ==
      l : L SF := [a * xCoord p, a * yCoord p, a * zCoord p, getColor p]
      point l

    p0 + p1 ==
      l : L SF := [xCoord p0 + xCoord p1, yCoord p0 + yCoord p1, _
                   zCoord p0 + zCoord p1, getColor2(p0, p1)]
      point l

    p0 - p1 ==
      l : L SF := [xCoord p0 - xCoord p1, yCoord p0 - yCoord p1, _
                   zCoord p0 - zCoord p1, getColor2(p0, p1)]
      point l

    dot(p0, p1) ==
      (xCoord p0 * xCoord p1) + (yCoord p0 * yCoord p1) +_
        (zCoord p0 * zCoord p1)

    cross(p0, p1) ==
      x0 := xCoord p0; y0 := yCoord p0; z0 := zCoord p0;
      x1 := xCoord p1; y1 := yCoord p1; z1 := zCoord p1;
      l : L SF := [y0 * z1 - y1 * z0, z0 * x1 - z1 * x0, _
                   x0 * y1 - x1 * y0, getColor2(p0, p1)]
      point l

    unitVector p == (inv sqrt dot(p, p)) * p

    cosSinInfo n ==
      ans : L L SF := []
      theta : SF := 2 * pi()/n
      for i in 1..(n-1) repeat             --!! make more efficient
        angle := i * theta
        ans := concat([cos angle, sin angle], ans)
      ans

    loopPoints(ctr, pNorm, bNorm, rad, cosSin) ==
      ans : L Pt := []
      while not(empty?(cosSin)) repeat
        cossin := first cosSin; cos := first cossin; sin := second cossin
        ans := cons(ctr + rad * (cos * pNorm + sin * bNorm), ans)
        cosSin := rest cosSin
      pt := ctr + rad * pNorm
      concat(pt, concat(ans, pt))

)abbrev package EXPRTUBE ExpressionTubePlot
++ Author: Clifton J. Williamson
++ Date Created: Bastille Day 1989
++ Keywords:
++ Examples:
++ Package for constructing tubes around 3-dimensional parametric curves.
-- Currently unused.
ExpressionTubePlot() : Exports == Implementation where
  B   ==> Boolean
  I   ==> Integer
  FE  ==> Expression Integer
  SY  ==> Symbol
  SF  ==> DoubleFloat
  L   ==> List
  S   ==> String
  SEG ==> Segment
  F2F ==> MakeFloatCompiledFunction(FE)
  Pt  ==> Point SF
  PLOT3 ==> Plot3D
  TUBE  ==> TubePlot Plot3D

  Exports ==> with
    constantToUnaryFunction : SF -> (SF -> SF)
      ++ constantToUnaryFunction(s) is a local function which takes the
      ++ value of s, which may be a function of a constant, and returns
      ++ a function which always returns the value \spadtype{DoubleFloat} s.
    tubePlot : (FE, FE, FE, SF -> SF, SEG SF, SF -> SF, I) -> TUBE
      ++ tubePlot(f, g, h, colorFcn, a..b, r, n) puts a tube of radius r(t) with
      ++ n points on each circle about the curve \spad{x = f(t)},
      ++ \spad{y = g(t)}, \spad{z = h(t)} for t in \spad{[a, b]}.
      ++ The tube is considered to be open.
    tubePlot : (FE, FE, FE, SF -> SF, SEG SF, SF -> SF, I, S) -> TUBE
      ++ tubePlot(f, g, h, colorFcn, a..b, r, n, s) puts a tube of radius \spad{r(t)}
      ++ with n points on each circle about the curve \spad{x = f(t)},
      ++ \spad{y = g(t)},
      ++ \spad{z = h(t)} for t in \spad{[a,b]}. If s = "closed", the tube is
      ++ considered to be closed; if s = "open", the tube is considered
      ++ to be open.
    tubePlot : (FE, FE, FE, SF -> SF, SEG SF, SF, I) -> TUBE
      ++ tubePlot(f, g, h, colorFcn, a..b, r, n) puts a tube of radius r with
      ++ n points on each circle about the curve \spad{x = f(t)},
      ++ \spad{y = g(t)}, \spad{z = h(t)} for t in \spad{[a, b]}.
      ++ The tube is considered to be open.
    tubePlot : (FE, FE, FE, SF -> SF, SEG SF, SF, I, S) -> TUBE
      ++ tubePlot(f, g, h, colorFcn, a..b, r, n, s) puts a tube of radius r with
      ++ n points on each circle about the curve \spad{x = f(t)},
      ++ \spad{y = g(t)}, \spad{z = h(t)} for t in \spad{[a, b]}.
      ++ If s = "closed", the tube is
      ++ considered to be closed; if s = "open", the tube is considered
      ++ to be open.

  Implementation ==> add
    import from Plot3D
    import from F2F
    import from TubePlotTools

--% variables

    getVariable : (FE, FE, FE) -> SY
    getVariable(x, y, z) ==
      varList1 := variables x
      varList2 := variables y
      varList3 := variables z
      (not (# varList1 <= 1)) or (not (# varList2 <= 1)) or _
       (not (# varList3 <= 1)) =>
        error "tubePlot: only one variable may be used"
      empty?(varList1) =>
        empty?(varList2) =>
          empty?(varList3) =>
            error "tubePlot: a variable must appear in functions"
          first varList3
        t2 := first varList2
        empty?(varList3) => t2
        not (first varList3 = t2) =>
          error "tubePlot: only one variable may be used"
      t1 := first varList1
      empty?(varList2) =>
        empty?(varList3) => t1
        not (first varList3 = t1) =>
          error "tubePlot: only one variable may be used"
        t1
      t2 := first varList2
      empty?(varList3) =>
        not (t1 = t2) =>
          error "tubePlot: only one variable may be used"
        t1
      not (first varList3 = t1) or not (t2 = t1) =>
        error "tubePlot: only one variable may be used"
      t1

--% tubes: variable radius

    tubePlot(x : FE, y : FE, z : FE, colorFcn : SF -> SF, _
             tRange : SEG SF, radFcn : SF -> SF, n : I, string : S) ==
      -- check value of n
      n < 3 => error "tubePlot: n should be at least 3"
      -- check string
      flag : B :=
        string = "closed" => true
        string = "open" => false
        error "tubePlot: last argument should be open or closed"
      -- check variables
      t := getVariable(x, y, z)
      -- coordinate functions
      xFunc := makeFloatFunction(x, t)
      yFunc := makeFloatFunction(y, t)
      zFunc := makeFloatFunction(z, t)
      -- derivatives of coordinate functions
      xp := differentiate(x, t)
      yp := differentiate(y, t)
      zp := differentiate(z, t)
      -- derivative of arc length
      sp := sqrt(xp ^ 2 + yp ^ 2 + zp ^ 2)
      -- coordinates of unit tangent vector
      Tx := xp/sp; Ty := yp/sp; Tz := zp/sp
      -- derivatives of coordinates of unit tangent vector
      Txp := differentiate(Tx, t)
      Typ := differentiate(Ty, t)
      Tzp := differentiate(Tz, t)
      -- K = curvature = length of curvature vector
      K := sqrt(Txp ^ 2 + Typ ^ 2 + Tzp ^ 2)
      -- coordinates of principal normal vector
      Nx := Txp / K; Ny := Typ / K; Nz := Tzp / K
      -- functions SF->SF giving coordinates of principal normal vector
      NxFunc := makeFloatFunction(Nx, t)
      NyFunc := makeFloatFunction(Ny, t)
      NzFunc := makeFloatFunction(Nz, t)
      -- coordinates of binormal vector
      Bx := Ty * Nz - Tz * Ny
      By := Tz * Nx - Tx * Nz
      Bz := Tx * Ny - Ty * Nx
      -- functions SF -> SF giving coordinates of binormal vector
      BxFunc := makeFloatFunction(Bx, t)
      ByFunc := makeFloatFunction(By, t)
      BzFunc := makeFloatFunction(Bz, t)
      -- create Plot3D
      parPlot := plot(xFunc, yFunc, zFunc, colorFcn, tRange)
      tvals := first tValues parPlot
      curvePts := first listBranches parPlot
      cosSin := cosSinInfo n
      loopList : L L Pt := []
      while not(empty?(tvals)) repeat
        -- note: tvals and curvePts have the same number of elements
        tval := first tvals; tvals := rest tvals
        ctr := first curvePts; curvePts := rest curvePts
        pNormList : L SF :=
          [NxFunc tval, NyFunc tval, NzFunc tval, colorFcn tval]
        pNorm : Pt := point pNormList
        bNormList : L SF :=
          [BxFunc tval, ByFunc tval, BzFunc tval, colorFcn tval]
        bNorm : Pt := point bNormList
        lps := loopPoints(ctr, pNorm, bNorm, radFcn tval, cosSin)
        loopList := cons(lps, loopList)
      tube(parPlot, reverse! loopList, flag)

    tubePlot(x : FE, y : FE, z : FE, colorFcn : SF -> SF, _
             tRange : SEG SF, radFcn : SF -> SF, n : I) ==
      tubePlot(x,y,z,colorFcn,tRange,radFcn,n,"open")

--% tubes: constant radius

    project : (SF, SF) -> SF
    project(x, y) == x

    constantToUnaryFunction x == s +-> project(x, s)

    tubePlot(x : FE, y : FE, z : FE, colorFcn : SF -> SF, _
             tRange : SEG SF, rad : SF, n : I, s : S) ==
      tubePlot(x, y, z, colorFcn, tRange, constantToUnaryFunction rad, n, s)

    tubePlot(x : FE, y : FE, z : FE, colorFcn : SF -> SF, _
             tRange : SEG SF, rad : SF, n : I) ==
      tubePlot(x,y,z,colorFcn,tRange,rad,n,"open")

)abbrev package NUMTUBE NumericTubePlot
++ Author: Clifton J. Williamson
++ Date Created: Bastille Day 1989
++ Keywords:
++ Examples:
++ Package for constructing tubes around 3-dimensional parametric curves.
NumericTubePlot(Curve) : Exports == Implementation where
  Curve : PlottableSpaceCurveCategory
  B   ==> Boolean
  I   ==> Integer
  SF  ==> DoubleFloat
  L   ==> List
  S   ==> String
  SEG ==> Segment
  Pt  ==> Point SF
  TUBE ==> TubePlot Curve
  Triad ==> Record(tang : Pt, norm : Pt, bin : Pt)

  Exports ==> with
    tube : (Curve, SF, I) -> TUBE
      ++ tube(c, r, n) creates a tube of radius r around the curve c.

  Implementation ==> add
    import from TubePlotTools

    LINMAX  := convert(0.995)@SF
    XHAT := point(1, 0, 0, 0)
    YHAT := point(0, 1, 0, 0)
    PREV0 := point(1, 1, 0, 0)
    PREV := PREV0

    colinearity : (Pt, Pt) -> SF
    colinearity(x, y) == dot(x, y)^2/(dot(x, x) * dot(y, y))

    orthog : (Pt, Pt) -> Pt
    orthog(x, y) ==
      if colinearity(x, y) > LINMAX then y := PREV
      if colinearity(x, y) > LINMAX then
        y := (colinearity(x, XHAT) < LINMAX => XHAT; YHAT)
      a := -dot(x, y)/dot(x, x)
      PREV := a*x + y

    poTriad : (Pt, Pt, Pt) -> Triad
    poTriad(pl, po, pr) ==
      -- use divided difference for t.
      t := unitVector(pr - pl)
      -- compute n as orthogonal to t in plane containing po.
      pol := pl - po
      n   := unitVector orthog(t, pol)
      [t, n, cross(t, n)]

    curveTriads : L Pt -> L Triad
    curveTriads l ==
      (k := #l) < 2 => error "Need at least 2 points to specify a curve"
      PREV := PREV0
      k = 2 =>
        t := unitVector(second l - first l)
        n := unitVector(t - XHAT)
        b := cross(t, n)
        triad : Triad := [t, n, b]
        [triad, triad]
      -- compute interior triads using divided differences
      midtriads : L Triad :=
        [poTriad(pl, po, pr) for pl in l for po in rest l _
               for pr in rest rest l]
      -- compute first triad using a forward difference
      x := first midtriads
      t := unitVector(second l - first l)
      n := unitVector orthog(t, x.norm)
      begtriad : Triad := [t, n, cross(t, n)]
      -- compute last triad using a backward difference
      x := last midtriads
      -- efficiency!!
      t := unitVector(l.k - l.(k-1))
      n := unitVector orthog(t, x.norm)
      endtriad : Triad := [t, n, cross(t, n)]
      concat(begtriad, concat(midtriads, endtriad))

    curveLoops : (L Pt, SF, I) -> L L Pt
    curveLoops(pts, r, nn) ==
      triads := curveTriads pts
      cosSin := cosSinInfo nn
      loops : L L Pt := []
      for pt in pts for triad in triads repeat
        n := triad.norm; b := triad.bin
        loops := concat(loopPoints(pt, n, b, r, cosSin), loops)
      reverse! loops

    tube(curve, r, n) ==
      n < 3 => error "tube: n should be at least 3"
      brans := listBranches curve
      loops : L L Pt := []
      for bran in brans repeat
        loops := concat(loops, curveLoops(bran, r, n))
      tube(curve, loops, false)

--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.