Subdivision in 2D

// Subdivision in 2D.scad
//
// Version 0.7
// November 13, 2024
// By: Stone Age Sculptor
// License: CC0 (Public Domain)
//
// This script is a test to explore subdivision.
// My goal was something that feels like NURBS.
// Easy:   The control points are just points
//         without angle.
// Smooth: The curve does not go through
//         the original points.
// Simple: Instead of calculating the final spline,
//         the functions go through the basic
//         iterations for smoothing.

$fn = 100;

ShowTests(15,"1");
ShowTests(-15,"weighted");
ShowPathTest();

// The functions work in 3D as well,
// but that is slow and can not be rendered
// at the same time with 2D.
*ShowTube();
*translate([14,0,5])
  rotate([0,180,0])
    ShowTube();

module ShowTube()
{
  tubepoints = [[-15,-2,5],[-3,-2,5],[0,-2,0],[10,-2,0],[18,-12,8],[18,12,8],[10,2,0],[0,2,0],[-3,2,5],[-15,2,5]];
  tubepath = Subdivision(tubepoints,4,method="weightedpath");
  tube_thickness = 3;

  // the tube in 3D
  for(i=[0:len(tubepath)-2])
  {
    hull()
    {
      translate(tubepath[i])
        sphere(d=tube_thickness);
      translate(tubepath[i+1])
        sphere(d=tube_thickness);
    }
  }
}

module ShowTests(y_offset, method)
{
  // The original shape to start with.
  shape = [[-20,-10],[-12,15],[8,7],[13,-3]];

  // The shape with points for explanation.
  translate([0,y_offset])
  {
    color("LightSteelBlue")
      translate([0,0,-2])
        polygon(shape);

    // Show the points of the shape in Blue.
    color("Blue")
      for(i=[0:len(shape)-1])
        translate(shape[i])
          circle(1);

    // Show the points in between in Purple.
    PurplePoints = [for(i=[0:len(shape)-1]) ((shape[i] + shape[i+1 > (len(shape)-1) ? 0 : i+1])/2)];
    color("Purple")
      for(i=[0:len(PurplePoints)-1])
        translate(PurplePoints[i])
          circle(1);

    if(method=="1")
    {
      // Method "1" is the 1,1 weighted subdivision

      // The text of the used method
      color("Black")
        translate([-65,0])
          text("\"1\"",size=6);

      // Show the new points between the Blue and the Purple points.
      OrangePoints = [ for(i=[0:len(shape)-1]) each [ (shape[i] + PurplePoints[i])/2, (PurplePoints[i] + shape[i+1 > (len(shape)-1) ? 0 : i+1])/2] ];
      color("Orange")
        for(i=[0:len(OrangePoints)-1])
          translate(OrangePoints[i])
            circle(1);
    }
    else
    {
      // The other method is the weighted method.

      // The text of the used method
      color("Black")
        translate([-65,0])
          text("\"weighted\"",size=6);

      // Show the average of two Purple points in Black
      black_point = (PurplePoints[0]+PurplePoints[1])/2;
      color("Black")
        translate(black_point)
          circle(1);

      // Show the line that the new points is on,
      // according to a weight variable.
      color("Gray")
        hull()
        {
          translate(black_point)
            circle(0.2);
          translate(shape[1])
            circle(0.2);
        }

      // Show the new points.
      RedPoints = [for(i=[0:len(shape)-1]) ((shape[i] + (PurplePoints[i] + PurplePoints[(i-1) < 0 ? (len(shape)-1) : i-1])/2)/2)];
      color("Red")
        for(i=[0:len(RedPoints)-1])
          translate(RedPoints[i])
            circle(1);

      // Show a Green triangle.
      color("LawnGreen")
        translate([0,0,-1])
          polygon([shape[1],PurplePoints[0],PurplePoints[1]]);
    }
  }

  // Example with the shape.
  translate([35,y_offset])
  {
    color("LightSteelBlue")
      polygon(shape);
    color("Red")
      translate([0,0,2])
        polygon(Subdivision(shape,4,method=method));
  }

  // Example with triangle.
  translate([65,y_offset])
  {
    triangle = [for(a=[0,120,240]) [15*sin(a),15*cos(a)] ];

    color("LightSteelBlue")
      polygon(triangle);

    color("Red")
      translate([0,0,1])
        polygon(Subdivision(triangle,4,method=method));
  }

