perspective.py

   1 #!/usr/bin/env python3
   2
   3 import numpy
   4 import numpy.linalg
   5
   6 # call with 2x4 vector (4 points of x, y)
   7 def calc_basis(p):
   8   # see https://math.stackexchange.com/questions/296794/finding-the-transform-matrix-from-4-projected-points-with-javascript
   9   A = numpy.ones((3, 4), numpy.double)
  10   A[:2, :] = p
  11   return A[:, :3] * numpy.linalg.solve(A[:, :3], A[:, 3])[numpy.newaxis, :]
  12
  13 def calc_transform(p, q):
  14   A = calc_basis(p)
  15   #y = A @ numpy.array(
  16   #  [
  17   #    [1., 0., 0., 1.],
  18   #    [0., 1., 0., 1.],
  19   #    [0., 0., 1., 1.]
  20   #  ],
  21   #  numpy.double
  22   #)
  23   #print ('xxx', y[:2, :] / y[2:, :], p)
  24
  25   B = calc_basis(q)
  26   #y = B @ numpy.array(
  27   #  [
  28   #    [1., 0., 0., 1.],
  29   #    [0., 1., 0., 1.],
  30   #    [0., 0., 1., 1.]
  31   #  ],
  32   #  numpy.double
  33   #)
  34   #print ('yyy', y[:2, :] / y[2:, :], q)
  35
  36   #return B @ numpy.linalg.inv(A)
  37   return numpy.linalg.solve(A.transpose(), B.transpose()).transpose()
  38
  39 # call with 2xN vector (N points of x, y), computes 2xN vector
  40 def apply_transform_multi(A, p):
  41   x = numpy.ones((3, p.shape[1]), numpy.double)
  42   x[:2, :] = p
  43   y = A @ x
  44   return y[:2, :] / y[2:, :]
  45
  46 # call with 2-vector (x, y), computes 2-vector
  47 def apply_transform(A, p):
  48   return apply_transform_multi(A, p[:, numpy.newaxis])[:, 0]
  49
  50 if __name__ == '__main__':
  51   x = numpy.array([[11., 52., 23., 74.], [15., 36., 27., 58.]], numpy.double)
  52   y = numpy.array([[31., 92., 73., 24.], [65., 26., 87., 18.]], numpy.double)
  53   A = calc_transform(x, y)
  54   print(apply_transform_multi(A, x), y)