Anyone play with 3D Alpha Shapes ?

Code:
-module(wpc_k_nearest).

%% This module attempts to implement concepts describe in ...

%%  http://repositorium.sdum.uminho.pt/bitstream/1822/6429/1/ConcaveHull_ACM_MYS.pdf



%% It appears to work on some base case ... test cases ... but also fails 

%%  Not fully tested yet. Will do at least 8K points fairly quickly (several seconds)b



-export([

    k_nearest/2, 

    k_nearest_to_polygon/2,

    jiggle/1,

    degrees_oriented/3,

    sample/0,

    bisect_segments/1,

    intersects_polyline/3,

    intersecting_segments/4

    ]).



-export([init/0]).





-include("e3d.hrl").





init() -> false.



%% we can pass parameters in more of an eclectic manner using records

-record(bounds, 

     {

     k=4,

     dir={-1.0,0.0,0.0},  % for sure will be left of all points after initial sort.

     pt=nil,              % don't have one yet ... set to nil ... because why not !

     pts=[],              

     first_k=[]  %% first k candidates for first point !

     }).





%% Goldfish !

sample() -> 

   XY = [

      { 0, 0},

      { 1, 1 },

      { 2, 2 },

      { 3, 3 },

      { 2, 4 },

      { 1, 5 },

      { 2, 6 },

      { 3, 7 },

      { 4, 8 },



      {  2, 8 },

      {  0, 8 },

      { -2, 8 },





      { -4, 8 },

      { -3, 7 },

      { -2, 6  },

      { -3, 5 },

      { -4, 4 },

      { -3, 3},

      { -2, 2 },

      { -1, 1 }

      ],

  [ {1.0*X, 1.0*Y, 0.0}   || {X,Y} <- XY ].  % shuffle to make it point-cloudier !





%% Use this to increase problem size on goldfish or sample !

%% Just cuts edges into two equal parts

bisect_segments([]) -> [];

bisect_segments([{X1,Y1,Z1}]) ->  [{X1,Y1,Z1}];

bisect_segments([{X1,Y1,Z1},{X2,Y2,Z2}]) ->  

    Center = e3d_vec:average([{X1,Y1,Z1},{X2,Y2,Z2}]),

    [{X1,Y1,Z1},Center,{X2,Y2,Z2}];

bisect_segments([{X1,Y1,Z1},{X2,Y2,Z2}|T]) ->  

    [A,B,_] = bisect_segments([{X1,Y1,Z1},{X2,Y2,Z2}]),

    List2 = bisect_segments([{X2,Y2,Z2}|T]),

    lists:append([A,B], List2).





jiggle([]) -> [];    % jiggle points in our sample cloud a little bit. 

jiggle([{X,Y,Z}|T]) -> 

   [{X + random:uniform()*0.01, Y + random:uniform()*0.01, Z}|jiggle(T)].





%% This module attempts to implement concepts describe in ...

%%  http://repositorium.sdum.uminho.pt/bitstream/1822/6429/1/ConcaveHull_ACM_MYS.pdf



%% Pts should be a fairly planar cloud of points. Very flat.

%% Also should be oriented so that the cloud is tangent to y=0

%  for now we will leave these responsibilties to the calling module functions. 

%% Also ... cloud should be laying in z=0 plane. 

