/*********************************************************************
 * Main for final assignment in Computational Geometry, June 1999.
 * ===============================================================
 * 
 * The main program reads two files:
 * 1) input-file, which defines the set of points and polygons.
 * 2) command-file, which includes insert, delete, and draw
 *      commands to the algorithm.
 * The program then uses the Algorithm class, defined in the file
 * "algorithm.h", to perform the tasks described in the exercise.
 * The program displays the solutions computed by the algorithm
 * (when there is a 'draw' command), and writes them to a file (if
 * requested). It also prints the total running time of the algorithm.
 *
 * Compilation: use the supplied 'makefile' (on valis).
 * Usage: "main radius input-file command-file [output-file]"
 *   radius        - the radius (an integer!)
 *   input-file    - file with points and polygons
 *   command-file  - file with commands
 *   output-file   - file to write solutions (optional)
 * 
 * You are expected to implement the Algorithm class - see
 * files "algorithm.h" and "algorithm.C".
 * Do NOT change this file.
 * In case of bugs/problems/requests, please contact Chaim Linhart
 * at email:  chaim@math.tau.ac.il .
 *
 *********************************************************************/

#include "algorithm.h"
#include <CGAL/IO/Window_stream.h>
#include <CGAL/Timer.h>

#define WINDOW_SIZE  400
#define M_LINE       100

#define COMM_INSERT_STR   "Insert"
#define COMM_DELETE_STR   "Delete"
#define COMM_DRAW_STR     "Draw"

#define MIN(x,y)  ((x) < (y) ? (x) : (y))
#define MAX(x,y)  ((x) > (y) ? (x) : (y))

#define CHECK_RANGE(i,m)                                      \
            do {                                              \
               if (((i) < 0) || ((i) >= (int) (m))) {         \
                  cerr << "ERROR: Index (" << i               \
                       << ") out of range (" << m << ")!\n";  \
                  exit (1);                                   \
               }                                              \
            } while (0)


//////////////////////////////////////////////////////////////
// Wait for button click in window
//////////////////////////////////////////////////////////////
void any_button(CGAL_Window_stream &W)
{
    double x, y;
    cout << "Click on window to continue...\n";
    W.read_mouse (x,y);
}


//////////////////////////////////////////////////////////////
// Draw the polygons
//////////////////////////////////////////////////////////////
void draw_polygons (CGAL_Window_stream& W, vector<Polygon> &polygons)
{
    EdgeIterator  e_begin, e_end, e;
    int i;
    
    W << CGAL_BLACK;

    for (i=0; i < (int) polygons.size(); i++) {

        e_begin = polygons[i].edges_begin();
        e_end   = polygons[i].edges_end();
        for (e=e_begin; e != e_end; e++)
            W << (*e);
        
    }
}

//////////////////////////////////////////////////////////////
// Draw the points
//////////////////////////////////////////////////////////////
void draw_points (CGAL_Window_stream& W, 
                  vector<Point> &points, vector<bool> &active)
{
    int i;
    
    for (i=0; i < (int) points.size(); i++) {
        if (active[i])
            W << CGAL_BLUE;
        else
            W << CGAL_RED;
        W << points[i];
    }
}

//////////////////////////////////////////////////////////////
// Clear the window, and draw the points and polygons
//////////////////////////////////////////////////////////////
void draw (CGAL_Window_stream& W, 
           vector<Point> &points, vector<bool> &active,
           vector<Polygon> &polygons)
{
    W << CGAL_WHITE;
    W.clear();
    draw_polygons (W, polygons);
    draw_points   (W, points, active);
}

//////////////////////////////////////////////////////////////
// Read points and polygons from file
//////////////////////////////////////////////////////////////
int read_input_file (char *filename, 
                     vector<Point> &points, 
                     vector<Polygon> &polygons)
{
    ifstream ifs(filename);
    int n_points, n_polygons, n_vertices;
    int i,j, x,y;

    // Read number of points
    if (ifs.eof()) return (0);
    ifs >> n_points;
    if (n_points < 1) return (0);

    // Read points
    for (i=0; i < n_points; i++) {
        if (ifs.eof()) return (0);
        ifs >> x >> y;
        Point p(x,y);
        points.push_back (p);
    }
        
    // Read number of polygons
    if (ifs.eof()) return (0);
    ifs >> n_polygons;
    if (n_polygons < 0) return (0);

    // Read polygons
    for (i=0; i < n_polygons; i++) {
        
        Polygon poly;
        
        // Read number of vertices
        if (ifs.eof()) return (0);
        ifs >> n_vertices;
        if (n_vertices < 3) return (0);
    
        // Read vertices
        for (j=0; j < n_vertices; j++) {
            if (ifs.eof()) return (0);
            ifs >> x >> y;
            Point p(x,y);
            poly.push_back (p);
        }
        
        polygons.push_back (poly);
        
    }
        
    ifs.close();
    return (1);
}


