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

Some ideas related to FullControlXYZ/fullcontrol#68

This is still very much WIP - just uploading so the progress can be seen, and parts can potentially be extracted from it for use already.

There are several additional features, functionalities, and quality of life improvements I'm planning to implement, and it's likely there are some bugs in the existing code.

Current status is at least capable of this:

Some progression ideas

• Polygon contains should be possible using Sasa Milenkovic's Point in Polygon code (GPLv3 licensed 👍)
  • can be optionally compiled and parallelised with numba as presented here. (Care should be taken around open/closed Polygon definitions
  • it would be useful to store a closed polygon in _boundary_points, and have a closed (or open) boolean flag option in the generate_boundary_points function)
• I want to check the pure python performance
  • given Google Colab comes with numba pre-installed it should at least be fast for most people who need it to
  • I'm also tangentially interested in whether the code can be made a bit more understandable, but that's not particularly important given it's being taken from someone else's project
• It's likely advisable to also implement a performant/exact intercept method
  • the default binary containment search will likely be sufficient as a starting point, in which case it makes more sense to focus on other features first
• AxisAlignedRectangle.distance_to_boundary method should be straightforward (just clamp x and y to within bounds and it should end up with the closest point on the boundary)
• AreaGroup.generate_boundary_points methods can filter out some child area points by ensuring they meet the containment requirements for the group type
• It should be possible to brute-force generate boundary points using containment checks and intercepts with a grid size smaller than the smallest expected feature difference
  • That's calculation heavy and may be a bit buggy, but could be a valuable starting point for people trying to use arbitrary functions that may not have nicely defined boundaries
  • It may be a bit complex to handle areas with multiple boundaries (e.g. holes, or multiple shapes)
• In a related vein, if a mathematically defined area is (very) complex to compute containment checks for, it may be worth inheriting from Polygon and generating a polygonal boundary once that can then be used for subsequent containment checks for cutting fill designs and the like
  • This probably isn't useful for a single layer, but could be for designs that will generate several layers with the same/similar boundaries
• It's likely worth adding an Area.generate_bounding_rectangle(angle) method to help with the initial infill generations (prior to cutting)
• Concavity likely needs to be handled by marking the movements over gaps as 'travel' moves, since the movements do still need to happen
  • Determining these may be somewhat expensive, as the movement line for every path cut would need to be checked for out of bounds points down to a distance tolerance (perhaps with repeated bisections?)
  • Marking these should be handled with an extra array, or a data-frame column or something
• Dynamic path thickness would be useful for controlling fill density
  • Relevant to FullControlXYZ/fullcontrol#23 when based on vector fields / programmatic functions
  • Relevant to some ideas from FullControlXYZ/fullcontrol#15 when based on raster images
• Specifying start and end paths to interpolate between could result in a very generic approach to creating CONVEX-like fills
  • Non-straight edges could be used as boundaries to follow along at each side
    
  • Number of length-wise lines could either be explicitly specified, or determined automatically from minimum and maximum thicknesses as bounds
  • It should be possible for many polygons to be automatically split into 2 pairs of opposing edge groups, where the shorter pair are used as boundaries
    • Highly non-convex shapes (e.g. a cross or T) may be best split into more than one fill region, but doing that automatically seems hard
  • Handling holes would likely be quite complex, but could be worth implementing
    • This could maybe be done as a post-processing push-aside or cut operation
• For FullControl integration it may be relevant to create a nice data-structure for a GCode object
  • annotations would need to be kept separately but with known indices, such that path modifiers (e.g. rotations, translations, etc) can act efficiently on all components involving physical positions

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