inscribe-hatch.py

#!/usr/bin/env python3
import itertools
import typing
import math

import cairo
import shapely.affinity
import shapely.geometry

DISTANCE_CUTOFF = 2  # max jump length to bridge with a curve
OUTLINE_DISTANCE = .5  # ideally half the line width
INSET_DISTANCE = .9  # distance from outside edge to begin hatching
HATCH_OVERLAP = .3  # amount to overlap hatch lines to account for slop; 0 = no overlap
CURVE_JOIN_DISTANCE = .25

BACK_AND_FORTH = True


class Segment (typing.NamedTuple):
    '''
    '''
    row: int
    x1: float
    y1: float
    x2: float
    y2: float
    
    def distance(self, other):
        ''' Distance from self end to other beginning
        '''
        return math.hypot(self.x2 - other.x1, self.y2 - other.y1)


class Fragment (typing.NamedTuple):
    '''
    '''
    p1: int
    p2: int
    p: int


class MultiPath (typing.NamedTuple):
    '''
    '''
    points: shapely.geometry.MultiPoint
    fragments: typing.List[typing.List[Fragment]]


def get_hatch_segment_xs(xmin, xmax):
    '''
    '''
    if not BACK_AND_FORTH:
        # All parallel lines point in the same direction and ends do not connect
        return itertools.cycle([(xmax, xmin)])

    return itertools.cycle([(xmin, xmax), (xmax, xmin)])


def fill_hatch_segments(polygon):
    '''
    '''
    assert polygon.type == 'Polygon'
    
    xmin, y, xmax, ymax = polygon.bounds
    rows, xs = itertools.count(), get_hatch_segment_xs(xmin, xmax)
    segments = []
    
    for (row, (x1, x2)) in zip(rows, xs):
        line = shapely.geometry.LineString([(x1, y), (x2, y)])
        
        if y > ymax:
            break
        elif not line.intersects(polygon):
            continue
        else:
            y += 1
        
        scan = line.intersection(polygon)
        segments.extend([
            Segment(
                row,
                part.coords[0][0],
                part.coords[0][1],
                part.coords[-1][0],
                part.coords[-1][1],
            )
            for part in getattr(scan, 'geoms', [scan])
            if part.type == 'LineString'
        ])
    
    return segments


def fill_hatch_sequences(polygon):
    '''
    '''
    assert polygon.type == 'Polygon'
    
    distance_cutoff = DISTANCE_CUTOFF if BACK_AND_FORTH else 0
    segments = fill_hatch_segments(polygon)
    sequences = []
    
    while segments:
        segment = segments.pop(0)
        sequence = [(segment.x1, segment.y1), (segment.x2, segment.y2)]
        
        while True:
            others = [
                (i, other)
                for (i, other) in enumerate(segments)
                if abs(other.row - segment.row) == 1
                and segment.distance(other) < distance_cutoff
            ]
        
            if not others:
                sequences.append(sequence)
                break

            others.sort(key=lambda o: segment.distance(o[1]))
            segment = segments.pop(others[0][0])
            
            prior_length = abs(sequence[-1][-2] - sequence[-2][-2])
            follow_length = abs(segment.x2 - segment.x1)
            
            prior_diff = min(CURVE_JOIN_DISTANCE, prior_length)
            follow_diff = min(CURVE_JOIN_DISTANCE, follow_length)

            if segment.row % 2 == 0:
                # pointing right
                prior_diff, follow_diff = -prior_diff, -follow_diff
            
            sequence[-1] = tuple(
                list(sequence[-1][:-2]) + [sequence[-1][-2] - prior_diff, sequence[-1][-1]]
            )
            
            x3 = segment.x1 - follow_diff
            y3 = segment.y1

            x1, y1 = sequence[-1][-2] + 2*prior_diff, sequence[-1][-1]
            x2, y2 = x3 + 2*follow_diff, y3
            
            sequence.extend([
                (x1, y1, x2, y2, x3, y3),
                (segment.x2, segment.y2),
            ])

    return sequences


def inscribed_hatch_multipath(polygon, scale, angle=None):
    '''
    '''
    assert polygon.type == 'Polygon'
    
    # Add a bit of overlap to the hatched lines when the Axidraw is sloppy
    _scale = scale / (1.0 + HATCH_OVERLAP)
    
    # Track the best option of the angle deltas we try
    _count, multipath = math.inf, None
    
    # Try several candidate angle deltas?
    angles = [angle] if angle is not None else (0, 26, 52, 77, 103, 129, 154)

    for _angle in angles:
        transformed_polygon = shapely.affinity.rotate(
            shapely.affinity.scale(polygon, 1/_scale, 1/_scale, origin=(0, 0)),
            _angle,
            origin=(0, 0),
        )
        transformed_multipath = native_hatch_multipath(transformed_polygon)
        untransformed_points = shapely.affinity.scale(
            shapely.affinity.rotate(transformed_multipath.points, -_angle, origin=(0, 0)),
            _scale,
            _scale,
            origin=(0, 0),
        )
        if len(transformed_multipath.fragments) < _count:
            # Choose a new, better multipath candidate
            _count = len(transformed_multipath.fragments)
            multipath = MultiPath(untransformed_points, transformed_multipath.fragments)
        if _count <= 3:
            # Eh, good enough
            break

    return multipath

    
