#!/usr/bin/env python
'''
Copyright (c) 2017 Bruce Hill <bruce@bruce-hill.com>
This work is free. You can redistribute it and/or modify it under the
terms of the Do What The Fuck You Want To Public License, Version 2,
as published by Sam Hocevar. See http://www.wtfpl.net/ for more details.
'''
import math
import collections

Triangle = collections.namedtuple("Triangle", "v1 v2 v3")
Edge = collections.namedtuple("Edge", "e1 e2")

def find_contours(elevation_function, xmin, xmax, ymin, ymax, spacing=1, elevation=0.5):
    '''
    Return a list of contour lines that have a constant elevation near the specified
    value. Lines will be truncated to the sepcified x/y range, and the elevation
    function will be sampled at even intervals determined by the spacing parameter.
    '''
    elevation_data = dict()
    for x in range(xmin, xmax+1, spacing):
        for y in range(ymin, ymax+1, spacing):
            elevation_data[(x,y)] = elevation_function(x,y)

    triangles = []
    for x in range(xmin, xmax, spacing):
        for y in range(ymin, ymax, spacing):
            triangles.append(Triangle((x,y), (x+spacing,y), (x,y+spacing)))
            triangles.append(Triangle((x+spacing,y), (x,y+spacing), (x+spacing,y+spacing)))

    contour_segments = []
    for triangle in triangles:
        below = [v for v in triangle if elevation_data[v] < elevation]
        above = [v for v in triangle if elevation_data[v] >= elevation]
        # All above or all below means no contour line here
        if len(below) == 0 or len(above) == 0:
            continue
        # We have a contour line, let's find it
        minority = above if len(above) < len(below) else below
        majority = above if len(above) > len(below) else below

        contour_points = []
        crossed_edges = (Edge(minority[0], majority[0]), Edge(minority[0], majority[1]))
        for triangle_edge in crossed_edges:
            # how_far is a number between 0 and 1 indicating what percent of the way
            # along the edge (e1,e2) the crossing point is
            how_far = ((elevation - elevation_data[triangle_edge.e2])
                     / (elevation_data[triangle_edge.e1] - elevation_data[triangle_edge.e2]))
            crossing_point = (
                how_far * triangle_edge.e1[0] + (1-how_far) * triangle_edge.e2[0],
                how_far * triangle_edge.e1[1] + (1-how_far) * triangle_edge.e2[1])
            contour_points.append(crossing_point)
        contour_segments.append(Edge(contour_points[0], contour_points[1]))

    unused_segments = set(contour_segments)
    segments_by_point = collections.defaultdict(set)
    for segment in contour_segments:
        segments_by_point[segment.e1].add(segment)
        segments_by_point[segment.e2].add(segment)

    contour_lines = []
    while unused_segments:
        # Start with a random segment
        line = collections.deque(unused_segments.pop())
        while True:
            tail_candidates = segments_by_point[line[-1]].intersection(unused_segments)
            if tail_candidates:
                tail = tail_candidates.pop()
                line.append(tail.e1 if tail.e2 == line[-1] else tail.e2)
                unused_segments.remove(tail)
            head_candidates = segments_by_point[line[0]].intersection(unused_segments)
            if head_candidates:
                head = head_candidates.pop()
                line.appendleft(head.e1 if head.e2 == line[0] else head.e2)
                unused_segments.remove(head)
            if not tail_candidates and not head_candidates:
                # There are no more segments touching this line, so we're done with it.
                contour_lines.append(list(line))
                break

    return contour_lines


if __name__ == "__main__":
    # Run a test case with some arbitrary function
    # that maps (x,y) to an elevation between 0 and 1
    def elevation_function(x,y):
        return 1/(2+math.sin(2*math.sqrt(x*x + y*y))) * (.75+.5*math.sin(x*2))
    print(find_contours(elevation_function, 0,100, 0,100, spacing=1, elevation=0.5))