In [48]:
import copy
import math

import numpy as np
from numpy.linalg import inv

from skimage import transform as tf
from skimage import io

import grass.script as gscript
from grass.script import array as garray

In [49]:
im = io.imread('/home/anna/Documents/Projects/Hipp_STC/pictures/camera_9955.jpg')


(480, 640, 3)

In [50]:
dst = np.array((
(458.018387, 151.694736), 
(49.573341, 136.87887),
(101.722019, 218.340527),
(621.549653, 243.89273),
(330.494088, 183.206247),
(581.225082, 140.086904),
(288.971758, 268.646427),
(40.071481, 186.059015),
(125.698245, 223.043346),
(550.20581, 245.02926)
))

dst[:, 1] = im.shape[0] - dst[:, 1]

src = np.array(((557017.666939, 5533510.90344),
                (556984.233666, 5533500.45554),
                (556985.446273, 5533533.8606),
                (557052.614491, 5533543.83213),
                (557004.272519, 5533522.77226),
                (557027.087377, 5533507.93463),
                (557006.506142, 5533563.93473),
                (556981.936295, 5533521.97454),
                (556987.360807, 5533537.76944),
                (557041.925012, 5533549.57573)
                ))

centerx, centery = np.min(src[:, 0]), np.min(src[:, 1])
src[:, 0] -= centerx
src[:, 1] -= centery
revers = 200
src[:, 1] = revers - src[:, 1]


[[ 458.018387  328.305264]
 [  49.573341  343.12113 ]
 [ 101.722019  261.659473]
 [ 621.549653  236.10727 ]
 [ 330.494088  296.793753]
 [ 581.225082  339.913096]
 [ 288.971758  211.353573]
 [  40.071481  293.940985]
 [ 125.698245  256.956654]
 [ 550.20581   234.97074 ]]

In [51]:
tform3 = tf.ProjectiveTransform()
tform3.estimate(src, dst)
warped = tf.warp(im, tform3)

In [52]:
gscript.run_command('g.region', w=centerx, e=centerx + im.shape[1],
                    n=centery + revers, s=centery - im.shape[0] + revers, res=1)

name = 'rectified'
for num, color in zip([0, 1, 2], 'rgb'):
    rectified = garray.array()
    for y in range(rectified.shape[0]):
        for x in range(rectified.shape[1]):
            rectified[y,x] = round(255 * warped[y, x, num])
    rectified.write(mapname=name + '_' + color, overwrite=True)
gscript.run_command('r.colors', map=[name + '_r', name + '_g', name + '_b'], color='grey')
gscript.run_command('r.composite', red=name + '_r', green=name + '_g', blue=name + '_b',
                    output=name, overwrite=True)


Out[52]:
0

In [66]:
H = inv(tform3.params)
new = []
for point in dst:
    Z = float(point[0]) * H[2, 0] + point[1] * H[2, 1] + H[2, 2]
    X = (point[0] * H[0, 0] + point[1] * H[0, 1] + H[0, 2]) / Z
    Y = (point[0] * H[1, 0] + point[1] * H[1, 1] + H[1, 2]) / Z
    X += centerx
    Y = revers - Y + centery
    new.append((X, Y))

gscript.write_command('v.in.ascii', input='-',
                      output='rectified_points', overwrite=True,
                      stdin='\n'.join(['{}|{}'.format(X, Y) for X, Y in new]))


Out[66]:
0