micelle_utility.py

"""
A utility file for all things related to Micelles.

@author Peter Zaffetti
@date 2018
"""
from enum import Enum
from logger import get_logger
from primatives.point import euclidean_distance_between_np_array_points
from quick_hull import qhull3d
# from scipy.spatial import ConvexHull
from scipy.cluster.vq import kmeans2
from primatives.triangle import generate_triangles_from_quick_hull_output, cull_pt_vs_bounding_box_plus_tolerance, point_triangle_distance, Triangle
from primatives.point import convert_to_raw_points, compute_point_set_difference, compute_point_set_intersection


import numpy as np
from micelle_plot import *
# import matplotlib.pyplot as plt
# # need the below line for the 3d plotting.
# from mpl_toolkits.mplot3d import axes3d, Axes3D


class Axis(Enum):
    X_AXIS = 0
    Y_AXIS = 1
    Z_AXIS = 2


# TODO: need to update these values, they are randomly picked at the moment.
class MicelleType(Enum):
    def __str__(self):
        return str(self.value)

    UNKNOWN = 6
    CYLINDER = 7
    SPHERE = 8
    ELLIPSOID = 9
    WORM = 10


def identify_micelle(micelle, boundary_points, floating_point_tolerance=0.001, hausdorff_limit=1.00000, plot_output=False):
    """
    Identifies the Micelle shape from the given micelle and boundary points.

    :param micelle: the micelle object that was created from the points in the micelle file provided.
    :param boundary_points: the set of boundary points that was selected. This is a list of Point objects.
    :param floating_point_tolerance: the allowable amount of error when doing floating point calculations. Consider this an epsilon.
    :param hausdorff_limit: the allowable distance between two points before they are considered too far apart.
    :param plot_output: whether or not to plot the resulting micelle.
    :return: the identifier of the micelle. Corresponds to the MicelleType class.
    """
    if micelle is None:
        raise RuntimeError("The Micelle object should not be none.")

    if boundary_points is None:
        raise RuntimeError("The boundary points of the Micelle must be supplied.")

    logger = get_logger(__name__)

    # get the Micelle points
    np_micelle_points = micelle.get_points_as_np_arrays(unique_points=False)
    raw_micelle_points = micelle.get_raw_points()

    raw_boundary_points = convert_to_raw_points(boundary_points)

    micelle_convex_hull_vert_list, micelle_convex_hull_triangles_indices = qhull3d(vertices=raw_micelle_points)
    micelle_convex_hull_triangles =\
        generate_triangles_from_quick_hull_output(point_list=micelle_convex_hull_vert_list,
                                                  triangle_indices=micelle_convex_hull_triangles_indices)

    boundary_convex_hull_vert_list, boundary_convex_hull_triangles_indices = qhull3d(vertices=raw_boundary_points)
    boundary_convex_hull_triangles =\
        generate_triangles_from_quick_hull_output(point_list=boundary_convex_hull_vert_list,
                                                  triangle_indices=boundary_convex_hull_triangles_indices)
    micelle_id = MicelleType.UNKNOWN

    non_convex_points, convex_points = is_worm_micelle(boundary_points=boundary_points,
                                        micelle_convex_hull_vert_list=micelle_convex_hull_vert_list,
                                        micelle_convex_hull_triangles=micelle_convex_hull_triangles, hausdorff_limit=hausdorff_limit,
                                        floating_point_tolerance=floating_point_tolerance)
    if len(non_convex_points) > 0:
        micelle_id = MicelleType.WORM
    else:
        # generate the centroid (center point) of the set.
        # kmeans_result = kmeans(obs=points, k_or_guess=1)
        kmeans2_result = kmeans2(data=np_micelle_points, k=1)
        cluster_centroid = kmeans2_result[0][0]
        logger.info("Cluster centroid: " + str(cluster_centroid))

        # compute the diameter of the convex hull
        diameter = compute_diameter_from_convex_hull(micelle_convex_hull_vert_list)
        logger.info("Diameter of convex hull: " + str(diameter))

        if is_spherical_convexity(convex_hull_vert_list=micelle_convex_hull_vert_list, centroid=cluster_centroid,
                                  allowable_error=floating_point_tolerance):
            micelle_id = MicelleType.SPHERE
        elif is_cylindrical_micelle(convex_hull_vert_list=micelle_convex_hull_vert_list, centroid=cluster_centroid,
                                    allowable_error=floating_point_tolerance):
            micelle_id = MicelleType.CYLINDER
        elif is_ellipsoidal_micelle(convex_hull_vert_list=micelle_convex_hull_vert_list, axis_of_alignment=Axis.Y_AXIS,
                                    convex_hull_triangles_indices=micelle_convex_hull_triangles_indices,
                                    allowable_error=floating_point_tolerance):
            micelle_id = MicelleType.ELLIPSOID

    logger.info("The Micelle shape was identified as: {}".format(micelle_id.name))

    dat = {'micelle_points' : np_micelle_points,
            'boundary_points' : convex_points, # boundary_points
            'convex_hull_points' : micelle_convex_hull_vert_list,
            'convex_hull_triangles_indices' : micelle_convex_hull_triangles_indices,
            'non_convex_points' : non_convex_points}

    if plot_output:
        # show_micelle_plot(**dat)
        prog = InteractiveMicellePlot(**dat)
        prog.cmdloop()
    # TODO: fix this default value
    return micelle_id