//////////////////////////////////////////////////////////////
// Read commands from file
//////////////////////////////////////////////////////////////
int read_commands_file (char *filename, int n_points,
                        vector<command_t> &commands) 
{
    ifstream ifs(filename);
    char line[M_LINE+1], comm_str[M_LINE+1];
    int  ind, rc;

    // Read commands
    while (! ifs.eof()) {
        
        command_t comm;
       
        // Analyze next line
        ifs.getline (line, M_LINE);
        ind = -1;
        rc = sscanf (line, "%s %d\n", comm_str, &ind);
        if ((rc <= 0) || (strlen (comm_str) <= 0)) break;

        if (strcmp (comm_str, COMM_INSERT_STR) == 0) {
            // 'Insert' command
            if (rc < 2) return (0);
            if ((ind < 0) || (ind >= n_points)) return (0);
            comm.type = COMM_INSERT;
        }
        else if (strcmp (comm_str, COMM_DELETE_STR) == 0) {
            // 'Delete' command
            if (rc < 2) return (0);
            if ((ind < 0) || (ind >= n_points)) return (0);
            comm.type = COMM_DELETE;
        }
        else if (strcmp (comm_str, COMM_DRAW_STR) == 0) {
            // 'Draw' command
            comm.type = COMM_DRAW;
        }
        else {
            cerr << "Unknown command: '" << comm_str << "'.\n";
            return (0);
        }

        comm.pt_index = ind;
        commands.push_back (comm);
    
    }
    
    ifs.close();
    return (1);
}


//////////////////////////////////////////////////////////////
// Find the bounding-box (for display)
//////////////////////////////////////////////////////////////
void find_bbox (vector<Point> &points, vector<Polygon> &polygons,
                int radius, int *x, int *y, int *width)
{
    int i, x_min,x_max,y_min,y_max;
    
    x_min = (int) (points[0].x().to_double());
    x_max = x_min;
    y_min = (int) (points[0].y().to_double());
    y_max = y_min;

    for (i=0; i < (int) points.size(); i++) {
        x_min = MIN( x_min, (int) (points[i].x().to_double()) );
        x_max = MAX( x_max, (int) (points[i].x().to_double()) );
        y_min = MIN( y_min, (int) (points[i].y().to_double()) );
        y_max = MAX( y_max, (int) (points[i].y().to_double()) );
    }
    
    for (i=0; i < (int) polygons.size(); i++) {
        x_min = MIN( x_min, 
                     (int) (polygons[i].left_vertex()->x().to_double()) );
        x_max = MAX( x_max, 
                     (int) (polygons[i].right_vertex()->x().to_double()) );
        y_min = MIN( y_min, 
                     (int) (polygons[i].bottom_vertex()->y().to_double()) );
        y_max = MAX( y_max, 
                     (int) (polygons[i].top_vertex()->y().to_double()) );
    }
    
    *x = x_min - 2*radius;
    *y = y_min - 2*radius;
    *width = MAX( (x_max-x_min), (y_max-y_min) ) + 4*radius;
}


