Hermite interpolation demo

Hermite.py

#!/usr/bin/env python3

from argparse import ArgumentParser
from cmath import sin
from random import random
import sys, math

def main():
    arguments = ArgumentParser()
    arguments.add_argument('--width', '-w', default=1024, type=int, help='svg canvas width')
    arguments.add_argument('--height', '-H', default=768, type=int, help='svg canvas height')
    arguments.add_argument('--num', '-n', default=10, type=int, help='number of points')
    options = arguments.parse_args(sys.argv[1:])
    margin = 10
    width = options.width
    height = options.height
    num = options.num
    axis = []
    while len(axis) < num:
        axis.append((width - 2*margin)*random() + margin)
    axis.sort()
    for n in range(len(axis)):
        if n > 0 and axis[n] - axis[n-1] < width/(2*num): 
            axis[n] = axis[n-1] + width/(2*num)
            if axis[n] > width: axis[n] = width - margin

    points = []
    posY = height/2
    sign = 1
    while len(points) < len(axis):
        sign *= -1
        posY += 200*random()*sign
        if posY < 0: posY = margin
        if posY > height: posY = height - margin
        points.append((axis[len(points)], posY))
    n = 1
    slop = [0]
    while n + 1 < len(points):
        p0 = points[n-1]
        p1 = points[n]
        p2 = points[n+1]
        i = (p1[1] - p0[1]) / (p1[0] - p0[0])
        o = (p2[1] - p1[1]) / (p2[0] - p1[0])
        slop.append((i+o)/2)
        n += 1
    slop.append(0)
    n = 0
    svg = open('Hermite.svg', 'w+')
    svg.write('<?xml version="1.0"?>\n')
    svg.write('<svg xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" width="{0}" height="{1}" >\n'.format(width, height))
    svg.write('<g style="fill:none;stroke:red">\n')
    while n + 1 < len(points):
        p0 = points[n]
        p1 = points[n+1]
        s0 = slop[n]
        s1 = slop[n+1]
        t1 = p1[0] - p0[0]
        t2 = t1 * t1
        t3 = t2 * t1
        #  2 -2  1  1
        # -3  3 -2 -1
        #  0  0  1  0
        #  1  0  0  0
        a =  2 * p0[1]/t3 - 2 * p1[1]/t3 +     s0/t2 + s1/t2
        b = -3 * p0[1]/t2 + 3 * p1[1]/t2 - 2 * s0/t1 - s1/t1
        c =                                    s0
        d =      p0[1]
        svg.write('<path style="stroke-width:2" d="M{:.4f} {:.4f}'.format(*p0))
        total = 50
        i = 0
        while total > 0 and i <= total:
            t = t1 * i / total
            y = t * (t * (t * a + b) + c) + d
            x = t + p0[0]
            svg.write(' L{:.4f} {:.4f}'.format(x, y))
            i += 1
        svg.write('" />\n')
        n += 1
    n = 0
    while n < len(points):
        p = points[n]
        s = slop[n]
        dx = 40 * math.cos(math.atan(s))
        p0 = (p[0] - dx, s * (-dx) + p[1])
        p1 = (p[0] + dx, s * (+dx) + p[1])
        svg.write('<line style="stroke-with:1;stroke-dasharray:8,2,1,2" x1="{:.4f}" y1="{:.4f}" x2="{:.4f}" y2="{:.4f}"/>\n'.format(*p0, *p1))
        n += 1
    svg.write('</g>\n')
    n = 0
    svg.write('<path style="fill:none;stroke:DodgerBlue;stroke-width:1;stroke-dasharray:8,2" d="')
    while n < len(points):
        p = points[n]
        if n == 0: svg.write(' M{:.4f} {:.4f}'.format(*p))
        else: svg.write(' L{:.4f} {:.4f}'.format(*p))
        n += 1
    svg.write('" />\n')
    n = 0
    while n < len(points):
        svg.write('<circle cx="{0:.4f}" cy="{1:.4f}" r="4" stroke="none" fill="red"/>\n'.format(*points[n]))
        n += 1
    svg.write('</svg>')
    svg.seek(0)
    print(svg.read())
    svg.close()

if __name__ == '__main__': main()