def is_worm_micelle(hausdorff_limit, boundary_points, micelle_convex_hull_vert_list,
                    micelle_convex_hull_triangles, floating_point_tolerance):
    """
    Determines whether the Micelle is a worm by looking at the Hausdorff distance. If the minimum distance between any
    of the points not in the convex hull measured to the convex hull is greater than the hausdorff_limit then it is a
    worm.

    :param hausdorff_limit: the maximum allowable distance before the shape is identified as a worm. (Name might need changing)
    :param boundary_points: the points that make up the boundary points of the micelle
    :param micelle_convex_hull_vert_list: the points that make up the convex hull of the micelle.
    :param micelle_convex_hull_triangles: indices into the micelle points list that each make up a triangle in the micelle convex hull.
    :param floating_point_tolerance: the amount of tolerance allowable when doing a floating point comparison before the comparison is considered false.
    :return: True if it is a worm, False otherwise.
    """
    if boundary_points is None:
        raise RuntimeError("The Micelle boundary points can't be None when identifying if the Micelle is a worm.")

    if micelle_convex_hull_vert_list is None:
        raise RuntimeError("The convex hull object can't be None when identifying if the Micelle is a worm.")

    if micelle_convex_hull_triangles is None:
        # raise RuntimeError("The convex hull triangles can't be None when identifying if the Micelle is a worm.")
        return [None]  # TODO: PDZ- just return a list of one element for now since it will count as a worm. If there is a micelle of 4 points then the identification, presently, marks it as a worm.

    points_common_to_boundary_and_micelle_convex_hull = compute_point_set_intersection(boundary_points, micelle_convex_hull_vert_list)
    points_in_boundary_but_not_micelle_convex_hull = compute_point_set_difference(boundary_points, points_common_to_boundary_and_micelle_convex_hull)

    # check to see if there are no points in the boundary that are not shared with the convex hull. If no points exist, then our assumption that the boundary will contain an interior point of the micelle failed.
    if len(points_in_boundary_but_not_micelle_convex_hull) == 0:
        raise RuntimeError("Boundary data fails key assumption of containing some interior point of the micelle's convex hull, CH(S), for non-convex sets. This is a special case and this algorithm might not be able to reliably decide convexity.")


    points_that_break_convexity = list()
    points_that_dont_break_convexity = list()
    # check to see if any of the boundary points that are not on the micelle's convex hull are too far from the micelle's convex hull.
    for point in points_in_boundary_but_not_micelle_convex_hull:
        count_far = 0  # counts the number of triangles which are farther than the hausdorff limit
        for triangle in micelle_convex_hull_triangles:
            # check to see if the point is outside of the bounding box of the triangle, if it is increase the count by one.
            if cull_pt_vs_bounding_box_plus_tolerance(boundary_point=point, triangle=triangle, comparison_value=hausdorff_limit, tolerance=floating_point_tolerance):
                count_far += 1
            # otherwise, if the bounding box is huge, check to see if the distance to the triangle is closer than the limit.
            else:
                distance, closest_point_on_triangle = point_triangle_distance(triangle=triangle.as_numpy_array(), point=point)
                if distance > hausdorff_limit:
                    count_far += 1

        # check to see if the point is too far from all the triangles, if so then it is a worm.
        if count_far == len(micelle_convex_hull_triangles):
            points_that_break_convexity.append(point)
        else:
            points_that_dont_break_convexity.append(point)

    return points_that_break_convexity, points_that_dont_break_convexity