//////////////////////////////////////////////////////////////
// Solve the equation: ax^2 + bx + c = 0.
//////////////////////////////////////////////////////////////
void solve_equ2 (NT &a, NT &b, NT &c, 
                 NT &min_sol, NT &max_sol)
{
    NT s1, s2;
    
    if (a == 0) {
        if (b == 0) {
            cerr << "ERROR: equation has no solution (a=b=0)!\n";
            exit (1);
        }
        s1 = - (c / b);
        s2 = s1;
    }
    else {
        NT det = b*b - 4*a*c;
        if (det < 0) {
            cerr << "ERROR: equation has no solution (det<0)!\n";
            exit (1);
        }
        NT  sdet = sqrt (det);
        s1 = (-b + sdet)/(2*a);
        s2 = (-b - sdet)/(2*a);
    }
    min_sol = MIN (s1, s2);
    max_sol = MAX (s1, s2);
}

    
//////////////////////////////////////////////////////////////
// Compute solution point from combinatorial data
//////////////////////////////////////////////////////////////
int compute_solution (vector<Point> &points, 
                      vector<Polygon> &polygons,
                      NT &radius,
                      solution_t &solution,
                      double *x, double *y)
{
    if (solution.type == SOL_NONE) {

        // No solution
        return (0);
        
    }
    else if (solution.type == SOL_POINT) {
        
        // Solution determined by one point
        CHECK_RANGE(solution.point1, points.size());
        *x = points[solution.point1].x().to_double() -
             radius.to_double();
        *y = points[solution.point1].y().to_double();
        
    }
    else if (solution.type == SOL_POINT_POINT) {
        
        // Solution determined by two points
        CHECK_RANGE(solution.point1, points.size());
        CHECK_RANGE(solution.point2, points.size());
        NT  x1 = points[solution.point1].x(),
            y1 = points[solution.point1].y(),
            x2 = points[solution.point2].x(),
            y2 = points[solution.point2].y();
        NT  C = x2 - x1,
            D = y1 - y2,
            A = (x1*x1 - x2*x2) - (y2*y2 - y1*y1);
        NT  a = 4 * (C*C + D*D),
            b = 4*A*C - 8*(C*D*y1 + D*D*x1),
            c = 4*D*D*(x1*x1+y1*y1-radius*radius) - 4*A*D*y1 + A*A;
        NT  s1, s2;
        solve_equ2 (a, b, c, s1, s2);
        NT sy;
        if (D != 0)
            sy = (C/D)*s1 + (A/(2*D));
        else {
            a = 1;
            b = -2 * y1;
            c = y1 * y1;
            NT sy2;
            solve_equ2 (a, b, c, sy, sy2);
        }
        *x = s1.to_double();
        *y = sy.to_double();
        
    }
    else if (solution.type == SOL_POINT_EDGE) {
        
        // Solution determined by a point and an edge
        int poly = solution.polygon;
        CHECK_RANGE(solution.point1, points.size());
        CHECK_RANGE(poly           , polygons.size());
        CHECK_RANGE(solution.edge  , polygons[poly].size());
        NT  x1 = points[solution.point1].x(),
            y1 = points[solution.point1].y();
        VertexIterator  v, v_b, v_e;
        int i;
        v_b = polygons[poly].vertices_begin();
        v_e = polygons[poly].vertices_end();
        for (v=v_b, i=0; i < solution.edge; i++)
            v++;
        NT  v1x = v->x(), v1y = v->y();
        NT  v2x, v2y;
        v++;
        if (v == v_e) {
            v2x = polygons[poly].vertices_begin()->x(); 
            v2y = polygons[poly].vertices_begin()->y();
        }
        else {
            v2x = v->x();
            v2y = v->y();
        }
        if (v1x == v2x) {
            // vertical edge
            NT  a = 1,
                b = -2 * y1,
                c = y1*y1 + (v1x-x1)*(v1x-x1) - radius*radius;
            NT  s1, s2;
            solve_equ2 (a, b, c, s1, s2);
            if ((s1 >= MIN(v1y,v2y)) && (s1 <= MAX(v1y,v2y)))
                *y = s1.to_double();
            else if ((s2 >= MIN(v1y,v2y)) && (s2 <= MAX(v1y,v2y)))
                *y = s2.to_double();
            else {
                cerr << "ERROR: Edge (" << solution.edge << ") "
                     << "and point (" << solution.point1 << ") "
                     << "don't define a solution!\n";
                exit (1);
            }
            *x = v1x.to_double();
        }
        else {
            // non-vertical edge
            NT  M = (v2y - v1y) / (v2x - v1x);
            NT  B = v1y - (M * v1x);
            NT  a = 1 + M*M,
                b = 2 * ( (B-y1)*M - x1 ),
                c = x1*x1 + y1*y1 - radius*radius + B*B - 2*B*y1;
            NT  s1, s2;
            solve_equ2 (a, b, c, s1, s2);
            if ((s1 >= MIN(v1x,v2x)) && (s1 <= MAX(v1x,v2x)))
                *x = s1.to_double();
            else if ((s2 >= MIN(v1x,v2x)) && (s2 <= MAX(v1x,v2x)))
                *x = s2.to_double();
            else {
                cerr << "ERROR: Edge (" << solution.edge << ") "
                     << "and point (" << solution.point1 << ") "
                     << "don't define a solution!\n";
                exit (1);
            }
            NT sy = M * (*x) + B;
            *y = sy.to_double();
        }
        
    }
    else if (solution.type == SOL_VERTEX) {
        
        // Solution determined by a vertex
        int poly = solution.polygon;
        CHECK_RANGE(poly           , polygons.size());
        CHECK_RANGE(solution.vertex, polygons[poly].size());
        VertexIterator  v;
        int i;
        v = polygons[poly].vertices_begin();
        for (i=0; i < solution.vertex; i++)
            v++;
        *x = v->x().to_double();
        *y = v->y().to_double();
        
    }
    else {
        
        cerr << "Unknown type of solution: " << solution.type << "\n";
        exit (1);
            
    }
        
    return (1);
}


