Some tools to (reasonably efficiently) fill shapes with infill patterns, as is common in 3D printing

point_in_polygon.py

import numpy as np

try:
    # If available, import numba to enable compiling the functions for extra speed
    from numba import njit, boolean

except ImportError:
    print("Numba compilation not available - falling back to pure Python.")

    def njit(func, *args, **kwargs):
        return func

    boolean = bool


@njit
def is_inside_sm(polygon, point):
    '''
    Determines whether the given point is inside or outside the specified polygon
    Unlike the standard ray-casting algorithm, this one works on edges! (with no performance cost)

    Algorithm from Sasa Milenkovic (GPL-v3.0 licensed):
       https://github.com/sasamil/PointInPolygon_Py/blob/master/pointInside.py#L134-L168

    arguments:
      Polygon - searched polygon (closed path - the last point should match the first one)
      Point - an arbitrary point that can be inside or outside the polygon

    return value:
      0 - the point is outside the polygon
      1 - the point is inside the polygon 
      2 - the point is on an edge (boundary)
    '''
    length = len(polygon)-1
    dy2 = point[1] - polygon[0][1]
    intersections = 0
    ii = 0
    jj = 1

    while ii<length:
        dy  = dy2
        dy2 = point[1] - polygon[jj][1]

        # consider only lines which are not completely above/bellow/right from the point
        if dy*dy2 <= 0.0 and (point[0] >= polygon[ii][0] or point[0] >= polygon[jj][0]):

            # non-horizontal line
            if dy<0 or dy2<0:
                F = dy*(polygon[jj][0] - polygon[ii][0])/(dy-dy2) + polygon[ii][0]

                if point[0] > F: # if line is left from the point - the ray moving towards left, will intersect it
                    intersections += 1
                elif point[0] == F: # point on line
                    return 2

            # point on upper peak (dy2=dx2=0) or horizontal line (dy=dy2=0 and dx*dx2<=0)
            elif dy2==0 and (point[0]==polygon[jj][0] or (dy==0 and (point[0]-polygon[ii][0])*(point[0]-polygon[jj][0])<=0)):
                return 2

            # there is another posibility: (dy=0 and dy2>0) or (dy>0 and dy2=0). It is skipped 
            # deliberately to prevent break-points intersections to be counted twice.
 
        ii = jj
        jj += 1

    # an odd number of intersections means the point is inside the polygon
    return intersections & 1


@njit
def points_in_polygon(points, polygon):
    ln = len(points)
    D = np.empty(ln, dtype=boolean)
    for i in range(ln):
        D[i] = is_inside_sm(polygon,points[i])
    return D



if __name__ == '__main__':
    from time import perf_counter
    
    polygon_sides = 500
    checked_points = 50_000
    x = np.linspace(0, 2*np.pi, polygon_sides)
    polygon = np.vstack([np.sin(x)+0.5, np.cos(x)+0.5]).T
    points = np.random.uniform(-1.5, 1.5, size=(checked_points, 2))
    
    start = perf_counter()
    points_in_polygon(points, polygon)
    time_elapsed = perf_counter() - start
    print(f'{polygon_sides = }\n{checked_points = :_}\n{time_elapsed = :.5}s')

shape_fill.py

from itertools import chain
import math

import numpy as np

from point_in_polygon import points_in_polygon


def rectangle_fill(length, width, spacing, max_feature_size):
    '''
    'length' is the lengthwise distance being filled.
    'width' is the transverse distance being filled.
    'spacing' is the desired spacing between rectilinear lines.
    'max_feature_size' is the maximum segment length beyond which segments get sub-sampled.
        If set to None, fill line segments are directly from bottom to top.
    '''
    lengthwise_samples = math.ceil(length / max_feature_size) if max_feature_size is not None else 2
    y = np.linspace(0, length, lengthwise_samples)
    x = np.arange(0, width, spacing)
    X, Y = np.meshgrid(x, y)
    Y[:,1::2] = y[::-1, None]  # flip every second column
    path = np.empty((X.size+1, 2))
    path[:-1, 0] = X.T.flatten()
    path[:-1, 1] = Y.T.flatten()
    path[-1] = [width, path[-2,1]]
    return path

