align curves test

align.py

import numpy as np
import matplotlib.pyplot as plt
import scipy

def offset_and_scale(curve, c, a, b):
  return np.stack([curve[:, 0] + c, curve[:, 1] * a + b], axis=1)


def align_curves(curve1, curve2):
  def distance(args):
    c, a, b = args
    curve_offset = offset_and_scale(curve1, c, 1 + a, b)
    dists = scipy.spatial.distance.cdist(curve2, curve_offset)
    return np.concatenate([dists.min(axis=1), dists.min(axis=0)])

  res = scipy.optimize.least_squares(distance, x0=[0, 0, 0])
  print(res)
  return offset_and_scale(curve1,res.x[0], 1 + res.x[1], res.x[2]), (res.x[0], res.x[1], res.x[2])

def smooth_curve(curve, window_size=11, repeat=3):
  
  for x in range(repeat):
    smoothed = np.zeros_like(curve)
    smoothed[:, 0] = np.convolve(curve[:, 0], np.ones((window_size,)) / window_size, mode='same')
    smoothed[:, 1] = np.convolve(curve[:, 1], np.ones((window_size,)) / window_size, mode='same')

    curve = smoothed
  return smoothed

  

while True:
  a, b, c = np.random.rand(3) * 0.4
  x = np.linspace(0, 20, 1000)
  y = np.sin(x)

  noise_level = 0.1 * np.random.rand()

  noise1 = smooth_curve(np.random.randn(1000, 2) * noise_level, window_size=11)
  noise2 = smooth_curve(np.random.randn(1000, 2) * noise_level, window_size=11) 

  print("noise level", noise_level)

  curve1 = np.stack([x, y], axis=1) + noise1
  curve2 = offset_and_scale(curve1, c, 1 + a, b) + noise2

  print("truth", c, a, b)
  
  curve2_aligned, aligned_params = align_curves(curve1, curve2)
  print("aligned", aligned_params)

  plt.plot(curve1[:, 0], curve1[:, 1], label='curve1')
  plt.plot(curve2[:, 0], curve2[:, 1], label='curve2')
  plt.plot(curve2_aligned[:, 0], curve2_aligned[:, 1], label='curve2_aligned')
  plt.legend()
  plt.show()