// Demo code for affine intersection path of two linear extruded shapes 
// 

// the two shapes. !! In case of instabilities rotate shape by a small angle
shape1 = Tx(7, circle(10, 120)); 
shape2 = Tx(10, Rz(10, circle(12.2, 3))); 

echo(len(shape1)); 
w1 = 60;  // angle shape1
w2 = -39;   // angle shape2

// display shapes
#rotate([w1, 0, 0])linear_extrude(100, center= true) polygon(vec2(shape1)); 
#rotate([w2, 0, 0])linear_extrude(100, center= true) polygon(vec2(shape2)); 
   
// calc intersection paths   
res1 = intersect( shape1, shape2, w1, w2); 
res2 = intersect( shape2, shape1, w2, w1); 

// display point the two clouds
if(res1!=[])color("blue")   Points(res1, r=.5); 
if(res2!=[])color("green") Points(res2, r=.6); 


/////// 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 circle(r, $fn=$fn?$fn:360/$fa)=  [for(i=[0:360/$fn:359.99]) r*[cos(i), 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))
let(N1 =  Rx(w1, [0,0,1]), N2 = Rx(w2, [0,0,1]))
let(L1 = len(Sh1), L2=len(Sh2))
let(sign = L1<L2)
let(result =
  [
    for(i = [0:L2-1], j = [0:L1-1])
    // define the current face of Sh2
    let(P1 = Sh2[i], P2=Sh2[(i+1)%L2]) 
    let(P3 = P1+N2, P4 = P2+N2)
    let(d = norm(P1-P2))
    
    // calc the intersection of all lines of Sh1 with current face of Sh2
    let(A = Sh1[j], B = A+N1)   
      let(I = line_plane_intersection(A, B, P1, P2, P3))  
      let(dist1 = distance_point_to_line(P1, P3, I))
      let(dist2 = distance_point_to_line(P2, P4, I))
    // filter points that closely touch the face
      if(dist1+dist2<=d*1.0001) 
      I
  ])
  result; 
;