image source1;
source1.open("input_photograph.jpg");

image target1;
target1.open("output_artwork.png");

textbuffer output;
output.open("output.txt");
output.clear();

// Select the four corners of the artwork in the photograph
// numbered from left-to-right from top-to-bottom.
point p1 = source1.getPoint(1);
point p2 = source1.getPoint(2);
point p3 = source1.getPoint(3);
point p4 = source1.getPoint(4);

int result_size = 500;

point center;
{
    point q1 = p1;
    point q2 = p4;
    point q3 = p2;
    point q4 = p3;
    double q21x = (q2.x-q1.x);
    double q21y = (q2.y-q1.y);
    double q43x = (q4.x-q3.x);
    double q43y = (q4.y-q3.y);
    double qdiv = q43x * q21y - q43y * q21x;
    bool qc = qdiv > 0.0001 || qdiv < -0.0001;
    if (qc)
    {
        double qs = ((q3.y-q1.y) * q21x - (q3.x-q1.x) * q21y) / qdiv;
        center.x = q3.x + qs * q43x;
        center.y = q3.y + qs * q43y;
    }
}

output.append(format("center %lf,%lf\n", center.x, center.y));

double rig_f = 0;
double lef_f = 0;
double top_f = 0;
double bot_f = 0;

point horz;
bool h_cross;
{
    point q1 = p1;
    point q2 = p2;
    point q3 = p3;
    point q4 = p4;
    double q21x = (q2.x-q1.x);
    double q21y = (q2.y-q1.y);
    double q43x = (q4.x-q3.x);
    double q43y = (q4.y-q3.y);
    double qdiv = q43x * q21y - q43y * q21x;
    h_cross = qdiv > 0.0001 || qdiv < -0.0001;
    if (h_cross)
    {
        double qs = ((q3.y-q1.y) * q21x - (q3.x-q1.x) * q21y) / qdiv;
        horz.x = q3.x + qs * q43x;
        horz.y = q3.y + qs * q43y;
    }
}

if (h_cross)
{
    output.append(format("horz %lf,%lf\n", horz.x, horz.y));
    {
        point q1 = center;
        point q2 = horz;
        point q3 = p1;
        point q4 = p3;
        double q21x = (q2.x-q1.x);
        double q21y = (q2.y-q1.y);
        double q43x = (q4.x-q3.x);
        double q43y = (q4.y-q3.y);
        double qdiv = q43x * q21y - q43y * q21x;
        if (qdiv > 0.0001 || qdiv < -0.0001)
        {
            double qs = ((q3.y-q1.y) * q21x - (q3.x-q1.x) * q21y) / qdiv;
            rig_f = (2 * qs - 1) / (1 - qs);
            output.append(format("s_rig %f rig_f %f\n", qs, rig_f));
        }
    }
    {
        point q1 = center;
        point q2 = horz;
        point q3 = p2;
        point q4 = p4;
        double q21x = (q2.x-q1.x);
        double q21y = (q2.y-q1.y);
        double q43x = (q4.x-q3.x);
        double q43y = (q4.y-q3.y);
        double qdiv = q43x * q21y - q43y * q21x;
        if (qdiv > 0.0001 || qdiv < -0.0001)
        {
            double qs = ((q3.y-q1.y) * q21x - (q3.x-q1.x) * q21y) / qdiv;
            lef_f = (2 * qs - 1) / (1 - qs);
            output.append(format("s_lef %f lef_f %f\n", qs, lef_f));
        }
    }
}
else
    output.append("horz parallel\n");

point vert;
bool v_cross;
{
    point q1 = p1;
    point q2 = p3;
    point q3 = p2;
    point q4 = p4;
    double q21x = (q2.x-q1.x);
    double q21y = (q2.y-q1.y);
    double q43x = (q4.x-q3.x);
    double q43y = (q4.y-q3.y);
    double qdiv = q43x * q21y - q43y * q21x;
    v_cross = qdiv > 0.0001 || qdiv < -0.0001;
    if (v_cross)
    {
        double qs = ((q3.y-q1.y) * q21x - (q3.x-q1.x) * q21y) / qdiv;
        vert.x = q3.x + qs * q43x;
        vert.y = q3.y + qs * q43y;
    }
}