def advanced_fill(): # TODO
    '''
    'angle' is an optional anti-clockwise angle (in degrees) to rotate the fill lines by.
    'transverse_offset' is an optional function offsetting long segment points in the transverse direction
        e.g. to create a sinusoid instead of a straight line
    '''
    ...

def cut(path, area, join_collinear_segments=False, remove_internal_segments=False, copy=False):
    ''' Assumes 'area' is strictly contained within the 'path' region. '''
    # make sure it's a numpy array - copy data if necessary/requested
    path = np.array(path, copy=copy)
    # determine path segments that cross into and out of the area
    internal_points = area.contains(path)
    intersecting_segments = np.diff(np.int8(internal_points))
    shape_entries = np.zeros(len(path), dtype=bool)
    shape_entries[:-1][intersecting_segments == 1] = 1
    shape_exits = np.zeros(len(path), dtype=bool)
    shape_exits[1:][intersecting_segments == -1] = 1
    # push the outer points from the intersecting segments onto the area intercepts
    path[shape_entries] = [
        area.intercept(internal, external)
        for internal, external
        in zip(path[1:][shape_entries[:-1]], path[shape_entries])
    ]
    path[shape_exits] = [
        area.intercept(internal, external)
        for internal, external
        in zip(path[:-1][shape_exits[1:]], path[shape_exits])
    ]
    # cut the path down to only the internal and boundary points
    useful_points = shape_entries | shape_exits
    if not remove_internal_segments:
        useful_points |= internal_points
    cut_path = path[useful_points]

    if join_collinear_segments:
        # collinear segments are those with no bending (0 second derivative)
        collinear_segments = np.empty(len(cut_path), dtype=bool)
        deltas = np.diff(cut_path, axis=0)
        deltas /= np.linalg.norm(deltas, axis=1, keepdims=True)
        collinear_segments[1:-1] = np.all(np.isclose(np.diff(deltas, axis=0), 0), axis=1)
        collinear_segments[[0,-1]] = 0  # make sure to include the first and last points of the path
        cut_path = cut_path[~collinear_segments]

    return cut_path

def advanced_cut(): # TODO
    '''
    pushes cut lines onto geometry
    deals with cutting non-convex shapes
    '''
    ...


class FillableRegion:
    def __init__(self, area, fill_density):
        ''' fill_density is a percentage or a vector field of percentages '''
        self._area = area
        self._fill_density = fill_density

    def next_point(self, current_point, target_direction):
        ...


class Area:
    def __contains__(self, pt):
        return self.contains(np.array([pt]))[0]

    def intercept(self, internal_point, external_point, tolerance=1e-4, binary_preferred=False):
        ''' NOTE: default implementations are numerical, and assume "distance_to_boundary" calculations
            (if available, otherwise containment checks) are reasonably cheap.
            More direct/optimal methods of intercept estimation/calculation should generally be preferred.
        '''
        internal_point = np.asarray(internal_point)
        external_point = np.asarray(external_point)
        if not binary_preferred and hasattr(self, 'distance_from_boundary'):
            return self._intercept_distance(internal_point, external_point, tolerance)
        return self._intercept_binary(internal_point, external_point, tolerance)

    def _intercept_binary(self, internal_point, external_point, tolerance):
        while "meeting tolerance":
            # binary search using containment checks
            candidate = (internal_point + external_point) / 2
            if candidate in self:
                internal_point = candidate
            else:
                external_point = candidate

            if np.linalg.norm(external_point - internal_point) <= tolerance:
                break # within tolerance - stop trying to improve

        return candidate

    def _intercept_distance(self, internal_point, query_point, tolerance):
        direction = internal_point - query_point # want to move towards internal point
        direction /= np.linalg.norm(direction) # normalise to unit vector
        while abs(distance := self.distance_from_boundary(query_point)) > tolerance:
            # move towards the internal point by the current distance to the boundary
            query_point += distance * direction
        return query_point

    def generate_boundary_points(self, n_min=3, max_dist=None, max_error=1e-4, closed=False):
        ''' max error requires knowledge of absolute distance from boundary '''
        ...

    def __invert__(self):
        return Exclusion(self)

    def __or__(self, area):
        return Union(self, area)

    def __xor__(self, area):
        return ExclusiveUnion(self, area)

    def __and__(self, area):
        return Intersection(self, area)


