areacon.h

Go to the documentation of this file.

/********************************************/
#ifndef __AreaCon__areacon__
#define __AreaCon__areacon__

#include <stdio.h>
#include <cmath>
#include <algorithm>
#include <iostream>
#include <fstream>
#include "clipper.hpp"




namespace AreaCon {
    // Point Class-------------------------------------------------------------------------------------------------------
    // ------------------------------------------------------------------------------------------------------------------
    class Point
    {
    public:
        double x;
        double y;
        static double Robustness_Constant;
        Point(const double x = INFINITY, const double y = INFINITY);
        

        double Norm() const;
        Point AddPoint(const Point Test) const;
        Point FindPerpDirection(const Point Test, const double Norm) const;
        void FlipDirection(void);
        double PerpDistanceToLine(const Point Test1, const Point Test2) const;
        bool AreCollinear(const Point Test1, const Point Test2)const;
        bool AreBetween(const Point Test1, const Point Test2) const;
        void Mult(const double factor);
        

        static bool IsEqual(const Point Test1, const Point Test2);
        static double Distance(const Point Test1, const Point Test2);
        static Point FindPointAlongLine(const Point Test1, const Point Test2, const double distance);
        static double Norm(const Point Test);
        static Point AddPoints(const Point Test1, const Point Test2);
        static Point FindPerpDirection(const Point Test1, const Point Test2, const double Norm);
        static double PerpDistanceToLine(const Point Test1, const Point Test2, const Point Test3);
        static bool AreCollinear(const Point Test1, const Point Test2, const Point Test3);
        static bool AreBetween(const Point Test1, const Point Test2, const Point Test3);
        static std::vector<Point> FindCollinearIntersection(const Point p1, const Point p2, const Point p3, const Point p4);
        
    };
    // Poly Class-------------------------------------------------------------------------------------------------------
    // ------------------------------------------------------------------------------------------------------------------
    class Poly
    {
    public:

        Poly(std::vector<Point> Vertices = {});

        void SetVertices(const std::vector<Point> Vertices, const bool GetExtrema = true);
        std::vector<Point> GetVertices(void) const;
        int GetNVertices(void) const;
        bool pnpoly(const Point Test) const;
        void GetExtrema(double &minx, double &miny, double &maxx, double &maxy) const;
        
    private:
        double minx;
        double miny;
        double maxx;
        double maxy;
        std::vector<Point> Vertices;
        int NPoly;
        void InitializePoly(void);
    };
    // Mult_Array Class--------------------------------------------------------------------------------------------------
    // ------------------------------------------------------------------------------------------------------------------
    class Mult_Array
    {
    public:

        Mult_Array(const int N);
        Mult_Array(const Mult_Array &obj);
        Mult_Array& operator=(const Mult_Array &obj);
        ~Mult_Array(void);
        double **Array;
        const int N;
    };
    // Delaunay Class--------------------------------------------------------------------------------------------------
    // ------------------------------------------------------------------------------------------------------------------
    class DelaunayGraph
    {
    public:

        DelaunayGraph(const int NRegions);
        DelaunayGraph(const DelaunayGraph &obj);
        DelaunayGraph& operator=(const DelaunayGraph &obj);
        ~DelaunayGraph();
        
        Point ***Graph;
        const int NRegions;
    };
    // Int_Params Class--------------------------------------------------------------------------------------------------
    // ------------------------------------------------------------------------------------------------------------------
    class Int_Params
    {
    public:
        std::vector<double> Coefficient_a, Coefficient_b, Coefficient_c, Coefficient_d;
        std::vector<double> Int, Intx, Inty;
        double Unweighted_Area; 
        
        /*
         * Default Constructor.
         * @param[in] Coefficient_a A coefficient used in quickly calculating area integrals
         * @param[in] Coefficient_b A coefficient used in quickly calculating area integrals
         * @param[in] Coefficient_c A coefficient used in quickly calculating area integrals
         * @param[in] Coefficient_d A coefficient used in quickly calculating area integrals
         * @param[in] Int Vector representing the total integral over individual grid squares
         * @param[in] Intx Vector representing the integral of the function x*f(x,y) over individual grid squares
         * @param[in] Inty Vector representing the integral of the function y*f(x,y) over individual grid squares
         * @param[in] Unweighted_Area The overall area of some polygonal region
         */
        Int_Params(std::vector<double> Coefficient_a = {}, std::vector<double> Coefficient_b= {}, std::vector<double> Coefficient_c= {}, std::vector<double> Coefficient_d= {}, std::vector<double> Int = {}, std::vector<double> Intx = {}, std::vector<double> Inty = {}, double UnweightedArea = 0);