  // Example with square.
  translate([90,y_offset])
  {
    rectangle = [[-10,-10],[10,-10],[10,10],[-10,10]];

    color("LightSteelBlue")
      polygon(rectangle);

    color("Red")
      translate([0,0,2])
        polygon(Subdivision(rectangle,5,method=method));

    color("Black")
      translate([0,0,4])
        difference()
        {
          circle(10.0);
          circle(9.0);
        }
  }

  // Example with a more complex shape
  translate([115,y_offset])
  {
    testshape = [[-10,-10],[-10,10],[0,-10],[10,10],[10,-10]];

    color("LightSteelBlue")
      polygon(testshape);

    color("Red",0.6)
      translate([0,0,2])
        polygon(Subdivision(testshape,5,method=method));
  }

  // Example with a heart shape
  translate([140,y_offset])
  {
    // Both sides are defined.
    // It is a single solid 2D object, therefor both
    // sides of the shape have to be defined.
    heartshape = [[0,6],[-0.5,10],[-10,9],[-10,0],[-1,-4],[0,-10],[0,-10],[1,-4],[10,0],[10,9],[0.5,10],[0,6]];

    color("LightSteelBlue",0.2)
      polygon(heartshape);

    color("Red",0.6)
      translate([0,0,2])
        polygon(Subdivision(heartshape,5,method=method));
  }
}

module ShowPathTest()
{
  // Examples with paths
  points = [[-20,10],[-20,-10],[0,-10],[-15,15],[5,15],[20,0]];

  translate([0,-50])
  {
    // Show the text.
    color("Red")
      translate([-65,5])
        text("\"1path\"",size=5);

    color("Blue")
      translate([-65,-5])
        text("\"weightedpath\"",size=5);

    for(i=[0,1])
    {
      x_offset = i==0 ? 5 : 50;
      divisions = i==0 ? 1 : 3;

      translate([x_offset,0])
      {
        // Show the points.
        color("Black")
          for(i=[0:len(points)-1])
            translate(points[i])
              circle(1);

        // Show the number of divisions.
        color("Black")
        {
          t = str("divisions=", divisions);
          translate([-3,0])
            text(t,size=2.5);
        }

        // Draw the path.
        path1 = Subdivision(points,divisions,method="1path");
        color("Red")
        {
          for(i=[0:len(path1)-2])
          {
            hull()
            {
              translate(path1[i])
                circle(0.2);
              translate(path1[i+1])
                circle(0.2);
            }
          }
        }

        path2 = Subdivision(points,divisions,method="weightedpath");
        color("Blue")
        {
          for(i=[0:len(path2)-2])
          {
            hull()
            {
              translate(path2[i])
                circle(0.2);
              translate(path2[i+1])
                circle(0.2);
            }
          }
        }
      }
    }

    // The shape of a horseshoe
    translate([90,-5])
    {
      horsepoints = [[-7.5,16],[-9,14],[-12,1],[0,-8],[12,1],[9,14],[7.5,16]];
      horsepath = Subdivision(horsepoints,4,method="weightedpath");
      color("Black")
        translate([0,0,1])
          for(i=[0:len(horsepoints)-1])
            translate(horsepoints[i])
              circle(0.5);

      color("Brown",0.5)
        offset(1.8)
          for(i=[0:len(horsepath)-2])
          {
            hull()
            {
              translate(horsepath[i])
                circle(0.1);
              translate(horsepath[i+1])
                circle(0.1);
            }
          }
    }

    // An arc
    translate([120,-10])
    {
      arcpoints = [[10,0],[10,12],[10,16],[1,21],[0,24]];
      arcpath = Subdivision(arcpoints,4,method="weightedpath");
      a2 = concat(MirrorList(arcpath),ReverseList(arcpath));

      color("Olive",0.6)
        difference()
        {
          polygon(a2);
          offset(-2)
            polygon(a2);
        }

      color("Black")
        translate([0,0,1])
        for(i=[0:len(arcpoints)-1])
          translate(arcpoints[i])
            circle(0.5);

      color("Gray")
        for(i=[0:len(arcpath)-2])
        {
          hull()
          {
            translate(arcpath[i])
              circle(0.1);
            translate(arcpath[i+1])
              circle(0.1);
          }
        }
    }
  }
}

