Transform (morph) image according to control points from Hugin

morph.py

import sys
import numpy as np
from sklearn import linear_model
import imageio

# control point file format:
# x1 tab y1 tab x2 tab y2
# x1 tab y1 tab x2 tab y2
# ...
def read_cp(cp_path):
    cps = []
    with open(cp_path) as fp:
        for line in fp:
            if line.startswith('#'):
                continue
            cps.append([float(x) for x in line.strip().split("\t")])
    return cps

def get_poly(x, y):
    x /= 1000 # Make numbers smaller. Not sure if it helps with floating-point error
    y /= 1000
    return [x, y, x*x, x*y, y*y,
            x*x*x, x*x*y, x*y*y, y*y*y,
            x*x*x*x, x*x*x*y, x*x*y*y, x*y*y*y, y*y*y*y]

def main():
    if len(sys.argv) != 5:
        print("Usage: python3 morph.py control-points.txt in-photo1.jpg in-photo2.jpg out.jpg")
        print("out.jpg will have the same width and height as in-photo2.jpg, but with in-photo1.jpg's content.")
        return
    cp_path = sys.argv[1]
    in1_path = sys.argv[2]
    in2_path = sys.argv[3]
    out_path = sys.argv[4]

    cps = read_cp(cp_path)
    reg_in = []
    reg_x1 = []
    reg_y1 = []
    for x1, y1, x2, y2 in cps:
        reg_x1.append(x1)
        reg_y1.append(y1)
        reg_in.append(get_poly(x2, y2))
    reg = linear_model.LinearRegression()
    reg.fit(reg_in, reg_x1)
    coef_x1 = reg.coef_.tolist().copy()
    icpt_x1 = reg.intercept_
    reg.fit(reg_in, reg_y1)
    coef_y1 = reg.coef_.tolist().copy()
    icpt_y1 = reg.intercept_

    # check error
    print("Error: error, x1, y1, x2, y2")
    for x1, y1, x2, y2 in cps:
        x1p = icpt_x1 + sum([a*b for a, b in zip(get_poly(x2, y2), coef_x1)])
        y1p = icpt_y1 + sum([a*b for a, b in zip(get_poly(x2, y2), coef_y1)])
        dist = (x1p-x1)*(x1p-x1) + (y1p-y1)*(y1p-y1)
        print(dist, x1, y1, x2, y2)

    img1 = imageio.imread(in1_path)
    print("image1:", img1.shape)
    img2 = imageio.imread(in2_path)
    print("image2:", img2.shape)
    img2.fill(0)

    # check if cp falls in image dimentation
    for x1, y1, x2, y2 in cps:
        if x1 < 0 or x1 >= img1.shape[1] or y1 < 0 or y1 >= img1.shape[0]:
            print(x1, y1, x2, y2, "falls outside", in1_path)
            return
        if x2 < 0 or x2 >= img2.shape[1] or y2 < 0 or y2 >= img2.shape[0]:
            print(x1, y1, x2, y2, "falls outside", in2_path)
            return

    print("Generating image. It can take some time.")
    for y2, x2 in np.ndindex(img2.shape[:2]):
        x1 = icpt_x1 + sum([a*b for a, b in zip(get_poly(x2, y2), coef_x1)])
        y1 = icpt_y1 + sum([a*b for a, b in zip(get_poly(x2, y2), coef_y1)])

        # nearest neighbour (no interpolation)
        #x1 = int(x1 + 0.5);
        #y1 = int(y1 + 0.5);
        #if x1 >= 0 and y1 >= 0 and x1 < img1.shape[1] and y1 < img1.shape[0]:
        #    img2[y2, x2, :] = img1[y1, x1, :]
        #else:
        #    img2[y2, x2, :] = 0

        # bilinear interpolation
        x1i = int(x1)
        y1i = int(y1)
        if x1i >= 0 and y1i >= 0 and x1i+1 < img1.shape[1] and y1i+1 < img1.shape[0]:
            img2[y2, x2, :] = img1[y1i, x1i, :] * (y1i+1-y1) * (x1i+1-x1) \
                    + img1[y1i, x1i+1, :] * (y1i+1-y1) * (x1-x1i) \
                    + img1[y1i+1, x1i, :] * (y1-y1i) * (x1i+1-x1) \
                    + img1[y1i+1, x1i+1, :] * (y1-y1i) * (x1-x1i)
        else:
            img2[y2, x2, :] = 0
    imageio.imwrite(out_path, img2, quality=95)

if __name__ == '__main__':
    main()