        void CheckParameters(void) const;
    };
    // Parameters Class--------------------------------------------------------------------------------------------------
    // ------------------------------------------------------------------------------------------------------------------
    class Parameters
    {
    public:

        Parameters(const double line_int_step = 0.1,const double weights_step = 0.1, const double centers_step = 1,const double volume_tolerance = 0.002, const double convergence_criterion = 0.02, const int max_iterations_volume = 200, const int max_iterations_centers = 500, const double Volume_Lower_Bound = 10e-6, const double Robustness_Constant = 10e-8);

        const double line_int_step;
        const double weights_step; 
        const double centers_step; 
        const double volume_tolerance; 
        const double convergence_criterion;
        const int max_iterations_volume;
        const int max_iterations_centers;
        const double Volume_Lower_Bound;
        const double Robustness_Constant;
    private:
        void CheckParameters(void);
    };
    
    // Density Class--------------------------------------------------------------------------------------------------
    // ------------------------------------------------------------------------------------------------------------------
    class Density
    {
    public:

        Density();
        Density(const Poly Region, const int Nx = 0, const int Ny = 0, const std::vector<double> Values = {});

        void SetNewRegion(const Poly Region, const int Nx = 0, const int Ny = 0, const std::vector<double> Values = {});
        void SetParameters(const int Nx,const int Ny,const std::vector<double> Values);
        int GetNx(void) const {return Nx;};
        int GetNy(void) const {return Ny;};;
        Poly GetRegion(void) const {return Region;};;
        std::vector<bool> GetGridInRegion(void) const {return GridInRegion;};
        Int_Params GetIntegral(void) const {return Integral;};
        void GetExtrema (double &minx, double &miny, double &maxx, double &maxy) const {minx = this->minx; miny = this->miny; maxx = this->maxx; maxy = this->maxy;}
        double LineIntegral(double spacing, const Point &p1, const Point &p2) const;
        double CalculateWeightedArea(const Poly Region) const;
        Point CalculateCentroid(const Poly Region, const double &Volume) const;
        void SetVolumeLowerBound(const double VolumeLowerBound);
        double GetVolumeLowerBound(void);
    private:
        Poly Region;
        int Nx;
        int Ny;
        double dx;
        double dy;
        double minx;
        double miny;
        double maxx;
        double maxy;
        double Volume_Lower_Bound;
        std::vector<double> Values;
        std::vector<bool> GridInRegion;
        Int_Params Integral;
        void CheckParameterSizes(void);
        void SetExtrema(void);
        void Setdxy(void);
        void PreprocessIntegral(void);
        void CreateIntegralCoefficients(void);
        double CreateIntegralVector(void);
        void NormalizeIntegralVector(const double &Total);
        
        double InterpolateValue(const Point &Test) const;
        
        Point ConvertIndextoWorld(const int ii) const;
    };
    // Partition Class--------------------------------------------------------------------------------------------------
    // ------------------------------------------------------------------------------------------------------------------
    class Partition
    {
    public:

        Partition(int NRegions = 0, Density Prior = Density(), std::vector<double> desired_area = {}, Parameters Alg_Params = Parameters());


        void SetPartitionVariables(int NRegions = 0, Density Prior = Density(), std::vector<double> desired_area = {});
        void InitializePartition(std::vector<Point> Centers = {},std::vector<double> Weights = {});
        std::vector<Poly> GetCovering(void){return Covering;}
        std::vector<Point> GetCenters(void){return Centers;}
        std::vector<double> GetWeights(void){return Weights;}
        void CalculatePartition(bool WriteToFile, std::string filename_partition = "", std::string filename_centers = "" );
    private:
        std::vector<Point> Centers; 
        std::vector<Poly> Covering; 
        std::vector<double> Weights; 

        const Parameters Alg_Params;
        std::vector<double> desired_area;
        Density Prior;
        int NRegions;


        void CheckParams(void);
        bool CreateDefaultCenters(const Poly Region, const double multiplier);
        void CreateDefaultCenters(const Poly Region, const double initial_multiplier, const int max_steps);
        bool CreatePowerDiagram(void);
        void CleanCovering(const double tolerance, const long int &mult);
        void CreateDelaunayGraph(DelaunayGraph &Delaunay) const;
        double GradientStepCenter(const std::vector<double> &volumes);
        void GradientStepCenter(const double &temp_step, const std::vector<double> &volumes);
        void GradientStepWeights(const std::vector<double> &volumes, const DelaunayGraph &SharedEdges);
        double CalculateError(const std::vector<double> &volumes);
        std::vector<double> CalculateVolumes(void);
    };
    
}//the AreaCon namespace
#endif /* defined(__AreaCon__areacon__) */