def is_spherical_convexity(convex_hull_vert_list=None, centroid=None, allowable_error=None):
    """
    Checks to see if the Micelle is a sphere.

    :param convex_hull_vert_list: the list of points (x,y,z) from the original Micelle point set that make up the
    convex hull.
    :param centroid: the center point (x,y,z) that makes up the
    :param allowable_error: the amount difference that is allowed when checking all the distances to the centroid
    against each other, past which the distances are not equal (epsilon, amount of error, etc.)
    :return: True if the Micelle is a sphere, False otherwise
    """
    if convex_hull_vert_list is None:
        raise RuntimeError("The convex hull vertices list cannot be None when identifying if the Micelle is a sphere.")
    elif centroid is None:
        raise RuntimeError("The centroid cannot be none when identifying if the Micelle is a sphere.")
    elif allowable_error is None:
        raise RuntimeError("The allowable error must be specified when identifying if the Micelle is sphere.")

    logger = get_logger(__name__)

    # find the distance between all points on the convex hull and the centroid.
    vertex_distance_list_for_sphere = list()
    for point in convex_hull_vert_list:
        vertex_distance_list_for_sphere.append(euclidean_distance_between_np_array_points(point, centroid))
        logger.info("Distance is: " + str(euclidean_distance_between_np_array_points(point, centroid)))

    for distance in vertex_distance_list_for_sphere:
        for comparing_distance in vertex_distance_list_for_sphere:
            # check to see if there is too great a difference between two distances to the centroid.
            # If so, it isn't a sphere.
            if abs(distance - comparing_distance) > allowable_error:
                return False

    return True


def is_cylindrical_micelle(convex_hull_vert_list=None, centroid=None, allowable_error=None):
    """
    Checks to see if the Micelle is cylindrical.

    Checks to see if there exists an axis (x,y, or z) which is equidistant from every point on the convex hull. If such
    an axis exists, then the Micelle is a cylinder.

    :param convex_hull_vert_list: the points in the convex hull.
    :param centroid: the centroid (center of the cluster) of the set of points.
    :param allowable_error: the amount of acceptable error in estimation of equivalent distances.
    :return: True if the Micelle is a cylinder, false otherwise.
    """
    logger = get_logger(__name__)

    if convex_hull_vert_list is None:
        raise RuntimeError("The convex hull vert list cannot be None when identifying if the Micelle is a cylinder.")
    if allowable_error is None:
        raise RuntimeError("The allowable error cannot be None when identifying if the Micelle is a cylinder.")

    x_axis = 0
    y_axis = 1
    z_axis = 2

    x_bounds = (min(x_point[x_axis] for x_point in convex_hull_vert_list), max(x_point[x_axis]
                                                                               for x_point in convex_hull_vert_list))
    y_bounds = (min(y_point[y_axis] for y_point in convex_hull_vert_list), max(y_point[y_axis]
                                                                               for y_point in convex_hull_vert_list))
    z_bounds = (min(z_point[z_axis] for z_point in convex_hull_vert_list), max(z_point[z_axis]
                                                                               for z_point in convex_hull_vert_list))

    num_axis_divisions = 10
    min_index = 0
    max_index = 1
    axis_info_index = 1

    # iterate through each axis
    for axis in [(x_axis, x_bounds), (y_axis, y_bounds), (z_axis, z_bounds)]:
        min_distance_to_axis = list()
        axis_info = axis[axis_info_index]
        axis_min = axis_info[min_index]
        axis_max = axis_info[max_index]
        axis_point_step_size = (axis_max - axis_min) / num_axis_divisions

        # break the axis down into 10 equally distant points
        for division_count in range(0, num_axis_divisions):
            axis_increment = axis_min + (division_count * axis_point_step_size)
            axis_point = [0, 0, 0]
            axis_point[axis[0]] = axis_increment
            # compare each axis point to each convex hull point and take the minimum distance.
            min_distance_to_axis.append(min(euclidean_distance_between_np_array_points(axis_point, point)
                                            for point in convex_hull_vert_list))

        distances_are_equal = True
        compare_dist = min_distance_to_axis[0]
        # check all distances to see if they are even (within tolerance).
        for distance in min_distance_to_axis:
            if abs(distance - compare_dist) > allowable_error:
                distances_are_equal = False
                break

        if distances_are_equal:
            height_axis = "X axis"
            if axis == y_axis:
                height_axis = "Y axis"
            elif axis == z_axis:
                height_axis = "Z axis"

            logger.info("The axis that all convex hull points are equidistant to is: %s" % height_axis)
            return True

    return False