class AreaGroup(Area):
    def __init__(self, *areas):
        assert len(areas) > 1, "Groups must contain multiple areas."
        self._areas = areas

    def intercept(self, internal_point, external_point, *args, **kwargs):
        points = [
            area.intercept(internal_point, external_point, *args, **kwargs)
            for area in self._areas
            if internal_point in area # only check relevant sub-regions
        ]
        if len(points) > 1:
            # return the furthest out point
            # TODO: check reasoning / behaviour with nested regions and venn diagram intersection
            return min((point for point in points), key=lambda pt: np.linalg.norm(external_point-pt))
        return points[0]


class Union(AreaGroup):
    ''' In any sub-region(s). '''
    def contains(self, points):
        in_any = self._areas[0].contains(points)

        for area in self._areas[1:]:
            outside = ~in_any
            # only check the points that are not yet confirmed to be in at least one sub-region
            in_any[outside] |= area.contains(points[outside])

        return in_any

    def generate_boundary_points(self, *args, **kwargs):
        boundary_points = []
        for index, area in enumerate(self._areas):
            local_points = area.generate_boundary_points(*args, **kwargs)
            other_areas = Union(self._areas[:index] + self._areas[index+1:]) if len(self._areas) > 2 else self._areas[~index]
            boundary_points.append(local_points[~other_areas.contains(local_points)])
            
        return np.vstack(boundary_points)


class Intersection(AreaGroup):
    ''' In all sub-regions. '''
    def contains(self, points):
        in_all = self._areas[0].contains(points)

        for area in self._areas[1:]:
            # only check the points that have been in all sub-regions so far
            in_all[in_all] &= area.contains(points[in_all])

        return in_all


class ExclusiveUnion(AreaGroup):
    ''' In only one sub-region. '''
    def contains(self, points):
        in_one = self._areas[0].contains(points)

        for area in self._areas[1:]:
            # only check the points that are only known to be within one sub-region
            in_one[in_one] ^= area.contains(points[in_one])

        return in_one


class Exclusion(Area):
    def __init__(self, area):
        self._area = area

    def contains(self, points):
        return ~self._area.contains(points)

    def intercept(self, internal_point, external_point, *args, **kwargs):
        return self._area.intercept(external_point, internal_point, *args, **kwargs)

    def generate_boundary_points(self, *args, **kwargs):
        # reversed, in case that's relevant later
        return self._area.generate_boundary_points(*args, **kwargs)[::-1]


class AxisAlignedRectangle(Area):
    def __init__(self, xmin, xmax, ymin, ymax):
        self.xmin = xmin
        self.xmax = xmax
        self.ymin = ymin
        self.ymax = ymax

    def contains(self, pts):
        pts = np.asarray(pts)
        x = pts[:,0]
        y = pts[:,1]
        return (self.xmin <= x) & (x <= self.xmax) & (self.ymin <= y) & (y <= self.ymax)

    def distance_to_boundary(self, point):
        nearest_boundary_point = np.max([
            np.min([point, [self.xmax, self.ymax]], axis=0), 
            [self.xmin, self.ymin]
        ], axis=0)
        return np.linalg.norm(point - nearest_boundary_point)

    def generate_boundary_points(self, n_min=4, max_sep=None, max_error=None, closed=False):
        xmin, xmax, ymin, ymax = self.xmin, self.xmax, self.ymin, self.ymax
        if max_sep is None: 
            if n_min <= 4:
                boundary_points = [[xmin, ymin], [xmax, ymin], [xmax, ymax], [xmin, ymax]]
                if closed:
                    boundary_points.append(boundary_points[0])
                return boundary_points
            max_sep = np.inf

        # calculate point numbers to meet maximum separation requirement
        xn = math.ceil((xmax - xmin) / max_sep) + 2
        yn = math.ceil((ymax - ymin) / max_sep) + 2

        # add extra points if required by n_min
        n_ratio = 2*(xn+yn-2) / n_min
        if n_ratio < 1:
            xn = math.ceil(xn / n_ratio)
            yn = math.ceil(yn / n_ratio)

        x_pts = np.linspace(xmin, xmax, xn)
        y_pts = np.linspace(ymin, ymax, yn)

        boundary_points = [
            *([x, ymin] for x in x_pts[:-1]),
            *([xmax, y] for y in y_pts[:-1]),
            *([x, ymax] for x in x_pts[:0:-1]),
            *([xmin, y] for y in y_pts[:0:-1])
        ]
        if closed:
            boundary_points.append(boundary_points[0])
        return np.array(boundary_points)


