06-16-2016, 04:24 PM
06-23-2016, 09:32 PM
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.
06-26-2016, 11:23 AM
Some more explanation of what you mean? A visual explanation maybe?
06-27-2016, 02:17 PM
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.
http://repositorium.sdum.uminho.pt/bitst...CM_MYS.pdf
There is a 3D approach that I'm also working on ... which uses TETGEN.
07-01-2016, 05:06 PM
07-01-2016, 10:39 PM
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.