def is_ellipsoidal_micelle(convex_hull_vert_list=None,  convex_hull_triangles=None, convex_hull_triangles_indices=None,
                           axis_of_alignment=Axis.Y_AXIS, allowable_error=.0001):
    # TODO: remove this return. This is only to avoid a current bug with the ellipsoid code.
    return True
    X_COORD = 0
    Y_COORD = 1
    Z_COORD = 2
    W_COORD = 3

    if convex_hull_vert_list is None:
        raise RuntimeError("The convex hull list cannot be None when checking if the Micelle is an ellipsoid.")
    if convex_hull_triangles_indices is None:
        raise RuntimeError("The convex hull triangle indices list cannot be None when checking if the Micelle is an ellipsoid.")

    num_convex_hull_points = len(convex_hull_vert_list)

    max_point_1, max_point_2 = None, None
    major_axis_len = 0
    num_equivalent_max_distance = 0  # counts the number of times the max distance is encountered.

    # compute the max distance over all pairs of boundary points
    for point_index in range(0, num_convex_hull_points):

        # check against every other point but don't check against the point itself (so add 1 to avoid the check)
        comparison_point_index = point_index + 1
        while comparison_point_index < num_convex_hull_points:
            point = convex_hull_vert_list[point_index]
            comparison_point = convex_hull_vert_list[comparison_point_index]
            distance = euclidean_distance_between_np_array_points(point, comparison_point)
            print("Distance is: {}".format(distance))

            # if we reencounter the same max distance (above the error amount) for a different point, increment the count
            if 0 <= (distance - major_axis_len) <= allowable_error:
                num_equivalent_max_distance += 1

            # if the distance is greater than the current max, set it to the max and save the point.
            elif distance > major_axis_len:
                num_equivalent_max_distance = 1
                major_axis_len = distance
                max_point_1, max_point_2 = point, comparison_point

            comparison_point_index += 1

    if major_axis_len < 1:
        # PDZ- This comes from the core-PDZ-6a.txt file where we need to make sure that the max distance is greater than 1
        raise RuntimeError("The max distance is too small.")

    transform_matrix = align_line_segment_midpoint_to_axis(point_1=max_point_1, point_2=max_point_2,
                                                           axis_to_align_to=axis_of_alignment)

    # translate and rotate all the points of the convex hull. It is important to keep the order the same so that the triangles can be regenerated properly.
    transformed_convex_hull = list()
    for point in convex_hull_vert_list:
        homogeneous_point = list(point)

        # homogenize the point with a 4th coordinate if it only has 3
        if len(homogeneous_point) == 3:
            homogeneous_point.append(1)
        transformed_convex_hull.append(np.dot(homogeneous_point, transform_matrix))

    transformed_convex_hull_triangles = generate_triangles_from_quick_hull_output(transformed_convex_hull,
                                                                                  convex_hull_triangles_indices)

    # e comes from the core-PDZ-6a.txt file
    e = .25 * major_axis_len
    positive_e_point = [0, e, 0]
    negative_e_point = [0, -e, 0]
    t = 0.00001

    triangle_set_for_positive_e = list()
    triangle_set_for_negative_e = list()

    # collect two sets of triangles: those close to e and those close to -e
    for triangle in transformed_convex_hull_triangles:
        point_1, point_2, point_3 = triangle.as_numpy_array()
        max_y = max(point_1[Y_COORD], point_2[Y_COORD], point_3[Y_COORD])
        min_y = min(point_1[Y_COORD], point_2[Y_COORD], point_3[Y_COORD])

        # check to see if the triangle is locally close for e
        if (max_y + t) < e or e < (min_y + t):
            pass
        else:
            triangle_set_for_positive_e.append(triangle)

        # check to see if the triangle is locally close for -e
        if (max_y + t) < -e or -e < (min_y + t):
            pass
        else:
            triangle_set_for_negative_e.append(triangle)

    # convert the triangles back to non-homogeneous points to use the pt triangle distance algo
    non_homogeneous_pos_e_triangles = list()
    non_homogeneous_neg_e_triangles = list()
    for triangle in triangle_set_for_positive_e:
        point_1, point_2, point_3 = triangle[0], triangle[1], triangle[2]
        non_homogeneous_pos_e_triangles.append(Triangle(point_1[:3], point_2[:3], point_3[:3]))

    for triangle in triangle_set_for_negative_e:
        point_1, point_2, point_3 = triangle[0], triangle[1], triangle[2]
        non_homogeneous_neg_e_triangles.append(Triangle(point_1[:3], point_2[:3], point_3[:3]))

    pos_e_dist_list = list()
    neg_e_dist_list = list()

    # determine the max distance from the triangle list to positive e
    for triangle in non_homogeneous_pos_e_triangles:
        distance, point = point_triangle_distance(triangle.as_numpy_array(), positive_e_point)
        pos_e_dist_list.append(distance)

    for triangle in non_homogeneous_neg_e_triangles:
        distance, point = point_triangle_distance(triangle.as_numpy_array(), negative_e_point)
        neg_e_dist_list.append(distance)

    max_pos_dist = max(pos_e_dist_list)
    max_neg_dist = max(neg_e_dist_list)

    if abs(max_pos_dist - max_neg_dist) <= allowable_error and max_pos_dist < (major_axis_len * .5):
        return True

    return False


