/*
 * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
 * Applied Mathematics, Norway.
 *
 * Contact information: E-mail: tor.dokken@sintef.no                      
 * SINTEF ICT, Department of Applied Mathematics,                         
 * P.O. Box 124 Blindern,                                                 
 * 0314 Oslo, Norway.                                                     
 *
 * This file is part of SISL.
 *
 * SISL is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Affero General Public License as
 * published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version. 
 *
 * SISL is distributed in the hope that it will be useful,        
 * but WITHOUT ANY WARRANTY; without even the implied warranty of         
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          
 * GNU Affero General Public License for more details.
 *
 * You should have received a copy of the GNU Affero General Public
 * License along with SISL. If not, see
 * <http://www.gnu.org/licenses/>.
 *
 * In accordance with Section 7(b) of the GNU Affero General Public
 * License, a covered work must retain the producer line in every data
 * file that is created or manipulated using SISL.
 *
 * Other Usage
 * You can be released from the requirements of the license by purchasing
 * a commercial license. Buying such a license is mandatory as soon as you
 * develop commercial activities involving the SISL library without
 * disclosing the source code of your own applications.
 *
 * This file may be used in accordance with the terms contained in a
 * written agreement between you and SINTEF ICT. 
 */

#include <limits>
#include <iostream>
#include <fstream>
#include <string>
#include <stdexcept>
#include <time.h>
#include <cstdlib>

#include "sisl.h"
#include "GoReadWrite.h"

using namespace std;

namespace {
    string OUT_FILE_POINTS = "example15_points.g2";
    string IN_FILE_SURFACE = "example10_surf.g2";

    string DESCRIPTION = 
    //==========================================================
    "This program will do a simple ray-tracing of a bspline \n"
    "surface.  From one point P in space, we will extend lines \n"
    "('rays') towards the surface, and calculate the point where \n"
    "the line intersects with the surface.  All points obtained \n"
    "in this way thus describe how the surface is 'visible' from \n"
    "P.  They will be stored to the file '" + OUT_FILE_POINTS + "'\n"
    "The surface used is the one generated by 'example10', namely \n"
    + IN_FILE_SURFACE + "'.\n\n"
    "Note that in this program, two ray-tracing routines are \n"
    "provided; (1) one slow and robust one, which would always \n"
    "work, and one quick but more fragile one (2), which bases \n"
    "its point estimates on previously obtained results, and \n"
    "thus is likely to fail when some parts of a surface shadows \n"
    "other parts.  The default method is (1); however the user \n"
    "can invoke the quick method by running the program with the\n"
    "additional argument 'q' (quick).  An examination of the \n"
    "obtained point cloud will in this case reveal erroneous \n"
    "behaviour where the examined surface 'folds over' itself.\n"
    "The user will notice, though, that the execution time is \n"
    "significantly shorter.  This method is therefore useful in\n"
    "some cases where surfaces are not self-obstructing.\n\n";
    //==========================================================
 
    double P[] = {0, 0, -30};

    // definition of view window
    const double WIN_LL[] = {-5, -5, -10}; // lower left corner
    const double WIN_LR[] = { 5, -5, -10}; // lower right corner
    const double WIN_UL[] = {-5,  5, -10}; // upper left corner
    const double WIN_UR[] = { 5,  5, -10}; // upper right corner
    
    const int HRES = 200;
    const int VRES = 200;

    const double EPSGE = 1.0e-5;  // geometrical tolerance
    const double EPSCO = 1.0e-15; // computational tolerance

    void raytracer_simple(SISLSurf* surf, vector<double>& points);
    void raytracer_advanced(SISLSurf* surf, vector<double>& points);
    void calculate_direction(int vscan, int hscan, double* dir);
    bool trace_ray_robust(SISLSurf* surf, double* param_val, double* dir);
    bool trace_ray_rapid(SISLSurf* surf, double* init_guess, double* dir) ;

    bool USE_LOCAL_METHOD = false;
}; // end anonymous namespace 

