#ifndef __STREAM_SURFACES__
#define __STREAM_SURFACES__

// ========================================================================== //
// Forward declarations
class Options;
class PDS;
class FrameField;

namespace CGLA {
    class Vec3d;
};

// Local libraries
#include "Utilities/BBoxMap.h"
#include "Pointcloud/Pointcloud.hpp"

// External libraries
#include <GEL/HMesh/Manifold.h>

// Standard libraries
#include <random>

// ========================================================================== //
// Start of the class

/** Stream surface class.
 *
 * This class is designed to handle the stream surface creation. This use
 * containers for frame fields, currently we support axis aligned grid of
 * frames.
 *
 * The Surface Creation pipeline is to run the following commands:
 *
 * - Construct with an Options file. Need to contain 'data_name' option.
 * - Creation surfaces through the function create_surfaces.
 * - Compute the probe activation through compute_indicator.
 * - Export the constructed system using the save function.
 *
 * Then the external Matlab system can be used to run a subselection
 * optimisation which will choose which surfaces should be included in the
 * optimal set.
 *
 * The Super Sampling pipeline can now be conducted:
 *
 * - Define the option super_sampling to be true.
 * - Use the import_selected to use the known set of surfaces.
 * - Run create_surfaces which will now super sample surfaces.
 * - Use the triangle_converter to obtain meshes if desired.
 * - Export the surfaces using the dedicated function export_surfaces.
 *
 * See the list of "Protected Attributes" to see modifiable Options.
 */
class StreamSurfaces {

  protected:
    // ====================================================================== //
    // User Options with default values, modified through options file

    // ...................................................................... //
    // Specification of input names and structure

    std::string data_name     = "";               // Name of the input data
    std::string data_location = "../../../Data/"; // Location of the input data

    // ...................................................................... //
    // Meta settings

    bool   debug           = false; // Show debugging information
    double tolerated_error = 0.2;   // The tolerated error (eps)
    double desired_spacing = 10;    // The desired spacing (Ns)

    // ...................................................................... //
    // Generation of surfaces

    double surf_pds_r      = 1.0;     // PDS radius fraction of diag
    int    surf_pds_k      = 30;      // Maximal PDS attempts
    int    surf_min_number = 0;       // Minimal number of surfaces
    int    surf_max_number = INT_MAX; // Maximal number of surfaces

    // ...................................................................... //
    // Toggles for additional features

    bool singularity_toggle = false; // Toggle for singularities
    bool include_hull       = false; // Include the hull

    // ...................................................................... //
    // Topology optimization settings

    double threshold_void  = -1e16; // Thickness clipping threshold
    double threshold_solid = 1e16;  // Density clipping threshold

    // ...................................................................... //
    // Super Sampling and Post-processing

    bool   super_sampling     = false; // Super sample toggle
    double super_sample_pds_r = 0.5;   // Super sample scaling value
    int    smooth_max_iter    = 0;     // Maximal number of smoothing steps

    // ...................................................................... //
    // Input specifications

    std::string      seed_method      = "pds";             // Seed placement
    std::vector<int> seed_grid_amount = {10, 10, 10};      // Grid size
    double           seed_pds_r       = 0.05;              // PDS radius
    std::string      seed_file_type   = "index";           // File type
    std::string      seed_file_name   = "seed_points.csv"; // Filename

    // ...................................................................... //
    // Computation of voxel indicator response

    std::string      probe_method    = "grid";             // Method
    std::string      probe_file_name = "probe_points.csv"; // Filename
    std::vector<int> probe_grid_max  = {100, 100, 100};    // Max grid

    // ...................................................................... //
    // Export Settings

    std::string export_name     = "result";  // Export base file name
    std::string export_location = "result/"; // Export file locations

  public:
    // ====================================================================== //
    // Public API

