+#!/usr/bin/env python3
+
+import gamma
+import numpy
+import perspective
+import scipy
+import scipy.ndimage
+import sys
+
+sys.stderr.write(
+ '''stdin: control blocks description
+xs0,ys0 xs1,ys1 xs2,ys2 xs3,ys3 : size (pixels)
+xf0,yf0 xf1,yf1 xf2,yf2 xf3,yf3 : freedom (fraction of size)
+in.jpg out.jpg x0,y0 x1,y1 x2,y2 x3,y3 : top left (pixels)
+...
+'''
+)
+
+size = numpy.array(
+ [
+ [int(j) for j in i.split(',')]
+ for i in sys.stdin.readline().split()
+ ],
+ numpy.int32
+).transpose()
+assert size.shape == (2, 4)
+freedom = numpy.array(
+ [
+ [float(j) for j in i.split(',')]
+ for i in sys.stdin.readline().split()
+ ],
+ numpy.double
+).transpose()
+assert freedom.shape == (2, 4)
+
+files = []
+line = sys.stdin.readline()
+while len(line):
+ fields = line.split()
+ in_jpg = fields[0]
+ out_jpg = fields[1]
+ top_left = numpy.array(
+ [
+ [int(j) for j in i.split(',')]
+ for i in fields[2:]
+ ],
+ numpy.int32
+ ).transpose()
+ assert top_left.shape == (2, 4)
+ files.append((in_jpg, out_jpg, top_left))
+ line = sys.stdin.readline()
+
+in_jpg, out_jpg, top_left0 = files[0]
+sys.stderr.write(f'read {in_jpg:s}\n')
+image0 = gamma.read_image(in_jpg)
+shape = image0.shape
+sys.stderr.write(f'shape {str(shape):s}\n')
+
+sys.stderr.write(f'write {out_jpg:s}\n')
+gamma.write_image(out_jpg, image0)
+
+A = numpy.identity(3, numpy.double)
+for i in range(len(files) - 1):
+ in_jpg, out_jpg, top_left1 = files[i + 1]
+ sys.stderr.write(f'read {in_jpg:s}\n')
+ image1 = gamma.read_image(in_jpg)
+ assert image1.shape == shape
+
+ p = []
+ q = []
+ for j in range(4):
+ xs, ys = size[:, j]
+ xf, yf = freedom[:, j]
+ xf = int(round(xf * xs))
+ yf = int(round(yf * ys))
+ x0, y0 = top_left0[:, j]
+ x1, y1 = top_left1[:, j]
+
+ block0 = image0[
+ max(y0, 0):max(y0 + ys, 0),
+ max(x0, 0):max(x0 + xs, 0),
+ :
+ ]
+ #gamma.write_image(f'diag{j:d}.jpg', block0)
+ mean0 = numpy.mean(numpy.mean(block0, 0), 0)
+ data0 = block0 - mean0[numpy.newaxis, numpy.newaxis, :]
+
+ best = None
+ for x in range(-xf, xf + 1):
+ for y in range(-yf, yf + 1):
+ block1 = image1[
+ max(y1 + y, 0):max(y1 + y + ys, 0),
+ max(x1 + x, 0):max(x1 + x + xs, 0),
+ :
+ ]
+ if block1.shape != block0.shape:
+ continue
+ mean1 = numpy.mean(numpy.mean(block1, 0), 0)
+ data1 = block1 - mean1[numpy.newaxis, numpy.newaxis, :]
+
+ score = numpy.sum(data0 * data1)
+ if best is None or score > best[0]:
+ best = (score, x, y)
+
+ score, x, y = best
+ sys.stderr.write(f'score {score:.3f} x {x / xs:.3f} y {y / ys:.3f}\n')
+ p.append(numpy.array([x0 + xs / 2, y0 + ys / 2], numpy.double))
+ q.append(numpy.array([x1 + x + xs / 2, y1 + y + ys / 2], numpy.double))
+
+ top_left1[:, j] = (x1 + x, y1 + y)
+ p = numpy.stack(p, 1)
+ q = numpy.stack(q, 1)
+ A = perspective.calc_transform(p, q) @ A # makes the transform cumulative
+
+ sys.stderr.write('remap\n')
+ ys, xs, cs = shape
+ #def mapping(p):
+ # #return tuple(
+ # # perspective.apply_transform(A, numpy.array(p, numpy.double)[::-1])[::-1]
+ # #)
+ # q = A[:, 0] * p[1] + A[:, 1] * p[0] + A[:, 2]
+ # return (q[1] / q[2], q[0] / q[2])
+ #out_image1 = numpy.stack(
+ # [
+ # scipy.ndimage.geometric_transform(image1[:, :, j], mapping)
+ # for j in range(cs)
+ # ],
+ # 2
+ #)
+ coords = numpy.zeros((2, ys, xs), numpy.double)
+ coords[0, :, :] = numpy.arange(ys, dtype = numpy.double)[:, numpy.newaxis]
+ coords[1, :, :] = numpy.arange(xs, dtype = numpy.double)[numpy.newaxis, :]
+ mapped_coords = perspective.apply_transform_multi(
+ A,
+ coords[::-1, :, :].reshape((2, ys * xs))
+ ).reshape(2, ys, xs)[::-1, :, :]
+ out_image1 = numpy.stack(
+ [
+ scipy.ndimage.map_coordinates(image1[:, :, j], mapped_coords)
+ for j in range(cs)
+ ],
+ 2
+ )
+
+ sys.stderr.write(f'write {out_jpg:s}\n')
+ gamma.write_image(out_jpg, out_image1)
+
+ image0 = image1
+ top_left0 = top_left1
+
+sys.stdout.write(
+ ' '.join(
+ [
+ ','.join([str(size[j, i]) for j in range(2)])
+ for i in range(4)
+ ]
+ ) + '\n'
+)
+sys.stdout.write(
+ ' '.join(
+ [
+ ','.join([str(freedom[j, i]) for j in range(2)])
+ for i in range(4)
+ ]
+ ) + '\n'
+)
+for in_jpg, out_jpg, top_left in files:
+ sys.stdout.write(
+ ' '.join(
+ [in_jpg, out_jpg] +
+ [
+ ','.join([str(top_left[j, i]) for j in range(2)])
+ for i in range(4)
+ ]
+ ) + '\n'
+ )