def graham(meshes, xs, ys, zs, expose_horizon=False):
    """
    convex_graham
    """

    distance_factor = 1.0  # TODO: make this an option (what does it even do)?

    # first lets handle the horizon
    visibilities, weights, horizon = only_horizon(meshes, xs, ys, zs)

    # Order the bodies from front to back.  We can do this through zs (which are
    # the z-coordinate positions of each body in the system.
    # Those indices are in the same order as the meshes['comp_no']
    # (needing to add 1 because comp_no starts at 1 not 0)

    # TODO: I don't think this whole comp_no thing is going to work with nested
    # triples.. may need to rethink the whole meshes.component_by_no thing

    front_to_back_comp_nos = np.argsort(zs)[::-1]+1
    nbodies = len(front_to_back_comp_nos)

    for i_front in range(0, nbodies-1):
        # for a binary, i_front will only be 0
        comp_no_front = front_to_back_comp_nos[i_front]
        comp_front = meshes.component_by_no(comp_no_front)
        mesh_front = meshes[comp_front]
        visibility_front = visibilities[comp_front]

        for i_back in range(i_front+1, nbodies):
            # for a binary, i_back will only be 1
            comp_no_back = front_to_back_comp_nos[i_back]
            comp_back = meshes.component_by_no(comp_no_back)
            mesh_back = meshes[comp_back]
            visibility_back = visibilities[comp_back]

            # If mesh_back is entirely hidden, then we can skip any checks and
            # leave it hidden.  Note that here we need to check visibility instead
            # of mesh['visibility'] or mesh_back['visibility'] since we are not
            # adjusting those values in memory.
            if np.all(visibility_back==0.0):
                continue

            # Determine a scale factor for the triangle
            min_size_back = mesh_back.areas.min()
            distance = distance_factor * 2.0/3**0.25*np.sqrt(min_size_back)

            # Select only those triangles that are not hidden
            tri_back_vis = mesh_back.vertices_per_triangle[visibility_back > 0.0].reshape(-1,9)
            tri_front_vis = mesh_front.vertices_per_triangle[visibility_front > 0.0].reshape(-1,9)

            back = np.vstack([tri_back_vis[:,0:2], tri_back_vis[:,3:5], tri_back_vis[:,6:8]])
            front = np.vstack([tri_front_vis[:,0:2], tri_front_vis[:,3:5], tri_front_vis[:,6:8]])

            # Star in front ---> star in back
            if not front.shape[0]:
                continue

            hull, inside = _graham_scan_inside_hull(front, back)
            hidden = inside.reshape(3,-1).all(axis=0)
            visible = ~(inside.reshape(3,-1).any(axis=0))

            # Triangles that are partially hidden are those that are not
            # completely hidden, but do have at least one vertex hidden
            partial = ~hidden & ~visible

            # These returned visibilities are only from mesh_back_vis
            # So to apply back to our master visibilities parameter, we need
            # to find the correct inds.
            visibility_back[visibility_back > 0.0] = 1.0*visible + 0.5*partial


            ###################################################
            ## TODO: port the following code to be included  ##
            ###################################################


            # It's possible that some triangles that do not have overlapping
            # vertices, still are partially covered (i.e. if the objects fall
            # in between the vertices)
            #visible1 = ~star1.mesh['hidden']
            #backc = star1.mesh['center'][visible1,0:2]
            #if not len(backc):
            #    continue
            #
            #subset_partial = closest_points(hull, backc, distance)
            #
            #if not len(subset_partial):
            #    continue
            #
            #arg = np.where(visible1)[0][subset_partial]
            #star1.mesh['visible'][arg] = False
            #star1.mesh['partial'][arg] = True


    return visibilities, None, None
convex_graham