//===========================================================================
int main(int varnum, char** vararg)
//===========================================================================
{
    cout << '\n' << vararg[0] << ":\n" << DESCRIPTION << endl;
    cout << "To proceed, press enter, or ^C to quit." << endl;
    getchar();

    try {
	// parsing command line argument
	if (varnum > 1) {
	    string arg1(vararg[1]);
	    if (arg1 == "q" || arg1 == "Q") {
		USE_LOCAL_METHOD = true;
	    } else {
		cerr << "Invalid command argument '" << arg1 << "' specified.\n";
		return 0;
	    }
	}

	ifstream is(IN_FILE_SURFACE.c_str());
	if (!is) {
	    throw runtime_error("Unable to open input file.");
	}
	ofstream os(OUT_FILE_POINTS.c_str());
	if (!os) {
	    throw runtime_error("Unable to open output file.");
	}

	// read surface to be raytraced
	SISLSurf* surf = readGoSurface(is);
	vector<double> points;

	// the complete raytracing is carried out here
	clock_t starttime = clock();
	if ( ! USE_LOCAL_METHOD ) {
	    raytracer_simple(surf, points);
	} else {
	    raytracer_advanced(surf, points);
	}
	clock_t endtime = clock();
	cout << "Time used: " << (endtime - starttime) / double(1000) << endl;

	// saving to file
	writeGoPoints(points.size() / 3, &points[0], os);

	// cleaning up
	freeSurf(surf);
	is.close();
	os.close();

    } catch (exception& e) {
	cerr << "Exception thrown: " << e.what() << endl;
	return 0;
    }

    return 1;
};