// Reverse a list of points.
function ReverseList (list) =
  let(n=len(list)-1)
  [for(i=[0:n]) list[n-i] ];

// Mirror a list of 2D points.
// The points on the positive x-axis
// are mirrored to the negative x-axis.
function MirrorList (list) =
  let(n=len(list)-1)
  [for(i=[0:n]) [-list[i].x, list[i].y] ];


// ------------------------------------------------
// Subdivision functions
// ------------------------------------------------
// Version 1
// November 13, 2024
// By: Stone Age Sculptor
// License: CC0 (Public Domain)

// Subdivision
// -----------
// Parameters:
//   list:      A list of coordinates in 2D or 3D.
//              But not a 3D surface.
//   divisions: The number of divisions.
//              0 for no smoothing, 5 is extra smooth.
//   method:    The subdivision method:
//              "1"        A basic 1,1 weighted subdivision
//                         for a closed shape.
//              "weighted" A weighted average subdivision
//                         for a closed shape.
//              "1path"    A basic 1,1 weighted subdivision
//                         for a path.
//              "weightedpath"
//                         A weighted average subdivision
//                         for a path.
// Return:
//   A new list with a smoother shape.
//   The new list is often used with the polygon() function.
function Subdivision(list,divisions,method) =
  divisions > 0 ?
    Subdivision(
      method=="1"            ? _Subdivision11(list) :
      method=="1path"        ? _SubdivisionPath11(list) :
      method=="weighted"     ? _SubdivisionWeighted(list) :
      method=="weightedpath" ? _SubdivisionPathWeighted(list) :
      assert(false,"Unknown method in Subdivision"),
      divisions - 1,method) :
      list;

// This is the most basic subdivision with 1,1 weighting.
// It will work in 2D and 3D.
//
// The average in OpenSCAD is: (current_point + next_point) / 2
// The index of the list wraps around for the next point.
//
// The 'list2' is the average points between the original points.
// The returned list is the new average between the average points
// and the original points.
function _Subdivision11(list) =
  let (n = len(list)-1)
  let (list2 = [ for(i=[0:n]) (list[i] + list[(i+1) > n ? 0 : i+1])/2 ])
  [ for(i=[0:n]) each [ (list[i] + list2[i])/2, (list2[i] + list[(i+1) > n ? 0 : i+1])/2 ]];

// My own attempt with variable weighting.
// The goal was a smoothing algoritme that feels
// like NURBS.
// Explanation:
//   The average points between the original
//   points are calculated. These are kept.
//   A second set of average points between
//   those average points are temporarely calculated.
//   A new point is created on the line between
//   an original point and the temporarely point.
//   The position on that line is defined by
//   a 'weight' variable.
//   The result is the combination of the
//   kept points and the new points.
//
// When the 'weight' variable is set to
// sqrt(2) - 1, then the result approximates
// a circle when the control points is a square.
// It is some kind of cubic B-spline, but I don't
// know if it matches with one of the known algoritmes.
function _SubdivisionWeighted(list) =
  let (weight = sqrt(2) - 1)
  let (n = len(list)-1)
  let (list2 = [ for(i=[0:n]) (list[i] + list[(i+1) > n ? 0 : i+1])/2 ])
  let (list3 = [ for(i=[0:n]) (weight*list[i] + (1-weight)/2*(list2[i] + list2[(i-1) < 0 ? n : i-1])) ])
  [ for(i=[0:n]) each [list3[i], list2[i]] ];

// The basic subdivision with 1,1 weighting.
// But now for a path with a open begin and end.
function _SubdivisionPath11(list) =
  let (n = len(list)-2)
  let (list2 = [ for(i=[0:n]) (list[i] + list[i+1])/2 ])
  [ list[0], for(i=[1:n]) each [ (list[i] + list2[i-1])/2, (list[i] + list2[i])/2 ], list[n+1]];

// My own attempt with variable weighting.
// But now for a path with a open begin and end.
function _SubdivisionPathWeighted(list) =
  let (weight = sqrt(2) - 1)
  let (n = len(list)-2)
  let (list2 = [ for(i=[0:n]) (list[i] + list[i+1])/2 ])
  let (list3 = [ list[0], for(i=[1:n]) (weight*list[i] + (1-weight)/2*(list2[i] + list2[i-1])), list[n+1] ])
  [ for(i=[0:n]) each [list3[i], list2[i]], list3[n+1] ];