Alternative sphere implementations in OpenSCAD userspace

alternative_spheres.scad

translate([-30,0,0]) poly3d(sphere(r=9,$fn=80));
translate([-10,0,0]) poly3d(normalized_cube(r=9,div_count=20));
translate([10,0,0]) poly3d(spherified_cube(9,div_count=20));
translate([30,0,0]) poly3d(icosahedron(9,n=4));

module poly3d(p) {
  polyhedron(points=p[0],faces=p[1]);
}

function normalize(v) = v / norm(v); // convert vector to unit vector
function flatten(l) = [ for (a = l) for (b = a) b ];
function range(a,b) = let (step = a > b ? -1 : 1) [ for (i = [a:step:b]) i ];

function translate(v, points) = 
  let(
    lp = len(points[0]), // 2d or 3d point data?
    lv = len(v),         // 2d or 3d translation vector?
    V = lp > lv ? [v.x,v.y,0] : v // allow adding a 2d vector to 3d points
  )
  (lv > lp) ? 
    [for (p = points) let(pz = (p.z == undef) ? 0 : p.z) [p.x,p.y,pz]+v] : // allow adding a 3d vector to 2d point data
    [for (p = points) p+v];

// rotate a list of 3d points by angle vector v
// vector v = [pitch,roll,yaw] in degrees
function rotate(v, points) = 
  let(
    lv = len(v),
    V = (lv == undef ? [0,0,v] : 
      (lv > 2 ? v : 
        (lv > 1 ? [v.x,v.y,0] :
          (lv > 0 ? [v.x,0,0] : [0,0,0])
        )
      )
    ),
    cosa = cos(V.z), sina = sin(V.z),
    cosb = cos(V.y), sinb = sin(V.y),
    cosc = cos(V.x), sinc = sin(V.x),
    Axx = cosa*cosb,
    Axy = cosa*sinb*sinc - sina*cosc,
    Axz = cosa*sinb*cosc + sina*sinc,
    Ayx = sina*cosb,
    Ayy = sina*sinb*sinc + cosa*cosc,
    Ayz = sina*sinb*cosc - cosa*sinc,
    Azx = -sinb,
    Azy = cosb*sinc,
    Azz = cosb*cosc
  )
  [for (p = points)
    let( pz = (p.z == undef) ? 0 : p.z )
    [Axx*p.x + Axy*p.y + Axz*pz, 
     Ayx*p.x + Ayy*p.y + Ayz*pz, 
     Azx*p.x + Azy*p.y + Azz*pz]
  ];

// based on get_fragments_from_r documented on wiki
// https://en.wikibooks.org/wiki/OpenSCAD_User_Manual/The_OpenSCAD_Language
function fragments(r=1) = ($fn > 0) ? 
  ($fn >= 3 ? $fn : 3) : 
  ceil(max(min(360.0 / $fa, r*2*PI / $fs), 5));

// Draw a circular arc with center c, radius r, etc.
// "center" parameter centers the sweep of the arc about the offsetAngle (half to each side of it)
// "internal" parameter does enables polyhole radius correction
function arcPath(r=1, angle=360, offsetAngle=0, c=[0,0], center=false, internal=false) = 
  let (
    fragments = ceil((abs(angle) / 360) * fragments(r,$fn)),
    step = angle / fragments,
    a = offsetAngle-(center ? angle/2 : 0),
    R = internal ? r / cos (180 / fragments) : r,
    last = (abs(angle) == 360 ? 1 : 0)
  )
  [ for (i = [0:fragments-last] ) let(a2=i*step+a) c+R*[cos(a2), sin(a2)] ];

function circle(r=1, c=[0,0], internal=false, d) = 
  arcPath(r = d == undef ? r : d/2,c=c,angle=-360,internal=internal);

// return [points,faces] for a sphere
// should produce identical geometry to builtin openscad sphere module
function sphere(r=1,d) = 
  let(
    R = d == undef ? r : d/2,
    fragments = fragments(R),
    rings = floor((fragments+1)/2),
    c1 = circle(r=R),
    points = flatten([for (i = [0:rings-1])
      let(angle =  (180 * (i + 0.5)) / rings)
      translate([0,0,R*cos(angle)], sin(angle)*c1) 
    ]),
    lp = len(points),
    faces = concat(
      flatten([
        for (i = [0:rings-2], j = [0:fragments-1])
          let (
            il = i * fragments,
            il1 = il + fragments,
            j1 = (j == fragments-1) ? 0 : j+1,
            i0 = il + j, 
            i1 = il + j1, 
            i2 = il1 + j, 
            i3 = il1 + j1
          )
          [[i0,i2,i1],[i1,i2,i3]]
      ]),
      [range(0,fragments-1)], // top ring face
      [range(lp-1,lp-fragments)] // bottom ring face
    )
  )
  [points,faces];


