Anyone play with 3D Alpha Shapes ?

[ggaliens]

Anyone play with 3D Alpha Shapes ?

Generalized method for hulling ?

[ggaliens]

I'm in the middle of adding support of importing point-clouds in PLY file format. After that happens ... and alpha shape will be used to approximate the boundary of the point cloud.

[Dimitri]

Some more explanation of what you mean? A visual explanation maybe?

[ggaliens]

Below link is 2D approach for concave hulls. I should have this implemented soon.

http://repositorium.sdum.uminho.pt/bitst...CM_MYS.pdf


There is a 3D approach that I'm also working on ... which uses TETGEN.

[ggaliens]

[ggaliens]

The code in case anyone is interested in algorithms ...

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.