def align_line_segment_midpoint_to_axis(point_1, point_2, axis_to_align_to):
    """
    Generate the matrix which translates and rotates the given line segment (formed between point_1 and point_2) to
    align with the given axis. The midpoint of the line segment will be calculated and aligned to the origin.

    :param point_1: One point that makes up a line segment.
    :param point_2: Another point that makes up a line segment.
    :param axis_to_align_to: the axis to align the line segment to.
    :return: the transformation matrix to align the line segment to the axis.
    """
    X_COORD = 0
    Y_COORD = 1
    Z_COORD = 2
    W_COORD = 3

    # error out if we get 2D points
    if len(point_1) < 3 or len(point_2) < 3:
        raise RuntimeError("Alignment of a line segment to an axis can only be done in 3D at present.")

    # convert the points to lists so we can easily append later if needed.
    point_1 = list(point_1)
    point_2 = list(point_2)

    # create a 3D homogeneous vector of the midpoint of the line segment
    midpoint = np.array([[(point_1[0] + point_2[0]) / 2],
                         [(point_1[1] + point_2[1]) / 2],
                         [(point_1[2] + point_2[2]) / 2],
                         [1]
                         ])

    # create the 3D homogeneous translation matrix to shift the midpoint to the origin
    origin_translation_matrix = np.array([[1, 0, 0, -midpoint[0, 0]],
                                          [0, 1, 0, -midpoint[1, 0]],
                                          [0, 0, 1, -midpoint[2, 0]],
                                          [0, 0, 0, 1]])

    # make point 1 and point 2 homogeneous if they are not. This should probably be a stronger check rather than the length.
    if len(point_1) == 3: # 3D points are only homogeneous if they have 4 coordinates
        point_1.append(1)
    if len(point_2) == 3:
        point_2.append(1)

    point_1_translated = np.dot(point_1, origin_translation_matrix)
    point_2_translated = np.dot(point_2, origin_translation_matrix)

    if axis_to_align_to == Axis.X_AXIS:
        raise NotImplemented("Not yet implemented for the X-Axis")
    elif axis_to_align_to == Axis.Y_AXIS:
        z_rotation_matrix = None
        x_rotation_matrix = None

        # rotate around the Z axis. Begin by zeroing out the z coordinate of the translated max Y point. This is done by only taking the 2 other coordinate points
        selected_point = point_1_translated if point_1_translated[Y_COORD] > point_2_translated[Y_COORD] else point_2_translated
        zeroed_z_point = [selected_point[X_COORD], selected_point[Y_COORD]]

        # create a right triangle with the y axis by projecting the zeroed Z- point to the Y-axis (by setting X = 0)
        y_axis_point = [0, zeroed_z_point[1]]

        # check to see if the line is aligned with the z axis, if so we need a 90 degree (pi/2 radian) rotation
        if y_axis_point[1] == 0:
            z_axis_angle_of_rotation = np.pi / 2.0
        else:
            # now that we have a triangle on the X-Y plane, get the angle between the line segment and the y-axis
            # the arctan of the opposite (x value) over the adjacent (y value) will give the angle
            z_axis_angle_of_rotation = np.arctan(zeroed_z_point[0] / y_axis_point[1])

        z_rotation_matrix = generate_z_axis_rotation_matrix(z_axis_angle_of_rotation)

        # rotate around the X axis.
        zeroed_x_point = [selected_point[Y_COORD], selected_point[Z_COORD]]

        # create a right triangle with the y axis by projecting the zeroed X point to the Y-Axis by setting Z = 0
        y_axis_point = [0, zeroed_x_point[0]]

        # check to see if the line is aligned with the x axis, if so we need a 90 degree (pi/2 radian) rotation
        if y_axis_point[1] == 0:
            x_axis_angle_of_rotation = np.pi / 2.0
        else:
            # now that we have a triangle on the Z-Y plane, get the angle between the line segment and the y-axis
            # the arctan of the oppposite (Z value) over the adjacent (y value) will give the angle
            x_axis_angle_of_rotation = np.arctan(zeroed_x_point[1] / y_axis_point[1])

        x_rotation_matrix = generate_x_axis_rotation_matrix(x_axis_angle_of_rotation)

        # combine the translation and rotation operations matrices to reduce the computational load per point operation.
        result_matrix = np.dot(origin_translation_matrix, x_rotation_matrix)
        result_matrix = np.dot(z_rotation_matrix, result_matrix)
        return result_matrix
    else:
        raise NotImplemented("Not yet implemented for the Z-Axis")

    return None


