Winding Number Algorithm - Openscad

windingnumber.scad

$fa = 12;                  // minimum angle
$fs = 2;                   // minimum size
$fn = $preview ? 24 : 64;  // number of fragments

micron = 0.001;      // for showing nice on screen OpenSCAD
print_margin = 0.3;  // print_margin for plastic expansion and some free room

// inside or outside: https://en.wikipedia.org/wiki/Point_in_polygon#Winding_number_algorithm
// radius vs distance to polygon: https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line

// Define the coordinates of the enclosed area
area_points = [ [ 10, 10 ], [ 100, 10 ], [ 80, 80 ], [ 40, 120 ], [ 10, 100 ], [ 10, 10 ] ];

diameter = 18;
save_radius = diameter / 2 + 0.55;
radius = (diameter  + 0.55)/2;

function distanceToLineSegment3(x, y, projection_x, projection_y) =
  sqrt((x - projection_x) * (x - projection_x) + (y - projection_y) * (y - projection_y));
                
function distanceToLineSegment2(x, y, x1, y1, x2, y2, dx, dy, vx, vy, length_squared, dot_product) =
  (dot_product <= 0) ? sqrt(vx*vx + vy*vy) : ((dot_product>=length_squared) ? sqrt((x-x2) * (x-x2) + (y-y2) * (y-y2)) : distanceToLineSegment3(x, y, x1 + dot_product / length_squared * dx, y1 + dot_product / length_squared * dy));

function distanceToLineSegment(x, y, x1, y1, x2, y2) =
  distanceToLineSegment2(x, y, x1, y1, x2, y2, x2-x1, y2-y1, x-x1, y-y1, (x2-x1)*(x2-x1)+(y2-y1)*(y2-y1), (x-x1)*(x2-x1)+(y-y1)*(y2-y1));

function windingNumber2(x, y, x1, y1, x2, y2) =
  (x < (x1 + ((y - y1) / (y2 - y1)) * (x2 - x1))) ? 1 : 0;

function windingNumber(x, y, x1, y1, x2, y2) = 
  (((y1 <= y) && (y2 > y)) || ((y1 > y) && (y2 <= y)))
  ? windingNumber2(x, y, x1, y1, x2, y2)
  : 0;

function PointInsideArea(x, y, area_points, i) = 
  (i == len(area_points) - 1)
  ? 0
  : PointInsideArea(x, y, area_points, i + 1) + windingNumber(x, y, area_points[i][0], area_points[i][1], area_points[i + 1][0], area_points[i + 1][1]);

function isPointInsideArea(x, y, area_points) = (PointInsideArea(x, y, area_points, 0) % 2 == 1) ? true : false;

function RadiusToLine(x, y, area_points, r, i) =
  (i == len(area_points)-1) 
  ? 999999 
  : min(RadiusToLine(x, y, area_points, r, i+1), distanceToLineSegment(x, y, area_points[i][0], area_points[i][1], area_points[i+1][0], area_points[i+1][1]) - r);

function isRadiusInsideArea(x, y, area_points, r) =
  RadiusToLine(x, y, area_points, r, 0);

difference() {
  color("black") linear_extrude(1) polygon(area_points);
  for(xi=[0.22: 1: 8]) {
    for (yi = [-0.05: 1: 8]) {
      x = 2*radius*sqrt(3/4)*xi;
      y = ((xi%2)*0.5+yi) * 2*radius;
    
      if ((isPointInsideArea(x, y, area_points) == true) && (isRadiusInsideArea(x, y, area_points, save_radius) > 0)) {
        translate([ x, y, -1 ])
        cylinder(h = 3, d = diameter);
      }
    }
  }
}