class Circle(Area):
    def __init__(self, center, radius):
        self.center = np.asarray(center)
        self.radius = radius

    def contains(self, pts):
        pts = np.asarray(pts)
        return np.sum((pts - self.center)**2, axis=1) <= self.radius**2

    def distance_to_boundary(self, point):
        return np.linalg.norm(point - self.center) - self.radius

    def generate_boundary_points(self, n_min=10, max_sep=None, max_error=None, closed=False, optimise_perimeter=False):
        max_sep = np.inf if max_sep is None else max_sep
        max_sep_n = math.ceil(2 * math.pi * self.radius / max_sep)
        if max_error is None:
            max_error_n = 0
        else:
            max_error_n = math.ceil(math.pi / math.acos(1 - max_error / self.radius))

        n = max(n_min, max_sep_n, max_error_n, 4) + closed
        thetas = np.linspace(0, 2*np.pi, n, endpoint=closed)
        pts = np.empty((n, 2))
        pts[:,0] = np.cos(thetas)
        pts[:,1] = np.sin(thetas)
        # TODO: generate points outside the circle to better approximate total line correctness
        return (pts * self.radius) + self.center


def create_donut(center, inner_radius, outer_radius):
    return Circle(center, outer_radius) & ~Circle(center, inner_radius)


class Ellipse(Area):
    def __init__(self, major, minor, center):
        self.major = np.asarray(major)
        self.minor = np.asarray(minor)
        self.center = np.asarray(center)

        # calculate defining vectors (for efficient boundary generation)
        self._f_a = self.major - self.center
        self._f_b = self.minor - self.center
        
        # calculate pins+string representation (for efficient containment checks)
        a = self.a = np.linalg.norm(self._f_a)
        b = self.b = np.linalg.norm(self._f_b)
        e = np.sqrt(1 - (b*b)/(a*a))
        df = e * self._f_a

        self._focus1 = self.center + df
        self._focus2 = self.center - df
        self._string_length = 2 * a
        
    def contains(self, pts):
        return np.linalg.norm(pts - self._focus1, axis=1) + np.linalg.norm(pts - self._focus2, axis=1) <= self._string_length

    def generate_boundary_points(self, n_min=12, max_sep=None, max_error=None, closed=False, optimise_perimeter=False):
        if max_error is None:
            max_error_n = 0
        else:
            max_error_n = NotImplemented  # TODO

        n = max(n_min, max_error_n, 8)

        if max_sep is not None:
            n += n % 4  # max_sep compensation assumes n is a multiple of 4
            # The longest segment is at the quarter turn approaching the minor axis
            while (max_sep_sq := (self.a * np.cos(np.pi/2 - 2*np.pi/n))**2 + (self.b * (1 - np.sin(np.pi/2 - 2*np.pi/n)))**2) > max_sep**2:
                n += 4
        
        thetas = np.linspace(0, 2*np.pi, n+closed, endpoint=closed)
        # TODO: generate points outside the ellipse to better approximate total line correctness
        return self._f_a * np.cos(thetas)[:,np.newaxis] + self._f_b * np.sin(thetas)[:,np.newaxis] + self.center

    @classmethod
    def from_extents(cls, a, b, center=(0,0), theta=0):
        unit_direction = np.array([np.cos(theta), np.sin(theta)])
        major = a * np.array([np.cos(theta), np.sin(theta)]) + center
        minor = b * np.array([-np.sin(theta), np.cos(theta)]) + center
        return cls(major, minor, center)


