dmj5
Class DMJ_PolyCurve

Object
  extended by dmj5:DMJ_PolyCurve

class 
Object:DMJ_PolyCurve

PolyCurve tool class. PolyCurve is a helper class that you can use to define, manage, and render a shape composed of connected lines and curves. Each segment may be one of these: - a straight line - a quadratic Bézier curve; these are used to create TrueType font shapes, because they are mathematically easy to render at any resolution - a cubic Bézier curve; these are used in high-end vector illustration programs (e.g. Adobe Illustrator) and in PostScript fonts, as they are more flexible than quadratic curves, but they are harder to render Cubic Bézier curves are so difficult to render that this class actually does not render them directly; instead, it converts cubic curves to a set of very small straight line segments, and renders from that. Usually this is not noticeable, but you should be aware of this if you use this class in unusual contexts. This class is not derived from anything as it does not have any user-exposed parameters. Your code must provide all the parameters you require. You must call Rasterize() before making any queries about the curve, or you will get erroneous results. Changing curve data after it has been rasterized will not affect any queries. This code is adapted from DeluxeClipping, which I wrote for Janet Parke's course on Masking with Ultra Fractal. The cubic Bézier support is new for UF5. The additional query modes are also new.


Ultra Fractal Source

Toggle UF Source Code Display

 class DMJ_PolyCurve {
   ; PolyCurve tool class.
   ;
   ; PolyCurve is a helper class that you can use to define,
   ; manage, and render a shape composed of connected lines
   ; and curves. Each segment may be one of these:
   ;
   ; - a straight line
   ; - a quadratic Bézier curve; these are used to create
   ;   TrueType font shapes, because they are mathematically
   ;   easy to render at any resolution
   ; - a cubic Bézier curve; these are used in high-end
   ;   vector illustration programs (e.g. Adobe Illustrator)
   ;   and in PostScript fonts, as they are more flexible
   ;   than quadratic curves, but they are harder to render
   ;
   ; Cubic Bézier curves are so difficult to render that this
   ; class actually does not render them directly; instead,
   ; it converts cubic curves to a set of very small straight
   ; line segments, and renders from that. Usually this is
   ; not noticeable, but you should be aware of this if you
   ; use this class in unusual contexts.
   ;
   ; This class is not derived from anything as it does not
   ; have any user-exposed parameters. Your code must provide
   ; all the parameters you require. You must call Rasterize()
   ; before making any queries about the curve, or you will
   ; get erroneous results. Changing curve data after it has
   ; been rasterized will not affect any queries.
   ;
   ; This code is adapted from DeluxeClipping, which I wrote
   ; for Janet Parke's course on Masking with Ultra Fractal.
   ; The cubic Bézier support is new for UF5. The additional
   ; query modes are also new.
   
 public:
   import "common.ulb"
 ;  $define debug
   
   func DMJ_PolyCurve()
     m_WindingMode = 0
     m_Roots = new ComplexArray(3)
   endfunc
 
   ; set the winding mode
   ; pmode = winding mode:
   ;  0 = fully self-intersecting (default) (like "invert overlaps" in DeluxeClipping)
   ;  1 = individual loops not self-intersecting, subsequent loops may intersect
   ;  2 = entire curve not self-intersecting (like "all inside" in DeluxeClipping)
   func SetWindingMode(int pmode)
     m_WindingMode = pmode
   endfunc
 
   ; set the endpoint extension mode
   ; pmode = extension mode:
   ;  0 = none (cap endpoints)
   ;  1 = extend lines
   ;  2 = extend curves (not implemented)
   func SetExtensionMode(int pmode)
     m_EndpointExtensionMode = pmode
   endfunc
 
   ; set the number of curve segments
   ; you must call this prior to setting any curve data
   func SetSegmentCount(int ppoints)
     setLength(m_Anchors, ppoints)
     setLength(m_Control1, ppoints)
     setLength(m_Control2, ppoints)
     setLength(m_Lengths, ppoints)
     setLength(m_Types, ppoints)
     setLength(m_Auto, ppoints)
     m_AnchorNext = 0
     m_AnchorHighest = -1
   endfunc
   
   ; set parameters for a segment
   ; note: do not manually specify a close point; use the
   ; auto-close flag instead. Otherwise, your shape may not
   ; render correctly (it will "leak").
   ; ppoint = slot number (use -1 for next slot)
   ; ptype = segment type: 0 = none (close), 1 = linear, 2 = quadratic (use pcontrol1), 3 = cubic (use pcontrol1 and 2)
   ; pauto = auto-compute controls: 1 = close, 2 = smooth, 4 = curvature smooth
   func SetSegment(int ppoint, int ptype, int pauto, complex panchor, complex pcontrol1, complex pcontrol2)
     if (ppoint == -1)
       ppoint = m_AnchorNext
       m_AnchorNext = m_AnchorNext + 1
     endif
     if (ppoint > m_AnchorHighest)
       m_AnchorHighest = ppoint
     endif
 
     print("segment ", ppoint, ": type ", ptype, " auto ", pauto, " ", panchor, " ", pcontrol1, " ", pcontrol2)
 
     m_Anchors[ppoint] = panchor
     m_Control1[ppoint] = pcontrol1
     m_Control2[ppoint] = pcontrol2
     m_Types[ppoint] = ptype
     m_Auto[ppoint] = pauto
   endfunc
   
   ; rasterize a set of curve segments that may contain
   ; cubic segments into lines and quadratic segments;
   ; also performs auto-compute at this time
   ; pflags = extra computation flags: 1 = compute segment lengths (needed for ClosestPoint, DistanceAlongCurve), 2 = compute normals (needed for ClosestCurve)
   func Rasterize(int pflags)
     ; count number of straight/quadratic segments we will need
     ; and perform auto computations
     int j = 0
     int k = 0
     int l = 0
     int m = 0
     int n = 0
     int segmentcount = 0
     int quadraticcount = 0
     int loopstart = 0
     int loopnext
     int loopprevious
     int looprealstart
     bool computenormals = (pflags % 4 >= 2)
     bool computelengths = (pflags % 2 == 1)
     float s
     float t
     float totallength = 0
     complex point
     complex point2
     print("Rasterize() trace")
     
     while (j <= m_AnchorHighest)
       print("segment ", j)
       ; determine next loop segment (auto-close handled here)
       loopnext = j + 1
       if (m_Auto[j] % 2 == 1)
         loopnext = loopstart
         segmentcount = segmentcount + 1
         print("auto-close to segment ", loopnext)
       endif
 
       ; count the segment
       if (m_Types[j] < 2)
         segmentcount = segmentcount + 1
       elseif (m_Types[j] == 2)
         segmentcount = segmentcount + 1
         quadraticcount = quadraticcount + 1
         if (m_Auto[j] % 4 >= 2)
           loopprevious = j - 1
           if (loopprevious < loopstart)  ; auto-smooth at the start of the curve; find the end
             k = j
             while (k < length(m_Anchors) && loopprevious < j)
               if (m_Auto[k] % 2 == 1 || m_Types[k] == 0)
                 loopprevious = k
               endif
               k = k + 1
             endwhile
             if (m_Types[loopprevious] < 2)
               m_Anchors[j] = m_Anchors[loopprevious]
             else
               m_Anchors[j] = (m_Control1[j] + m_Control1[loopprevious]) / 2
             endif
           else
             m_Anchors[j] = (m_Control1[j] + m_Control1[loopprevious]) / 2
           endif
           print("smoothed to ", loopprevious, " ", m_Anchors[j])
         endif
       else
         segmentcount = segmentcount + 1024  ; **** this number should be dynamic
         ; **** perform auto here
       endif
 
       if (m_Auto[j] % 2 == 1 || m_Types[j] == 0)
         loopstart = j + 1
       endif
 
       j = j + 1
     endwhile
     print("total expanded segments: ", segmentcount)
     
     ; create updated arrays and copy over segments
     setLength(m_RealIsStartPoint, segmentcount)
     setLength(m_RealIsEndPoint, segmentcount)
     setLength(m_RealAnchors, segmentcount)
     setLength(m_RealControls, segmentcount)
     setLength(m_RealAnchorNormals, segmentcount)
     setLength(m_RealControlNormals, segmentcount)
     setLength(m_RealExitNormals, segmentcount)
     setLength(m_RealNormals, segmentcount)
     setLength(m_RealLengths, segmentcount)
     setLength(m_RealSummedLengths, segmentcount)
     setLength(m_RealTypes, segmentcount)
     setLength(m_RealA, segmentcount)
     setLength(m_RealB, segmentcount)
 
     if (computelengths)
       setLength(m_RealSubLengths, quadraticcount*256)
       setLength(m_RealIndex, segmentcount)
     endif
 
     j = 0
     k = 0
     l = 0
     loopstart = 0
     loopnext = 0
     looprealstart = 0
     while (j <= m_AnchorHighest)
       ; set up endpoint
       if (looprealstart == k)
         m_RealIsStartPoint[k] = true
       else
         m_RealIsStartPoint[k] = false
       endif
       m_RealIsEndPoint[k] = false
     
       ; determine next loop segment (auto-close handled here)
       loopnext = j + 1
       if (m_Auto[j] % 2 == 1 || m_Types[j] == 0)
         loopnext = loopstart
         loopstart = j + 1
       endif
       
       ; copy over segment
       if (m_Types[j] == 0)      ; explicit point (normally these are automatic)
         print("close point ", m_Anchors[j])
         m_RealIsEndPoint[k] = true  ; this marks an endpoint
         m_RealTypes[k] = 0
         m_RealAnchors[k] = m_Anchors[j]
         m_RealSummedLengths[k] = totallength
         m_RealLengths[k] = 0
         k = k + 1
         looprealstart = k      ; next segment will be a starting point
       else
         if (m_Types[j] == 1)
           print("line segment ", m_Anchors[j], " ", m_Anchors[loopnext])
           m_RealTypes[k] = 1
           m_RealAnchors[k] = m_Anchors[j]
           ; precompute length squared used to determine distance
           m_RealLengths[k] = |m_Anchors[loopnext] - m_Anchors[j]|
           m_RealSummedLengths[k] = totallength
           totallength = totallength + m_RealLengths[k]
 ;          print("length: ", m_RealLengths[k])
 ;          print("summed: ", m_RealSummedLengths[k])
 ;          print("total: ", totallength)
           k = k + 1
         elseif (m_Types[j] == 2)
           print("quadratic segment ", m_Anchors[j], " ", m_Control1[j], " ", m_Anchors[loopnext])
           m_RealTypes[k] = 2
           m_RealAnchors[k] = m_Anchors[j]
           m_RealControls[k] = m_Control1[j]
           ; precompute values used to determine distance
           m_RealA[k] = m_Anchors[j] - 2*m_Control1[j] + m_Anchors[loopnext]
           m_RealB[k] = -2*m_Anchors[j]+2*m_Control1[j]
           t = |m_RealA[k]|
           if (t < 1e-25)      ; curve segment is too straight; call it a line! faster and avoids divide-by-zero in ClosestPoint()
             print("straight quadratic (m = ", t, "), replacing with line")
             m_RealTypes[k] = 1
             m_RealLengths[k] = |m_Anchors[loopnext] - m_Anchors[j]|
             m_RealSummedLengths[k] = totallength
             totallength = totallength + m_RealLengths[k]
           else
             if (computelengths)
               m_RealIndex[k] = l*256
               m = m_RealIndex[k]
               n = m + 256
               s = 0
               t = 0
               point = m_Anchors[j]
 ;              $ifdef debug
 ;                print("filling slot ", l, " from ", m, " to ", n)
 ;              $endif
               while (m < n)
                 t = t + 0.00390625      ; 1/256
                 point2 = InterpolateQuadratic(t, m_Anchors[j], m_Control1[j], m_Anchors[loopnext])
                 m_RealSubLengths[m] = cabs( point2-point )
                 point = point2
                 s = s + m_RealSubLengths[m]
 ;                $ifdef debug
 ;                  if (m % 16 == 0)
 ;                    print("sublength: ", m_RealSubLengths[m])
 ;                    print("segment total: ", s)
 ;                  endif
 ;                $endif
                 m = m + 1
               endwhile
               l = l + 1
               m_RealLengths[k] = s      ; total accumulated length
               m_RealSummedLengths[k] = totallength
               totallength = totallength + m_RealLengths[k]
 ;              print("length: ", m_RealLengths[k])
 ;              print("summed: ", m_RealSummedLengths[k])
 ;              print("total: ", totallength)
             endif
           endif
           k = k + 1
         elseif (m_Types[j] == 3)
           print("cubic segment ", m_Anchors[j], " ", m_Control1[j], " ", m_Control2[j], " ", m_Anchors[loopnext])
           ; **** rasterize cubic curve here
         endif
         if (loopstart > j)      ; this segment closed a loop
           m_RealIsStartPoint[looprealstart] = false  ; make sure start/end points are cleared of endpoint status
           m_RealIsEndPoint[k] = false
           m_RealTypes[k] = 0
           m_RealAnchors[k] = m_Anchors[loopnext]
           k = k + 1
           looprealstart = k    ; next segment will be a starting point
         endif
       endif
 
       j = j + 1
     endwhile
     
     m_RealAnchorCount = k
     
     if (computenormals)
       ComputeNormals()
     endif
   endfunc
 
   ; Compute normals for a curve.
   ; For each segment, we compute the normal at the beginning, the
   ; control point, and the end. We then interpolate a "point normal"
   ; for each endpoint. Since computing normals is expensive, we don't
   ; automatically compute normals for every curve.
   func ComputeNormals()
     int segmentcount = m_RealAnchorCount
     int j = 0
     complex exitnormal = 0
     while (j < segmentcount)
       if (m_RealTypes[j] == 0)
         m_RealAnchorNormals[j] = exitnormal
         m_RealControlNormals[j] = exitnormal
         m_RealExitNormals[j] = exitnormal
         m_RealNormals[j] = exitnormal
       elseif (m_RealTypes[j] == 1)
         m_RealAnchorNormals[j] = m_RealAnchors[j+1] - m_RealAnchors[j]
         m_RealAnchorNormals[j] = m_RealAnchorNormals[j] / cabs(m_RealAnchorNormals[j])
         exitnormal = m_RealAnchorNormals[j]
         m_RealControlNormals[j] = exitnormal
         m_RealExitNormals[j] = exitnormal
       elseif (m_RealTypes[j] == 2)
         m_RealAnchorNormals[j] = m_RealControls[j] - m_RealAnchors[j]
         m_RealAnchorNormals[j] = m_RealAnchorNormals[j] / cabs(m_RealAnchorNormals[j])
         m_RealControlNormals[j] = m_RealAnchors[j+1] - m_RealAnchors[j]
         m_RealControlNormals[j] = m_RealControlNormals[j] / cabs(m_RealControlNormals[j])
         m_RealExitNormals[j] = m_RealAnchors[j+1] - m_RealControls[j]
         m_RealExitNormals[j] = m_RealExitNormals[j] / cabs(m_RealExitNormals[j])
         exitnormal = m_RealExitNormals[j]
       endif
       if (m_RealTypes[j] > 0)
         if (j > 0 && m_RealTypes[j-1] > 0)  ; not the first point, and previous segment has exit normal (not a point)
           m_RealNormals[j] = m_RealAnchorNormals[j] + m_RealExitNormals[j-1]
           m_RealNormals[j] = m_RealNormals[j] / cabs(m_RealNormals[j])
         else                ; first point in a sequence; use starting normal for this segment
           m_RealNormals[j] = m_RealAnchorNormals[j]
         endif
       endif
       j = j + 1
     endwhile
   endfunc
 
   ; query: is a point to the left of a line segment?
   ; positive value = no, negative value = yes, 0 = on line
   static float func IsLeftLine(complex pz, complex panchor1, complex panchor2)
     return -( (real(panchor2) - real(panchor1)) * (imag(pz) - imag(panchor1)) - \
           (real(pz) - real(panchor1)) * (imag(panchor2) - imag(panchor1)) )
   endfunc
 
   ; query: is a point to the left of a parabolic segment?
   ; note that if the point is outside the triangle enclosing
   ; the segment, we treat this as the same as if the point
   ; is to the left of the lines bounding the parabolic arc
   ; (we extend the arc by straight lines)
   static float func IsLeftQuadratic(complex pz, complex panchor1, complex pcontrol, complex panchor2)
     ; see if it's within the triangle
     float isleft1 = IsLeftLine(pz, panchor1, panchor2)
     float isleft2 = IsLeftLine(pz, panchor2, pcontrol)
     float isleft3 = IsLeftLine(pz, pcontrol, panchor1)
 
     if ((isleft1 <= 0 && isleft2 < 0 && isleft3 < 0) || (isleft1 >= 0 && isleft2 > 0 && isleft3 > 0))
       ; point is inside the triangle; determine if/where the ray
       ; intersects the parabola.
       return IsLeftQuadraticCore(pz, panchor1, pcontrol, panchor2)
     else
       float isleft5 = IsLeftLine(pcontrol, panchor1, panchor2)
       if (isleft5 > 0)
         if (isleft2 > 0 && isleft3 > 0)
           return -1.0
         else
           return 1.0
         endif
       elseif (isleft5 < 0)
         if (isleft2 < 0 && isleft3 < 0)
           return 1.0
         else
           return -1.0
         endif
       else
         return isleft1
       endif
     endif
   endfunc
   
   ; query: is a point to the left of a parabolic segment?
   ; this does not care whether the point is inside the enclosing triangle
   static float func IsLeftQuadraticCore(complex pz, complex panchor1, complex pcontrol, complex panchor2)
     float isleft4 = 0    ; parabola flag
     complex c        ; work variables for quadratic segments
     complex q
     complex s
     complex t
     float d
     float tx
     float ty
     float isy
 
     q = 2*conj(panchor2-panchor1) / |panchor2-panchor1|  ; rotation and scaling vector
     c = (panchor1 + panchor2) * 0.5
     t = (pz - c) * q          ; transform pixel so line segment is rotated to X-axis and centered at origin
     s = (pcontrol - c) * q        ; transform control point the same way
     isy = 1/imag(s)            ; precalc
     tx = real(t) - imag(t)*real(s)*isy  ; shear
     ty = imag(t)*isy          ; squash
 
     d = 0.5-0.5*sqr(tx)          ; height of parabola at this point
     isleft4 = ty - d          ; parabola isleft flag
     if (imag(s) < 0)
       isleft4 = -isleft4
     endif
 
     return -isleft4
   endfunc
 
   ; query: is a point to the left of a particular curve segment?
   ; this determines the type of segment automatically
   float func IsLeft(complex pz, int psegment)
     while (psegment >= 0)
       if (m_RealTypes[psegment] == 0)
         ; point; if we have a previous segment, use it
         psegment = psegment - 1
       elseif (m_RealTypes[psegment] == 1)
         return IsLeftLine(pz, m_RealAnchors[psegment], m_RealAnchors[psegment+1])
       elseif (m_RealTypes[psegment] == 2)
         return IsLeftQuadratic(pz, m_RealAnchors[psegment], m_RealControls[psegment], m_RealAnchors[psegment+1])
       endif
     endwhile
     ; no valid segment, just points; answer "no" (not a valid answer)
     return 0
   endfunc
 
   ; given a complex number pz, return whether it is inside
   ; the curve or not; this is the fastest query
   bool func IsInside(complex pz)
     int segmentcount = m_RealAnchorCount
     int j = 0        ; indices for current/next
     int k = 1
     int fullwinding = 0    ; winding numbers
     int segmentwinding = 0
     float isleft1 = 0    ; triangle flags  
     float isleft2 = 0
     float isleft3 = 0
     float isleft4 = 0    ; parabola flag
 
     $ifdef debug
       int testx = 500
       int testy = 480
       if (#x == testx && #y == testy)
         print("IsInside() trace")
         print("total anchor count: ", m_RealAnchorCount)
         print("test point: ", testx, " ", testy, " = ", pz)
       endif
     $endif
     
     while (j < segmentcount)
       if (m_RealTypes[j] == 0)
         ; close point; deal with winding numbers
         if (m_WindingMode == 0)
           fullwinding = fullwinding + segmentwinding
         elseif (m_WindingMode == 1)
           if (segmentwinding != 0)
             fullwinding = fullwinding + 1
           endif
         else
           if (segmentwinding != 0)
             fullwinding = 1
           endif
         endif
         $ifdef debug
           if (#x == testx && #y == testy)
             print("segment ", j, " segmentwinding ", segmentwinding, " fullwinding ", fullwinding)
           endif
         $endif
         segmentwinding = 0
 
       elseif (m_RealTypes[j] == 1)
         ; line segment; update winding number based solely on line segment
 
         ; The winding number method counts all line segments
         ; that cross the vertical position of the pixel to
         ; the RIGHT of it. Segments that cross up increment
         ; the count, segments that cross down decrement it.
         ; With a closed shape, this will equal zero if the
         ; point is not inside (regardless of the clockwise/
         ; counter-clockwise direction of points). Very
         ; clever algorithm. Google it.
         
         if (imag(m_RealAnchors[j]) <= imag(pz))  ; line segment imag <= point imag
           if (imag(m_RealAnchors[k]) > imag(pz))  ; upward crossing
             $ifdef debug
               if (#x == testx && #y == testy)
                 print("segment ", j, " upward crossing")
               endif
             $endif
               
             if ( (real(m_RealAnchors[k]) - real(m_RealAnchors[j])) * (imag(pz) - imag(m_RealAnchors[j])) - \
                (real(pz) - real(m_RealAnchors[j])) * (imag(m_RealAnchors[k]) - imag(m_RealAnchors[j])) > 0 )
               segmentwinding = segmentwinding + 1
               $ifdef debug
                 if (#x == testx && #y == testy)
                   print("winding + 1")
                 endif
               $endif
             endif
           endif
         else
           if (imag(m_RealAnchors[k]) <= imag(pz))  ; downward crossing
             $ifdef debug
               if (#x == testx && #y == testy)
                 print("segment ", j, " downward crossing")
               endif
             $endif
             if ( (real(m_RealAnchors[k]) - real(m_RealAnchors[j])) * (imag(pz) - imag(m_RealAnchors[j])) - \
                (real(pz) - real(m_RealAnchors[j])) * (imag(m_RealAnchors[k]) - imag(m_RealAnchors[j])) < 0 )
               segmentwinding = segmentwinding - 1
               $ifdef debug
                 if (#x == testx && #y == testy)
                   print("winding - 1")
                 endif
               $endif
             endif
           endif
         endif
         
       elseif (m_RealTypes[j] == 2)
         ; quadratic segment; update winding number based on a parabola
 
         ; Just so you know, I worked this out myself. No doubt there's
         ; source code I could have ripped into a formula, but I wanted
         ; the satisfaction of knowing I could do it. I am quite sure
         ; this is not the most optimal method, but it does work. Pay
         ; particular attention to the use of < vs. <= as if you make a
         ; mistake, you will have horizontal line segment errors in the
         ; drawn shape.
 
         ; This is essentially the same algorithm as arbitrary polygon,
         ; but extended so that each segment can be a quadratic Bézier
         ; curve rather than just a straight line. Each curve is
         ; described by a triangle, and a parabola inscribed within the
         ; triangle such that two sides of the triangle are tangents
         ; to the parabola. Points that do not lie within the triangle
         ; are treated the same as the arbitrary polygon (only the chord
         ; cutting across the parabola, the third side of the triangle,
         ; matters) but for points inside the triangle a determination
         ; must be made as to whether the point is to the left of the
         ; parabola or not. There are a few edge cases where the
         ; parabola loops back on itself and the winding number can be
         ; changed multiple times for each segment.
 
         ; three cases: pixel is below triangle containing curve, above it, or between top and bottom of it
         ; if either of the first two, don't examine this curve further (it is irrelevant to the winding number)
         if ((imag(pz) >= imag(m_RealAnchors[j]) && \
            (imag(pz) < imag(m_RealAnchors[k]) || imag(pz) < imag(m_RealControls[j]))) || \
           (imag(pz) < imag(m_RealAnchors[j]) && \
            (imag(pz) >= imag(m_RealAnchors[k]) || imag(pz) >= imag(m_RealControls[j]))))
           $ifdef debug
             if (#x == testx && #y == testy)
               print(j, " considered for curve segment")
             endif
           $endif
 
           ; pixel is between top and bottom of curve
           ; see if it's within the triangle
           isleft1 = (real(m_RealAnchors[k]) - real(m_RealAnchors[j])) * (imag(pz) - imag(m_RealAnchors[j])) - \
                 (real(pz) - real(m_RealAnchors[j])) * (imag(m_RealAnchors[k]) - imag(m_RealAnchors[j]))
           isleft2 = (real(m_RealControls[j]) - real(m_RealAnchors[k])) * (imag(pz) - imag(m_RealAnchors[k])) - \
                 (real(pz) - real(m_RealAnchors[k])) * (imag(m_RealControls[j]) - imag(m_RealAnchors[k]))
           isleft3 = (real(m_RealAnchors[j]) - real(m_RealControls[j])) * (imag(pz) - imag(m_RealControls[j])) - \
                 (real(pz) - real(m_RealControls[j])) * (imag(m_RealAnchors[j]) - imag(m_RealControls[j]))
 
           if ((isleft1 <= 0 && isleft2 < 0 && isleft3 < 0) || (isleft1 >= 0 && isleft2 > 0 && isleft3 > 0))
             ; point is inside the triangle; determine if/where the ray
             ; intersects the parabola. start by normalizing the parabola
             ; and the ray
             $ifdef debug          
               if (#x == testx && #y == testy)
                 print(j, " considered inside triangle")
               endif
             $endif
 
             isleft4 = -IsLeftQuadraticCore(pz, m_RealAnchors[j], m_RealControls[j], m_RealAnchors[k])
             
             if (imag(m_RealAnchors[k]) > imag(m_RealAnchors[j]))  ; upward crossing
               if (imag(pz) >= imag(m_RealAnchors[j]) && \
                 imag(pz) < imag(m_RealAnchors[k]))        ; within line segment
                 if (isleft4 > 0)
                   segmentwinding = segmentwinding + 1      ; confirmed upward crossing
                 endif
               else                        ; outside line segment
                 if (isleft1 < 0)                ; loop back occurs to the left
                   if (isleft4 > 0)
                     segmentwinding = segmentwinding + 1
                   endif
                 else                      ; loop back occurs to the right
                   if (isleft4 < 0)
                     segmentwinding = segmentwinding - 1    ; loop back; downward crossing
                   endif
                 endif
               endif
 
             else                          ; downward crossing
               if (imag(pz) < imag(m_RealAnchors[j]) && \
                 imag(pz) >= imag(m_RealAnchors[k]))        ; within line segment
                 if (isleft4 < 0)
                   segmentwinding = segmentwinding - 1      ; confirmed downward crossing
                 endif
               else                        ; outside line segment
                 if (isleft1 > 0)                ; loop back occurs to the left
                   if (isleft4 < 0)
                     segmentwinding = segmentwinding - 1
                   endif
                 else                      ; loop back occurs to the right
                   if (isleft4 > 0)
                     segmentwinding = segmentwinding + 1    ; loop back; upward crossing
                   endif
                 endif
               endif
 
             endif    ; upward/downward
 
           else
             $ifdef debug
               if (#x == testx && #y == testy)
                 print(j, " tested against line segment only")
               endif
             $endif
 
             ; point is outside the triangle; all that matters is its relation to the line segment
             if (imag(m_RealAnchors[j]) <= imag(pz))          ; line segment imag < point imag
               if (imag(m_RealAnchors[k]) > imag(pz))        ; upward crossing
                 if ( isleft1 > 0 )
                   segmentwinding = segmentwinding + 1
                 endif
               endif
             else
               if (imag(m_RealAnchors[k]) <= imag(pz))        ; downward crossing
                 if ( isleft1 < 0 )
                   segmentwinding = segmentwinding - 1
                 endif
               endif
             endif
 
           endif
         endif
       
       endif
 
       j = j + 1
       k = k + 1
     endwhile
 
     $ifdef debug
       if (#x == testx && #y == testy)
         print("final winding ", fullwinding)
       endif
     $endif
     
     ; set solid flag with results
     return (abs(fullwinding) % 2 != 0)
   endfunc
 
   ; given a complex number pz, locate the closest
   ; point on the curve, the segment number it is on,
   ; and the distance squared to that point; note
   ; that distance can never be negative (use IsLeft()
   ; to determine sign if you need it)
   complex func ClosestPoint(complex pz, float &psegment, float &psegmentoffset, float &pdistancesquared)
     int segmentcount = m_RealAnchorCount
     int j = 0        ; indices for current/next
     int k = 1
     float t
     float t2
     complex c        ; third term of quadratic curve function (x and y parts)
     float m          ; coefficients of distance-squared function (a cubic polynomial)
     float n
     float r
     float s
     int rootcount      ; count of roots to cubic equation
     complex point
     complex pointmin = m_RealAnchors[0]  ; closest point on curve; start with first point
     float ds = 0
     float dsmin = |pz-m_RealAnchors[0]|  ; closest distance squared; start with simple distance to first point
     float segment
     float segmentoffset
     float segmentmin = 0
     float segmentminoffset = 0
     complex vl
     complex vp
     float vls
     float vps
     
     $ifdef debug
       int testx = 700
       int testy = 400
       if (#x == testx && #y == testy)
         print("ClosestPoint() trace")
         print("total anchor count: ", m_RealAnchorCount)
         print("test point: ", testx, " ", testy, " = ", pz)
         print("starting dsmin: ", dsmin)
       endif
     $endif
     
     while (j < segmentcount)
       if (m_RealTypes[j] == 1)
         ; compute distance from line segment to point
         ; start by finding closest point along the line
         vl = m_RealAnchors[k] - m_RealAnchors[j]  ; vector for line
         vp = pz - m_RealAnchors[j]          ; vector from line start to point
         t = real(vl)*real(vp) + imag(vl)*imag(vp)  ; dot product, |vl| |vp| cos theta
         t2 = t / m_RealLengths[j]          ; normalize = cos theta * |vp|/|vl|
 
         $ifdef debug
           if (#x == testx && #y == testy)
             print("segment ", j, " type 1 (line) ", m_RealAnchors[j], " ", m_RealAnchors[k])
             print("start: ", m_RealIsStartPoint[j], ", end: ", m_RealIsEndPoint[k])
             print("length: ", m_RealLengths[j])
             print("line vector: ", vl)
             print("point vector: ", vp)
             print("dot product: ", t)
             print("normalized: ", t2)
           endif
         $endif
 
         ; cap position along line if we're not extending it in either direction
         if (t2 <= 0 && (m_EndpointExtensionMode == 0 || !m_RealIsStartPoint[j]))
           point = m_RealAnchors[j]
           ds = |pz-point|
           segment = j
         elseif (t2 >= 1 && (m_EndpointExtensionMode == 0 || !m_RealIsEndPoint[k]))
           point = m_RealAnchors[k]
           ds = |pz-point|
           segment = k
         else
           point = m_RealAnchors[j] + vl*t2
           ds = |pz-point|
           if (t2 < 0 || t2 > 1)  ; this result is not really "on" the curve; log it as a segment + offset
             segment = j
             segmentoffset = t2
           else
             segment = j+t2
             segmentoffset = 0
           endif
         endif
 
         $ifdef debug
           if (#x == testx && #y == testy)
             print("ds: ", ds)
             print("point: ", point)
           endif
         $endif
         
         ; update closest point, if this one is closer
         if (ds < dsmin)
           dsmin = ds
           pointmin = point
           segmentmin = segment
           segmentminoffset = segmentoffset
 
           $ifdef debug
             if (#x == testx && #y == testy)
               print("updated: dsmin = ", dsmin, ", pointmin = ", pointmin, ", segmentmin = ", segmentmin, ", segmentminoffset = ", segmentminoffset)
             endif
           $endif
         endif
 
       elseif (m_RealTypes[j] == 2)
         ; compute distance from parametric parabola to point
         ; we define D(t) to be the distance from P(t) to pz
         ; and look for minimum values of D(t) by taking its
         ; derivative D'(t) and solving that; this gives us
         ; local minima and maxima, so we simply compute the
         ; distances at these points and choose the lowest
         ;
         ; this was a pain for me to work out, mainly because
         ; of a stupid mistake in my equation scribbling that
         ; I kept making over and over again...
 
         c = m_RealAnchors[j] - pz        ; remaining coefficients (depend on point we're finding distance to)
 
         m = |m_RealA[j]|            ; get coefficients to cubic equation
         n = real(m_RealA[j])*real(m_RealB[j]) + imag(m_RealA[j])*imag(m_RealB[j])
         r = 2 * ( real(m_RealA[j])*real(c) + imag(m_RealA[j])*imag(c) ) + |m_RealB[j]|
         s = real(m_RealB[j])*real(c) + imag(m_RealB[j])*imag(c)
         
         ; D'(t) = 4mt^3 + 6nt^2 + 2rt + 2s
         ; note that only r and s vary with c
         rootcount = Math.SolveRealCubicPolynomial(2*m, 3*n, r, s, 3, m_Roots)
         
         $ifdef debug
           if (#x == testx && #y == testy)
             print("segment ", j, " type 2 (quadratic) ", m_RealAnchors[j], " ", m_RealControls[j], " ", m_RealAnchors[k])
             print("start: ", m_RealIsStartPoint[j], ", end: ", m_RealIsEndPoint[k])
             print("a: ", m_RealA[j])
             print("b: ", m_RealB[j])
             print("c: ", c)
             print("m: ", m)
             print("n: ", n)
             print("r: ", r)
             print("s: ", s)
             print("rootcount: ", rootcount)
             print("root 0: ", m_Roots.m_Elements[0])
             print("root 1: ", m_Roots.m_Elements[1])
             print("root 2: ", m_Roots.m_Elements[2])
           endif
         $endif
         
         while (rootcount > 0)
           rootcount = rootcount - 1
           
           ; examine this root
           t = real(m_Roots.m_Elements[rootcount])
           
           ; determine point at t
           ; we may need to compute the curve at the t value if it's within range
           ; also, if we are extending beyond endpoints, we may need to recompute
           ; t as well
           if (t <= 0)
             if (m_EndpointExtensionMode == 0 || !m_RealIsStartPoint[j])
               point = m_RealAnchors[j]
               t = 0
             else
               vl = m_RealControls[j] - m_RealAnchors[j]  ; vector for line
               vp = pz - m_RealAnchors[j]          ; vector from line start to point
               t2 = real(vl)*real(vp) + imag(vl)*imag(vp)  ; dot product, |vl| |vp| cos theta
               vls = |vl|
               vps = cabs(vp)
               t = t2 / vls                ; normalize = cos theta
               point = m_RealAnchors[j] + vl*t
               t = t * sqrt(vls)
               $ifdef debug
                 if (#x == testx && #y == testy)
                   print("extended start line")
                   print("line vector: ", vl)
                   print("point vector: ", vp)
                   print("line vector length sqd: ", vls)
                   print("point vector length sqd: ", vps)
                   print("dot product: ", t2)
                   print("normalized: ", t)
                   print("point: ", point)
                 endif
               $endif
             endif
           elseif (t >= 1)
             if (m_EndpointExtensionMode == 0 || !m_RealIsEndPoint[k])
               point = m_RealAnchors[k]
               t = 1
             else
               vl = m_RealAnchors[k] - m_RealControls[j]  ; vector for line
               vp = pz - m_RealControls[j]          ; vector from line start to point
               t2 = real(vl)*real(vp) + imag(vl)*imag(vp)  ; dot product, |vl| |vp| cos theta
               vls = |vl|
               vps = cabs(vp)
               t = t2 / vls                ; normalize = cos theta
               point = m_RealControls[j] + vl*t
               t2 = (t-1) * sqrt(vls)
               $ifdef debug
                 if (#x == testx && #y == testy)
                   print("extended end line")
                   print("line vector: ", vl)
                   print("point vector: ", vp)
                   print("line vector length sqd: ", vls)
                   print("point vector length sqd: ", vps)
                   print("dot product: ", t2)
                   print("normalized: ", t)
                   print("point: ", point)
                 endif
               $endif
             endif
           else
             point = InterpolateQuadratic(t, m_RealAnchors[j], m_RealControls[j], m_RealAnchors[k])
           endif
           
           ; update closest point
           ds = |pz-point|
           $ifdef debug
             if (#x == testx && #y == testy)
               print("ds: ", ds)
             endif
           $endif
           if (ds < dsmin)
             dsmin = ds
             pointmin = point
             if (t < 0)
               segmentmin = j
               segmentminoffset = t
             elseif (t > 1)
               segmentmin = j+1
               segmentminoffset = t2
             else
               segmentmin = j+t
               segmentminoffset = 0
             endif
           endif
         endwhile
 
       endif  ; segment type
       
       j = j + 1
       k = k + 1
     endwhile
     
     $ifdef debug
       if (#x == testx && #y == testy)
         print("result t: ", segmentmin)
         print("result to: ", segmentminoffset)
         print("result B(t): ", pointmin)
         print("result ds: ", dsmin)
       endif
     $endif
     psegment = segmentmin
     psegmentoffset = segmentminoffset
     pdistancesquared = dsmin
     return pointmin
   endfunc
 
   ; Given a complex number pz, determine the closest
   ; similar curve to it. That is, assume the endpoints
   ; and control point can slide along their respective
   ; normals in unison; find the appropriate scaling
   ; factor that will allow them to create a curve that
   ; passes through our test point. This is useful for
   ; creating smooth mappings around curves without the
   ; discontinuities on the concave side that using
   ; ClosestPoint() would.
   complex func ClosestCurve(complex pz, float &psegment, float &pdistancesquared)
     return 0
   endfunc
 
   ; given a fractional segment number, determine the
   ; distance from the curve start to that point
   ; if ptotal is set, the distance will be cumulative
   ; from all closed loops, not just from the start of
   ; the closest closed loop
   float func DistanceAlongCurve(float psegment, bool ptotal)
     int segment = floor(psegment)
     float s = 0
     float t = psegment - segment
     int m
     int n
     
     ; to start with, use all the summed lengths of previous segments
     s = m_RealSummedLengths[segment]
     
     ; now determine a partial length within this segment
     if (t > 0)
       if (m_RealTypes[segment] == 1)
         ; linear segment
         s = s + t*m_RealLengths[segment]
         
       elseif (m_RealTypes[segment] == 2)
         ; quadratic segment
         ; determine number of sub-segments to use
         t = t * 256
         m = m_RealIndex[segment]
         n = m + floor(t)
         t = t - floor(t)
         while (m < n)
           s = s + m_RealSubLengths[m]
           m = m + 1
         endwhile
         if (t > 0)
           s = s + t*m_RealSubLengths[m]
         endif
         
       endif
     endif
     
     ; if we're not using the total distance, find the loop start
     ; and subtract all the distance accumulated before it
     
     return s
   endfunc
 
   ; given an interpolant and two endpoints, compute a linear curve
   static complex func InterpolateLinear(float pt, complex panchor1, complex panchor2)
     return panchor1 + pt * (panchor2-panchor1)
   endfunc
   
   ; given an interpolant, two endpoints, and a control point, compute a parabolic curve
   static complex func InterpolateQuadratic(float pt, complex panchor1, complex pcontrol, complex panchor2)
     return sqr(1-pt)*panchor1 + 2*pt*(1-pt)*pcontrol + sqr(pt)*panchor2
   endfunc
 
   ; given an interpolant, two endpoints, and two control points, compute a cubic curve
   static complex func InterpolateCubic(float pt, complex panchor1, complex pcontrol1, complex pcontrol2, complex panchor2)
     return (1-pt)^3*panchor1 + 3*pt*(1-pt)^2*pcontrol1 + 3*pt^2*(1-pt)*pcontrol2 + pt^3*panchor2
   endfunc
 
   ; public variables
   ; you can set these yourself, but if you create invalid
   ; curve data it may not render correctly
   complex m_Anchors[]
   complex m_Control1[]
   complex m_Control2[]
   float m_Lengths[]
   int m_Types[]
   int m_Auto[]
   
 protected:
   int m_WindingMode
   int m_EndpointExtensionMode
   int m_AnchorNext
   int m_AnchorHighest
 
   bool m_RealIsStartPoint[]  
   bool m_RealIsEndPoint[]
   complex m_RealAnchors[]
   complex m_RealControls[]
   complex m_RealAnchorNormals[]
   complex m_RealControlNormals[]
   complex m_RealExitNormals[]
   complex m_RealNormals[]
   float m_RealLengths[]
   float m_RealSummedLengths[]
   float m_RealSubLengths[]
   int m_RealIndex[]
   int m_RealTypes[]
   complex m_RealA[]
   complex m_RealB[]
   int m_RealAnchorCount
 
   ComplexArray m_Roots
 }
 


Constructor Summary
DMJ_PolyCurve()
          $define debug
 
Method Summary
 complex ClosestCurve(complex pz, float psegment, float pdistancesquared)
          Given a complex number pz, determine the closest similar curve to it.
 complex ClosestPoint(complex pz, float psegment, float psegmentoffset, float pdistancesquared)
          given a complex number pz, locate the closest point on the curve, the segment number it is on, and the distance squared to that point; note that distance can never be negative (use IsLeft() to determine sign if you need it)
 void ComputeNormals()
          Compute normals for a curve.
 float DistanceAlongCurve(float psegment, boolean ptotal)
          given a fractional segment number, determine the distance from the curve start to that point if ptotal is set, the distance will be cumulative from all closed loops, not just from the start of the closest closed loop
static complex InterpolateCubic(float pt, complex panchor1, complex pcontrol1, complex pcontrol2, complex panchor2)
          given an interpolant, two endpoints, and two control points, compute a cubic curve
static complex InterpolateLinear(float pt, complex panchor1, complex panchor2)
          given an interpolant and two endpoints, compute a linear curve
static complex InterpolateQuadratic(float pt, complex panchor1, complex pcontrol, complex panchor2)
          given an interpolant, two endpoints, and a control point, compute a parabolic curve
 boolean IsInside(complex pz)
          given a complex number pz, return whether it is inside the curve or not; this is the fastest query
 float IsLeft(complex pz, int psegment)
          query: is a point to the left of a particular curve segment? this determines the type of segment automatically
static float IsLeftLine(complex pz, complex panchor1, complex panchor2)
          query: is a point to the left of a line segment? positive value = no, negative value = yes, 0 = on line
static float IsLeftQuadratic(complex pz, complex panchor1, complex pcontrol, complex panchor2)
          query: is a point to the left of a parabolic segment? note that if the point is outside the triangle enclosing the segment, we treat this as the same as if the point is to the left of the lines bounding the parabolic arc (we extend the arc by straight lines)
static float IsLeftQuadraticCore(complex pz, complex panchor1, complex pcontrol, complex panchor2)
          query: is a point to the left of a parabolic segment? this does not care whether the point is inside the enclosing triangle
 void Rasterize(int pflags)
          rasterize a set of curve segments that may contain cubic segments into lines and quadratic segments; also performs auto-compute at this time pflags = extra computation flags: 1 = compute segment lengths (needed for ClosestPoint, DistanceAlongCurve), 2 = compute normals (needed for ClosestCurve)
 void SetExtensionMode(int pmode)
          set the endpoint extension mode pmode = extension mode: 0 = none (cap endpoints) 1 = extend lines 2 = extend curves (not implemented)
 void SetSegment(int ppoint, int ptype, int pauto, complex panchor, complex pcontrol1, complex pcontrol2)
          set parameters for a segment note: do not manually specify a close point; use the auto-close flag instead.
 void SetSegmentCount(int ppoints)
          set the number of curve segments you must call this prior to setting any curve data
 void SetWindingMode(int pmode)
          set the winding mode pmode = winding mode: 0 = fully self-intersecting (default) (like "invert overlaps" in DeluxeClipping) 1 = individual loops not self-intersecting, subsequent loops may intersect 2 = entire curve not self-intersecting (like "all inside" in DeluxeClipping)
 
Methods inherited from class Object
 

Constructor Detail

DMJ_PolyCurve

public DMJ_PolyCurve()
$define debug

Method Detail

SetWindingMode

public void SetWindingMode(int pmode)
set the winding mode pmode = winding mode: 0 = fully self-intersecting (default) (like "invert overlaps" in DeluxeClipping) 1 = individual loops not self-intersecting, subsequent loops may intersect 2 = entire curve not self-intersecting (like "all inside" in DeluxeClipping)


SetExtensionMode

public void SetExtensionMode(int pmode)
set the endpoint extension mode pmode = extension mode: 0 = none (cap endpoints) 1 = extend lines 2 = extend curves (not implemented)


SetSegmentCount

public void SetSegmentCount(int ppoints)
set the number of curve segments you must call this prior to setting any curve data


SetSegment

public void SetSegment(int ppoint,
                       int ptype,
                       int pauto,
                       complex panchor,
                       complex pcontrol1,
                       complex pcontrol2)
set parameters for a segment note: do not manually specify a close point; use the auto-close flag instead. Otherwise, your shape may not render correctly (it will "leak"). ppoint = slot number (use -1 for next slot) ptype = segment type: 0 = none (close), 1 = linear, 2 = quadratic (use pcontrol1), 3 = cubic (use pcontrol1 and 2) pauto = auto-compute controls: 1 = close, 2 = smooth, 4 = curvature smooth


Rasterize

public void Rasterize(int pflags)
rasterize a set of curve segments that may contain cubic segments into lines and quadratic segments; also performs auto-compute at this time pflags = extra computation flags: 1 = compute segment lengths (needed for ClosestPoint, DistanceAlongCurve), 2 = compute normals (needed for ClosestCurve)


ComputeNormals

public void ComputeNormals()
Compute normals for a curve. For each segment, we compute the normal at the beginning, the control point, and the end. We then interpolate a "point normal" for each endpoint. Since computing normals is expensive, we don't automatically compute normals for every curve.


IsLeftLine

public static float IsLeftLine(complex pz,
                               complex panchor1,
                               complex panchor2)
query: is a point to the left of a line segment? positive value = no, negative value = yes, 0 = on line


IsLeftQuadratic

public static float IsLeftQuadratic(complex pz,
                                    complex panchor1,
                                    complex pcontrol,
                                    complex panchor2)
query: is a point to the left of a parabolic segment? note that if the point is outside the triangle enclosing the segment, we treat this as the same as if the point is to the left of the lines bounding the parabolic arc (we extend the arc by straight lines)


IsLeftQuadraticCore

public static float IsLeftQuadraticCore(complex pz,
                                        complex panchor1,
                                        complex pcontrol,
                                        complex panchor2)
query: is a point to the left of a parabolic segment? this does not care whether the point is inside the enclosing triangle


IsLeft

public float IsLeft(complex pz,
                    int psegment)
query: is a point to the left of a particular curve segment? this determines the type of segment automatically


IsInside

public boolean IsInside(complex pz)
given a complex number pz, return whether it is inside the curve or not; this is the fastest query


ClosestPoint

public complex ClosestPoint(complex pz,
                            float psegment,
                            float psegmentoffset,
                            float pdistancesquared)
given a complex number pz, locate the closest point on the curve, the segment number it is on, and the distance squared to that point; note that distance can never be negative (use IsLeft() to determine sign if you need it)


ClosestCurve

public complex ClosestCurve(complex pz,
                            float psegment,
                            float pdistancesquared)
Given a complex number pz, determine the closest similar curve to it. That is, assume the endpoints and control point can slide along their respective normals in unison; find the appropriate scaling factor that will allow them to create a curve that passes through our test point. This is useful for creating smooth mappings around curves without the discontinuities on the concave side that using ClosestPoint() would.


DistanceAlongCurve

public float DistanceAlongCurve(float psegment,
                                boolean ptotal)
given a fractional segment number, determine the distance from the curve start to that point if ptotal is set, the distance will be cumulative from all closed loops, not just from the start of the closest closed loop


InterpolateLinear

public static complex InterpolateLinear(float pt,
                                        complex panchor1,
                                        complex panchor2)
given an interpolant and two endpoints, compute a linear curve


InterpolateQuadratic

public static complex InterpolateQuadratic(float pt,
                                           complex panchor1,
                                           complex pcontrol,
                                           complex panchor2)
given an interpolant, two endpoints, and a control point, compute a parabolic curve


InterpolateCubic

public static complex InterpolateCubic(float pt,
                                       complex panchor1,
                                       complex pcontrol1,
                                       complex pcontrol2,
                                       complex panchor2)
given an interpolant, two endpoints, and two control points, compute a cubic curve