Circle spline 2D in Nim

circlespline2d.nim

import vmath
import math

type
  ColinearSegmentError* = object of Defect
  
  Input2D = object
    dir: Vec2        # tangent vector at first input knot
    knot: seq[Vec2]
  
  CircleSegment = object
    center: Vec2
    radius: float
    startAngle: float
    endAngle: float
    startPoint: Vec2
    endPoint: Vec2
    arcLength: float
    outTangent: Vec2  # Store outgoing tangent

  CircleSpline = seq[CircleSegment]


proc wedgeProd*[T](a, b: GVec2[T]): float =
  ## 2D wedge product or exterior product
  return a.x * b.y - a.y * b.x


proc isCollinear(p1, p2, tangent: Vec2): bool =
  ## Check if points p1, p2 and tangent vector are collinear
  let chord = (p2 - p1).normalize()
  let tangentNorm = tangent.normalize()
  let crossProduct = abs(wedgeProd(chord, tangentNorm))
  crossProduct < 1e-6


proc calcCenterRadius(p1, p2, tangent: Vec2): tuple[center: Vec2, radius: float] =
  ## Calculate center and radius of circle segment given two points and tangent
  if isCollinear(p1, p2, tangent):
    raise newException(ColinearSegmentError, "Points and tangent are collinear")    
  let chord = p2 - p1
  let chordLength = chord.length
  let chordDir = chord / chordLength
  let tangentNorm = tangent.normalize 
  # Calculate angle between tangent and chord
  let cosTheta = dot(tangentNorm, chordDir)
  let theta = arccos(cosTheta)  
  # Calculate radius using the chord length and angle
  let radius = chordLength / (2.0 * sin(theta)) 
  # Calculate direction to center: rotate tangent by ±90° based on desired curve direction
  let normalDir = vec2(-tangentNorm.y, tangentNorm.x)
  let toCenter = if wedgeProd(tangentNorm, chordDir) > 0: normalDir else: -normalDir 
  # Calculate distance from p1 to center using right triangle formed by radius and half-chord
  let distToCenter = abs(radius) 
  # Calculate center point
  let center = p1 + toCenter * distToCenter 
  return (center: center, radius: abs(radius))


proc initCircleSegment(startPoint, endPoint, tangent: Vec2): CircleSegment =
  ## Initialize a circle segment from two points and tangent vector
  let (center, radius) = calcCenterRadius(startPoint, endPoint, tangent)  
  # Calculate start and end angles
  let startVec = (startPoint - center).normalize
  let endVec = (endPoint - center).normalize  
  var startAngle = arctan2(startVec.y, startVec.x)
  var endAngle = arctan2(endVec.y, endVec.x)  
  # Make sure angles follow the correct direction based on tangent
  let crossStart = wedgeProd(startVec, tangent)
  if crossStart < 0:
    while endAngle <= startAngle:
      endAngle += 2 * PI
  else:
    while endAngle >= startAngle:
      endAngle -= 2 * PI      
  let arcLength = abs(endAngle - startAngle) * radius
  # Calculate outgoing tangent at end point
  let toEnd = endPoint - center
  let outTangent = if crossStart < 0:
    vec2(toEnd.y, -toEnd.x).normalize
  else:
    vec2(-toEnd.y, toEnd.x).normalize  
  return CircleSegment(
    center: center,
    radius: radius,
    startAngle: startAngle,
    endAngle: endAngle,
    startPoint: startPoint,
    endPoint: endPoint,
    arcLength: arcLength,
    outTangent: outTangent
  )


proc initCircleSpline*(input: Input2D): CircleSpline =
  ## Generate a sequence of circle segments from input points and initial tangent
  var cspl = CircleSpline(newSeq[CircleSegment]())
  if input.knot.len < 2:
    return    
  var currentTangent = input.dir.normalize  
  # Generate segments between consecutive points
  for i in 0 ..< input.knot.len - 1:
    let seg = initCircleSegment(input.knot[i], input.knot[i+1], currentTangent)
    cspl.add(seg)
    currentTangent = seg.outTangent  # Use the calculated outgoing tangent
  return cspl


proc addCircleSegment*(spline: var CircleSpline, knot: Vec2)=
  let startKnot = spline[^1].endPoint
  let tangent = spline[^1].outTangent
  let seg = initCircleSegment(startKnot, knot, tangent)
  spline.add(seg)



proc getTotalLength*(spline: CircleSpline): float =
  ## Calculate total length of spline
  var splLen = 0.0
  for seg in spline:
    splLen += seg.arcLength
  return splLen


proc getPoint*(spline: CircleSpline, t: float): Vec2 =
  ## Get point on spline at parameter t ∈ [0,1]
  if spline.len == 0:
    return vec2(0.0, 0.0)    
  let totalLength = spline.getTotalLength()
  var remainingLength = t * totalLength  
  # Find segment containing point
  for seg in spline:
    if remainingLength <= seg.arcLength:
      let segT = remainingLength / seg.arcLength
      let angle = seg.startAngle + segT * (seg.endAngle - seg.startAngle)
      return seg.center + vec2(cos(angle), sin(angle)) * seg.radius
    remainingLength -= seg.arcLength    
  # return last point
  return spline[^1].endPoint


when is_main_module:
  import pixie
  let
    width = 1000
    height = 1000
    image = newImage(width, height)
    ctx = newContext(image)
    

    input = Input2D(
      dir: vec2(0.5, -1), # vary (small steps)
      knot: @[
        vec2( -0.0,  -0.0), 
        vec2( 50.0,  20.0),
        vec2(100.0, 100.0), 
        vec2( 90.0, 200.0), 
        vec2(100.0, 300.0), 
      ]
    )
    resolution = 500

  var spline = initCircleSpline(input)  
  spline.addCircleSegment(vec2(250.0, 350.0))
  spline.addCircleSegment(vec2(150.0, 300.0))


  image.fill(rgba(255, 255, 255, 255)) 
  ctx.fillStyle = rgba(255, 0, 0, 255) 
  let halfw = width.float / 2.0
  let halfh = height.float / 2.0
  ctx.fillCircle(circle( vec2( halfw, halfh), 4))   
  ctx.fillStyle = rgba(0, 0, 255, 255) 
  

  for i in 0 ..< resolution:
    let t = i / (resolution - 1)
    let point = getPoint(spline, t)
    let p = vec2( point.x + halfw,   - point.y + halfh) # transform for image, centr (0,0), Y = up
    ctx.fillCircle(circle(p, 2))

  image.writeFile("cspl.png")