function normalized_cube(r=1,div_count=12,d) = 
  let(
    R = d == undef ? r : d/2,
    origin=[0,0,1],
    right = [1,0,0],
    up = [0,1,0],
    div2 = div_count/2,
    face_points = [
      for (j = [0:div_count], 
           i = [0:div_count])
        let(
          face_point = origin + 2.0 * (right * (i-div2) + up * (div2-j)) / div_count
        )
        R*normalize(face_point)
    ],
    lface = len(face_points),
    points = concat(
      face_points,
      rotate([90,0,0],face_points),
      rotate([-90,0,0],face_points),
      rotate([180,0,0],face_points),
      rotate([0,-90,0],face_points),
      rotate([0,90,0],face_points)
    ), 
    faces = [
      for (f = [0:6-1],
           j = [0:div_count-1], 
           i = [0:div_count-1]) 
        let (
          i0=f*lface+j*(div_count+1)+i,
          i1=i0+1,
          i2=i0+(div_count+1),
          i3=i2+1
          //tri1 = 
        )
        [[i0,i1,i3],
         [i0,i3,i2]]
    ] // groups of two need to be concat'd
  )
  [points,flatten(faces)];


function spherified_cube(r=1,origin=[0,0,1],div_count=12,d) = 
  let(
    R = d == undef ? r : d/2,
    right = [1,0,0],
    up = [0,1,0],
    div2 = div_count/2,
    face_points = [
      for (j = [0:div_count], 
           i = [0:div_count])
        let(
          p = origin + 2.0 * (right * (i-div2) + up * (div2 - j)) / div_count,
          p2 = [p.x*p.x,p.y*p.y,p.z*p.z],
          rx = p.x * sqrt(1.0 - 0.5 * (p2.y + p2.z) + p2.y*p2.z/3.0),
          ry = p.y * sqrt(1.0 - 0.5 * (p2.z + p2.x) + p2.z*p2.x/3.0),
          rz = p.z * sqrt(1.0 - 0.5 * (p2.x + p2.y) + p2.x*p2.y/3.0)
        )
        R*[rx,ry,rz]
    ],
    lface = len(face_points),
    points = concat(
      face_points,
      rotate([90,0,0],face_points),
      rotate([-90,0,0],face_points),
      rotate([180,0,0],face_points),
      rotate([0,-90,0],face_points),
      rotate([0,90,0],face_points)
    ), 
    faces = [
      for (f = [0:6-1],
           j = [0:div_count-1], 
           i = [0:div_count-1]) 
        let (
          i0=f*lface+j*(div_count+1)+i,
          i1=i0+1,
          i2=i0+(div_count+1),
          i3=i2+1
          //tri1 = 
        )
        [[i0,i1,i3],
         [i0,i3,i2]]
    ] // groups of two need to be concat'd
  )
  [points,flatten(faces)];


// n = subdivision iterations (number of faces = 20 * 4^n)
// points constructed from golden rectangles as shown here https://commons.wikimedia.org/wiki/File:Icosahedron-golden-rectangles.svg
function icosahedron(r=1,d,n=0) =
  let(
    R = d == undef ? r : d/2,
    phi = (1+sqrt(5))/2,
    magnitude = norm([1,phi]),
    a = 1/magnitude,
    b = phi/magnitude,
    points = [
      [-a,b,0],[a,b,0],[a,-b,0],[-a,-b,0], // golden rectangle in XY plane
      [-b,0,a],[b,0,a],[b,0,-a],[-b,0,-a], // golden rectangle in XZ plane
      [0,-a,b],[0,a,b],[0,a,-b],[0,-a,-b]  // golden rectangle in YZ plane
    ],
    faces = [
      [0,9,4],[4,9,8],[8,9,5],[5,9,1],[1,9,0],
      [2,11,3],[3,11,7],[7,11,10],[10,11,6],[6,11,2],
      [8,5,2],[2,5,6],[6,5,1],[6,1,10],[10,1,0],[10,0,7],[7,0,4],[7,4,3],[3,4,8],[8,2,3]
    ],
    poly = n==0 ? [points,faces] : subdivide_faces([points,faces], n=n)
  )
  [[for (p = poly[0]) R*normalize(p)],poly[1]];

// faces better be triangles... or else!
function subdivide_faces(poly3d, n=1) = 
  let(
    points = poly3d[0],
    tris = poly3d[1],
    newpoints = flatten([for (t = tris) 
      let(
        p0 = points[t[0]],
        p1 = points[t[1]],
        p2 = points[t[2]]
      )
      [(p0+p1)/2, (p1+p2)/2, (p2+p0)/2]
    ]),
    lp = len(points),
    allpoints = concat(points, newpoints),
    faces = flatten([for (i = [0:len(tris) - 1])
      let(
        f = tris[i],
        p0 = f[0], p1 = f[1], p2 = f[2],
        m0 = lp + i * 3, m1 = m0 + 1, m2 = m1 + 1
      )
      [[p0,m0,m2],
      [m0,p1,m1],
      [m0,m1,m2],
      [m2,m1,p2]]
    ])
  )
  n > 1 ? subdivide_faces([allpoints,faces],n-1) : [allpoints,faces];