PAEK_reverse_engineering.py

from math import *

# a performant solution would store a prefix sum of line lengths to
# a sidechannel and then use that to do a bsearch; on the GPU,
# you'd do a sum tree / histopyramid as a preprocessing step
def find_point(points, d):
    d = d
    for i in range(1, len(points)):
        x0,y0 = points[i-1]
        x1,y1 = points[i]
        dx = (x1 - x0)
        dy = (y1 - y0)
        l = (dx*dx + dy*dy)**0.5
        if d < l:
            return (i - 1), float(d) / l
        d = d - l
    return len(points) - 2, 1.0

def sample_point(points, d):
    result = find_point(points, d)
    k,t = result
    x0,y0 = points[k]
    x1,y1 = points[k+1]
    x = x0*(1-t) + x1*t
    y = y0*(1-t) + y1*t
    return x,y

def draw_points(points):
    for i in range(1, len(points)):
        x0,y0 = points[i-1]
        x1,y1 = points[i]
        line(x0,y0,x1,y1)

points = [
(674,948),
(570,838),
(586,663),
(116,714),
(22,645),
(196,574),
(108,182),
(176,50),
(428,142),
(574,51),
(549,333),
(262,207),
(255,412),
(430,414),
(506,558),
(550,472),
(705,573),
(658,842),
(836,948),
]

size(1000,1000)

# compute the gaussian weight for a given distance from the center
rho = 1.0
def weight(d):
    return pow(e, -(d**2 / (2.0 * rho * rho))) / sqrt(2.0 * pi * rho * rho)

draw_points(points)
for i in range(1000):
    stp = 10
    k = -stp
    sumx,sumy = 0,0
    wsum = 0
    while k <= stp:
        x,y = sample_point(points, (i+k)*10)
        w = weight(abs(k)/float(stp))
        sumx += x * w
        sumy += y * w
        wsum += w      
        k += 1
    x = sumx / wsum
    y = sumy / wsum    
    ellipse(x,y,10,10)