    /** Construct a new Stream Surfaces object.
     *
     * The Options class is used to specify input data and settings.
     *
     */
    StreamSurfaces(Options &options);
    /** Construct a new Stream Surfaces object.
     *
     * Construct a StreamSurfaces object with Options specified directly from a
     * file.
     *
     *
     * @param options_filename default is "options.cfg".
     */
    StreamSurfaces(std::string options_filename = "options.cfg");
    /** Destructor of the class. */
    ~StreamSurfaces(){};

    /** Import surfaces which have been selected by optimisation.
     *
     * This function is used to import the StreamSurfaces selected by the
     * Optimiser step of the pipeline.
     *
     */
    void import_selected();

    /**
     *
     */
    void prepare_fixed();

    /** Main surface creation.
     *
     * This function will start the generation process of StreamSurfaces.
     */
    void create_surfaces();

    /** Smooth all created surfaces.
     *
     * This function will run through and smooth all the generated surfaces
     */
    void smooth();

    /** PointCloud cleanup.
     *
     * @deprecated old function not maintained atm.
     */
    void clean_surfaces();

    /** Compute indicator response.
     *
     * This function will compute the indicator response. Each probe is
     * activated if a surface point is nearby. The frame index is then defined
     * by looking up the local frame at the probe point, and matching the
     * surface normal to this frame.
     */
    void compute_indicator();

    /** Convert to triangle meshes.
     *
     * Go through all generated surfaces and run a simple meshing algorithm to
     * convert these to triangle meshes.
     *
     * The algorithm is based on local 2D Delaunay Triangulation and a point
     * based voting scheme to eleminate problems at the boundary of local
     * patches.
     */
    void triangle_converter();

    /** Run all the exports.
     *
     * This is the main function which run all exports:
     *
     * - Export of full surface list.
     * - Export of seed point list.
     * - Export of probe point list.
     * - Export of indicator response.
     *
     */
    void save();

    /** Export all generated surfaces.
     *
     * This function will export to the folder designated through the options
     * framework.
     *
     * Default folder: result/
     *
     */
    void export_surf();
    /** Export all generated surfaces.
     *
     * This function will export all surfaces to a userdefined folder specified
     * directly.
     *
     * @param folder name.
     */
    void export_surf(std::string folder);

    // ====================================================================== //
    // Public Printing functions

    /** Print list of surface sizes. */
    void print_surfaces();

  private:
    // ====================================================================== //
    // Internal data fields

    /** Structure designed to represent a stream surface. */
    struct StreamSurface {
        PointCloud      pcl  = PointCloud();
        HMesh::Manifold mesh = HMesh::Manifold();

        StreamSurface(){};
        StreamSurface(PointCloud input_pcl) : pcl(input_pcl){};

        size_t size() { return pcl.size(); };
    };

    // Design domain variables
    int    dimX, dimY, dimZ; // Number of elements in each dimension.
    int    dimS;             // Number of frame vectors.
    BBox   bbox;             // Bounding Box structure.
    double bbox_diag;        // Size of the bounding box.
    int    desired_surfaces; // Desired number of surfaces.

    // Frame variables
    FrameField *frame;

    // Probe point variables
    std::vector<CGLA::Vec3d> probe_points;       // List of probes
    std::vector<int>         indicator_response; // Indicator response list

    // Pointcloud variables
    std::vector<StreamSurface> surfaces; // Listing the surfaces

    std::vector<int> fixed; // List of indices which are fixed.

  private:
    // ====================================================================== //
    // Private IO functions

    /** Function designed for internal construction.
     *
     * This is the construction function used by all constructors.
     *
     * @param options object.
     */
    void constructor(Options &options);

    /** Evaluate options from the object. */
    void read_options(Options &options);
    /** Read information from data file. */
    void read_information();
    /** Construct seeds based on the options. */
    void read_input();
    /** Construct probe points based on options. */
    void read_probe();

    /** Import all mesh matching the specified pattern. */
    void import_surf(std::string pattern);