//////////////////////////////////////////////////////////////
// Print combinatorial solution to output stream
//////////////////////////////////////////////////////////////
ostream& operator<< (ostream& os, solution_t &solution)
{
    if (solution.type == SOL_NONE)
        os << "NONE";
    else if (solution.type == SOL_POINT)
        os << "Point " << solution.point1;
    else if (solution.type == SOL_POINT_POINT)
        os << "Point " << solution.point1 
           << ", Point " << solution.point2;
    else if (solution.type == SOL_POINT_EDGE)
        os << "Point " << solution.point1 
           << ", Edge " << solution.edge
           << " of polygon " << solution.polygon;
    else if (solution.type == SOL_VERTEX)
        os << "Vertex " << solution.vertex 
           << " of polygon " << solution.polygon;
    else
        os << "Unknown!";
    
    return os;
}


//////////////////////////////////////////////////////////////
// Main
//////////////////////////////////////////////////////////////
int main(int argc, char* argv[])
{  
    vector<command_t>   commands;
    vector<Point>       points;
    vector<Polygon>     polygons;
    vector<bool>        active;
    vector<solution_t>  solutions;
    solution_t          last_solution;
    Algorithm           alg;
    int    i, x_min, y_min, width;
    int    radius_i;
    NT     radius;  
    double x, y;
    CGAL_Timer tm;

  
    // Analyze command-line parameters
    if ((argc < 4) || (argc > 5)) {
        cout << "Usage: " << argv[0] 
             << " radius input-file command-file [output-file]\n";
        exit(-1);
    }
    
    radius_i = atoi (argv[1]);
    if (radius_i <= 0) {
        cerr << "Radius must be a positive integer.\n";
        exit (1);
    }
    radius = NT(radius_i);
    
    cout << "Radius is " << radius_i << ".\n";
                                                   
    // Read points and polygons from input file
    if (! read_input_file (argv[2], points, polygons)) {
        cerr << "Problems reading input file.\n";
        exit (1);
    }
    
    cout << "Read " << points.size() << " points and " 
         << polygons.size() << " polygons from input file.\n";
    
    // Initialize all points as "active"
    for (i=0; i < (int) points.size(); i++)
        active.push_back (true);
    
    // Read commands from commands file
    if (! read_commands_file (argv[3], points.size(), commands)) {
        cerr << "Problems reading commands file.\n";
        exit (1);
    }
    
    cout << "Read " << commands.size() << " commands.\n";

    // Find bounding-box (for display purposes)
    find_bbox (points, polygons, radius_i, &x_min, &y_min, &width);
    
    // Draw points and polygons
    CGAL_Window_stream  W(WINDOW_SIZE, WINDOW_SIZE);
    W.init (x_min, x_min+width, y_min);
    W.set_node_width (1);
    W.display ();
    draw (W, points, active, polygons);

    // Part 1
    // ======
    // Initialize algorithm and compute the solution for
    // all the points
    tm.start();
    alg.init (points, polygons, radius, last_solution);
    tm.stop();

    cout << "Initialization time: " << tm.time() << "\n";
    solutions.push_back (last_solution);

    // Part 2
    // ======
    // Execute the commands, one by one
    tm.start();
    for (i=0; i < (int) commands.size(); i++) {
        
        if (commands[i].type == COMM_DRAW) {
            // 'Draw' command - draw the last output point
            tm.stop();
            draw (W, points, active, polygons);
            if (compute_solution (points, polygons, radius,
                                  last_solution, &x, &y)) {
                W << CGAL_GREEN << Point(x,y);
                W.draw_circle (x, y, radius_i);
            }
            any_button (W);
            tm.start();
        }
        else {
            // 'Insert'/'Delete' command - send to algorithm
            alg.update (commands[i], last_solution);
    
            // Update active points vector
            if (commands[i].type == COMM_INSERT)
                active[commands[i].pt_index] = true;
            else
                active[commands[i].pt_index] = false;
            
            solutions.push_back (last_solution);
        }
    }
    tm.stop();
    
    cout << "Total time: " << tm.time() << "\n";

    
    // Write solutions to output file, if requested
    if (argc > 4) {
        ofstream ofs(argv[4]);
        for (i=0; i < (int) solutions.size(); i++)
            ofs << i << ": " << solutions[i] << "\n";
        ofs.close();
    }
    
    return (0);
}