perspective.py

   1 #!/usr/bin/env python3
   2
   3 import numpy
   4 import numpy.linalg
   5 import scipy
   6 import scipy.ndimage
   7
   8 # call with 2x4 vector (4 points of x, y)
   9 def calc_basis(p):
  10   # see https://math.stackexchange.com/questions/296794/finding-the-transform-matrix-from-4-projected-points-with-javascript
  11   A = numpy.ones((3, 4), numpy.double)
  12   A[:2, :] = p
  13   return A[:, :3] * numpy.linalg.solve(A[:, :3], A[:, 3])[numpy.newaxis, :]
  14
  15 def calc_transform(p, q):
  16   A = calc_basis(p)
  17   #y = A @ numpy.array(
  18   #  [
  19   #    [1., 0., 0., 1.],
  20   #    [0., 1., 0., 1.],
  21   #    [0., 0., 1., 1.]
  22   #  ],
  23   #  numpy.double
  24   #)
  25   #print ('xxx', y[:2, :] / y[2:, :], p)
  26
  27   B = calc_basis(q)
  28   #y = B @ numpy.array(
  29   #  [
  30   #    [1., 0., 0., 1.],
  31   #    [0., 1., 0., 1.],
  32   #    [0., 0., 1., 1.]
  33   #  ],
  34   #  numpy.double
  35   #)
  36   #print ('yyy', y[:2, :] / y[2:, :], q)
  37
  38   #return B @ numpy.linalg.inv(A)
  39   return numpy.linalg.solve(A.transpose(), B.transpose()).transpose()
  40
  41 # call with 2xN vector (N points of x, y), computes 2xN vector
  42 def apply_transform_multi(A, p):
  43   x = numpy.ones((3, p.shape[1]), numpy.double)
  44   x[:2, :] = p
  45   y = A @ x
  46   return y[:2, :] / y[2:, :]
  47
  48 # call with 2-vector (x, y), computes 2-vector
  49 def apply_transform(A, p):
  50   return apply_transform_multi(A, p[:, numpy.newaxis])[:, 0]
  51
  52 def remap_image(A, image):
  53   ys, xs, cs = image.shape
  54   coords = numpy.zeros((2, ys, xs), numpy.double)
  55   coords[0, :, :] = numpy.arange(ys, dtype = numpy.double)[:, numpy.newaxis]
  56   coords[1, :, :] = numpy.arange(xs, dtype = numpy.double)[numpy.newaxis, :]
  57   mapped_coords = apply_transform_multi(
  58     A,
  59     coords[::-1, :, :].reshape((2, ys * xs))
  60   ).reshape(2, ys, xs)[::-1, :, :]
  61   return numpy.stack(
  62     [
  63       scipy.ndimage.map_coordinates(image[:, :, j], mapped_coords)
  64       for j in range(cs)
  65     ],
  66     2
  67   )
  68
  69 if __name__ == '__main__':
  70   x = numpy.array([[11., 52., 23., 74.], [15., 36., 27., 58.]], numpy.double)
  71   y = numpy.array([[31., 92., 73., 24.], [65., 26., 87., 18.]], numpy.double)
  72   A = calc_transform(x, y)
  73   print(apply_transform_multi(A, x), y)