    /** Export of a single surface.
     *
     * This function export a given surface to a specified position. Each part
     * the structure is exported as a ply PointCloud, wavefront mesh and an
     * internal specification of the implicit.
     *
     * @param index of the stream surface which should be exported.
     * @param folder of where the export should be placed.
     */
    void export_surf(int index, std::string folder);
    /** Ply file manipulation.
     *
     * @param fileName
     * @param mesh
     */
    void obj_save(std::string fileName, HMesh::Manifold &mesh);

    /** Export the list of activated probes in a sparse format.*/
    void export_adjacency();
    /** Export the list of probe locations.*/
    void export_probes();
    /** Export list of the seed locations and associated directions.*/
    void export_inputs();

    // ====================================================================== //
    // Single Surface manipulation

    /** Creation of a single surface.
     *
     * @param surf
     * @param pds_r
     * @param rand_dev
     */
    void single_surface(StreamSurface &surf, const double pds_r,
                        std::mt19937 &rand_dev);
    /** Smoothing of single surface.
     *
     * This internal function will run the smoothing of a single surface.
     * Smoothing is done by re-estimating all point positions from the
     * underlying frame field. This mean we propergate information in all
     * direction instead of just an expanding front.
     *
     * @param surf
     * @param pds
     */
    void smooth_single(StreamSurface &surf, PDS &pds);
    /** Smoothing of single surface.
     *
     * This internal function will run the smoothing of a single surface.
     * Smoothing is done by re-estimating all point positions from the
     * underlying frame field. This mean we propergate information in all
     * direction instead of just an expanding front.
     *
     * @param surf
     * @param pds_r
     */
    void smooth_single(StreamSurface &surf, const double pds_r);

    /** Computation of the indicator response for a single surface.
     *
     * This function will compute the indicator response of a single stream
     * surface. A KDTree is build for all point in the surface. For each probe
     * point we locate the nearest surface point. If it is close enough we
     * interpolate the field at the probe point, and match the surface normal to
     * the field. The index of the closest aligned field is stored in the
     * indicator response list.
     *
     * @param surf StreamSurface in question.
     * @param id of the current stream surface.
     */
    void indicator(StreamSurface &surf, const int id);

    // ====================================================================== //
    // Pre processor

    /** Extract the outer hull for the topology optimization examples.
     * @return HMesh::Manifold
     */
    HMesh::Manifold isocontour_hull();

    // ====================================================================== //
    // Mathematical operations

    /**
     *
     * @param p0
     * @param step_size
     * @param d0
     * @param norm0
     * @return CGLA::Vec3d
     */
    CGLA::Vec3d runge_kutta_4(CGLA::Vec3d &p0, double step_size,
                              CGLA::Vec3d &d0, CGLA::Vec3d &norm0);
    /**
     *
     * @param p0
     * @param pds
     * @param surf
     * @param R
     * @param step_size
     * @return CGLA::Vec3d
     */
    CGLA::Vec3d estimate_position(CGLA::Vec3d &p0, PDS &pds,
                                  StreamSurface &surf, double R,
                                  double step_size);
    /**
     *
     * @param p0
     * @param pds
     * @param surf
     * @param R
     * @param step_size
     * @return CGLA::Vec3d
     */
    CGLA::Vec3d estimate_normal(CGLA::Vec3d &p0, PDS &pds, StreamSurface &surf,
                                double R, double step_size);

    // ====================================================================== //
    // Boolean Operations

    /**
     *
     * @param p
     * @param norm
     * @param pds
     * @param surf
     * @param dist
     * @return true
     * @return false
     */
    bool point_is_valid(CGLA::Vec3d &p, CGLA::Vec3d &norm, PDS &pds,
                        StreamSurface &surf, double dist);

    /**
     *
     * @param idx
     * @return true if
     * @return false if
     */
    bool is_fixed(int idx);
    /**
     *
     * @param p
     * @return true
     * @return false
     */
    bool inside_domain(CGLA::Vec3d &p);

    HMesh::Manifold triangulate(PointCloud &pcl);
};

#endif /* __STREAM_SURFACES__ */