if (v_cross)
{
    output.append(format("vert %lf,%lf\n", vert.x, vert.y));
    {
        point q1 = center;
        point q2 = vert;
        point q3 = p1;
        point q4 = p2;
        double q21x = (q2.x-q1.x);
        double q21y = (q2.y-q1.y);
        double q43x = (q4.x-q3.x);
        double q43y = (q4.y-q3.y);
        double qdiv = q43x * q21y - q43y * q21x;
        if (qdiv > 0.0001 || qdiv < -0.0001)
        {
            double qs = ((q3.y-q1.y) * q21x - (q3.x-q1.x) * q21y) / qdiv;
            top_f = (2 * qs - 1) / (1 - qs);
            output.append(format("s_top %f top_f %f\n", qs, top_f));
        }
    }
    {
        point q1 = center;
        point q2 = vert;
        point q3 = p3;
        point q4 = p4;
        double q21x = (q2.x-q1.x);
        double q21y = (q2.y-q1.y);
        double q43x = (q4.x-q3.x);
        double q43y = (q4.y-q3.y);
        double qdiv = q43x * q21y - q43y * q21x;
        if (qdiv > 0.0001 || qdiv < -0.0001)
        {
            double qs = ((q3.y-q1.y) * q21x - (q3.x-q1.x) * q21y) / qdiv;
            bot_f = (2 * qs - 1) / (1 - qs);
            output.append(format("s_bot %f bot_f %f\n", qs, bot_f));
        }
    }
}
else
    output.append("vert parallel\n");

output.append(format("%lf ", atan2(2, 1)));
output.append(format("%lf ", atan2(0.5, 1)));
output.append("\n");

point mid = p1 + p2 + p3 + p4;
double corner_distance = distance(mid, p1);
double r_x = distance(p1, p3) / distance(p2, p4);
double x_cor = r_x > 1.0 ? sqrt(1 + (r_x-1)*(r_x-1)) : sqrt(1 + (1/r_x-1)*((1/r_x-1)));
double r_y = distance(p1, p2) / distance(p3, p4);
double y_cor = r_y > 1.0 ? sqrt(1 + (r_y-1)*(r_y-1)) : sqrt(1 + (1/r_y-1)*((1/r_y-1)));
output.append(format("r_x = %lf (%lf), r_y = %lf (%lf)\n", r_x, x_cor, r_y, y_cor));
double sel_height = (distance(p1, p3) + distance(p2, p4))/2;
double sel_width = (distance(p1, p2) + distance(p3, p4))/2;
output.append(format("sel_height = %lf, sel_width = %lf\n", sel_height, sel_width));
double sel_height = sel_height * y_cor;
double sel_width = sel_width * x_cor;
output.append(format("sel_height = %lf, sel_width = %lf\n", sel_height, sel_width));
double ratio = sel_height / sel_width;
double sqrt_ratio = sqrt(sel_height / sel_width);
output.append(format("ratio_sqrt = %lf\n", sqrt_ratio));
// If the ratio of the art work is known, set sqrt_ratio to sqrt(height/width) here
ratio = sqrt_ratio * sqrt_ratio;
int vert_steps = round(result_size / sqrt_ratio);
int horz_steps = round(result_size * sqrt_ratio);
output.append(format("imgratio=%5f\n", ratio, sqrt_ratio));
if (sqrt_ratio > 1.0)
    output.append(format("WIDTH=\"%d\" HEIGHT=\"125\"\n", round(100/ratio)));
else
    output.append(format("WIDTH=\"125\" HEIGHT=\"%d\"\n", round(100*ratio)));

target1.new(vert_steps, horz_steps);

int loaded_width = source1.width;
int loaded_height = source1.height;

int sample_rate = round(sqrt(sel_height / horz_steps * sel_width / vert_steps));
if (sample_rate < 3)
    sample_rate = 3;

output.append(format("sample_rate = %d\n", sample_rate));

for x from 0 to vert_steps
{
    double d_x = 1.0 / vert_steps;
    double ds_x = d_x / sample_rate;
    double dsh_x = ds_x / 2;
    double v_x = dsh_x + x * d_x;
    for y from 0 to horz_steps
    {
        double d_y = 1.0 / horz_steps;
        double ds_y = d_y / sample_rate;
        double dsh_y = ds_y / 2;
        double v_y = dsh_y + y * d_y;
        color col;
        for xi from 0 to sample_rate
        {
            double t_x = v_x + xi * ds_x;
            double top_x = t_x * (1 + top_f)/(1 + top_f * t_x);
            double bot_x = t_x * (1 + bot_f)/(1 + bot_f * t_x);
            double dif_x = bot_x - top_x;
            for yi from 0 to sample_rate
            {
                double t_y = v_y + yi * ds_y;
                double rig_y = t_y * (1 + rig_f)/(1 + rig_f * t_y);
                double lef_y = t_y * (1 + lef_f)/(1 + lef_f * t_y);
                double dif_y = lef_y - rig_y;
                double div = 1 - dif_x * dif_y;
                double f_y = (rig_y + top_x * dif_y) / div;
                double g_y = 1.0 - f_y;
                double f_x = (top_x + rig_y * dif_x) / div;
                double g_x = 1.0 - f_x;
                point p = g_x*(g_y*p1 + f_y*p3) + f_x*(g_y*p2 + f_y*p4);
                col += source1.getInterpolatedColor(p);
            }
        }
        
        target1.setColor(x, y, col);
    }
}