class Polygon(Area):
    def __init__(self, boundary_points, closed=False):
        if closed:
            self._boundary_points = np.asarray(boundary_points)
        else:
            self._boundary_points = np.vstack([boundary_points, [boundary_points[0]]])
        self._N = len(self._boundary_points) - 1

    def contains(self, pts):
        return points_in_polygon(pts, self._boundary_points)

    def generate_boundary_points(self, n_min=3, max_sep=None, max_error=None, closed=False):
        if max_sep is None and n_min <= self._N:
            return self._boundary_points if closed else self._boundary_points[:-1]

        line_lengths = np.linalg.norm(np.diff(self._boundary_points, axis=0), axis=1)
        n_segment_points = np.empty(self._N, dtype=int)
        if max_sep is None:
            n_segment_points.fill(2)
        else:
            n_segment_points[:] = np.ceil(line_lengths / max_sep) + 1

        # Try to equally scale up the points per segment if there aren't enough
        n_ratio = (n_segment_points.sum() - (self._N + 1)) / n_min  # Don't double-count the endpoints
        if n_ratio > 1:
            ratiod_point_counts = n_segment_points * n_ratio
            if (proposed_segment_points := np.round(ratiod_point_counts)).sum() >= n_min:
                n_segment_points[:] = proposed_segment_points
            else:
                n_segment_points[:] = np.ceil(ratiod_point_counts)

        boundary_points = [
            *chain(
                *(np.linspace(p1, p2, n)[:-1]
                  for (p1, p2, n) in 
                  zip(self._boundary_points[:-1], np.roll(self._boundary_points[:-1], -1), n_segment_points)
                )
            )
        ]
        if closed:
            boundary_points.append(boundary_points[0])

        return np.array(boundary_points)


if __name__ == '__main__':
    import matplotlib.pyplot as plt
    
    pts = []

    outer_path = rectangle_fill(5, 15, 0.3, 0.4)
    areas = [
        AxisAlignedRectangle(0.5, 1.5, 1, 2.4),
        Circle((4, 3), 1.5),
        Circle((7, 3), 1) | AxisAlignedRectangle(7, 9, 1, 4),
        Ellipse.from_extents(2, 0.5, (12, 4), np.pi/8),
        Polygon([[10,0],[10.5,2.5],[12,2],[14,3.5],[13,2],[15,1]]),
    ]

    cuts = [
        cut(outer_path, area, *args, copy=True)
        for area, args in zip(areas, [
            (False, True),
            (True, False),
            (False, True),
            (False, True),
            (False, True),
        ])
    ]

    fig, ax = plt.subplots()
    ax.set_aspect('equal')
    x, y = outer_path.T
    ax.plot(x, y, marker='x', color='blue', linewidth=1, markersize=3)
    print('area:')
    for index, area in enumerate(areas):
        print(f'\b', index, end='')
        x, y = area.generate_boundary_points(50, closed=True).T
        ax.plot(x, y, linestyle='--', color='black', linewidth=3)

    print('\ncut:')
    for index, c in enumerate(cuts):
        print('\b', index, end='') 
        x, y = c.T
        ax.plot(x, y, marker='o', color='orange', linewidth=2, markersize=6)
    print()
    plt.show()

    my_rect = AxisAlignedRectangle(-1,1,-2,2)
    print(f'{[0,0] in my_rect = }')
    print(f'{[-2,1] in my_rect = }')
    my_donut = create_donut((1,1), 1, 2)
    print(f'{my_donut.contains(np.array([[1,1],[2.5,1], [1, -0.5], [5,5]])) = }')
    print(f'{my_donut.intercept((2.5,1),(1,1), binary_preferred=True) = }')
    print(f'{my_donut.intercept((2.5,1),(2.5,5)) = }')

    # Plan:
    ... # rectilinear infill of an axis-aligned rectangle, defined with boundaries
    ... # rectilinear infill of a concave polygon
    ... # rectilinear infill of a shape defined programmatically/mathematically
    ... # gyroid infill of a polygon
    ... # rectilinear polygon infill generated with thickness dependent on a voxel grid (or image)
    ... # gyroid infill of programmatic shape with thickness dependent on a vector field

Update (28/Feb/2024)

• Completed Ellipse implementation
• Added Polygon support
  • It's now possible to input arbitrary shapes from a list of points, instead of needing to define your own area containment functions

Usage Notes:

• Concavity is not yet handled, but if the shape is defined such that the infill path doesn't need to pass through successive regions of the shape then it should still work fine
• Boolean operations work for containment checks, but don't yet have decent handling of boundary generation (which is a nice-to-have for visualisation purposes)
• Sharp/thin features can be missed by the fill cutting algorithm if the fill's max_feature_size means that no points are generated inside the thin region