namespace {

void raytracer_simple(SISLSurf* surf, vector<double>& points)
{
    // a simple, slow and robust raytracer.  It scans the view
    // line by line, and find all intersections with the
    // current line and the surface.  When several intersections are found,
    // the one closest to the viewpoint (P) is chosen.
    points.clear();
    double dir[3];
    double param_val[2];
    for (int vscan = 0; vscan < VRES; ++vscan) {
	for (int hscan = 0; hscan < HRES; ++hscan) {
	    calculate_direction(vscan, hscan, dir);
	    bool found = trace_ray_robust(surf, param_val, dir);
	    if (found) {
		points.insert(points.end(), &dir[0], &dir[3]);
	    }
	}
    }
}

void raytracer_advanced(SISLSurf* surf, vector<double>& points)
{
    // A relatively quick raytracer, scanning the view line by line,
    // and using previously found points as estimates for detecting new
    // ones, thus speeding up the location of intersections.  However,
    // this is likely to give rise to points on hidden parts of the surface
    // if the surface 'hides' parts of itself as seen from the vieweing
    // point.  To avoid this problem, either the 'raytracer_simple' routine
    // should be used, or this routine should be elaborated to take into
    // account non-local properties of the surface.
    points.clear();
    double dir[3];
    double param_val[2];
    param_val[0] = surf->et1[surf->ik1 - 1] - 1;
    param_val[1] = surf->et2[surf->ik2 - 1] - 1;
    bool hscan_forwards = true;
    int start, stop, inc;
    bool found_last = false;

    for (int vscan = 0; vscan < VRES; ++vscan) {
	start = hscan_forwards ? 0 : HRES;
	stop  = hscan_forwards ? HRES - 1 : -1;
	inc   = hscan_forwards ? 1 : -1;
	hscan_forwards = !hscan_forwards;

	for (int hscan = start; hscan != stop; hscan += inc) {
	    calculate_direction(vscan, hscan, dir);
	    if (found_last) {
		found_last = trace_ray_rapid(surf, param_val, dir);
	    } else {
		found_last = trace_ray_robust(surf, param_val, dir);
	    }
 	    if (found_last) {
 		points.insert(points.end(), &dir[0], &dir[3]);
	    } 
	}
    }
}

void calculate_direction(int vscan, int hscan, double* dir)
{
    // find the direction of the 'line-of-sight', as seen from the point 
    // P through the 'window' specified by WIN_LL, WIN_LR, WIN_UL and WIN_UR.
    double v_ratio = double(vscan) / double(VRES - 1);
    double h_ratio = double(hscan) / double(HRES - 1);
    for (int d = 0; d < 3; ++d) {
	dir[d] = (1-h_ratio) * ((1-v_ratio) * WIN_LL[d] + v_ratio * WIN_LR[d]) +
	    h_ratio * ((1-v_ratio) * WIN_UL[d] + v_ratio * WIN_UR[d]) - P[d];
    }
}

bool trace_ray_rapid(SISLSurf* surf, double* init_guess, double* dir) 
{
    // Rapidly find the intersection of a ray with the surface, provided 
    // we have a good initial guess for the parameters in that point ('init_guess').
    // If this is not the case, the routine will invoke the slower
    // 'trace_ray_robust' routine.

    // caching of parameter domain in order to save time
    static double start_dom[2];
    static double end_dom[2];
    static SISLSurf* cached_surf = 0;
    if (surf != cached_surf) {
	// caching parameter domain
	start_dom[0] = surf->et1[surf->ik1 - 1];
	start_dom[1] = surf->et2[surf->ik2 - 1];
	end_dom[0] = surf->et1[surf->in1];
	end_dom[1] = surf->et2[surf->in2];
	cached_surf = surf;
    }

    if (init_guess[0] < start_dom[0]) {
	// there is no (valid) initial guess available
	return trace_ray_robust(surf, init_guess, dir);
    }

    double isect_uv[2];
    int jstat = 0;

    // Using local properties of the surface to converge to the solution point, 
    // provided we have a 'good' initial guess.

    s1518(surf,        // the surface to 'raytrace'
	  P,           // observer's position
	  dir,         // direction of ray 
	  EPSGE,       // geometric resolution
	  start_dom,   // lower limit of search rectangle (umin, vmin)
	  end_dom,     // upper limit of search rectangle (umax, vmax)
	  init_guess,  // initial guess - we can use the last found value
	  isect_uv,    // the intersection point - if found!
	  &jstat);     // status
    
    if (jstat < 0) {
	cout << jstat << endl;
	throw runtime_error("Error occurred in SISL routine s1518.");
    } else if (jstat != 1) {
	// no point found.  Resorting to robust technique
	return trace_ray_robust(surf, init_guess, dir);
    } 
    // we found a point.  Checking quality of result.
    
    //  do not want a point that is in reality hidden behind another, closer point.
    // If our initial guess was good, this would not happen all too often, but once
    // in a while, we can have 'marched across' a visible silhouette.  In order to check
    // if we might have done such a thing, we will verify that the scalar product of the 
    // normal in the previous point with the direction 'dir' has the same sign as the
    // scalar product of the normal in the calculated point with 'dir'.  If this is 
    // NOT the case, we will resort to the robust method
    
    // calculating normal for the initial guess point
    int temp1, temp2;
    double point[9]; // point and tangents in u and v parameter will be stored here
    double n[3]; // normal will be stored here
    
    s1421(surf, // the surface in question
	  1,    // we need at least one derivative in order to calculate normal
	  init_guess, // parameter value
	  &temp1,   // uninteresting for our purposes - parameter interval
	  &temp2,   // uninteresting for our purposes - parameter interval
	  point,   // storage for the point and its derivatives
	  n,       // storage for the normal
	  &jstat);

    if (jstat < 0) {
	throw runtime_error("Error occurred in SISL routine s1421.");
    } else if (jstat > 0) {
	    cerr << "WARNING: warning occured inside call to SISL routine s1857. \n" 
		 << endl;
    }
    double scalar_product_1 = n[0] * dir[0] + n[1] * dir[1] + n[2] * dir[2];
    
    // calculating normal for the found point
    s1421(surf, // the surface in question
	  1,    // we need at least one derivative in order to calculate normal
	  isect_uv, // parameter value
	  &temp1,   // uninteresting for our purposes - parameter interval
	  &temp2,   // uninteresting for our purposes - parameter interval
	  point,   // storage for the point and its derivatives
	  n,       // storage for the normal
	  &jstat);
    if (jstat < 0) {
	throw runtime_error("Error occurred in SISL routine s1421.");
    } else if (jstat > 0) {
	    cerr << "WARNING: warning occured inside call to SISL routine s1857. \n" 
		 << endl;
    }
    double scalar_product_2 = n[0] * dir[0] + n[1] * dir[1] + n[2] * dir[2];
    if (scalar_product_1 * scalar_product_2 < 0) {
	// the normals are pointing in distinct directions
	return trace_ray_robust(surf, init_guess, dir);
    }
    // if we got here, the point found in the 'rapid' way should be completely OK
    copy(&point[0], &point[3], dir);
    copy(&isect_uv[0], &isect_uv[2], init_guess);
    return true;
}



bool trace_ray_robust(SISLSurf* surf, double* param_val, double* dir) 
{
    // Find the intersection of a ray with a surface that is closest to the 
    // viewpoint P.  We do not depend on any initial guess for this point; the
    // routine does a global search.  If several intersection points are found
    // for the same ray, the one closest to the viewpoint P is chosen.

    int num_intpoints = 0;
    double* pointpar = 0;
    int num_intcurves;
    int jstat;
    SISLIntcurve** intcurves = 0;

    // intersect a surface with a line
    s1856(surf,    // the surface to intersect with
	  P,       // a point on the line (in our case: the viewpoint)
	  dir,     // the directional vector of the line
	  3,       // dimension of the Euclidean space
	  EPSCO,   // machine (computational) tolerance
	  EPSGE,   // geometric tolerance
	  &num_intpoints,  // reports the number of intersection points found
	  &pointpar,       // pointer to array containing the surface parameters for the 
	                   // found intersection points
	  &num_intcurves,  // returns the number of intersection curves found (usually none!)
	  &intcurves,      // pointer to an array containing the intersection curves found
	  &jstat);
    if (jstat < 0 ) {
	throw runtime_error("Error occured in SISL routine s1856.");
    } else if (jstat > 0) {
	cerr << "WARNING: warning occured inside call to SISL routine s1856."
	     << endl;
    }
    
    if (num_intpoints == 0) {
	return false; // no intersection with ray
    }
    // If we got here, there was at least one intersection with the ray.  We want to find
    // the closest one.
    double dist = numeric_limits<double>::max();
    int i;
    for (i = 0; i < num_intpoints; ++i) {
	int temp1, temp2;
	double temp_pos[3];
	// calculating 3D point
	s1424(surf,             // input surface
	      0,                // evaluate position only (no derivatives in u-param)
	      0,                // evaluate position only (no derivatives in v-param)
	      &pointpar[2 * i], // parameter values
	      &temp1,           // unused for our purposes (returns used u-interval)
	      &temp2,           // unused for our purposes (returns used v-interval)
	      temp_pos,         // the calculated position
	      &jstat);          // status variable 
	if (jstat < 0 ) {
	    throw runtime_error("Error occured in SISL routine s1424.");
	} else if (jstat > 0) {
	    cerr << "WARNING: warning occured inside call to SISL routine s1424."
		 << endl;
	}
	
	// calculating distance of 3D point
	double cur_dist = 0;
	for (int j = 0; j < 3; ++j) {
	    double tmp = (temp_pos[j] - P[j]);
	    tmp *= tmp;
	    cur_dist += tmp;
	}
	if (cur_dist < dist) {
	    dist = cur_dist;
	    param_val[0] = pointpar[2 * i];
	    param_val[1] = pointpar[2 * i + 1];
	    for (int k = 0; k < 3; ++k) {
		dir[k] = temp_pos[k];
	    }
	}
    }
    free(pointpar);
    for (i = 0; i < num_intcurves; ++i) {
	free(intcurves[i]);
    }
    free(intcurves);
    return true;
}


}; // end anonymous namespace