/*
 # This file is part of the Astrometry.net suite.
 # Licensed under a 3-clause BSD style license - see LICENSE
 */
#include <unistd.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <assert.h>

/**

 wget "http://antwrp.gsfc.nasa.gov/apod/image/0403/cmsky_cortner_full.jpg"
 solve-field --backend-config backend.cfg -v --keep-xylist %s.xy --continue --scale-low 10 --scale-units degwidth cmsky_cortner_full.xy --no-tweak
 cp cmsky_cortner_full.xy 1.xy
 cp cmsky_cortner_full.rdls 1.rd
 cp cmsky_cortner_full.wcs 1.wcs
 cp cmsky_cortner_full.jpg 1.jpg

 tweak -w 1.wcs -x 1.xy -r 1.rd -v

 **/

#include "starutil.h"
#include "mathutil.h"
#include "bl.h"
#include "matchobj.h"
#include "xylist.h"
#include "rdlist.h"
#include "ioutils.h"
#include "starkd.h"
#include "boilerplate.h"
#include "sip.h"
#include "sip_qfits.h"
#include "log.h"
#include "fitsioutils.h"
#include "fit-wcs.h"
#include "verify.h"
#include "histogram2d.h"

#include "plotstuff.h"
#include "plotimage.h"
#include "cairoutils.h"

static const char* OPTIONS = "hx:w:r:vj:I:";

void print_help(char* progname) {
    BOILERPLATE_HELP_HEADER(stdout);
    printf("\nUsage: %s\n"
           "   -w <WCS input file>\n"
           "   -x <xyls input file>\n"
           "   -r <rdls input file>\n"
           "   [-I <background-image>]: background for plots.\n"
           "   [-v]: verbose\n"
           "   [-j <pixel-jitter>]: set pixel jitter (default 1.0)\n"
           "\n", progname);
}