def native_hatch_multipath(polygon):
    '''
    '''
    assert polygon.type == 'Polygon'
    
    points, fragments = [], []
    
    outline = polygon.buffer(-OUTLINE_DISTANCE).boundary
    for line in getattr(outline, 'geoms', [outline]):
        if line.type == 'LineString' and True:
            fragments.append([
                Fragment(None, None, i)
                for i in range(len(points), len(points) + len(line.coords))
            ])
            points.extend(line.coords)
    
    inset_multipolygon = polygon.buffer(-INSET_DISTANCE)

    if inset_multipolygon.type not in ('Polygon', 'MultiPolygon') or inset_multipolygon.area == 0:
        return MultiPath(shapely.geometry.MultiPoint(points), fragments)
    
    for inset_polygon in getattr(inset_multipolygon, 'geoms', [inset_multipolygon]):
        for sequence in fill_hatch_sequences(inset_polygon):
            head, tail = sequence[0], sequence[1:]

            assert len(head) == 2

            fragments.append([Fragment(None, None, len(points))])
            points.append(head[-2:])

            for part in tail:
                assert len(part) in (2, 6)
                if len(part) == 2:
                    fragments[-1].append(Fragment(None, None, len(points)))
                    points.append(part)
                else:
                    fragments[-1].append(Fragment(len(points), len(points)+1, len(points)+2))
                    points.extend([part[:2], part[2:4], part[4:6]])
    
    return MultiPath(shapely.geometry.MultiPoint(points), fragments)


def draw_shapely_line_geom(ctx, geom):
    if getattr(geom, 'geoms', False):
        parts = [part for part in geom.geoms]
    else:
        parts = [geom]

    for part in parts:
        try:
            coords = part.coords
        except NotImplementedError:
            continue
        if not coords:
            continue
        ctx.move_to(*coords[0])
        for coord in coords[1:]:
            ctx.line_to(*coord)


def draw_multipath_strokes(ctx, multipath):
    '''
    '''
    for fragments in multipath.fragments:
        head, tail = fragments[0], fragments[1:]
        p = multipath.points.geoms[head.p]
        ctx.move_to(p.x, p.y)
        for fragment in tail:
            if fragment.p1 is None and fragment.p2 is None:
                p = multipath.points.geoms[fragment.p]
                ctx.line_to(p.x, p.y)
            else:
                p1 = multipath.points.geoms[fragment.p1]
                p2 = multipath.points.geoms[fragment.p2]
                p = multipath.points.geoms[fragment.p]
                ctx.curve_to(p1.x, p1.y, p2.x, p2.y, p.x, p.y)
        ctx.stroke()


def calculate_hatch_angle(shape):
    '''
    '''
    rectangle = shape.minimum_rotated_rectangle.exterior
    (x1, y1), (x2, y2), (x3, y3) = rectangle.coords[:3]
    
    if math.hypot(x3 - x2, y3 - y2) < math.hypot(x2 - x1, y2 - y1):
        return 180 * math.atan2(y2 - y3, x3 - x2) / math.pi
    else:
        return 180 * math.atan2(y1 - y2, x2 - x1) / math.pi


if __name__ == '__main__':
    font_svg = __import__('font-svg')
    color_cube = __import__('color-cube')

    FONT = font_svg.load_font('fonts/EMSReadability.svg')

    with cairo.SVGSurface('out.svg', 8.5*72, 5.5*72) as surface:
        context = cairo.Context(surface)
        
        hotdog1 = shapely.geometry.LineString([(72, 72), (144, 144)]).buffer(36, 5)
        
        draw_multipath_strokes(
            context,
            inscribed_hatch_multipath(hotdog1, 9, None),
        )
        
        hotdog2 = shapely.geometry.LineString([(216, 72), (252, 144)]).buffer(36, 5)
        
        draw_multipath_strokes(
            context,
            inscribed_hatch_multipath(hotdog2, 6, None),
        )
        
        hotdog3 = shapely.geometry.LineString([(288, 72), (360, 108)]).buffer(36, 5)
        
        draw_multipath_strokes(
            context,
            inscribed_hatch_multipath(hotdog3, 3, None),
        )
        
        context.stroke()

out.svg