k_nearest(K,Pts) when is_integer(K) ->

    Pts2 = lists:usort(fun({_,Y1,_},{_,Y2,_}) -> Y1 < Y2 end, Pts),

    [Pt|_] = Pts2,

    k_nearest(#bounds{k=K,pt=Pt}, Pts);

k_nearest(#bounds{pts=Boundary,pt=Pt}=B,Pts0) ->

    %%io:format("------------------------------------------------------\n", []),

    case length(Boundary) rem 25 of 

      0 -> io:format("Boundary Len = ~p\n", [ length(Boundary) ]);

      _ -> ok

    end, 

    [{_,Pt2}|_] = Candidates = candidates(#bounds{}=B, Pts0 -- [Pt|Boundary]),

    FK = 

    case B#bounds.first_k of 

      [] -> [Pt0||{_Ang,Pt0}<-Candidates];

      [_|_]=FK0 -> FK0

    end,

    %%io:format("First K = ~p\n", [ FK ]),





    %%io:format("Boundary = ~p\n", [ Boundary ]),

    %%io:format("Point Cur = ~p\n", [ Pt ]),

    %%io:format("Dir Cur = ~p\n", [ Dir0 ]),

    



    %%io:format("Point Next = ~p\n", [ Pt2 ]),

    



    %% upon looping around the sampled border and when it is closing down ... 

    %% first boundary point (last onto stack) will be a member of the first k candidates looked at.

    %% for some reason ... using the first pt is not loose enough logic. 

    DirNext = e3d_vec:sub(Pt2,Pt),  

    case Boundary of  

        [B1,_|_] ->

               B1_MEM = gb_sets:is_member(B1, gb_sets:from_list(FK)),

               case B1_MEM of

                  false ->  

                      k_nearest(B#bounds{pts=[Pt2|Boundary],dir=DirNext,pt=Pt2,first_k=FK},Pts0);

                  true -> 

                      Boundary

               end;

        _      ->  k_nearest(B#bounds{pts=[Pt2|Boundary],dir=DirNext,pt=Pt2,first_k=FK},Pts0)  

    end.





k_nearest_to_polygon(K, Vs) ->

    Pts = k_nearest(K,Vs),

    polygon_to_mesh(Pts).



    

%% Read about 2d alpha hulls. We want new points that form edges to right (most-clockwise) from 

%% last boundary vector ... even if angle is negative (left turn). 

%% We want to return K closest points and sort then my angle of turn ascending.

%% Lastly ... don't allow any new segments that would intersect the old boundary in a non

%% trivial way. Tangent segments are OK ? 

candidates(#bounds{k=K,pts=Boundary,dir=Dir0,pt=Pt0},Pts0a) ->

    SortByDistance = fun ({X,Y,_Z}, Pts0) ->

        List0 = 

           [{(X-X2)*(X-X2)+(Y-Y2)*(Y-Y2), {X2,Y2,Z2}} || {X2,Y2,Z2} <- Pts0 ],

        [_Pt0||{_,_Pt0}<-lists:sort(List0)]

    end,

    Pts0 = SortByDistance(Pt0, Pts0a),

    Normal = {0.0,0.0,1.0}, 

    MyAcc = fun({X0,Y0,Z0}=Pt2,Acc) -> 

       Y = Dir0,

       Z = e3d_vec:sub({X0,Y0,Z0}, Pt0),

       Ang = degrees_oriented(Normal,Y,Z),

       [{Ang,Pt2}|Acc]

    end, 

    List0 = lists:foldl(MyAcc, [], lists:sublist(Pts0,K)),

    %%io:format("Angles Dict = ~p\n", [ List0 ]),

    Candidates = lists:reverse(lists:sort(List0)),

    lists:filter(fun ({_Ang,Pt2}) -> not(intersects_polyline(Pt0,Pt2,Boundary))  end, Candidates).





%%  Google This :  atan2(norm(cross(a,b)),dot(a,b))

%%  The norm( ) shown above is length (e3d_vec:len)

%%  Test:  TESTING THIS IS EXTREMELY IMPORTANT !  We want -pi to +pi  returned here ... pretty sure (MEW).

%%  wpc_k_nearest:degrees_oriented({0.0,0.0,1.0}, {-1.0,0.0,0.0}, {1.0,1.0,0.0} ). 

-spec degrees_oriented(e3d_vector(), e3d_vector(), e3d_vector()) -> float().

%%   Oriented is more generic than signed. So we shall just call it oriented.

%%   Nrm ... the vector perpendicular to plane in which A and B vectors lay

%%   There is an obvious advantage to asserting the direction Nrm ... so use this 

%%   approach as often as possible.

degrees_oriented({_,_,_}=X, {_,_,_}=Y, {_,_,_}=Z) ->

    C = e3d_vec:cross(Y,Z),

    Radians = -1.0*signum(e3d_vec:dot(X,C))*math:atan2(e3d_vec:len(C),e3d_vec:dot(Y,Z)),

    180.0 * Radians / math:pi().

signum(X) when X < 0.0 ->  -1.0;

signum(X) when X > 0.0 ->  1.0;

signum(_X) ->  1.0.    %% grabbing at straws here. I think we want a SIGN one way or another !





%% Tests : 

%%  True case : wpc_k_nearest:intersecting_segments({0.0,0.0,0.0},{1.0,1.0,0.0}, {1.0,0.0,0.0},{0.0,1.0,0.0}).

%%  False case : wpc_k_nearest:intersecting_segments({0.0,0.0,0.0},{1.0,1.0,0.0}, {1.0,9.0,0.0},{0.0,9.0,0.0}).

intersects_polyline({_X1,_Y1,_Z1},{_X2,_Y2,_Z2},[]) -> false;

intersects_polyline({_X1,_Y1,_Z1},{_X2,_Y2,_Z2},[{_,_,_}]) -> false;

intersects_polyline({X1,Y1,Z1},{X2,Y2,Z2},[{_,_,_}=A,{_,_,_}=B|T]) -> 

    case intersecting_segments({X1,Y1,Z1},{X2,Y2,Z2},A,B) of 

        true   -> true;

        false  ->  intersects_polyline({X1,Y1,Z1},{X2,Y2,Z2},[B|T])

    end.





segments_tangent({X1,Y1,Z1},{X2,Y2,Z2},{X3,Y3,Z3},{X4,Y4,Z4}) ->

    D1 = e3d_vec:line_dist({X1,Y1,Z1},{{X3,Y3,Z3},{X4,Y4,Z4}}),

    D2 = e3d_vec:line_dist({X2,Y2,Z2},{{X3,Y3,Z3},{X4,Y4,Z4}}),

    D3 = e3d_vec:line_dist({X3,Y3,Z3},{{X1,Y1,Z1},{X2,Y2,Z2}}),

    D4 = e3d_vec:line_dist({X4,Y4,Z4},{{X1,Y1,Z1},{X2,Y2,Z2}}),

    lists:any(fun(D) -> abs(D)  < 0.001 end, [D1,D2,D3,D4]).





%% The main function that returns true if line segment 'p1q1'

%% and 'p2q2' intersect.

intersecting_segments({X1,Y1,Z1},{X2,Y2,Z2},{X3,Y3,Z3},{X4,Y4,Z4}) ->

    Tangent = segments_tangent(     {X1,Y1,Z1},{X2,Y2,Z2},{X3,Y3,Z3},{X4,Y4,Z4}),

    HasOne = intersecting_segments0({X1,Y1,Z1},{X2,Y2,Z2},{X3,Y3,Z3},{X4,Y4,Z4}),

    (HasOne andalso not(Tangent)).

intersecting_segments0(P1, Q1, P2, Q2) ->

    O1 = ccw(P1, Q1, P2),

    O2 = ccw(P1, Q1, Q2),

    O3 = ccw(P2, Q2, P1),

    O4 = ccw(P2, Q2, Q1),

    (O1 /= O2) andalso (O3 /= O4).

ccw({Ax,Ay,_Az},{Bx,By,_Bz},{Cx,Cy,_Cz}) ->

    (Cy-Ay) * (Bx-Ax) > (By-Ay) * (Cx-Ax).







polygon_to_mesh(Pts) ->

    Face1 = #e3d_face{vs=lists:seq(0,length(Pts)-1)},

    Face2 = Face1#e3d_face{

         vs=lists:reverse(Face1#e3d_face.vs)

         },

    Mesh1 = #e3d_mesh{},

    Mesh2 = Mesh1#e3d_mesh{

         vs=Pts,

         fs=[Face1,Face2]

         },

    %%Mesh3 = e3d_mesh:clean_faces(Mesh2),

    %%e3d_mesh:merge_vertices(Mesh3).

    Mesh2.