include <BOSL2/std.scad>
include <BOSL2/rounding.scad>
/**
    Calculate an ambiguous cylinder. This is a cylinder that looks different 
    from one side than the other. This was idea invented by Sugihara.
    
    https://www.cs.jhu.edu/~misha/ReadingSeminar/Papers/Sugihara18.pdf
    
    This method is inspired by https://www.youtube.com/watch?v=sQ-Fsv8vVKo
*/
$fn=1000;

sz      = 15;    
height  = 30;
w1      = 30; 
shape1  = ngon(2*sz+0.001, 40, ratio=cos(w1)); 
w2      = -20; 
shape2  = ngon(2*sz, 4, ratio=cos(w2)); 

// note to self: make sure the shapes have the same width
   
intersection = intersect( shape1, shape2, w1, w2); // Rudolf's code
intersection2d = path2d(intersection);
cleaned = hull2d_path(intersection2d, true);
top = [ for ( i=cleaned ) intersection[i] + [0,0,height] ];

path = smooth_path(
    resample_path(
        path_merge_collinear(top,closed=true)
        ,200,closed=true), 
    method="corners", closed=true);
    
path_proj=path2d(path);
inside=offset(path_proj,delta=1.5, same_length=true);
inside_top = hstack(inside, column(path,2));
vnf = vnf_vertex_array(
    [path, path3d(path_proj), path3d(inside), inside_top],
    col_wrap=true,
    row_wrap=true);

vnf_polyhedron(vnf);



// debug
*union()  {
    %zmove(height) {
        rotate([w1, 0, 0])linear_extrude(10*sz, center= true) polygon(vec2(shape1)); 
        rotate([w2, 0, 0])linear_extrude(10*sz, center= true) polygon(vec2(shape2)); 
   }

    *skin([bot,top], slices=10);
//    color("red") Points(bot,r=1);
//    color("red") Points(top,r=1);
}   


// --------------------------------------------------------
// Rudolf's code, modified circle and called it ngon


/////// code part /////////////////////////////

// intersection point of line AB and plane P1P2P3
function line_plane_intersection(A, B, P1, P2, P3) =
    let(
        N = cross(P2 - P1, P3 - P1), // normal of plane
        d = N*P1,                    // plane equation N * P = d
        AB = B - A,                  // line direction vector
        t_numer = d - N*A,
        t_denom = N*AB,
        t = t_numer/ t_denom
    )
    //assert(t_denom != 0, "N,AB are parallel") 
    (A + t * AB); 
    

// distance between point P and line L1,L2
function distance_point_to_line(L1, L2, P) =
  let(
      v = L2 - L1, 
      w = P - L1, 
      proj_length = (w*v) / (v*v), // Projektion von w auf v
      proj_point = L1 + proj_length * v, // Lotfußpunkt auf der Gerade
      distance = norm(P - proj_point) // Lotrechter Abstand
  ) distance;

// construct a circle shape
function ngon(r, $fn=$fn?$fn:360/$fa, ratio=0.7)=  [for(i=[0:360/$fn:359.99]) r*[cos(i), ratio*sin(i), 0]]; 

// rotate vector or list (of list of...) vectors around x, y, or z
function Rx(x, A) = A[0][0]!=undef?[for(i=A) Rx(x, i)]:
    A*[[1, 0, 0], [0, cos(x), sin(x)], [0, -sin(x), cos(x)]]; 
function Ry(y, A) = A[0][0]!=undef?[for(i=A) Ry(y, i)]:
    A*[[cos(y), 0, sin(y)], [0, 1, 0], [-sin(y), 0, cos(y)]]; 
function Rz(z, A) = A[0][0]!=undef? [for(i=A) Rz(z, i)]: 
    A*[[cos(z), sin(z), 0], [-sin(z), cos(z), 0], [0, 0, 1]]; 
    
// translate vector or list (of list of...) vectors with x
function Tx(x=0, v=undef) = v[0][0]!=undef?[for (i=v) Tx(x,i)]:v+[x,0,0]; 
function Ty(y=0, v=undef) = v[0][0]!=undef?[for (i=v) Ty(y,i)]:v+[0,y,0]; 

// reduce 3D point list to 2D
function vec2(v) = v[0][0]==undef?[v[0], v[1]]:[for(a=v) vec2(a)]; 

// show list of points as point cloud 
module Points(L, r)  assert(is_list(L) && is_list(L[0])) 
   for(p=L) translate(p) sphere(r); 

// intersect two 2D shapes at an angle w
function intersect(Sh1, Sh2, w1=w1, w2=w2) = 
    let(
        Sh1 = Rx(w1, Sh1), 
        Sh2 = Rx(w2, Sh2),
        N1  =  Rx(w1, [0,0,1]), 
        N2  = Rx(w2, [0,0,1]),
        L1  = len(Sh1), 
        L2  = len(Sh2),
    )
    [
        for(i = [0:L2-1], j = [0:L1-1])
        // define the current face of Sh2
            let(
                P1 = Sh2[i], 
                P2 = Sh2[(i+1)%L2],
                P3 = P1+N2, 
                P4 = P2+N2,
                d = norm(P1-P2),
                // calc the intersection of all lines of Sh1 with current face of Sh2
                A = Sh1[j], 
                B = A+N1,   
                I = line_plane_intersection(A, B, P1, P2, P3),  
                dist1 = distance_point_to_line(P1, P3, I),
                dist2 = distance_point_to_line(P2, P4, I),
            )
            // filter points that closely touch the face
            if(dist1+dist2<=d*1.0001) 
                    I
    ];