int main(int argc, char** args) {
    int c;
    char* xylsfn = NULL;
    char* wcsfn = NULL;
    char* rdlsfn = NULL;

    xylist_t* xyls = NULL;
    rdlist_t* rdls = NULL;
    sip_t sip;
    int i, j;
    int W, H;
    //double xyzcenter[3];
    //double fieldrad2;
    double pixeljitter = 1.0;
    int loglvl = LOG_MSG;
    double wcsscale;

    char* bgfn = NULL;

    //double nsigma = 3.0;

    fits_use_error_system();

    while ((c = getopt(argc, args, OPTIONS)) != -1) {
        switch (c) {
        case 'I':
            bgfn = optarg;
            break;
        case 'j':
            pixeljitter = atof(optarg);
            break;
        case 'h':
            print_help(args[0]);
            exit(0);
        case 'r':
            rdlsfn = optarg;
            break;
        case 'x':
            xylsfn = optarg;
            break;
        case 'w':
            wcsfn = optarg;
            break;
        case 'v':
            loglvl++;
            break;
        }
    }
    if (optind != argc) {
        print_help(args[0]);
        exit(-1);
    }
    if (!xylsfn || !wcsfn || !rdlsfn) {
        print_help(args[0]);
        exit(-1);
    }
    log_init(loglvl);

    // read WCS.
    logmsg("Trying to parse SIP header from %s...\n", wcsfn);
    if (!sip_read_header_file(wcsfn, &sip)) {
        logmsg("Failed to parse SIP header from %s.\n", wcsfn);
    }
    // image W, H
    W = sip.wcstan.imagew;
    H = sip.wcstan.imageh;
    if ((W == 0.0) || (H == 0.0)) {
        logmsg("WCS file %s didn't contain IMAGEW and IMAGEH headers.\n", wcsfn);
        // FIXME - use bounds of xylist?
        exit(-1);
    }
    wcsscale = sip_pixel_scale(&sip);
    logmsg("WCS scale: %g arcsec/pixel\n", wcsscale);

    // read XYLS.
    xyls = xylist_open(xylsfn);
    if (!xyls) {
        logmsg("Failed to read an xylist from file %s.\n", xylsfn);
        exit(-1);
    }

    // read RDLS.
    rdls = rdlist_open(rdlsfn);
    if (!rdls) {
        logmsg("Failed to read an rdlist from file %s.\n", rdlsfn);
        exit(-1);
    }

    // Find field center and radius.
    /*
     sip_pixelxy2xyzarr(&sip, W/2, H/2, xyzcenter);
     fieldrad2 = arcsec2distsq(sip_pixel_scale(&sip) * hypot(W/2, H/2));
     */

    {
        // (x,y) positions of field stars.
        double* fieldpix;
        int Nfield;
        double* indexpix;
        starxy_t* xy;
        rd_t* rd;
        int Nindex;

        xy = xylist_read_field(xyls, NULL);
        if (!xy) {
            logmsg("Failed to read xyls entries.\n");
            exit(-1);
        }
        Nfield = starxy_n(xy);
        fieldpix = starxy_to_xy_array(xy, NULL);
        logmsg("Found %i field objects\n", Nfield);

        // Project RDLS into pixel space.
        rd = rdlist_read_field(rdls, NULL);
        if (!rd) {
            logmsg("Failed to read rdls entries.\n");
            exit(-1);
        }
        Nindex = rd_n(rd);
        logmsg("Found %i indx objects\n", Nindex);
        indexpix = malloc(2 * Nindex * sizeof(double));
        for (i=0; i<Nindex; i++) {
            anbool ok;
            double ra = rd_getra(rd, i);
            double dec = rd_getdec(rd, i);
            ok = sip_radec2pixelxy(&sip, ra, dec, indexpix + i*2, indexpix + i*2 + 1);
            assert(ok);
        }

        logmsg("CRPIX is (%g,%g)\n", sip.wcstan.crpix[0], sip.wcstan.crpix[1]);

        /*

         // ??
         // Look for index-field pairs that are (a) close together; and (b) close to CRPIX.

         // Split the image into 3x3, 5x5 or so, and in each, look for a
         // (small) rotation and log(scale), then (bigger) shift, using histogram
         // cross-correlation.

         // Are the rotations and scales really going to be big enough that this
         // is required, or can we get away with doing shift first, then fine-tuning
         // rotation and scale?

         {
         // NxN blocks
         int NB = 3;
         int b;
         // HACK - use histogram2d machinery to split image into blocks.
         histogram2d* blockhist = histogram2d_new_nbins(0, W, NB, 0, H, NB);
         int* fieldi = malloc(Nfield * sizeof(int));
         int* indexi = malloc(Nindex * sizeof(int));
         // rotation bins
         int NR = 100;
         // scale bins (ie, log(radius) bins)
         double minrad = 1.0;
         double maxrad = 200.0;
         int NS = 100;
         histogram2d* rsfield = histogram2d_new_nbins(-M_PI, M_PI, NR,
         log(minrad), log(maxrad), NS);
         histogram2d* rsindex = histogram2d_new_nbins(-M_PI, M_PI, NR,
         log(minrad), log(maxrad), NS);
         histogram2d_set_y_edges(rsfield, HIST2D_DISCARD);
         histogram2d_set_y_edges(rsindex, HIST2D_DISCARD);

         for (b=0; b<(NB*NB); b++) {
         int bin;
         int NF, NI;
         double dx, dy;
         NF = NI = 0;
         for (i=0; i<Nfield; i++) {
         bin = histogram2d_add(blockhist, fieldpix[2*i], fieldpix[2*i+1]);
         if (bin != b)
         continue;
         fieldi[NF] = i;
         NF++;
         }

         for (i=0; i<Nindex; i++) {
         bin = histogram2d_add(blockhist, indexpix[2*i], indexpix[2*i+1]);
         if (bin != b)
         continue;
         indexi[NI] = i;
         NI++;
         }
         logmsg("bin %i has %i field and %i index stars.\n", b, NF, NI);

         logmsg("histogramming field rotation/scale\n");
         for (i=0; i<NF; i++) {
         for (j=0; j<i; j++) {
         dx = fieldpix[2*fieldi[i]] - fieldpix[2*fieldi[j]];
         dy = fieldpix[2*fieldi[i]+1] - fieldpix[2*fieldi[j]+1];
         histogram2d_add(rsfield, atan2(dy, dx), log(sqrt(dx*dx + dy*dy)));
         }
         }
         logmsg("histogramming index rotation/scale\n");
         for (i=0; i<NI; i++) {
         for (j=0; j<i; j++) {
         dx = indexpix[2*indexi[i]] - fieldpix[2*indexi[j]];
         dy = indexpix[2*indexi[i]+1] - fieldpix[2*indexi[j]+1];
         histogram2d_add(rsindex, atan2(dy, dx), log(sqrt(dx*dx + dy*dy)));
         }
         }


         }
         histogram2d_free(rsfield);
         histogram2d_free(rsindex);
         free(fieldi);
         free(indexi);
         histogram2d_free(blockhist);
         }
         */

        {
            double* fieldsigma2s = malloc(Nfield * sizeof(double));
            int besti;
            int* theta;
            double logodds;
            double Q2, R2;
            double qc[2];
            double gamma;

            // HACK -- quad radius-squared
            Q2 = square(100.0);
            qc[0] = sip.wcstan.crpix[0];
            qc[1] = sip.wcstan.crpix[1];
            // HACK -- variance growth rate wrt radius.
            gamma = 1.0;

            for (i=0; i<Nfield; i++) {
                R2 = distsq(qc, fieldpix + 2*i, 2);
                fieldsigma2s[i] = square(pixeljitter) * (1.0 + gamma * R2/Q2);
            }

            logodds = verify_star_lists(indexpix, Nindex,
                                        fieldpix, fieldsigma2s, Nfield,
                                        W*H,
                                        0.25,
                                        log(1e-100),
                                        log(1e100),
                                        &besti, NULL, &theta, NULL, NULL);

            logmsg("Logodds: %g\n", logodds);

            if (bgfn) {
                plot_args_t pargs;
                plotimage_t* img;
                cairo_t* cairo;
                char outfn[32];

                j = 0;
				
                plotstuff_init(&pargs);
                pargs.outformat = PLOTSTUFF_FORMAT_PNG;
                sprintf(outfn, "tweak-%03i.png", j);
                pargs.outfn = outfn;
                img = plotstuff_get_config(&pargs, "image");
                //img->format = PLOTSTUFF_FORMAT_JPG; // guess
                plot_image_set_filename(img, bgfn);
                plot_image_setsize(&pargs, img);
                plotstuff_run_command(&pargs, "image");
                cairo = pargs.cairo;
                // red circles around every field star.
                cairo_set_color(cairo, "red");
                for (i=0; i<Nfield; i++) {
                    cairoutils_draw_marker(cairo, CAIROUTIL_MARKER_CIRCLE,
                                           fieldpix[2*i+0], fieldpix[2*i+1],
                                           2.0 * sqrt(fieldsigma2s[i]));
                    cairo_stroke(cairo);
                }
                // green crosshairs at every index star.
                cairo_set_color(cairo, "green");
                for (i=0; i<Nindex; i++) {
                    cairoutils_draw_marker(cairo, CAIROUTIL_MARKER_XCROSSHAIR,
                                           indexpix[2*i+0], indexpix[2*i+1],
                                           3);
                    cairo_stroke(cairo);
                }

                // thick white circles for corresponding field stars.
                cairo_set_line_width(cairo, 2);
                for (i=0; i<Nfield; i++) {
                    if (theta[i] < 0)
                        continue;
                    cairo_set_color(cairo, "white");
                    cairoutils_draw_marker(cairo, CAIROUTIL_MARKER_CIRCLE,
                                           fieldpix[2*i+0], fieldpix[2*i+1],
                                           2.0 * sqrt(fieldsigma2s[i]));
                    cairo_stroke(cairo);
                    // thick cyan crosshairs for corresponding index stars.
                    cairo_set_color(cairo, "cyan");
                    cairoutils_draw_marker(cairo, CAIROUTIL_MARKER_XCROSSHAIR,
                                           indexpix[2*theta[i]+0],
                                           indexpix[2*theta[i]+1],
                                           3);
                    cairo_stroke(cairo);
					
                }

                plotstuff_output(&pargs);
            }


            free(theta);
            free(fieldsigma2s);
        }


        free(fieldpix);
        free(indexpix);
    }



    if (xylist_close(xyls)) {
        logmsg("Failed to close XYLS file.\n");
    }
    return 0;
}