def generate_x_axis_rotation_matrix(angle_of_rotation):
    """
    Creates a matrix that can be used to rotate a 3D point a certain angle around the x axis.

    :param angle_of_rotation: the angle of rotation around the point.
    :return: a matrix which can be used to rotate a 3D point.
    """
    return np.array([[1, 0, 0, 0],
                     [0, np.cos(angle_of_rotation), -np.sin(angle_of_rotation), 0],
                     [0, np.sin(angle_of_rotation), np.cos(angle_of_rotation), 0],
                     [0, 0, 0, 1]
                     ])


def generate_y_axis_rotation_matrix(angle_of_rotation):
    """
    Creates a matrix that can be used to rotate a 3D point a certain angle around the y axis.

    :param angle_of_rotation: the angle of rotation around the point.
    :return: a matrix which can be used to rotate a 3D point.
    """
    return np.array([[np.cos(angle_of_rotation), 0, np.sin(angle_of_rotation), 0],
                     [0, 1, 0, 0],
                     [-np.sin(angle_of_rotation), 0, np.cos(angle_of_rotation), 0],
                     [0, 0, 0, 1]
                     ])


def generate_z_axis_rotation_matrix(angle_of_rotation):
    """
    Creates a matrix that can be used to rotate a 3D point a certain angle around the z axis.

    :param angle_of_rotation: the angle of rotation around the point.
    :return: a matrix which can be used to rotate a 3D point.
    """
    return np.array([[np.cos(angle_of_rotation), -np.sin(angle_of_rotation), 0, 0],
                     [np.sin(angle_of_rotation), np.cos(angle_of_rotation), 0, 0],
                     [0, 0, 1, 0],
                     [0, 0, 0, 1]
                     ])


