import numpy as np
from scipy.spatial import ConvexHull


class HistogramBuilder:

    def __init__(self, center, number_of_spheres=None, max_radius=None, number_of_bins=None):

        # Define attributes of spherical histogram
        # These attributes will be passed in by user(?)
        self.center = center # Coordinate at center of histogram

        self.number_of_spheres = number_of_spheres # Number of concentric spheres defining layers of the histograms
        if number_of_spheres is None:
            self.number_of_spheres = 2

        self.max_radius = max_radius # Maximum size of outer sphere
        if max_radius is None:
            self.max_radius = 20.0

        self.number_of_bins = number_of_bins  # Number of bins PER concentric sphere
        if number_of_bins is None:
            self.number_of_bins = 16.0

        error = f"Center of type {type(center)} must be a list or numpy array with dimension 1x3."
        if isinstance(center, list):
            if not len(center) == 3:
                raise TypeError(error)

        elif isinstance(center, np.ndarray):
            if not center.shape == (3, 1) and not center.shape == (3,):
                raise TypeError(error)

        else:
            raise TypeError(error)

        # Define attributes calculated with passed parameters
        # V - E + F = 2, F is facets (num of bins) and 2E (edges) = 3F
        self.number_of_angles = (self.number_of_bins / 2) + 2 # Number of angles computed needed to obtain the number of bins
        if self.number_of_spheres == 1: # Define radius for case where only one sphere is present, otherwise geomspace returns one instead of going up to the max radius
          rad = [self.max_radius]
        else:
          rad = np.geomspace(1, self.max_radius, self.number_of_spheres) # Increment radii logarithmically so bins in sphere close to center are smaller
        self.radii = np.linspace(0, self.max_radius, self.number_of_spheres+1) # Use linear spacing for radii
        
    def get_angles(self):
        # Obtain angles for points distributed evenly across sphere's surface
        # Fibonacci spheres method approximates this by incrementing points by golden ratio
        # Distrubtion is approximately even because the golden ratio is the most irrational number
        # Source: http://arxiv.org/pdf/0912.4540
        idx = np.arange(0, self.number_of_angles, dtype=float)
        golden_ratio = (1 + 5**0.5) / 2 # 
        phi = np.arccos(1 - 2 * idx / self.number_of_angles)
        theta = 2 * np.pi *idx / golden_ratio
        return phi, theta

    def get_points(self, angles):
        # Convert polar coordinates to cartesian
        # Theta azimuthal angle (xy plane)
        # Phi is zenith angle (from z-axis)
        phi, theta = angles
        x, y, z = np.cos(theta) * np.sin(phi), np.sin(theta) * np.sin(phi), np.cos(phi)
        points = np.vstack((x,y,z)).T
        return points

    def get_bins(self, points, lrf=None):
        # Define hexahedral bins based on faces defined by connecting vertices
        hex_bins = []
        poly3d = []
        
        # print(points.shape)
        # Use built in function to construct convex hull(s) using points on sphere's surface
        if lrf is not None:
            points = points @ lrf.T
            # print(points.shape)
        hull = ConvexHull(points)
        
        # print(f'ConvexHull simplices length is {len(hull.simplices)}')
        for k in range(len(self.radii)-1):
            for i in range(len(hull.simplices)):
                # Evaluate set of points on outer sphere
                new_points = points[hull.simplices[i]]
                first_face = new_points*self.radii[k] + self.center
                second_face = new_points*self.radii[k+1] + self.center

                # Get vertices for bins
                plot_points = np.vstack((first_face, second_face))
                hex_bins.extend(plot_points)

                # Get plotting information
                vertices = [[0, 1, 2], [2, 5, 4, 1], [1, 0, 3, 4], [4, 5, 3], [3, 0, 2, 5]]
                poly3d.append([[plot_points[v] for v in vertex] for vertex in vertices])
        return hex_bins, poly3d

    @classmethod
    def get_histogram_builder(cls, center, number_of_spheres=None, max_radius=None, number_of_bins=None):
        histogram_builder = cls(center, number_of_spheres, max_radius, number_of_bins)
        return histogram_builder
