void homography_from_4pt(const float *x, const float *y, const float *z, const float *w, float h[8])
{
    float t1 = x[0];   float t2 = z[0];  float t4 = y[1];  float t5 = t1 * t2 * t4;
    float t6 = w[1];  float t7 = t1 * t6;  float t8 = t2 * t7;  float t9 = z[1];
    float t10 = t1 * t9;  float t11 = y[0];  float t14 = x[1];  float t15 = w[0];
    float t16 = t14 * t15;  float t18 = t16 * t11;  float t20 = t15 * t11 * t9;
    float t21 = t15 * t4;  float t24 = t15 * t9;  float t25 = t2 * t4;
    float t26 = t6 * t2;  float t27 = t6 * t11;  float t28 = t9 * t11;
    float t30 = 0.1e1 / (-t24 + t21 - t25 + t26 - t27 + t28);
    float t32 = t1 * t15;  float t35 = t14 * t11;  float t41 = t4 * t1;
    float t42 = t6 * t41;  float t43 = t14 * t2;  float t46 = t16 * t9;
    float t48 = t14 * t9 * t11;  float t51 = t4 * t6 * t2;  float t55 = t6 * t14;

    h[0] = -(-t5 + t8 + t10 * t11 - t11 * t7 - t16 * t2 + t18 - t20 + t21 * t2) * t30;
    h[1] = (t5 - t8 - t32 * t4 + t32 * t9 + t18 - t2 * t35 + t27 * t2 - t20) * t30;
    h[2] = t1;
    h[3] = (-t9 * t7 + t42 + t43 * t4 - t16 * t4 + t46 - t48 + t27 * t9 - t51) * t30;
    h[4] = (-t42 + t41 * t9 - t55 * t2 + t46 - t48 + t55 * t11 + t51 - t21 * t9) * t30;
    h[5] = t14;
    h[6] = (-t10 + t41 + t43 - t35 + t24 - t21 - t26 + t27) * t30;
    h[7] = (-t7 + t10 + t16 - t43 + t27 - t28 - t21 + t25) * t30;
}

void homography_from_4corresp(
    const float *p1, const float *p2, const float *p3, const float *p4,
    const float *q1, const float *q2, const float *q3, const float *q4, float H[3][3])
{
    float Hr[3][3], Hl[3][3];

    homography_from_4pt(p1, p2, p3, p4, &Hr[0][0]);
    homography_from_4pt(q1, q2, q3, q4, &Hl[0][0]);

    // the following code computes R = Hl * inverse Hr
    float t2 = Hr[1][1] - Hr[2][1] * Hr[1][2];
    float t4 = Hr[0][0] * Hr[1][1];
    float t5 = Hr[0][0] * Hr[1][2];
    float t7 = Hr[1][0] * Hr[0][1];
    float t8 = Hr[0][2] * Hr[1][0];
    float t10 = Hr[0][1] * Hr[2][0];
    float t12 = Hr[0][2] * Hr[2][0];
    float t15 = 1 / (t4 - t5*Hr[2][1] - t7 + t8*Hr[2][1] + t10*Hr[1][2] - t12*Hr[1][1]);
    float t18 = -Hr[1][0] + Hr[1][2] * Hr[2][0];
    float t23 = -Hr[1][0] * Hr[2][1] + Hr[1][1] * Hr[2][0];
    float t28 = -Hr[0][1] + Hr[0][2] * Hr[2][1];
    float t31 = Hr[0][0] - t12;
    float t35 = Hr[0][0] * Hr[2][1] - t10;
    float t41 = -Hr[0][1] * Hr[1][2] + Hr[0][2] * Hr[1][1];
    float t44 = t5 - t8;
    float t47 = t4 - t7;
    float t48 = t2*t15;
    float t49 = t28*t15;
    float t50 = t41*t15;
    H[0][0] = Hl[0][0] * t48 + Hl[0][1] * (t18*t15) - Hl[0][2] * (t23*t15);
    H[0][1] = Hl[0][0] * t49 + Hl[0][1] * (t31*t15) - Hl[0][2] * (t35*t15);
    H[0][2] = -Hl[0][0] * t50 - Hl[0][1] * (t44*t15) + Hl[0][2] * (t47*t15);
    H[1][0] = Hl[1][0] * t48 + Hl[1][1] * (t18*t15) - Hl[1][2] * (t23*t15);
    H[1][1] = Hl[1][0] * t49 + Hl[1][1] * (t31*t15) - Hl[1][2] * (t35*t15);
    H[1][2] = -Hl[1][0] * t50 - Hl[1][1] * (t44*t15) + Hl[1][2] * (t47*t15);
    H[2][0] = Hl[2][0] * t48 + Hl[2][1] * (t18*t15) - t23*t15;
    H[2][1] = Hl[2][0] * t49 + Hl[2][1] * (t31*t15) - t35*t15;
    H[2][2] = -Hl[2][0] * t50 - Hl[2][1] * (t44*t15) + t47*t15;
}