def show_micelle_plot(micelle_points=None, boundary_points=None, convex_hull_points=None, convex_hull_triangles_indices=None, non_convex_points=list()):
    """
    Draws the Micelle and its convex hull.

    :param micelle_points: the set of points that make up the Micelle.
    :param boundary_points: the points that make up the boundary.
    :param convex_hull_points: the convex hull points to be drawn.
    :param convex_hull_triangles_indices: the triangles that make up the convex hull. This is a list of indices into the convex_hull_points
    :param non_convex_points: the non-convex points, if any, that identify the worm.
    """
    if micelle_points is None:
        raise RuntimeError("The Micelle points cannot be None if a plot is desired.")

    undecimated_fig = plt.figure()
    undecimated_axis = undecimated_fig.add_subplot(111, projection="3d")
    for point in micelle_points:
        undecimated_axis.scatter(point[0], point[1], point[2], c="#000000")

    undecimated_fig.savefig("undecimated.png")
    undecimated_fig.show()

    fig = plt.figure()
    ax = fig.add_subplot(111, projection="3d")
    for point in boundary_points:
        ax.scatter(point[0], point[1], point[2], c="b")

    if convex_hull_points:
        for point in convex_hull_points:
            ax.scatter(point[0], point[1], point[2], c="g")

    if convex_hull_triangles_indices:
        # you need the vertices to be able to show the triangles.
        if not convex_hull_points:
            pass
        # for each triangle, get the three indexes which make up the triangle
        for (first_index, second_index, third_index) in convex_hull_triangles_indices:
            # convert each index into a hull point so the lines of the triangle can be connected.
            (first_pt_x, first_pt_y, first_pt_z) = convex_hull_points[first_index]
            (second_pt_x, second_pt_y, second_pt_z) = convex_hull_points[second_index]
            (third_pt_x, third_pt_y, third_pt_z) = convex_hull_points[third_index]

            # connect the first to the second point to form the first line.
            ax.plot([first_pt_x, second_pt_x],
                    [first_pt_y, second_pt_y],
                    [first_pt_z, second_pt_z], c="r")

            # connect the second point to the first point to form the second line.
            ax.plot([second_pt_x, third_pt_x],
                    [second_pt_y, third_pt_y],
                    [second_pt_z, third_pt_z], c="r")

            # connect the third point back to the first point to finish and close up the triangle.
            ax.plot([third_pt_x, first_pt_x],
                    [third_pt_y, first_pt_y],
                    [third_pt_z, first_pt_z], c="r")

    for point in non_convex_points:
        ax.scatter(point[0], point[1], point[2], c="#FF00FF")
        # value = str(point[0]) + ", " + str(point[1]) + ", " + str(point[2])
        # ax.text(point[0], point[1], point[2], s=value)
    plt.savefig("decimated.png")
    # plt.show()


def compute_diameter_from_convex_hull(convex_hull_vert_list=None):
    """
    Computes the diameter of a convex hull.

    :param convex_hull_vert_list: the set of vertices that make up a convex hull
    :return: the diameter across the convex hull.
    """
    if convex_hull_vert_list is None:
        raise RuntimeError("The convex hull should not be none when calculating the diameter.")

    # logger = get_logger(__name__)

    max_distance = 0
    for point in range(len(convex_hull_vert_list)):
        comparing_point = convex_hull_vert_list[point]

        # Copy the set of points on the convex hull.
        # Remove the points adjacent since a diameter can't be directly next to the selected point
        comparing_set = list(convex_hull_vert_list)
        del comparing_set[point + 1 if point + 1 < len(convex_hull_vert_list) else 0]
        del comparing_set[point - 1 if point - 1 >= 0 else -1]

        for remaining_point in comparing_set:
            max_distance = euclidean_distance_between_np_array_points(comparing_point, remaining_point) \
                if euclidean_distance_between_np_array_points(comparing_point, remaining_point) > max_distance \
                else max_distance

    return max_distance